format: text

text:
# Towards the Final Fate of an Unstable Black String

## I Introduction

The stability of four-dimensional black holes is a well known, and fundamental result of relativity theory chandra . The picture in higher dimensional spacetimes was shown to be quite different by Gregory and Laflamme gregorylaflame1 ; gregorylaflame2 , who demonstrated the existence of unstable long wavelength modes of the black string in perturbation theory. This finding, coupled to arguments based on entropy considerations, led to their conjecture that black strings might bifurcate into higher dimensional analogues of spherical black holes. Cosmic censorship would be violated were this the case, since a singularity must be encountered by a bifurcating black hole horizon, essentially as a consequence of the principle of equivalence hawkingellis .

The existence of the Gregory-Laflamme instability has been assumed in many subsequent studies of higher dimensional gravity theory, including the classical limit of string theory (see, for example, hms96 ; imsy98 ; agmoo00 ; hv00 ), some of which have also assumed the validity of the bifurcation conjecture. However, a linearized analysis can say little, if anything, concerning the nature of the full non-linear evolution of an unstable string, and the final end-state of such a configuration remained to be established.

Recently, Horowitz and Maeda were able to prove, under some assumptions, that black strings cannot bifurcate in finite time horowitzmaeda . Furthermore, they conjectured that the system is likely to approach a new stationary solution which is not translationally invariant along the string direction. However, even if the assumptions involved in the proof are sufficiently generic, their analysis cannot identify the final end-state of evolution. Partial answers can be sought via perturbation analysis as done by Gubser gubser . By assuming the Horowitz-Maeda conjecture and linearizing the solution at the critical length to first order (and a partial extension to second order), Gubser argued that the transition to the final solution must be of second order type (i.e. discontinuous). Despite these developments, it seems clear that a convincing answer to the question at hand can only be obtained by solving the full equations governing the problem. A step in this direction would be to search for special solutions, such as stationary ones, and compare the physical content of the obtained configurations with the black string solutions. This has recently been carried out by Wiseman wiseman , who numerically solves the equations resulting from a static ansatz. Interestingly, he finds non-uniform solutions with mass larger than that of the black string for a given compactification radius. Thus, he concludes that the solutions he finds cannot be the end-states conjectured by Horowitz and Maeda (also see related work by Kol kol ).

Additional work by Unruh and Wald unruhwald  studies the dynamics of a uniform cylindrical configuration of matter in Newtonian gravity. They observe that a perturbation of the density gives raise to a Jeans instability responsible for the collapse of the system along the cylinder’s length. They then argue that if the main features of this model are robust in the passage to the general relativistic system, one possible end-state for the perturbed black string would be collapse in the string direction, resulting in singularity formation. Note that this collapse need not lead to violations of cosmic censorship, as the final singularity could still be hidden by an event horizon geddes .

Clearly there are several distinct viable possibilities for the final end-state of a perturbed black string, with remarkably different consequences associated with the range of options. Current conjectures range from “nothing interesting happens”, to violations of cosmic censorship, to the arguably more extreme case of a complete collapse of the spacetime. In order to completely settle the issue, the full dynamics of the perturbed black string needs to be addressed. At least in principle, this will allow us to identify which of the above possibilities (if any!) is actually realized. In this paper we report on preliminary work in this direction—a program to simulate the dynamics of the black string through numerical solution of the Einstein equations. At this stage of the project, we cannot yet provide an answer to the question of the end-state; however, we have tantalizing results that show the spacetime going through a configuration resembling 5-dimensional spherical black holes connected by thin black strings that expand along the string dimension. Our simulations eventually crash while the spacetime is still fairly dynamical, and so we cannot determine whether what we see is near the end-state, or merely an intermediate configuration in a more complicated evolution. Underlying the current failure of our simulations is the fact that the coordinate system employed is not well adapted to the solution that unfolds at late times, wherein fatally steep gradients develop in metric functions 111This divergence of metric gradients does not happen earlier as resolution is increased, and is not accompanied by divergence of curvature invariants such as the Kretschmann scalar. This suggests that the code is evolving to a coordinate, rather than a geometric, singularity..

The outline of the remainder of the paper is as follows. In Sec. II, we begin by describing the equations of motion, our coordinate choices, generation of initial data, as well as our numerical solution scheme. Additionally, we also briefly mention the tools we employ to monitor the solution, deferring details to the Appendices. In Sec. III we discuss the results obtained with this code, and conclude in Sec. IV by mentioning directions for future work that may allow us to more definitively answer questions regarding the end-state of the Gregory-Laflamme instability.

## II Equations, Black Strings and Numerics

We wish to solve the vacuum Einstein equations in higher dimensional settings. For simplicity, and without loss of generality in studying the Gregory-Laflamme instability, we only consider the 5-dimensional case, and restrict attention to spherical symmetry within the 4-dimensional subspace tangent to the “extra” dimension. We also use the natural generalization of the ADM decomposition to derive the system of equations that we then solve numerically. Choosing units in which $`G=c=1`$, and adopting MTW MTW  conventions, our starting point is thus a metric element given by

$`ds^2=(\alpha ^2+\gamma _{AB}\beta ^A\beta ^B)dt^2+2\gamma _{AB}\beta ^Adx^Bdt+\gamma _{AB}dx^Adx^B+\gamma _\mathrm{\Omega }d\mathrm{\Omega }^2`$ (1)
where $`x^A=(r,z)`$, and $`d\mathrm{\Omega }^2`$ is the 2-spherical line element with coordinates chosen orthogonal to the $`t=`$ constant congruences (hence there is no shift corresponding to angular directions). All metric components defined via (1) depend upon $`(t,r,z)`$: $`t`$ is a time-like coordinate, $`r`$ is a radial coordinate, and $`z`$ is the coordinate along the length of the string. To further expedite the numerical implementation, we make $`z`$ a periodic coordinate by identifying $`z=0`$ and $`z=L`$. Then, following the results of Gregory and Laflamme, we expect black strings to be unstable only for $`L`$ greater than some critical length $`L_c`$.

The vacuum Einstein equations, written in ADM form MTW ; york\_78 , are 1) the Hamiltonian constraint

$$H^{(4)}R+K^2K_{ab}K^{ab}=0,$$
(2)
2) the momentum constraints

$$M_aK_{a|b}^bK_{|a}=0,$$
(3)
3) the evolution equations for the $`\gamma _{ab}`$

$$\frac{\gamma _{ab}}{t}=2\alpha K_{ab}+\beta _{a|b}+\beta _{b|a},$$
(4)
that follow from the definition of the extrinsic curvature $`K_{ab}`$ associated with $`t=\mathrm{const}.`$ slices, and 4) the evolution equations for the extrinsic curvature

$`{\displaystyle \frac{K_{ab}}{t}}`$ $`=`$ $`\alpha ({}_{}{}^{(4)}R_{ab}^{}+KK_{ab})2\alpha K_{ac}K^c{}_{b}{}^{}\alpha _{|ab}`$ (5)
$`+\beta ^c{}_{|a}{}^{}K_{cb}^{}+\beta ^c{}_{|b}{}^{}K_{ca}^{}+\beta ^cK_{ab|c}+\alpha F_a^cF_b^d\gamma _{cd}H.`$
In the above, $`a,b,\mathrm{}`$ are four-dimensional (spatial) indices, $`{}_{}{}^{(4)}R_{ab}^{}`$ and $`{}_{}{}^{(4)}R`$ are, respectively, the Ricci tensor and Ricci scalar intrinsic to the four-dimensional spatial hypersurfaces, $`\alpha `$ is the lapse function, $`\beta ^c`$ is the shift vector, the vertical bar | denotes covariant differentiation in the spatial hypersurfaces (compatible with $`\gamma _{ab}`$), and $`F_a^c=2\delta _r^c\delta _a^r`$. We note that the term proportional to $`H`$ in (5) has been added as a result of stability considerations; see for instance, the discussions in lehnerreview ; shinkai ; kelly . To simplify the final set of equations solved numerically, as well as to regularize certain terms that otherwise diverge at spatial infinity (see the discussion of our coordinate system in the next sub-section), we define the following variables

$$g_{rr}\gamma _{rr}g_{rz}\gamma _{rz}g_{zz}\gamma _{zz}$$

$$g_{\theta \theta }\gamma _\mathrm{\Omega }/r^2g_{\varphi \varphi }\gamma _\mathrm{\Omega }/(r^2\mathrm{sin}^2\theta )$$

$$k_{rr}r^2K_{rr}/\alpha k_{rz}K_{rz}k_{zz}K_{zz}$$

$$k_{\theta \theta }K_{\theta \theta }/\alpha k_{\varphi \varphi }K_{\theta \theta }/(\alpha \mathrm{sin}^2\theta )$$
(6)
and use them as the fundamental dynamical quantities in our numerical code. As discussed in the following sub-sections, to complete the prescription of the evolution problem we need to choose a suitable lapse and shift, specify initial and boundary conditions, and then implement these choices and specifications consistently.

### II.1 Boundary and Coordinate Conditions

A particular concern here is that “standard” outer boundary conditions lehnerreview , often imposed during numerical evolution of Einstein’s equations, might not be well suited for studying the string instability. In particular, we must be able to evolve for very long times while absolutely minimizing spurious influences from the outer boundary of the computational domain. In addition, in the present case we cannot assume, a priori, that any given initial configuration will settle down to some stationary solution; thus, boundary conditions predicated on such assumptions (such as a $`1/r`$ fall-off condition in a metric component), when imposed at a finite proper distance from the black string, could very well adversely affect the numerical results 222In fact, in an earlier version of the code that did not use a radially compactified coordinate system, we did encounter such problems, in that some artificial stationary non-homogeneous solution was apparently entirely “sourced” via an outer boundary located at a finite distance from the string.. To ensure minimal boundary influence we therefore extend the domain of integration to $`i^o`$ by radially compactifying the spacelike hypersurfaces via the introduction of a new coordinate, $`x`$, defined by

$$x\frac{r}{1+r}.$$
(7)
As might be expected, this transformation causes computational problems of its own—most notably decreased spatial resolution at large distances—but, as discussed in Sec. II.3, we can deal with these difficulties using numerical dissipation. Having introduced the new compactifying coordinate, we can directly impose boundary conditions derived from the demand of asymptotic flatness at spatial infinity, which lies at $`x=1`$.

We employ singularity excision techniques unruh  to allow us to evolve the entire perturbed black string spacetime exterior to the apparent horizon (plus a certain “buffer zone” that lies within the horizon). Hence, we do not need to impose inner boundary conditions as long as the $`t=\mathrm{const}.`$ hypersurfaces penetrate the horizon, and that all characteristics of the evolution equations are in-going on the boundary. Ensuring that this is the case involves choosing “good” coordinate conditions (choice of lapse and shift), which, for generic string evolutions, remains an open problem. As a preliminary step, we have based our coordinate choices on those that yield the ingoing Eddington-Finkelstein form of the unperturbed black string metric

$$ds_{\mathrm{BS}}^2=(12M/r)dt^2+4M/rdrdt+(1+2M/r)dr^2+dz^2+r^2d\mathrm{\Omega }^2.$$
(8)
Comparison with the general 5-dimensional ADM form provides the identifications

$`\alpha _{\mathrm{BS}}`$ $`=`$ $`(1+2M/r)^{(1/2)},`$ (9)
$`\beta _{\mathrm{BS}}^A`$ $`=`$ $`(2M/(r+2M))\delta _r^A.`$ (10)
For reference we also list the two non-trivial components of the extrinsic curvature of a $`t=\mathrm{const}.`$ slice defined by (8):

$`K_{rr}`$ $`=`$ $`2M{\displaystyle \frac{(r+M)}{r^3}}\sqrt{{\displaystyle \frac{r}{r+2M}}},`$
$`K_{\theta \theta }`$ $`=`$ $`2M\sqrt{{\displaystyle \frac{r}{r+2M}}}.`$ (11)
In generalizing (9) and (10) to the dynamical case, we have chosen $`\alpha =\alpha _{\mathrm{BS}}`$ and $`\beta ^z=\beta _{\mathrm{BS}}^z=0`$. In a preliminary version of our code, we also required that $`\beta ^r=\beta _{\mathrm{BS}}^r`$. This, however, caused a coordinate pathology to develop at late times during the evolution of unstable strings—specifically, some regions of the horizons approached a zero coordinate-radius, while maintaining finite proper radius. In our current efforts, we choose $`\beta ^r`$ such that $`g_{\theta \theta }`$ remains constant during evolution seidelexcision ; lehnercadm , by requiring that

$$\beta ^r=\frac{2\alpha K_{\theta \theta }}{\gamma _{\theta \theta }^{}{}_{,r}{}^{}}.$$
(12)
This shift condition performs reasonably well, as will be seen in Sec. III. However, our current simulations still suffer from “grid-stretching” problems in the $`z`$ direction at late times, suggesting that a more dynamical gauge condition for $`\beta ^z`$ (and possibly for $`\alpha `$ and $`\beta ^r`$) could be useful. This issue is discussed in more detail in Sec. IV.

### II.2 Initial Data

As anticipated, we observe that even numerical truncation errors, if non-uniform in the $`z`$ direction, are enough to trigger the Gregory-Laflamme instability in our simulations. However, to reduce the computational effort required to reach the “interesting” (i.e. non-perturbative) stages of evolution, we adopt initial configurations whose departure from the black string solution can be arbitrarily tuned. In order to find such data, we must solve the Hamiltonian constraint, and the $`r`$ and $`z`$ components of the momentum constraint (the other components of the momentum constraint are trivially satisfied because our coordinate system is adapted to spherical symmetry). The deviation—not necessarily small—from the black string solution, is introduced via $`g_{\theta \theta }`$, and takes the following form:

$$g_{\theta \theta }(0,r,z)=1+A\mathrm{sin}\left(z\frac{2\pi q}{L}\right)e^{(rr_o)^2/\delta _r^2}.$$
(13)
Here, $`A`$ is used to control the overall strength of the “perturbation”, while $`q`$ is an integer that controls the spatial frequency in the $`z`$ direction. For the results presented below, $`A=0.1`$, $`q=1`$, $`r_0=2.5`$ and $`\delta _r=0.5`$, and we perturb about a unit mass ($`M=1`$) black string solution. As described in more detail in Appendix A, $`g_{rr}`$, $`k_{rr}`$, $`k_{\theta \theta }`$ are then calculated by solving the constraint equations, with the remainder of the metric and extrinsic curvature variables set to the values they would take for an unperturbed black string (see (8) and (II.1)).

### II.3 Numerical Evolution

To numerically evolve the initial data sets described above, we discretize the evolution equations (4)-(5) using second-order accurate finite difference techniques that include Crank-Nicholson treatment of the temporal and spatial derivatives. We use a uniform distribution of grid points in $`z`$ and $`x`$ (recall that $`xr/(1+r)`$). The resulting implicit system of algebraic equations is solved iteratively. We initially implemented a serial version of the algorithm, and later coded a parallel version using the CACTUS Computational Toolkit cactus , wherein the equations of motion, monitoring tools and I/O were handled by our own routines, suitably interfaced to CACTUS.

Black hole excision is handled as follows. We periodically find the apparent horizon, as discussed in the following sub-section and Appendix B. We then define the surface along which we excise to be a certain number of “buffer” points inside the apparent horizon (typically $`10`$$`30`$ buffer points are used) 333In several tests, we also adopted an excision region given by the global minimum $`r`$ value of the apparent horizon and compared the results with those obtained when the excision region was defined by the apparent horizon. The agreement obtained gives extra indication that the excision implementation is consistent.. During each Crank-Nicholson iteration, all the evolution difference equations are applied up to the excision surface, and any function values referenced by finite difference stencils interior to this surface are defined via fourth order extrapolation. When the apparent horizon location, and hence excision surface changes during evolution, function values at all repopulated points (i.e. those that moved from inside to outside the excision surface during the time step) are computed via the same fourth order extrapolation routine. The one exception to this procedure is for the grid values of $`g_{\theta \theta }`$, which we specify a priori on the entire computational domain, and that remain fixed due to our gauge choice (12). Moreover, we have found it useful to choose a functional form for $`g_{\theta \theta }`$ that tends to zero at some positive value of $`r`$ (though inside the original apparent horizon location and outside of the limits of integration of the initial data). This causes the “pinching-off” of the unstable black string to be less severe in coordinate space, i.e. we approach zero areal radius at a finite coordinate $`r`$. In turn, this slightly reduces the virulence of the coordinate problems we observe at late times, and also provides better load-balancing of the parallel code, given the method CACTUS uses to distribute grids among processors.

For the evolution, we choose a time step $`dt=\lambda _{\mathrm{CFL}}\mathrm{min}(dr,dz)`$, where the constant $`\lambda _{\mathrm{CFL}}`$ must be set less than $`1/\sqrt{2}`$ in order to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition that results from our iterative solution of the Crank-Nicholson scheme (typically we use $`\lambda _{\mathrm{CFL}}=0.25`$). Note that this restriction on $`\lambda _{\mathrm{CFL}}`$ is based upon flat-space light speeds within our coordinate system ((8) with $`M=0`$), which, for the solutions presented here, are always greater than or equal to the actual coordinate light speeds. The function $`\mathrm{min}(dr,dz)`$ is calculated by only considering mesh spacings within the non-excised portion of the coordinate domain. Thus, as the excision surface moves, $`dt`$ changes with time since our grid is uniform in $`x`$, and hence non-uniform in $`r`$.

Crucially, we add Kreiss-Oliger-style kreissoliger  numerical dissipation to the evolution equations to control unphysical high-frequency solution components (“noise”) that would otherwise arise during the simulations. This is particularly helpful at the excision surface, and near $`i^0`$, where the radial compactification of points causes all outgoing wave-like components of variables to eventually become poorly resolved. Smoothing of the high frequency components via the Kreiss-Oliger dissipation—which only targets wavelengths of size on the order of the mesh spacing—prevents them from inducing numerical instabilities near the outer boundary.

#### II.3.1 Monitoring the Evolution

To elucidate the nature of our computed spacetimes, we monitor the following quantities: 1) the location of the apparent horizon (which is also used for excision as explained earlier), 2) the trajectories of null geodesics that, to a certain extent, should trace the event horizon, and 3) the Kretschmann invariant $`I`$ (the square of the Riemann tensor):

$`I=R_{\alpha \beta \gamma \delta }R^{\alpha \beta \gamma \delta }.`$ (14)
If cosmic censorship holds—and results from our current simulations provide no evidence to the contrary—then any apparent horizon found will always be inside an event horizon. As is well known, although the apparent horizon can often be used as a reasonable approximation to the event horizon, the two do not, in general, coincide 444Indeed, depending on the slicing an apparent horizon need not exist at all waldah .. Clearly, the event horizon is the quantity of interest in studying the Gregory-Laflamme instability, and therefore we would like to locate it, or at least a reasonably good approximation to it, in our simulation results. Such an approximation can be obtained by looking for the boundary of the causal past of some $`r=\mathrm{const}.`$ surface that is sufficiently far outside the apparent horizon that it is certain not to be inside the event horizon, yet close enough to the apparent horizon that its causal past, tracing backwards from the end-time of the simulation, probes the region of interest of the spacetime. We use a method to find the approximate event horizon discussed by Libson et al. eh . The approach is based on radial outgoing null geodesics; as explained in eh , the stable direction for the integration of null rays that emanate from the vicinity of an event horizon is backwards in time $`t`$. Thus, once we have the complete data from the entire evolution, we start with data from the latest time step available, and trace the null rays backwards in time.

Monitoring curvature scalars is useful in obtaining coordinate independent information about a numerical solution of the Einstein equations. In particular, $`I`$, as defined by (14) evaluates to $`I_{\mathrm{BS}}=48M^2/\gamma _{\theta \theta }^3`$ for the unperturbed black string solution, and as $`\gamma _{\theta \theta }`$ is an invariant in spherical symmetry (i.e. it is proportional to the area of an $`r=\mathrm{const}.`$ 2-sphere), we can compare $`I_{\mathrm{BS}}`$ to the values of $`I`$ computed from a numerical solution to get some indication of how close the computed solution is to a black string spacetime. Furthermore, we can examine $`I`$ to see whether curvature singularities (other than the central $`r=0`$ singularity) may be forming prior to the demise of the simulation that invariably occurs when sufficiently steep metric gradients develop. We note, however, that if $`I`$ does not diverge, it does not necessarily follow that the geometry is remaining non-singular; we would need to examine a larger set of curvature scalars to be certain that the solutions are remaining free of physical singularities.

Appendix B contains an explanation of the method we used to find apparent horizons, while details of the integration techniques aimed at approximately locating event horizons can be found in Appendix C.

## III Results

In this section we present results from our preliminary study of the black string instability. After briefly showing that we recover some of the key Gregory-Laflamme results in the next sub-section, we present a detailed analysis of a typical unstable case in Sec. III.2. In the following, we will use the value of $`L_c14.3M`$ (with $`M=1`$) found by Gubser gubser , which is more accurate than the value we can estimate from the zero-crossing of the (positive mode) interpolating curve presented in gregorylaflame1 .

### III.1 Recovery of Gregory-Laflamme Results

We ran a variety of simulations of black strings that were perturbed according to the prescription discussed in section II.2. We concentrated on cases with $`L`$ ranging from $`0.6L_c`$ to $`1.8L_c`$ and defined $`g_{\theta \theta }`$ via equation (13). In general, we observed the expected instability for $`L>L_c`$, though for the maximum resolution at which we performed this survey ($`800`$ grid points in $`r`$ and $`200`$ points in $`z`$), we could only confirm $`L_c`$ to within about $`2\%`$ of the expected value. In this regard we note that as $`L`$ approaches $`L_c`$ from above, the growth rate of the instability goes to $`0`$, requiring longer time integrations to identify the instability, which, in turn, demands ever increasing resolution to counter the effects of accumulating numerical errors. Furthermore, the initial configurations we have adopted contain energy in the form of gravitational waves, and some of this energy falls into the string early on during the evolution. The increase in the mass of the string (based upon the increase in area of the apparent horizon) is typically around $`0.30.5\%`$, and we would have to take this into account were we to attempt to determine $`L_c`$ from our simulations to a higher degree of accuracy.

As a demonstration of the ability of our code to “bracket” the instability, and following the notation of gubser , Fig. 1 shows a plot of, $`\lambda `$, defined by

$$\lambda =\frac{1}{2}\left(\frac{R_{\mathrm{max}}}{R_{\mathrm{min}}}1\right)$$
(15)
for $`L=0.975L_c`$ and $`L=1.03L_c`$. In the above, $`R_{\mathrm{max}}`$ and $`R_{\mathrm{min}}`$ are the maximum and minimum areal radii, respectively, of the apparent horizon at some $`t=\mathrm{const}.`$ slice of the spacetime. In particular, we have $`\lambda =0`$ for the static black string spacetime.

### III.2 Beyond the Linear Regime

We now present more detailed results from the simulation of an unstable black string evolution. Specifically, we take $`L=1.4L_c`$, since it is expected that this particular range for the $`z`$ coordinate will yield something close to the fastest growth rate for the shortest wavelength instability gregorylaflame1 . Because we are now probing uncharted territory with our computations, we rely on convergence tests to provide an intrinsic measure of the level of error in our calculations. To that end, we ran the simulation at several resolutions ($`n_r\times n_z`$): $`200\times 50`$ (grid spacing $`h`$), $`400\times 100`$ ($`h/2`$), $`800\times 200`$ ($`h/4`$), and $`1600\times 400`$ ($`h/8`$). Due to our use of a compactified radial coordinate, the lowest resolution calculation cannot adequately resolve the late-time behavior of the solution. However, for the “medium resolution” computation with mesh spacing $`h/4`$, we are apparently within the convergent regime—see Fig. 2 below for plots of the maximum and minimum areal radii of the apparent horizon as a function of time, as well as the quantity $`\lambda `$ defined by (15), and Fig. 3 for plots of the norm of the Hamiltonian constraint as a function of resolution. Therefore, unless otherwise noted, all the results shown below are taken from the $`h/4`$ simulation.

Fig. 4 shows embedding diagrams of the apparent horizon at several times during evolution of the string, and Fig. 5 shows the proper length of one period of the apparent horizon (suppressing the angular coordinates) versus time. (Our embedding uses the vertical axis to represent the areal radius of the apparent horizon—the horizontal axis is then uniquely determined by requiring that the length of the curve be equal to the proper length of the horizon). The simulation crashes shortly after the last time frame shown, apparently due to the coordinate pathologies that have been discussed previously. The embedding diagrams suggest that, at least in the vicinity of the apparent horizon, the solution is tending towards a spacetime that can be described as a sequence of spherical black holes connected by thin black strings. Additional, quantitative, evidence for this conjecture can be obtained through a computation of the curvature invariant, $`I`$, on the apparent horizon. For an exact black string solution, this quantity, which we denote $`I_{\mathrm{BS}}^0`$, is

$$I_{\mathrm{BS}}^0=\frac{12}{R_{\mathrm{AH}}^4},$$
(16)
while for the 5-dimensional spherical black hole, the equivalent quantity, $`I_{\mathrm{BH}}^0`$, is

$$I_{\mathrm{BH}}^0=\frac{72}{R_{\mathrm{AH}}^4},$$
(17)
where, in both of the above expressions, $`R_{\mathrm{AH}}`$ is the areal radius of the horizon. In Fig. 6 we plot the normalized quantity

$$I^{0N}\frac{I}{I_{\mathrm{BS}}^0}=\frac{IR_{\mathrm{AH}}^4}{12}$$
(18)
evaluated on the apparent horizon of our numerical solution of the unstable spacetime—$`I^{0N}`$ is $`1`$ for a black string, and $`6`$ for a black hole. The figure shows that, as judged by $`I^{0N}`$, the part of the apparent horizon that is forming a long neck always resembles a black string—the part that is forming a bulge, however, has a value of $`I^{0N}`$ tending towards that corresponding to a black hole. At $`R_{\mathrm{max}}`$, $`I^{0N}`$ has only reached $`5`$ by the time the simulation ends; however, the behavior of $`R_{\mathrm{max}}`$ seen in Fig. 2 suggests that the growth in the normalized curvature invariant, though slowing down, should continue. Fig. 6 also demonstrates the grid-stretching problems that we surmise are causing the code to eventually crash—in that plot we use the coordinate $`z`$ as the horizontal axis, and observe that the relatively small region where $`I^{0N}1`$ corresponds to the long neck in Fig. 4. In particular, in the vicinity of the “neck”, $`g_{zz}`$ becomes quite large, as do its derivatives.

Finally, in Fig. 7 we show plots of the approximate event horizon (as described in Sec. II.3.1), together with the apparent horizon for the simulation. The results shown in the plot suggest that our computed apparent horizon is an excellent approximation to the event horizon, at least at early times (not much can be said regarding the late time behavior of the event horizon, as the spacetime has not settled down to a stationary state when the simulation ends).

## IV Conclusions

We have performed a preliminary numerical study of the instability of the 5-dimensional black string. Coordinate pathologies prevent us from definitively identifying the final end-state(s) of an unstable black string. This claim is supported by the fact that the code crashes at very nearly the same time at varying resolution, and that curvature invariants remain well behaved throughout the evolution. The former suggests that a numerical instability is not responsible for the crash, while the latter indicates that a physical singularity is probably also not to blame. Despite the premature termination of the simulation, we find evidence that the spacetime evolves towards a configuration that looks like a sequence of black holes connected by thin black strings, and characterized by an expansion of the string direction. Since the spacetime is still fairly dynamical at the time our simulations end, we cannot deduce how close this state is to a final configuration. Nevertheless, the dynamical behavior observed is sufficiently robust for some comments to be made. For instance, the results are not inconsistent with Gregory and Laflamme’s conjecture that the solution bifurcates into a sequence of black holes—indeed, we can suggest at least two mechanisms by which this could occur:

1. Via a thinning neck that eventually vanishes, if the trend seen in the simulation continues. Note that this would require a) that the proper length of the string continue to grow, in order not to violate area theorems, and b) that the thin string be non-uniform, for otherwise it would be subject to a further Gregory-Laflamme-like instability.
2. Via a sequence of Gregory-Laflamme instabilities, if the thin neck stays “close” to a uniform black string, since the neck’s length is beyond the critical one for a string with an effective mass computed from the radius of the apparent horizon. In this case, one could envision a “cascade” of instabilities leading to the bifurcation.

We note that either scenario would not necessarily be inconsistent with Horowitz and Maeda’s results, should the vanishing of the neck take infinite affine time as measured by local null generators of the horizon. At the same time, a continuation of the observed trend would argue against achieving a stationary solution with a mild dependence on the string dimension (i.e. a small value of $`\lambda `$), as found in perturbative calculations gubser . For then, the rather extreme thinning/bulging that we see must be transient behavior that is “further” from the end-state than the perturbed black string was.

A more complete exposition of the nature of unstable black string evolution would appear to require coordinate conditions able to adapt to solution features that develop at late times—that is, in a manner that does not introduce severe metric gradients that are not correlated with large gradients in physical quantities. For example, it may help to replace the fixed-lapse slicing with maximal slicing, which enforces that the divergence of the local, spatial volume element be zero. Another, perhaps even more crucial option, would be to introduce a $`z`$ component to the shift vector that keeps $`g_{zz}`$ close to (or exactly) unity throughout the evolution. These options are currently under investigation.

Additionally, it would be interesting to explore a wider range of initial conditions describing “perturbed” black strings than that considered here. For example, the imposed $`z`$-periodicity implies that the equivalent, uncompactified space time consists of identical spherical-black-hole/black-string segments at late (intermediate) times. It would be instructive to see what happens should we break this symmetry, by making $`LL_c`$, and then introducing some higher-wavelength perturbation similar to $`q>1`$ in (13), but with more asymmetry in the initial data (note that this would be more computationally demanding). Finally, it would be very interesting to study the evolution of the solutions recently found by Wiseman wiseman , and mentioned in the introduction. These configurations actually correspond to stationary solutions, and their perturbative stability, or otherwise, is currently not known. Since Wiseman shows that his solutions cannot be the end-states conjectured by Horowitz and Maeda, it is important to understand their behavior, since if they are stable they may well represent physically meaningful states, while if unstable, they may be difficult to attain via dynamical evolution. An interesting observation from our simulations, to the length so far achieved, is that they do not display a conical “waist” like those presented in wiseman  and further analyzed in wisemankol .

## V Acknowledgments

We gratefully acknowledge support from the following agencies, institutes and grants: NSERC, NSF PHY-0099568, The Canadian Institute for Advanced Research, The Canadian Institute for Theoretical Astrophysics, The Pacific Institute for Mathematical Sciences, The Government of the Basque Country, The Izaak Walton Killam Fund and Caltech’s Richard Chase Tolman Fund. Computations were performed on (i) the vn.physics.ubc.ca cluster which was funded by the Canadian Foundation for Innovation (CFI) and the BC Knowledge Development Fund; (ii) LosLobos at Albuquerque High Performance Computing Center (iii) The high-performance computing facilities within LSU’s Center for Applied Information Technology and Learning, which is funded through Louisiana legislative appropriations, and (iv) The MACI cluster at the University of Calgary, which is funded by the Universities of Alberta, Calgary, Lethbridge and Manitoba, and by C3.ca, the Netera Alliance, CANARIE, the Alberta Science and Research Authority, and the CFI. We would like to thank G. Horowitz, W.G. Unruh, R. Wald, R. Myers, T. Wiseman and B. Kol for stimulating discussions.

## Appendix A Initial data solver

We solve the set of coupled constraint equations (2)-(3) via an iterative procedure, where at each sub-step of the iteration we solve a single equation for one of $`g_{rr}`$, $`k_{\theta \theta }`$ or $`k_{rr}`$, assuming that the values of the other variables are known. We iterate this process until the residuals of all the equations are simultaneously below a certain tolerance—a typical value is $`10^5`$. The overall iteration is initialized using values corresponding to an unperturbed black string solution. We now provide a few more details concerning the solution of each of the constraint equations.

The equations are discretized on a uniform grid of points ($`x_i,z_j`$) with $`i=1,\mathrm{},N_x`$, and $`j=1,\mathrm{},N_z`$ (recall that, from (7), the radial coordinate, $`r`$, is related to $`x`$ by $`r=x/(1x)`$). We first consider the Hamiltonian constraint (2) which, in the coordinate system we have adopted, has the following form:

$$F_1\frac{g_{rr}}{x}+F_2g_{rr}\frac{^2g_{rr}}{z^2}+F_3g_{rr}\frac{g_{rr}}{z}+F_4\left(\frac{g_{rr}}{z}\right)^2+F_5\left(g_{rr}\right)^2+F_6g_{rr}=0.$$
(19)
Here, the $`F_m,m=1,\mathrm{},6`$, are functions that generally depend on all the metric coefficient and their derivatives except $`g_{rr}`$ (and its derivatives). We discretize this equation to second order in the mesh spacing using a difference approximation centered at the points $`(x_{i+1/2},z_j)`$. Because of the discretization used and the form of (19), we can solve the resulting set of equations “line-by-line” in $`x`$, starting at the inner boundary, $`i=1`$, which is chosen well within the horizon of the string (typically at $`r=M`$), and where the boundary values, $`[g_{rr}]_{1,j},j=1,\mathrm{}N_z`$, are those corresponding to an unperturbed black string. As we integrate outwards in $`x`$, the determination of the $`i`$-th line of unknowns, $`[g_{rr}]_{i,j},j=1,\mathrm{}N_z`$, involves the solution of a $`N_z`$-dimensional, non-linear, cyclic (because of the $`z`$-periodicity), tridiagonal system that we solve using Newton’s method and a cyclic tridiagonal linear solver nr .

We now direct attention to the $`r`$ momentum constraint, which, viewed as an equation for $`k_{\theta \theta }`$, has the form:

$$G_1\frac{k_{\theta \theta }}{x}+G_2k_{\theta \theta }+G_3=0,$$
(20)
where the $`G_m,m=1,2,3`$ do not depend on $`k_{\theta \theta }`$ or its derivatives. We note that there is complete decoupling in the $`z`$-direction in this case; in effect, we have to solve an ODE along each $`z=\mathrm{const}.`$ line. We again discretize using second-order finite difference techniques, fix the boundary values $`[k_{\theta \theta }]_{N_x,j},j=1,\mathrm{}N_z`$, at $`x=1`$ ($`i_o`$) using the unperturbed black string solution, then solve for the remaining unknowns, marching inwards in $`x`$.

Finally, the $`z`$ component of the momentum constraint, which fixes $`k_{rr}`$, has the form:

$$H_1\frac{k_{rr}}{z}+H_2k_{rr}+H_3=0;$$
(21)
where the $`H_m,m=1,2,3`$ are independent of $`k_{rr}`$ and its derivatives. This equation is solved analogously to the $`r`$ momentum constraint, but now using a discretization that is centered at points ($`x_i,z_{j+1/2}`$). “Boundary values”, $`[k_{rr}]_{i,1},i=1,\mathrm{}N_x`$, are specified along the line $`z=z_{\mathrm{min}}`$, again using corresponding values from the black string solution, and the integration proceeds for $`j=2,3,\mathrm{}N_z`$.

## Appendix B Finding apparent horizons

We use a flow, or level-set method to search for apparent horizons within $`t=\mathrm{const}.`$ spatial slices of the spacetime. We restrict our search to simply connected apparent horizons that are periodic in z. Such an apparent horizon can be described by a curve in the $`(r,z)`$ plane, which we define to be the level surface $`F=0`$ of the function

$$F(r,z)=rR(z).$$
(22)
In other words, the apparent horizon will be given by the curve $`r=R(z)`$. The apparent horizon is the outermost, marginally trapped surface; hence, we want to find an equation for $`R(z)`$ such that the outward null expansion, normal to the corresponding surface $`F=0`$, is zero. To this end, we first construct the unit spatial vector $`s^a`$, normal to $`F=\mathrm{const}.`$:

$$s^a=\frac{g^{ab}F_{,b}}{\sqrt{g^{cd}F_{,c}F_{,d}}}.$$
(23)
Then, using $`s^a`$ and the $`t=\mathrm{const}.`$ hypersurface normal vector $`n^a`$, we can construct future-pointing outgoing($`+`$) and ingoing($``$) null vectors

$$\mathrm{}^a\pm =n^a\pm s^a.$$
(24)
The normalization of the null vectors is (arbitrarily) $`\mathrm{}_+^a\mathrm{}_a=2`$. The outward null expansion $`\theta _+`$ is then the divergence of $`\mathrm{}_+^a`$ projected onto an $`F=\mathrm{const}.`$ surface:

$$\theta _+=\left(g^{ab}s^as^b\right)_b\mathrm{}_{+a}.$$
(25)
Substituting expressions (22) and (23) into equation (25), with $`\theta _+=0`$, provides us with an ordinary differential equation for $`R(z)`$.

Initially we solved equation (25) via a “shooting” method—given a guess for $`R(0)`$, and assuming that $`R^{}(0)=0`$ (where a prime denotes differentiation with respect to $`z`$), we integrate the equation to $`z=z_L`$, and repeat the process until we find a solution where $`R(z_L)=R(z)`$ and $`R^{}(z_L)=0`$. An efficient sequence of guesses can be generated using a bisection search, as the qualitative behavior of the solution is different depending upon whether the initial guess for $`R(0)`$ is inside or outside the apparent horizon.

Since the shooting method is difficult to extend to a parallel implementation in an efficient way, we opted to use the following point-wise relaxation method (or flow method) to determine $`R(z)`$. We supply an initial guess, $`R_0(z)`$, for the entire function $`R(z)`$, and then iterate the following equation until the norm of the expansion $`\theta _+(z)`$ of $`F(0)`$ is below some desired threshold (in our runs we have typically set it to $`10^3h`$, where $`h`$ is the basic scale of discretization):

$$\mathrm{\Delta }R_n(z)R_{n+1}(z)R_n(z)=\theta _+(z)\mathrm{\Delta }\tau .$$
(26)
Here, $`R_n(z)`$ is the solution after the $`n^{th}`$ iteration. Equation (26) can be viewed as an explicit discretization of a parabolic evolution equation for $`R(z,\tau )`$, where $`\tau `$ is “time” and $`\mathrm{\Delta }\tau `$ is the time-step for the evolution (the parabolic nature of the equation is evident when $`\theta _+`$ is expanded in terms of $`R(z)`$ via (25)—for brevity we do not give the explicit form of the equation here). Thus, for stability $`\mathrm{\Delta }\tau `$ must be chosen to be less than $`\mathrm{\Delta }z^2`$.

Given a “reasonable” initial guess $`R_0(z)`$, one can see how iteration of equation (26) will cause $`R_n(z)`$ to “flow” to the apparent horizon: if $`R_n(z)`$ is outside of the apparent horizon, then typically the expansion $`\theta _+(z)`$ will be positive there, causing $`R_{n+1}(z)`$ to decrease towards the apparent horizon, and vice-versa if $`R_n(z)`$ is inside the apparent horizon. We use $`R(z)=2`$ as the initial guess at $`t=0`$; after $`t=0`$ we search for the apparent horizon every $`N`$ time steps (where $`N`$ is typically in the range $`1030`$), and use the shape found at the previous search as the initial guess for the next search. Usually, on the order of tens to thousands of iterations of (26) are required to solve for $`R(z)`$ to within a level of accuracy such that the approximate solution is roughly within a mesh spacing of the exact solution (as estimated in a few specific calculations by solving (26) close to machine precision). A single iteration of (26) can be computed very rapidly relative to the time taken to compute a metric-evolution step; however in a parallel environment, if thousands of iterations are needed on a regular basis (which is so at late times during the evolution of an unstable black string), the apparent horizon finder becomes a slight speed bottleneck in the code, due to the time it takes to communicate the results of each iteration amongst the processors involved.

## Appendix C Finding (Approximate) Event Horizons

Here we describe one method we use, following eh , to locate approximations to event horizons. This method involves locating the boundary of the causal past of some $`r=\mathrm{const}.`$ surface of the spacetime by following radial null geodesics.

We write the geodesic equation in Lagrangian form:

$$=g_{tt}\left(t^{}\right)^2+2g_{tr}t^{}r^{}+g_{rr}\left(r^{}\right)^2+2g_{rz}r^{}z^{}+g_{zz}\left(z^{}\right)^2,$$
(27)
where $`\lambda `$ is the affine parameter along the geodesic, and a prime denotes differentiation with respect to $`\lambda `$. Since we are interested in null trajectories, we set $`=0`$. For radial, ingoing geodesics, we have $`\theta ^{}=\varphi ^{}=0`$, and thus the geodesic equations reduce to the set:

$`\dot{r}`$ $`=`$ $`{\displaystyle \frac{\alpha }{\sqrt{g_{rr}}}}\beta ,`$ (28)
$`\dot{\lambda }`$ $`=`$ $`{\displaystyle \frac{2\alpha \sqrt{g_{rr}}}{\mathrm{\Pi }_r}},`$
where the dot denotes a derivative with respect to coordinate time, and $`\mathrm{\Pi }_r=/r^{}`$. Then, starting at a certain value of $`r=r_0`$, and for each grid point along the $`z`$ direction, equations (28) are integrated backwards in time using a second order Runge-Kutta scheme.

Following null geodesics along $`z=\mathrm{const}.`$ lines does not guarantee that we are tracing the causal past of $`r=r_0`$, though for the spacetimes studied here, and the coordinate system used, this should offer a good approximation. Furthermore, although radial geodesics might not be the best choice, since the event horizon is an attractor, they will trace it accurately. (For related discussions of approximate event horizon location in the axisymmetric four-dimensional case, see shapiroteukolskyEH ; seidelEVENTHORIZONS ; caveny .)


text:
# Physical wavelets: Lorentz covariant, singularity-free, finite energy, zero action, localized solutions to the wave equation

## 1 Motivation

While the particle physics community has for some time made extensive use of extended field configurations such as solitons, instantons, and sphalerons, no direct use has yet been made of the quite extensive literature on “localized wave” configurations developed by the engineering, optics, and mathematics communities. (For selected references see .) These localized waves are classical solutions of the wave equation that are partially localized in space or time, this localization generally coming at a cost such as infinite total energy and/or instability (leading to dispersion or diffraction). The catalogue of known localized waves is large and growing , but most of the known examples are not in a form that would be easy to apply to particle physics problems.

In this article I will exhibit a particularly simple “physical wavelet” that is more promising from a particle physics standpoint. It satisfies the properties that:

* It is a localized wave that solves the wave equation.
* It is a Lorentz covariant classical field configuration that lives in physical Minkowski space.
* The field is everywhere finite and nonsingular, and has quadratic falloff in both space and time.
* The total energy is finite, depending on the peak field and the width of the pulse.
* The total action is zero.

These physical wavelets can be constructed for both complex and real scalar fields. Extending these idea to the Maxwell and Yang-Mills fields is straightforward. The simplest case is that of the complex scalar field and it is to that case that I first turn.

## 2 Complex scalar field

Let $`\eta _{ab}=\text{diag}[+1,1,1,1]`$ be the Minkowski metric \[particle physics signature\], let $`x_0`$ be an arbitrary 4-vector, and let $`\zeta ^a`$ be arbitrary timelike 4-vector, then

$$\varphi (x)=\frac{\varphi _0(\eta _{ab}\zeta ^a\zeta ^b)}{\eta _{ab}[x^ax_0^ai\zeta ^a][x^bx_0^bi\zeta ^b]}$$
(1)
is a Lorentz covariant, finite energy, zero-action solution of the d’Alembertian wave equation $`\mathrm{\Delta }\varphi =0`$. The “center” of the pulse is at $`x_0`$ and its “width” is $`a=\sqrt{\eta _{ab}\zeta ^a\zeta ^b}`$. The field is everywhere finite and in fact

$$|\varphi (x)||\varphi _0|.$$
(2)
To see this, use the fact that $`\zeta `$ is timelike. Then, using the manifest Lorentz covariance of the field configuration, we can without loss of generality first translate $`x_00`$, and then go into the zero-momentum frame where

$$\zeta ^a=(a,0,0,0).$$
(3)
Then the field configuration is

$$\varphi (x)=\frac{\varphi _0a^2}{[tia]^2x^2y^2z^2}.$$
(4)
That is

$$\varphi (x)=\frac{\varphi _0a^2}{r^2t^2+a^2+2iat}.$$
(5)
Once written in this form it is a straightforward exercise to verify that the wave equation is satisfied. To see that the field is everywhere bounded note

$`|\varphi |^2`$ $`=`$ $`{\displaystyle \frac{|\varphi _0|^2a^4}{(r^2t^2+a^2)^2+4a^2t^2}}={\displaystyle \frac{|\varphi _0|^2a^4}{(r^2+t^2+a^2)^24r^2t^2}}`$ (6)
$``$ $`{\displaystyle \frac{|\varphi _0|^2a^4}{(r^2+t^2+a^2)^2(r^2+t^2)^2}}={\displaystyle \frac{|\varphi _0|^2a^4}{a^4+2a^2(r^2+t^2)}}|\varphi _0|^2.`$
From the penultimate inequality we also derive

$$|\varphi |^2\frac{1}{2}|\varphi _0|^2\frac{a^2}{r^2+t^2},$$
(7)
demonstrating the promised quadratic falloff in both space and time. Indeed for fixed $`t`$ the magnitude of the field is maximized when

$$r^2=\mathrm{max}\{t^2a^2,0\},$$
(8)
showing that the configuration disperses to spatial infinity at both $`t\pm \mathrm{}`$.

To calculate the 4-momentum, we remind the reader that the stress-energy tensor for a massless complex scalar is

$$T_{ab}=\frac{1}{2}\left[\varphi _a^{}\varphi _b+\varphi _a\varphi _b^{}\right]\frac{1}{2}\eta _{ab}|\varphi |^2.$$
(9)
Then

$`_aT^{ab}`$ $``$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }\varphi ^{}_b\varphi +{\displaystyle \frac{1}{2}}^a\varphi ^{}_a_b\varphi +{\displaystyle \frac{1}{2}}\mathrm{\Delta }\varphi _b\varphi ^{}+{\displaystyle \frac{1}{2}}^a\varphi _a_b\varphi ^{}`$ (10)
$`{\displaystyle \frac{1}{2}}^b\varphi ^{}_a_b\varphi ^{}{\displaystyle \frac{1}{2}}^a\varphi ^{}_a_b\varphi ,`$
$``$ $`{\displaystyle \frac{1}{2}}\left[\mathrm{\Delta }\varphi ^{}_b\varphi +\mathrm{\Delta }\varphi _b\varphi ^{}\right],`$ (11)
which vanishes by the equations of motion. But this means that

$$P^\mu =T^{ab}d\mathrm{\Sigma }_b$$
(12)
is a conserved quantity, the 4-momentum of the configuration, which is independent of the particular spacelike hypersurface $`\mathrm{\Sigma }`$ chosen to do the integration. By simple dimensional analysis

$$P^a=C|\varphi _0|^2\zeta ^a,$$
(13)
where $`C`$ is a dimensionless number to be calculated. \[Note that $`\zeta ^a`$ has the dimensions of a position vector — a distance.\] The energy density is

$$\rho =\frac{1}{2}\left[|_t\varphi |^2+|_r\varphi |^2\right],$$
(14)
and in the zero-momentum frame is easily calculated to be

$$\rho =\frac{2a^4|\varphi _0|^2(r^2+t^2+a^2)}{(r^2t^2+a^2+2iat)^2(r^2t^2+a^22iat)^2.}$$
(15)
For arbitrary $`t`$ this integrates to

$$=\mathrm{d}^3r\rho =_0^{\mathrm{}}4\pi r^2\rho =\frac{1}{2}\pi ^2|\varphi _0|^2a.$$
(16)
(This is independent of $`t`$ as it should be.) This is the invariant mass of the field configuration. By spherical symmetry, the total momentum is zero. Thus, for any timelike $`\zeta ^a`$

$$P^a=\frac{1}{2}\pi ^2|\varphi _0|^2\zeta ^a.$$
(17)
Furthermore the Lagrangian is

$$=\frac{1}{2}\left[|_t\varphi |^2|_r\varphi |^2\right],$$
(18)
which evaluates (in the zero-momentum frame) to

$$=\frac{2a^4|\varphi _0|^2(t^2+a^2r^2)}{(r^2t^2+a^2+2iat)^2(r^2t^2+a^22iat)^2}.$$
(19)
It is easy to check that

$$\mathrm{d}^4x=0,$$
(20)
so that the configuration is zero action.

In summary, what we have is a Lorentz covariant, singularity-free, finite energy, zero action, exact localized solution to the d’Alembertian equation. In many ways this configuration has more right to be called an “instanton” than do the instantons of QFT; those instantons live in Euclidean signature. This field configuration lives in real physical time.

Now the fact that there are finite energy solutions to the wave equation is not a surprise; that these finite energy solutions can coalesce, bounce, and disperse without producing field singularities is more interesting. One way of guessing that the field configuration in equation (1) is worth investigating is the following: It is easy to convince oneself that in 4 Euclidean dimensions the solution to Laplace’s equation with a delta function source at the origin is

$$\varphi (x)\frac{1}{x^2+y^2+z^2+t^2}.$$
(21)
Thus in (3+1) Lorentzian dimensions the \[singular\] solution to the wave equation with a delta function source at the origin is

$$\varphi (x)\frac{1}{x^2+y^2+z^2t^2}.$$
(22)
If the source is now moved to a real position $`x_0^a`$ we have

$$\varphi (x)\frac{1}{(xx_0)^2+(yy_0)^2+(zz_0)^2(tt_0)^2}.$$
(23)
which is still a singular field configuration. Finally, move the source away from physical Minkowski space to the complex position $`x_0^ai\zeta ^a`$, then

$$\varphi (x)\frac{1}{(xx_0+i\zeta ^1)^2+(yy_0+i\zeta ^2)^2+(zz_0+i\zeta ^3)^2(tt_0+i\zeta ^0)^2},$$
(24)
which is essentially equation (1) above. This style of approach has been particularly advocated by Kaiser . As we have just seen, if $`\zeta `$ is timelike the resulting field configuration is singularity free. However, for null and spacelike $`\zeta ^a`$, while the field is still a solution of the wave equation, the field is not bounded. Because of the singularities the energy and action integrals then diverge. Details are deferred for now and will be presented in sections 6 and 7 below.

One should also note that in the optics and engineering literature the most commonly used notations are not manifestly Lorentz covariant. Thus it is common to see expressions such as

$$\varphi \frac{1}{x^2+y^2+[b_1i(z+t)][b_2+i(zt)]}$$
(25)
(see, for instance,  ) whose Lorentz transformation properties are less than obvious — in fact this field configuration is equivalent to equation (1) with the identification

$$\zeta ^a=(\frac{b_1+b_2}{2},0,0,\frac{b_1b_2}{2});\zeta =b_1b_2;$$
(26)
## 3 Real scalar field

By taking real and imaginary parts of the complex solution above we can write down two solutions for the real scalar field. Namely

$$\varphi _1=\frac{\varphi _0(\eta _{ab}\zeta ^a\zeta ^b)\left\{\eta _{ab}[x^ax_0^a][x^bx_0^b]\eta _{ab}\zeta ^a\zeta ^b\right\}}{(\eta _{ab}[x^ax_0^a][x^bx_0^b]\eta _{ab}\zeta ^a\zeta ^b)^2+4(\eta _{ab}[x^ax_0^a]\zeta ^b)^2};$$
(27)
$$\varphi _2=\frac{2\varphi _0(\eta _{ab}\zeta ^a\zeta ^b)\left\{\eta _{ab}[x^ax_0^a]\zeta ^b\right\}}{(\eta _{ab}[x^ax_0^a][x^bx_0^b]\eta _{ab}\zeta ^a\zeta ^b)^2+4(\eta _{ab}[x^ax_0^a]\zeta ^b)^2}.$$
(28)
As previously, we can without loss of generality translate $`x_00`$ and go to the zero-momentum frame $`\zeta ^a=(a,0,0,0)`$, then

$$\varphi _1=\frac{\varphi _0a^2\left\{t^2r^2a^2\right\}}{(t^2r^2a^2)^2+4a^2t^2};$$
(29)
$$\varphi _2=\frac{\varphi _0a^2\mathrm{\hspace{0.33em}2}at}{(t^2r^2a^2)^2+4a^2t^2}.$$
(30)
The stress-energy for a real scalar field simplifies

$$T_{ab}=\varphi _a\varphi _b\frac{1}{2}\eta _{ab}|\varphi |^2.$$
(31)
Then

$`_aT^{ab}`$ $``$ $`\mathrm{\Delta }\varphi _b\varphi +^a\varphi _a_b\varphi ^b\varphi _a_b\varphi `$ (32)
$``$ $`\mathrm{\Delta }\varphi _b\varphi ,`$ (33)
which vanishes by the equations of motion. The calculation for the energy-momentum 4-vector now yields:

$$P_1^a=\frac{1}{4}\pi ^2|\varphi _0|^2\zeta ^a=P_2^a.$$
(34)
The action integral for both of these field configurations is still zero. \[$`\mathrm{d}^3xdt(\varphi _{1,2})=0`$.\]

## 4 Maxwell field

A similar construction can be performed for the Maxwell field. There are a number of choices one could make, and I will simply pick one that leads to a relatively simple field configuration. Start by picking a timelike 4-velocity $`V^a`$ and a spacelike unit vector $`m^a`$ orthogonal to it. Now adopt the ansatz

$$A^a=\left\{V^am^bm^aV^b\right\}_b\psi .$$
(35)
Here one is automatically in Lorenz gauge , $`_aA^a=0`$, and the Maxwell equations reduce to the wave equation $`\mathrm{\Delta }\psi =0`$ for the scalar potential $`\psi `$. The field configuration is manifestly Lorentz covariant and we can without loss of generality go to the inertial frame where

$$V^a=(1,0,0,0);m^a=(0;\stackrel{}{m}).$$
(36)
In this inertial frame

$$A_a=(\phi ,\stackrel{}{A})=([\stackrel{}{m}\stackrel{}{}]\psi ,\stackrel{}{m}\dot{\psi }),$$
(37)
where $`\stackrel{}{m}`$ is a constant unit vector. We shall soon see that this is tantamount to working in the zero-momentum frame of the field configuration. The electric field is

$$\stackrel{}{E}=\stackrel{}{}\left([\stackrel{}{m}\stackrel{}{}]\psi \right)+\stackrel{}{m}\ddot{\psi }=\stackrel{}{}\left([\stackrel{}{m}\stackrel{}{}]\psi \right)+\stackrel{}{m}^2\psi =\stackrel{}{}\times (\stackrel{}{}\times [\stackrel{}{m}\psi ]),$$
(38)
where we have used the wave equation for $`\psi `$. The magnetic field is

$$\stackrel{}{B}=\stackrel{}{}\times (\stackrel{}{m}\dot{\psi })=\stackrel{}{m}\times \stackrel{}{}\dot{\psi }.$$
(39)
The energy density and momentum flux (Poynting vector) are

$$\rho =\frac{1}{2}[\stackrel{}{E}^2+\stackrel{}{B}^2];\stackrel{}{S}=\stackrel{}{E}\times \stackrel{}{B}.$$
(40)
To evaluate total energy and momentum a useful integration by parts (subject to suitable falloff at spatial infinity) is

$`{\displaystyle \mathrm{d}^3x(\stackrel{}{}\times \stackrel{}{X}_1)(\stackrel{}{}\times \stackrel{}{X}_2)}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x(\stackrel{}{X}_1\times \stackrel{}{})(\stackrel{}{}\times \stackrel{}{X}_2)}`$ (41)
$`=`$ $`{\displaystyle \mathrm{d}^3x\stackrel{}{X}_1\left[\stackrel{}{}\times (\stackrel{}{}\times \stackrel{}{X}_2)\right]}`$ (42)
$`=`$ $`{\displaystyle \mathrm{d}^3x\stackrel{}{X}_1\left[\stackrel{}{}^2\stackrel{}{X}_2+\stackrel{}{}(\stackrel{}{}X_2)\right]}`$ (43)
$`=`$ $`{\displaystyle \mathrm{d}^3x\left\{\stackrel{}{X}_1\stackrel{}{}^2\stackrel{}{X}_2+(\stackrel{}{}\stackrel{}{X}_1)(\stackrel{}{}\stackrel{}{X}_2)\right\}}.`$ (44)
Consequently

$`{\displaystyle \mathrm{d}^3x\stackrel{}{E}^2}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x\left\{[\stackrel{}{}\times (\stackrel{}{m}\psi )]^2[\stackrel{}{}\times (\stackrel{}{m}\psi )]\right\}}`$ (45)
$`=`$ $`+{\displaystyle \mathrm{d}^3x\left\{(\stackrel{}{m}\psi )^2^2(\stackrel{}{m}\psi )+(\stackrel{}{}[\stackrel{}{m}\psi ])^2(\stackrel{}{}[\stackrel{}{m}\psi ])\right\}}`$ (46)
$`=`$ $`+{\displaystyle \mathrm{d}^3x\left\{\psi ^2^2\psi +(\stackrel{}{m}\stackrel{}{}\psi )^2(\stackrel{}{m}\stackrel{}{}\psi )\right\}}.`$ (47)
Now let us assume the potential $`\psi `$ is spherically symmetric $`\psi (r,t)`$. Then averaging over angular variables is the same as averaging over orientations of the unit vector $`\stackrel{}{m}`$ and under the angular integral we can effectively replace

$$m_im_j\frac{1}{3}\delta _{ij}.$$
(48)
Consequently

$`{\displaystyle \mathrm{d}^3x\stackrel{}{E}^2}`$ $`=`$ $`+{\displaystyle \mathrm{d}^3x\left\{^2\psi ^2\psi +\frac{1}{3}(\stackrel{}{}\psi )^2(\stackrel{}{}\psi )\right\}}`$ (49)
$`=`$ $`+{\displaystyle \mathrm{d}^3x\left\{(^2\psi )^2\frac{1}{3}(^2\psi )^2\right\}}`$ (50)
$`=`$ $`+{\displaystyle \frac{2}{3}}{\displaystyle \mathrm{d}^3x\left\{(^2\psi )^2\right\}}.`$ (51)
Similarly

$`{\displaystyle \mathrm{d}^3x\stackrel{}{B}^2}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x\left\{(\stackrel{}{m}\dot{\psi })^2(\stackrel{}{m}\dot{\psi })+(\stackrel{}{m}\stackrel{}{}\dot{\psi })(\stackrel{}{m}\stackrel{}{}\dot{\psi })\right\}}.`$ (52)
Again invoking spherical symmetry for $`\psi `$ we have

$`{\displaystyle \mathrm{d}^3x\stackrel{}{B}^2}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x\left\{\dot{\psi }^2\dot{\psi }+\frac{1}{3}(\stackrel{}{}\dot{\psi })(\stackrel{}{}\dot{\psi })\right\}}`$ (53)
$`=`$ $`+{\displaystyle \frac{2}{3}}{\displaystyle \mathrm{d}^3x\left\{(\stackrel{}{}\dot{\psi })(\stackrel{}{}\dot{\psi })\right\}}.`$ (54)
Therefore

$$=\mathrm{d}^3x\rho =\frac{1}{3}\mathrm{d}^3x\left\{(^2\psi )^2+(\stackrel{}{}\dot{\psi })(\stackrel{}{}\dot{\psi })\right\},$$
(55)
so the total energy is sensibly positive. Before making further choices regarding the potential $`\psi `$, consider the total momentum

$$\stackrel{}{\mathrm{}}=\mathrm{d}^3x\stackrel{}{S}=\mathrm{d}^3x\left[\stackrel{}{}\times (\stackrel{}{}\times [\stackrel{}{m}\psi ])\right]\times \left[\stackrel{}{}\times (\stackrel{}{m}\dot{\psi })\right].$$
(56)
Integration by parts implies

$$\stackrel{}{\mathrm{}}=\mathrm{d}^3x\left\{(\stackrel{}{m}\dot{\psi })^2(\stackrel{}{}\times [\stackrel{}{m}\psi ])\right\}=\mathrm{d}^3x\left\{\dot{\psi }\stackrel{}{m}[\stackrel{}{m}\times \stackrel{}{}(^2\psi )]\right\}=0,$$
(57)
so that the net momentum is zero. For the action, an integration by parts together with the equations of motion yields

$$\mathrm{d}^4x\frac{1}{2}\left[\stackrel{}{E}^2\stackrel{}{B}^2\right]=\frac{1}{3}\mathrm{d}^3x\left\{(^2\psi )^2(\stackrel{}{}\dot{\psi })(\stackrel{}{}\dot{\psi })\right\}=0.$$
(58)
So we still have a zero action solution.

Up to now the potential $`\psi `$ has only needed to be spherically symmetric and to satisfy the wave equation (plus some falloff constraint at spatial infinity to allow the integration by parts). The general solution to the wave equation in spherical symmetry is

$$\psi (r,t)=\frac{1}{r}\left[f(r+t)f(rt)\right].$$
(59)
Returning to the energy $``$ as given in equation (55), a further integration by parts, together with the wave equation for $`\psi `$ yields the computationally convenient form

$$=\frac{1}{3}\mathrm{d}^3x\left\{_t^2\psi _t^2\psi _t\psi _t^3\psi \right\}.$$
(60)
Whence

$``$ $`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}\mathrm{d}r\{(_t^2[f(r+t)f(rt)])^2`$
$`(_t[f(r+t)f(rt)])(_t^3[f(r+t)f(rt)])\}`$
$`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}\mathrm{d}r\{(_r^2[f(r+t)f(rt)])^2`$
$`(_r[f(r+t)+f(rt)])(_r^3[f(r+t)+f(rt)])\}`$
$`=`$ $`4\pi {\displaystyle _0^{\mathrm{}}}\mathrm{d}r\{(_r^2[f(r+t)f(rt)])^2`$
$`+(_r^2[f(r+t)+f(rt)])^2\}`$
$`=`$ $`8\pi {\displaystyle _0^{\mathrm{}}}dr\left\{[_r^2f(r+t)]^2+[_r^2f(rt)]^2\right\}`$ (64)
$`=`$ $`8\pi {\displaystyle _{\mathrm{}}^+\mathrm{}}ds[_s^2f(s)]^2.`$ (65)
Which verifies that the energy is constant in a model-independent manner. Furthermore if $`f(s)`$ is smooth and satisfies suitable falloff conditions at $`s\pm \mathrm{}`$ then the wavelet will be nonsingular and of finite energy.

To obtain a specific example it only remains to finish the complete specification of the potential $`\psi (r,t)`$. One particularly simple choice is to take one of the real scalar wavelets of the previous section $`\psi \varphi _{1,2}`$, with $`\zeta ^a=\zeta V^a=aV^a`$. Since the electric and magnetic fields are now specified in terms of derivatives of a smooth bounded function, the electromagnetic field is similarly smooth and bounded. For the energy, write it in the form

$$=\frac{1}{3}\mathrm{d}^3x\left\{\psi _t^4\psi _t\psi _t^3\psi \right\}.$$
(66)
Whence, for either $`\psi \varphi _{1,2}`$, an integration carried out at arbitrary $`t`$ yields the time-independent quantity

$$_{1,2}=\frac{1}{2}\frac{\varphi _0^2}{a}.$$
(67)
The calculation has for convenience been carried out in the zero-momentum frame of the wavelet. In a general Lorentz frame we would have

$$P_{1,2}^a=\frac{1}{2}\frac{\varphi _0^2}{a}V^a=\frac{1}{2}\frac{\varphi _0^2}{a^2}\zeta ^a.$$
(68)
Thus this field configuration, as for the scalar case, is Lorentz covariant, bounded, finite energy and zero action. The vector nature of the Maxwell field has added technical complications, but there is no real change in basic principles.

Aside: In closing this section, I should mention one other wavelet that is particularly attractive. If one makes use of the pseudo-differential operator $`(^2)^{1/2}`$ one could write

$$\psi =(^2)^{1/2}\varphi _{1,2}=\frac{1}{\mathrm{\Gamma }(1/2)}_0^{\mathrm{}}\frac{\mathrm{d}t}{t^{3/2}}\mathrm{exp}[t^2]\varphi _{1,2}.$$
(69)
With this choice of $`\psi `$ the energy integral simplifies

$$_{1,2}=\frac{1}{3}\mathrm{d}^3x\left\{(_t\varphi _{1,2})^2+(\varphi _{1,2})^2\right\}=\frac{1}{6}\pi ^2|\varphi _0^2|a.$$
(70)
The price paid for making the energy look simple is that the electric and magnetic fields are much more complicated to calculate.

Comment: Because the optics and engineering literature generally does not use manifestly Lorentz covariant notation, it can be very time consuming to calculate the 4-momentum of a specific pulse. Indeed only very recently  has Lekner provided specific and explicit computations of both $``$ and $`\stackrel{}{\mathrm{}}`$ (as well as the angular momentum) for a pulse similar to that considered above. As expected (once one has the covariant perspective advocated in this article) $`\stackrel{}{\mathrm{}}<`$, indicating the existence of a zero-momentum frame for the pulses of this type .

## 5 Yang-Mills field

Once we have the Maxwell wavelet above, a Yang–Mills wavelet is straightforward, indeed trivial. Let $`𝚲`$ be any constant matrix in the center of the gauge group and set $`𝐀^a=A^a𝚲`$. The construction is so simple that there is really nothing extra beyond the Maxwell wavelet considered above.

## 6 Null $`\zeta `$

Let us now return to the original scalar complex wavelet. Suppose the 4-vector $`\zeta `$ is null. Then because the numerator vanishes identically the original definition above gives $`\varphi 0`$. We should at a minimum change our field definition to read

$$\varphi (x)=\frac{\psi _0}{\eta _{ab}[x^ax_0^ai\zeta ^a][x^bx_0^bi\zeta ^b]}.$$
(71)
Then without loss of generality we go into the frame

$$\zeta ^a=(a,0,0,a),$$
(72)
and then

$$\varphi (x)=\frac{\psi _0}{[tia]^2x^2y^2[zia]^2}.$$
(73)
That is

$$\varphi (x)=\frac{\psi _0}{r^2t^2+2ia(tz)}.$$
(74)
Note

$$|\varphi (x)|=\frac{\psi _0}{\sqrt{(r^2t^2)^2+4a^2(tz)^2}}.$$
(75)
The denominator now vanishes when $`zt`$ and $`x=y=0`$, so that the field is divergent on the beam axis. Suppose we write $`R=\sqrt{x^2+y^2}`$ then

$$\varphi (x)=\frac{\psi _0}{R^2+z^2t^2+2ia(tz)}.$$
(76)
So we see that the field drops off as $`1/R^2`$ as we move away from the beam axis. \[And more critically, the field blows up as $`1/R^2`$ as we approach the beam axis.\] Attempts at calculating the energy and action now lead to divergent integrals. In other words, despite the fact that it still solves the wave equation, for null $`\zeta `$ this is not a particularly useful field configuration.

## 7 Spacelike $`\zeta `$

Foe the complex scalar wavelet, suppose the 4-vector $`\zeta `$ is spacelike. Then without loss of generality we go into the infinite velocity frame where

$$\zeta ^a=(0,0,0,a),$$
(77)
and then

$$\varphi (x)=\frac{\varphi _0a^2}{t^2x^2y^2[zia]^2}.$$
(78)
That is

$$\varphi (x)=\frac{\varphi _0a^2}{r^2t^2a^22iaz}.$$
(79)
Note

$$|\varphi (x)|=\frac{\varphi _0a^2}{\sqrt{(r^2t^2a^2)^2+4a^2z^2}}.$$
(80)
The denominator now vanishes when $`z=0`$ and $`x^2+y^2=a^2+t^2`$. That is, the field is divergent on a time-dependent circle orthogonal to the beam axis. There is a short distance singularity as one approaches this circle, and the energy and action integrals diverge. In other words, despite the fact that it still solves the wave equation, for spacelike $`\zeta `$ this is not a useful field configuration.

## 8 Discussion

The physical wavelet discussed in this article is important because it represents a qualitatively different extended field configuration of a type not normally encountered in particle physics. The wavelet is neither a soliton, nor an instanton, nor a sphaleron though it shares properties with all three of these extended objects:

* Like the soliton it lives in physical time (Minkowski space, not Euclidean space), and possesses a well-defined 4-velocity.
* Like the instanton it “dies away” in the infinite past and future.
* Like the instanton it possesses a continuously adjustable scale parameter.
* Like the sphaleron it is unstable to dispersal.

Because the wavelet fields are bounded and finite energy, wavelet configurations *will* be classically excited at any finite temperature. Because the wavelet configuration has zero action, arbitrarily complicated combinations of these physical wavelets can be added to the field configurations appearing in Feynman’s path integral without modifying the phase — quantum mechanically there is no “cost” in adding these configurations to the Lorentzian path integral and they *will* contribute.

Other “localized waves” might be interesting in specific applications but the particular example discussed in this article is important because of its extreme simplicity and pleasant behaviour.

###### Acknowledgments.

I wish to thank John Lekner for stimulating my interest in these issues. I also wish to thank Damien Martin for bringing the whole Lorenz/Lorentz issue to my attention.