text:
# A New Way to Make Waves

## I Introduction

I will describe a new way to “make waves”. Wave phenomena is treated mathematically by hyperbolic equations. Einstein’s theory of general relativity is the supreme example. It treats gravity as a distortion in the geometry of space-time. Gravitational waves are ripples in the geometry which travel through space changing the shape and size of objects in their path. At this historic time, Einstein’s theory is providing the basis for a new way to observe the distant universe. A worldwide network of gravitational wave observatories is under construction. It is designed to detect the gravitational waves produced throughout the cosmos by collisions between black holes.

The new way to make waves which I will describe was developed to simulate the production of gravitational waves by using a computer to solve Einstein’s equations. But the approach applies to any wave phenomena. In order to get the fundamental ideas across, I present them in terms of very simple systems. Only near the end will I give in to the urge to tell you about the application to black holes. However, even in simple examples, the ideas are easiest to understand in terms of space-time language and pictures.

The simplest hyperbolic equation is the 1-dimensional advection equation

$$(_t+c_x)\mathrm{\Phi }=0.$$
(1)
It has the general solution $`\mathrm{\Phi }(t,x)=F(xct)`$, which describes a wave traveling with velocity $`c`$ with initial waveform $`F(x)`$ at $`t=0`$. The curves $`x=x_0+ct`$ are called the characteristics of the equation. They are the paths in space-time (worldlines) along which you must move to “ride the wave” in the sense that $`\mathrm{\Phi }=const`$ along a characteristic worldline. Since $`F`$ is arbitrary, we can choose $`F=0`$ for $`x<0`$ and $`F=1`$ for $`x>0`$ to obtain a wave with a shock front along the characteristic $`x=ct`$. In modern parlance, this shock wave carries one bit of information. It illustrates the general way information is propagated by waves along characteristics in any hyperbolic system.

The standard way to solve a hyperbolic equation such as Eq. (1) is by means of the Cauchy problem. One poses initial Cauchy data $`\mathrm{\Phi }(t=0,x)`$ in some spatial region $`\mathrm{\Omega }`$ and then evolves the solution in time. The characteristics causally relate the initial data to a unique evolution throughout some domain of dependence $`D(\mathrm{\Omega })`$ in space-time. As a result, shock fronts, which represent a sudden signal, can only occur across characteristics. Hyperbolic equations lead to ordinary differential equations along the characteristics which govern the propagation of shock discontinuities. That makes it important for the purpose of numerical simulation to enforce the propagation along characteristics as extensively as possible. In complicated cases, where the wave velocity has functional dependence $`c(\mathrm{\Phi },x)`$ and an analytic solution is not possible, there still exists a well established theory of the general properties of hyperbolic systems based upon characteristics .

The method of characteristics  is a computational algorithm for evolving Cauchy data based upon this theory. As applied to the advection equation (1), the characteristic velocities are calculated at a given discretized time $`t_n=n\mathrm{\Delta }t`$ and the algorithm used to determine the field at time $`t_{n+1}`$ by requiring that $`\mathrm{\Phi }`$ be constant along the characteristic curve from a space-time point $`P`$ at time level $`n`$ to point $`Q`$ at time level $`n+1`$, i.e.

$$\mathrm{\Phi }_Q=\mathrm{\Phi }_P,$$
(2)
where $`x_Q=x_P+c\mathrm{\Delta }t`$. When the wave velocity is not constant in either space or time, the points $`P`$ and $`Q`$ along the characteristics do not in general lie exactly on spatial grid points, so that interpolations are necessary. But, because there are no sudden changes along characteristics (only across them), such interpolations give excellent accuracy by riding the wave like a surfer.

The method of characteristics for the Cauchy problem extends to the generalization of Eq. (1) to any symmetric hyperbolic system of first differential order equations in a multi-dimensional space with many evolution variables and with source terms for the creation of waves . However, there is another classification of hyperbolic systems, which is more familiar to physicists, based upon second order differential equations whose principal part has the form

$$g^{\alpha \beta }_\alpha _\beta \mathrm{\Phi }=0,$$
(3)
where $`\alpha `$ is a space-time index, i.e. $`x^\alpha =(t,x,y,z)`$ for 4-dimensional space-time, and a sum over the values of the repeated indices is implied. Equation (3) is classified as hyperbolic if the symmetric matrix $`g^{\alpha \beta }`$ has one negative eigenvalue and its remaining eigenvalues are positive. Displacements $`x^\alpha x^\alpha +dx^\alpha `$ along characteristic worldlines satisfy

$$g_{\alpha \beta }dx^\alpha dx^\beta =0,$$
(4)
where $`g_{\alpha \beta }`$, called the metric tensor, is the inverse matrix to $`g^{\alpha \beta }`$. The simplest example is the wave equation in one spatial dimension,

$$[\frac{1}{c^2}_t^2_x^2]\mathrm{\Phi }=0,$$
(5)
with the displacement along the characteristics satisfying $`c^2dt^2dx^2=0`$. By rewriting Eq. (5) in the form

$$(\frac{1}{c}_t_x)(\frac{1}{c}_t+_x)\mathrm{\Phi }=0,$$
(6)
a comparison with Eq. (1) shows that the general solution is $`\mathrm{\Phi }(t,x)=F(xct)+G(x+ct)`$. There are two characteristics $`x=x_0\pm ct`$ through each spatial point $`x_0`$. Information is propagated along these characteristics which now criss-cross the space-time. The conventional method of characteristics could be applied to Eq. (6) by the standard technique of reducing it to a coupled system of first order differential equations. But it is simpler to integrate Eq.(6) over a parallelogram $`\mathrm{\Sigma }`$ in space time whose sides are characteristics meeting at the corners $`P`$, $`Q`$, $`R`$ and $`S`$. This gives the 2-dimensional version of Eq. (2)

$$\mathrm{\Phi }_Q\mathrm{\Phi }_P+\mathrm{\Phi }_R\mathrm{\Phi }_S=0,$$
(7)
as illustrated in Fig. 1 for the case $`c=1`$ where the characteristics $`x=x_0\pm t`$ have $`45^o`$ slope in space-time. By using Eq. (7), the method of characteristics can be implemented as a computational algorithm for the Cauchy evolution of Eq. (5) by positioning the point $`Q`$ at time level $`n+1`$, the points $`P`$ and $`S`$ at time level $`n`$ and the point $`R`$ at time level $`n1`$. On a $`(t,x)`$ computational grid satisfying $`\mathrm{\Delta }x=c\mathrm{\Delta }t`$, the characteristics pass through diagonal grid points. The Cauchy data, for the wave equation, which consists of the initial values $`\mathrm{\Phi }(t=0,x)`$ and $`_t\mathrm{\Phi }(t=0,x)`$, is used to initialize an iterative evolution scheme at time levels $`n=0`$ and $`n=1`$. An exceptional feature is that the one-dimensional wave equation with constant wave velocity can be evolved without error by using Eq. (7) as a finite difference equation on such a grid. For a variable velocity, the points $`P`$, $`Q`$, $`R`$ and $`S`$ cannot be placed exactly on the grid and interpolations are necessary.

For the wave equation in two or more spatial dimensions, the manner in which characteristics determine domains of dependence and lead to propagation equations is qualitatively the same. The major difference is that an infinite number of characteristics now pass through each point. For the 3-dimensional wave equation,

$$[\frac{1}{c^2}_t^2_x^2_y^2_z^2]\mathrm{\Phi }=0$$
(8)
the characteristics which pass through the point $`(x_0,y_0,z_0)`$ at time $`t_0`$ are straight lines which trace out an expanding spherical wavefront of radius $`c(tt_0)`$ at a given time $`t`$; or, from a space-time point of view, these characteristics trace out the 3-dimensional characteristic cone

$$(xx_0)^2+(yy_0)^2+(zz_0)^2c^2(tt_0)^2=0.$$
(9)
In electrodynamics or relativity, the characteristics are light rays and this is called the light cone. I will often use that language here. The analogue in hydrodynamics is the Mach cone.

The future (past) light cone consists of the radially outward (inward) characteristics parameterized by $`t>t_0`$ ($`t<t_0`$). There is a 2-parameter set of characteristics through each point corresponding to the sphere of angular directions $`(\theta ,\varphi )`$ at that point. This leads to some arbitrariness in formulating an evolution algorithm for Cauchy data based upon the method of characteristics. There are an infinite number of characteristics and associated propagation equations which can be used to evolve Cauchy data from time $`t_0`$ to time $`t_0+\mathrm{\Delta }t`$.

For a practical numerical scheme, it is thus necessary either to average these propagation equations appropriately over the sphere of characteristic directions at or select out some finite number of characteristics whose propagation equations comprise a nonredundant set. The latter approach has been successfully carried out by Butler . For the case of plane flow of an inviscid fluid (a problem in two spatial dimensions), he formulates an algorithm based upon four “preferred” characteristics. In this problem, the geometry is further complicated because the characteristics are dynamically dependent upon the fluid variables, in contrast to the essentially time independent characteristic cones of Eq. (9). In the numerical scheme, the characteristics must themselves be determined by some finite difference approximation.

That summarizes the Cauchy problem and its solution by the method of characteristics.

## II Characteristic evolution

I now present a procedure which avoids the arbitrariness and awkwardness of the method of characteristics in a 3-dimensional Cauchy evolution. It is based upon a characteristic initial value approach rather than a Cauchy approach. In order to understand the distinction, it is essential to view space-time as 4-dimensional, with initial data given on a three-dimensional hypersurface. The Cauchy approach is based upon evolution of initial data given on a spatial hypersurface in space-time, i.e. points of space at the same instant of time $`t=t_0`$, to data at a later time $`t=t_0+\mathrm{\Delta }t`$. In the characteristic approach, data is initially posed on an outgoing characteristic cone (light cone), which defines a hypersurface at constant retarded time $`u=u_0`$. Characteristic evolution then proceeds iteratively to characteristic cones $`u_n=u_0+n\mathrm{\Delta }u`$. The numerical grid is intrinsically based upon these outgoing characteristic hypersurfaces. The algorithm for evolving from retarded time $`u_n`$ to $`u_{n+1}`$ is based upon characteristics which are uniquely and intrinsically picked out by the geometry of the retarded time hypersurfaces.

This characteristic evolution procedure radically differs from Cauchy evolution. It uses concepts developed in the 1960’s for studies of general relativity . These were prompted by the inability of the major mathematical tools such as Green’s functions and Fourier analysis to overcome the difficulties posed by nonlinearity of the equations and the ambiguity of the general coordinate freedom in the theory. Its success later motivated a more extensive mathematical treatment of the characteristic initial value problem . These new techniques were especially designed for investigation of gravitational waves. A computational algorithm based on this approach has been the most successful in simulating the production of gravitational waves from black holes.

The approach has two novel ingredients:

* the use of characteristic hypersurfaces to formulate a characteristic initial value problem (CIVP) and
* the use of compactification methods to describe on a finite grid the waves propagating to infinity.

Characteristic hypersurfaces are 3-dimensional sets traced out in space-time by characteristic worldlines or, equivalently, by their wavefronts. They provide a natural coordinate system to describe waves . It is very fruitful to use an initial value scheme which describes time evolution by means of the retarded time coordinate defined by characteristic hypersurfaces as a substitute for the familiar Cauchy scheme based upon constant time hypersurfaces. As will be shown, this approach uses a completely different form of the mathematical equations and the free initial data.

Compactification methods provide a rigorous description of the radiation field of a source observed in the asymptotic limit of going to infinity along a characteristic worldline.. The key idea is to introduce a new coordinate which ranges over values from 0 to 1 as the actual distance from the source ranges from 0 to $`\mathrm{}`$. The hyperbolic equations are rewritten in terms of these new coordinates. Asymptotic behavior at the “points at infinity” can then be studied in terms of the new coordinate which ranges over finite values. In this way, the concept of the radiation zone as an asymptotic limit at infinity is given rigorous meaning. Characteristic hypersurfaces are important here since waves travel to infinity along characteristic worldlines, not along Cauchy hypersurfaces of constant time. Even for field equations as complicated as those of general relativity, this procedure provides a finite geometrical description of waves travelling to infinity. The limit points at infinity form a boundary to the compactified space-time, which I will refer to as radiative infinity. At a given retarded time, this boundary has the topology of a sphere, representing a sphere of observers at infinity.

It should be emphasized that radiative infinity differs drastically from spatial infinity (the limit of going to infinity holding time constant). Early considerations of compactifying infinite space for computational purposes were discarded because they were based upon spatial infinity . An attempt to cover infinite space this way by a finite grid at constant time fails in a hyperbolic problem because there is necessarily an infinite physical distance between a grid point at infinity and its neighbors. This makes it impossible on a spatial grid at fixed time to resolve radiation with finite wave length which propagates to infinity. In contrast, for a grid constructed on a characteristic hypersurface, the grid points ride the wave without noticing its finite wavelength in the approach to radiative infinity.

In terms of characteristic coordinates $`u=ctx`$ and $`v=ct+x`$, Eq. (6) becomes

$$_u_v\mathrm{\Phi }=0,$$
(10)
Here $`\mathrm{\Psi }=_v\mathrm{\Phi }`$ satisfies the propagation equation

$$_u\mathrm{\Psi }=0$$
(11)
along the characteristics in the $`u`$-direction. This is the essence of how the use of characteristic coordinates simplifies the treatment of waves.

These ideas provide the physical basis for a new computational algorithm . I will illustrate how it applies to a nonlinear version of the wave equation in 3 spatial dimensions. However, this evolution algorithm can be taken over intact to other hyperbolic physical systems, including electromagnetic fields, the Yang-Mills gauge fields of elementary particle physics, as well as general relativity, because of the common mathematical structure of these theories as second differential order hyperbolic equations. Although it has not been explored how this approach might be implemented in the case of a first differential order hyperbolic system, such as hydrodynamics, the general ideas should be applicable.

## III Characteristic initial value problem for nonlinear waves

As a simple illustration of the characteristic initial value problem, consider the nonlinear scalar wave equation (SWE) in 3-spatial dimensions, which we write in spherical polar coordinates $`(t,r,\theta ,\varphi )`$ as

$`[{\displaystyle \frac{1}{c^2}}_t^2_x^2_y^2_z^2]\mathrm{\Phi }={\displaystyle \frac{1}{c^2}}_t^2\mathrm{\Phi }{\displaystyle \frac{1}{r}}_r^2(r\mathrm{\Phi })+{\displaystyle \frac{L^2\mathrm{\Phi }}{r^2}}=S(\mathrm{\Phi })`$ (12)
$`r=\sqrt{(}x^2+y^2+z^2),`$ (13)
where $`c`$ is the wave velocity, $`L^2`$ denotes the standard angular momentum operator

$$L^2\mathrm{\Phi }=\frac{_\theta (\mathrm{sin}\theta _\theta \mathrm{\Phi })}{\mathrm{sin}\theta }\frac{_\varphi ^2\mathrm{\Phi }}{\mathrm{sin}^2\theta }.$$
(14)
and $`S(\mathrm{\Phi })`$ represents a nonlinear source term. Rather than using ordinary time $`t`$, characteristic evolution uses the “retarded time” coordinate $`u=ctr`$. The outgoing radial characteristics are the curves of constant $`u`$, $`\theta `$ and $`\varphi `$. These are the curves in the $`r`$-direction, holding $`u=const`$, shown in the space-time Fig. 1. In $`(u,r,\theta ,\varphi )`$ coordinates, the SWE (12) takes the form

$$2_u_rg=_r^2g\frac{L^2g}{r^2}+rS.$$
(15)
where $`g=r\mathrm{\Phi }`$.

The striking feature about Eq.(15) is that it is only first differential order in retarded time $`u`$, unlike the more familiar form of the SWE (12) which is of second differential order in time $`t`$. When data is given on a characteristic initial hypersurface $`u=u_0`$, we need only specify the initial value of the field $`\mathrm{\Phi }`$, and then use Eq. (15) to calculate its retarded time derivative $`_u\mathrm{\Phi }`$ in order to evolve the initial data. This is in contrast to the conventional Cauchy scheme, where both $`\mathrm{\Phi }`$ and $`_t\mathrm{\Phi }`$ must be supplied at initial time $`t=t_0`$ and Eq. (12) is the used to compute the second time derivative $`_t^2\mathrm{\Phi }`$ in order to evolve the data.

In a computational implementation of the CIVP, rather than finite differencing Eq. (15) directly, it is advantageous to first convert it into an integral equation which is subsequently discretized. Using both the outgoing characteristic coordinate $`u`$ and the ingoing characteristic coordinate $`v=ct+r`$, the SWE (15) takes the form

$$4_u_vg=\frac{L^2g}{r^2}+rS(\frac{g}{r})$$
(16)
where $`g=r\mathrm{\Phi }`$. In the $`u`$-$`v`$ plane formed by fixing the angular coordinates $`(\theta ,\varphi )`$, we construct a parallelogram $`\mathrm{\Sigma }`$ made up of incoming and outgoing radial characteristics which intersect at vertices $`P,Q,R,S`$ as depicted in Fig. 1. By integrating Eq. (16) over the area $`\mathrm{\Sigma }`$ bounded by these vertices, we may establish the identity

$$g_Q=g_P+g_Sg_R+\frac{1}{2}_\mathrm{\Sigma }𝑑u𝑑r\left[\frac{L^2g}{r^2}+rS(\frac{g}{r})\right].$$
(17)
This simple identity is the 3-dimensional analogue of Eq. (7), adapted to include a source term. It is the starting point for an evolution algorithm which incorporates the essential role that characteristics play in the SWE.

In order to study the far field wave behavior, we transform this equation to the new radial coordinate

$$x=r/(1+r),\mathrm{\hspace{0.17em}\hspace{0.17em}0}x1.$$
(18)
This serves to map an infinite radial domain into a finite coordinate region, and assigns infinitely distant radial points to the edge of the coordinate patch (radiative infinity) at $`x=1`$, where the radiation signal can be identified.

### A Numerical Algorithm

To develop a discrete evolution algorithm, we work on the lattice of points

$`u_n`$ $`=`$ $`n\mathrm{\Delta }u`$ (19)
$`x_i`$ $`=`$ $`i\mathrm{\Delta }x`$ (20)
$`\zeta _{j,k}`$ $`=`$ $`(j+ik)\mathrm{\Delta }\varphi `$ (21)
where a complex stereographic coordinate $`\zeta `$ is used to cover the the sphere in two patches (North and South) in order to avoid the polar singularities of the spherical coordinates $`(\theta ,\varphi )`$. We denote the field at these sites by

$$g_{ijk}^n=g(u_n,x_i,\zeta _{j,k}).$$
(22)
(We will generally suppress the angular indices $`j`$ and $`k`$.)

With respect to the $`(u,x)`$ coordinate grid, it is not possible to place the corners P, Q, R and S at grid points since the slope of the characteristics in the compactified $`x`$ coordinate depends upon location. As a consequence, the field $`g`$ at these points must be interpolated from neighboring grid points. The essential feature of the placement of the parallelogram on the grid is that the sides formed by the ingoing characteristics intersect adjacent $`u`$-hypersurfaces at equal but opposite $`x`$-displacements from the neighboring grid points, as illustrated in Fig. 2. The field values at the vertices of the parallelogram are obtained by quadratic interpolation. Cancellations between the interpolation errors at the four vertices yields the accuracy

$$g_Qg_Pg_S+g_R=G_QG_PG_S+G_R+O((\mathrm{\Delta }x)^3\mathrm{\Delta }u),$$
(23)
where $`G`$ represents the exact analytic solution.

The integral in Eq. (17) can be evaluated by treating the integrand as a constant over the parallelogram, with value at the center. The radial coordinate of the point at the center is $`r_c=(r_P+r_S)/2`$. To compute the nonlinear term, the value of $`g`$ at $`r_c`$ is taken as the average $`g_c=(g_P+g_S)/2`$, with $`g_P`$ and $`g_S`$ evaluated from second-order linear interpolations over adjacent points on the grid. The angular derivatives in Eq. (17) are replaced with standard second-order-accurate finite difference approximations. $`L^2g`$ is calculated on the grid points, and the same interpolation procedure is used to obtain the value of $`L^2g_c`$. The integral term is then approximated by

$`{\displaystyle _\mathrm{\Sigma }}[L^2g+r^3S(g/r)]𝑑u𝑑r/r^2=[L^2g_c+r_c^3S(g_c/r_c)]{\displaystyle _\mathrm{\Sigma }}𝑑u𝑑r/r^2`$ (24)
$`=2\mathrm{log}\left({\displaystyle \frac{r_Qr_R}{r_Pr_S}}\right)\left[L^2g_c+r_c^3S({\displaystyle \frac{g_c}{r_c}})\right],`$ (25)
where the integrand is accurate to second order in $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }\varphi `$. The resulting finite difference equation

$`g_i^{n+1}(x_Q+x_Px_{i2}x_{i1})=`$ (30)
$`2g_{i1}^{n+1}(x_Q+x_Px_{i2}x_i)g_{i2}^{n+1}(x_Q+x_Px_{i1}x_i)`$
$`+\{g_{i+1}^n(x_S+x_Rx_{i1}x_i)2g_i^n(x_S+x_Rx_{i1}x_{i+1})`$
$`+g_{i1}^n(x_S+x_Rx_ix_{i+1})\}{\displaystyle \frac{(x_Sx_R)}{(x_Qx_P)}}`$
$`+{\displaystyle _\mathrm{\Sigma }}{\displaystyle \frac{drdu}{r^2}}\left[L^2(g_c)+r_c^3S({\displaystyle \frac{g_c}{r_c}})\right]{\displaystyle \frac{(\mathrm{\Delta }x)^2}{(x_Qx_P)}}.`$
relates values of $`g_i^{n+1}`$ with values at neighboring grid points which are either earlier in retarded time ($`g_{i1}^n,g_i^n,g_{i+1}^n`$), or else contemporary but located at smaller radius ($`g_{i2}^{n+1},g_{i1}^{n+1}`$). Consequently, it is possible to move through the grid, computing $`g_i^{n+1}`$ explicitly by an orderly march. This is achieved by starting at the origin at time $`u_{n+1}`$. Field values of $`g=r\mathrm{\Phi }`$ vanish there. Step outward to the next radial point, using Eq. (30) for all angular sites on the grid, and iterate this march out to radiative infinity thus updating the characteristic cone at $`u_{n+1}`$ and completing one retarded time step. This march is then iterated in retarded time.

The algorithm steps $`g`$ radially outward one cell with a local error of fourth order in grid size. This leads to second order global accuracy which is confirmed by convergence tests using known analytic solutions . A complete specification of the algorithm would require a description of how the startup procedure at the origin is handled and how stereographic coordinates are used to compute angular derivatives in a smooth way. Details of the use of North and South stereographic grids are given in Ref. . Physical behavior of nonlinear waves is treated in Ref. . Construction of an exact nonspherical solution for an $`S=\mathrm{\Phi }^3`$ self-interaction allows calibration of the algorithm in the nonlinear case where physical singularities form. The predicted second order accuracy is confirmed right up to the formation of the singularity. Other choices of nonlinear potential allow simulation of solitons.

The Courant-Friedrichs-Lewy (CFL) condition that the numerical domain of dependence contain the physical domain of dependence (determined by the characteristics) is a necessary condition for convergence of a finite difference algorithm. For a grid point at $`(u,r,\theta )`$, there are three critical grid points, at $`(u\mathrm{\Delta }u,r+\mathrm{\Delta }r,\theta )`$ and $`(u\mathrm{\Delta }u,r\mathrm{\Delta }r,\theta \pm \mathrm{\Delta }\theta )`$, which must lie inside its past characteristic cone. These gives rise to the inequalities $`\mathrm{\Delta }u<2\mathrm{\Delta }r`$ and $`\mathrm{\Delta }u<\mathrm{\Delta }r+(\mathrm{\Delta }r^2+r^2\mathrm{\Delta }\theta ^2)^{1/2}`$. At large $`r`$, the second inequality becomes $`\mathrm{\Delta }u<r\mathrm{\Delta }\theta `$ and the Courant limit on the time step is essentially the same as for a Cauchy evolution algorithm. However, near the vertex of the cone, the second inequality gives a stricter condition

$$\mathrm{\Delta }u<K\mathrm{\Delta }r\mathrm{\Delta }\theta ^2,$$
(31)
where the value of $`K`$ depends upon the start up procedure at the vertex. For the scalar wave equation, these stability limits were confirmed by numerical experiments  and it was found that $`K4`$.

Local von Neumann stability analysis leads to no constraints on the algorithm. This may seem surprising because no analogue of a CFL condition on the time step arises. It can be understood in the following vein. The local structure of the code is implicit, since it involves 3 points at the upper time level. Implicit algorithms do not necessarily lead to a CFL condition. However, the algorithm is explicit in the way that the evolution starts up as an outward radial march from the origin. It is this startup procedure that introduces a CFL condition.

Operating within the CFL limit, the algorithm gives a stable, globally second order accurate evolution on a compactified grid . (In some nonlinear applications, artificial dissipation is necessary for stability .) Radiative infinity behaves as a perfectly transmitting boundary so that no radiation is reflected back into the system. Numerical evolution satisfies a conservation law relating the loss of energy to the radiation flux at infinity.

That summarizes the characteristic initial value problem and its implementation as a new computational approach to simulate waves.

## IV Cauchy-characteristic matching

Characteristic evolution has many advantages over Cauchy evolution. Its one disadvantage is caused by the existence of either (i) caustics where neighboring characteristics focus or (ii), a milder version of this, cross-over points where two distinct characteristics collide. The vertex of the characteristic cone is a highly symmetric example of a point caustic where a complete sphere of characteristics focus. I have already discussed how a point focus gives rise to a strong limitation imposed by the CFL condition.

Cauchy-characteristic matching (CCM) is a way to avoid such limitations by combining the strong points of characteristic and Cauchy evolution in formulating a global evolution.  Here I illustrate the application of CCM to the nonlinear wave equation (12). This problem requires boundary conditions at infinity which ensure that the total energy and the energy loss due to radiation are both finite. In a 3-dimensional problem, these are the conditions responsible for the proper $`1/r`$ asymptotic decay of the radiation fields. However, for practical purposes, in the computational treatment of such a system by the Cauchy problem, an outer boundary is artificially established at some large but finite distance. Some condition is then imposed upon this boundary in an attempt to approximate the proper asymptotic behavior at infinity. Such an artificial boundary condition (ABC) typically causes partial reflection of the outgoing wave back into the system , which contaminates the accuracy of the evolution and the radiated signal. Furthermore, nonlinear wave equations often display backscattering so that it may not be correct to try to entirely eliminate incoming radiation from the numerical solution. The errors introduced by ABC’s are of an analytic origin, essentially independent of the computational discretization. In general, a systematic reduction of the error can only be achieved by simultaneously refining the grid and moving the computational boundary to a larger radius, which is computationally very expensive for three-dimensional simulations. CCM provides a global solution which does not introduce error at the analytic level.

For linear wave problems, a variety of ABC’s have been proposed. For recent reviews, see Ref’s. . During the last two decades, local ABC’s in differential form have been extensively employed by several authors  with varying success. Some local ABC’s have been derived for the linear wave equation by considering the asymptotic behavior of outgoing solutions ; this approach may be regarded as a generalization of the Sommerfeld outgoing radiation condition. Although such ABC’s are relatively simple to implement and have a low computational cost, their final accuracy is often limited because their simplifying assumptions are rarely met in practice . Systematic improvement of the accuracy of local ABC’s can only be achieved by moving the computational boundary to a larger radius.

The disadvantages of local ABC’s have led to implementation of nonlocal ABC’s based on integral representations of the infinite domain problem . Even for problems where the Green’s function is known and easily computed, such approaches were initially dismissed as impractical ; however, the rapid increase in computer power has made it possible to implement nonlocal ABC’s for the linear wave equation even in 3 space dimensions . For a linear problem, this can yield numerical solutions which converge to the exact infinite domain problem as the grid is refined, keeping the artificial boundary at a fixed distance. However, due to nonlocality, the computational cost per time step usually grows at a higher power of grid size ($`O(N^4)`$ per time step in a 3-dimensional problem with $`O(N^3)`$ spatial grid points) than in a local approach , which is demanding even for today’s supercomputers. Further, the applicability of current nonlocal ABC’s is restricted to problems where nonlinearity may be neglected near the grid boundary .

To my knowledge, only a few works have been devoted to the development of ABC’s for strongly nonlinear problems . In practice, nonlinear problems are often treated by linearizing the governing equations in the far field , using either local or nonlocal linear ABC’s . Besides introducing an approximation at the analytical level, this procedure requires that the artificial boundary be placed sufficiently far from the strong-field region, which sharply increases the computational cost in multidimensional simulations. There seems to be no currently available ABC which is able to produce numerical solutions which converge (as the discretization is refined) to the infinite domain exact solution of a strongly nonlinear 3-dimensional wave problem, keeping the artificial boundary at a fixed location.

For such nonlinear problems, CCM produces an accurate solution out to radiative infinity with effort $`O(N^3)`$ per time-step. CCM increases the total computational cost only by a factor $`2`$ with respect to a pure Cauchy algorithm with a local ABC. The use of numerical methods based upon matching a characteristic initial-value formulation and a Cauchy formulation can effectively remove the above difficulties associated with a finite computational boundary. There is no need to truncate space-time at a finite distance from the sources, since compactification of the radial coordinate makes it possible to cover the whole space-time with a finite grid. In this way, the true radiation zone signal may be computed. Although the characteristic formulation has stability limitations in interior region where the characteristic hypersurfaces can develop caustics, it proves to be both accurate and computationally efficient in the treatment of the exterior, caustic-free region.

CCM is a new approach to global numerical evolution which is free of error at the analytic level. The characteristic algorithm provides the outer boundary condition for the interior Cauchy evolution, while the Cauchy algorithm supplies the inner boundary condition for the characteristic evolution. Since CCM consists of discretizing an exact analytic treatment of the radiation from source to radiative infinity, it generates numerical solutions which converge to the exact analytic solution of the radiating system even in the presence of strong nonlinearity. Thus, any desired accuracy can be achieved by refining the grid, without moving the matching boundary. In practice, the method performs extremely well even at moderate resolutions.

### A CCM for nonlinear waves

As an illustration of the computational implementation of CCM , consider again the nonlinear 3-dimensional wave equation (12). For simplicity, set the velocity $`c=1`$. In the standard computational implementation of the Cauchy problem for (12), initial data $`\mathrm{\Phi }(t_0,x,y,z)`$ and $`_t\mathrm{\Phi }(t_0,x,y,z)`$ are assigned and evolved in a bounded spatial region, with some ABC imposed at the computational boundary. In a characteristic initial-value formulation the SWE (12) is reexpressed in the form of Eq. (15) in terms of $`g=r\mathrm{\Phi }`$, using standard spherical coordinates and a retarded time coordinate $`u=tr`$. The initial data $`g(u_0,r,\theta ,\varphi )`$, on an initial outgoing characteristic cone $`u=u_0`$ is then evolved globally out to radiative infinity.

In CCM, (12) is solved in an interior region $`rR_m`$ using a Cauchy algorithm, while a characteristic algorithm solves the retarded coordinate version (15) for $`rR_m`$. The matching procedures ensure that, in the continuum limit, $`\mathrm{\Phi }`$ and its gradient are continuous across the interface $`r=R_m`$. This is a requirement for any consistent matching algorithm, since a discontinuity in the field or its gradient could act as a spurious boundary source, contaminating both the interior and exterior evolutions.

For technical simplicity, I illustrate the details of the method here for spherically symmetric waves but it has been successfully implemented in fully nonlinear 3-dimensional problems without symmetry. In a 3-dimensional problem without symmetry, the characteristic evolution is carried out on an exterior spherical grid, while the Cauchy evolution uses a Cartesian grid covering the interior spherical region. Although a spherical grid could also be used in the interior, a Cartesian grid avoids the necessity of cumbersome numerical procedures to handle the singularity of spherical coordinates at the origin. A Cartesian discretization in the interior and a spherical discretization in the exterior are the coordinates natural to the geometries of the two regions. However, this makes the treatment of the interface somewhat involved; in particular, guaranteeing the stability of the matching algorithm requires careful attention to the details of the inter-grid matching. Nevertheless, there is a reasonably broad range of discretization parameters for which CCM is stable .

With the substitution $`G=r\varphi `$ and the use of spherical coordinates, the spherically symmetric version of the SWE (12) reduces to the 1-dimensional wave equation

$$_{tt}G=_{rr}G+rS.$$
(32)
The initial Cauchy data is $`G(t_0,r)`$ and $`_tG(t_0,r)`$ in the region $`0rR_m`$. Together with the regularity condition $`G(t,0)=0`$, these data determine a unique solution in the domain of dependence $`D_1^{}`$ indicated in Fig. 3. The outer boundary of the domain of dependence is the ingoing radial characteristic $`C_1^{}`$ described by $`r=R_mt+t_0`$. The solution cannot be constructed throughout the complete interior region $`rR_m`$ without additional information, which can be furnished by giving the value of $`G`$ on the outgoing characteristic $`C_{0^+}`$ described by $`r=R_m+tt_0`$ (see Fig. 3). In terms of the coordinates $`u=tr`$ and $`r`$, $`C_{0^+}`$ is described by $`u=t_0R_m`$ In these coordinates, expressing $`g(u,r)=G(u+r,r)`$, the spherically symmetric version of the wave equation (15) is

$$2_{ur}g=_{rr}g+rS.$$
(33)
A unique solution of (33) is determined by characteristic initial data consisting of the value of $`g`$ on the initial outgoing characteristic $`C_{0^+}`$ and on the ingoing characteristic $`C_1^{}`$. These data determine the solution uniquely throughout the future of $`C_1^{}`$ and $`C_{0^+}`$, i.e. the region $`D_{1^+}`$ in Fig. 3.

The matching scheme proceeds as shown in Fig. 3. First, initial Cauchy data are evolved from $`t_0`$ to $`t_1`$ throughout the region $`D_1^{}`$, which is in its domain of dependence. Next, this induces characteristic data on $`C_1^{}`$ which combined with the initial characteristic data on $`C_{0^+}`$ allows a characteristic evolution throughout the region $`D_{1^+}`$, bounded in the future by the characteristic $`C_{1^+}`$. The solution determined from this initial stage induces Cauchy data at time $`t_1`$ in the region $`rR_m`$, inside the matching boundary. This process can then be iterated to carry out the entire future evolution of the system.

CCM is a discretized version of this scheme in which the criss-cross pattern of characteristics inside the radius $`R_m`$ is at the scale of a grid spacing. The discretized evolution algorithm consists of the following steps (see Fig. 4):

Step 1. Cauchy evolution. The interior integration scheme is implemented on a uniform spatial grid $`r_i=i\mathrm{\Delta }r`$ ($`0iM`$) with outer radius $`R_B=M\mathrm{\Delta }r`$. We discretize (32) using the standard second-order finite difference scheme

$$\frac{G_i^{n+1}2G_i^n+G_i^{n1}}{(\mathrm{\Delta }t)^2}=\frac{G_{i+1}^n2G_i^n+G_{i1}^n}{(\mathrm{\Delta }r)^2}+r_iS_i^n,$$
(34)
where $`G_i^n=G(t_n,r_i)`$, $`S_i^n=S(t_n,r_i)`$, and $`t_n=t_0+n\mathrm{\Delta }t`$. The interior evolution is initialized by evaluating $`G_i^0`$ and $`G_i^1`$ $`(0iM)`$ to second order accuracy from the Cauchy initial data. In the $`n`$th time step, (34) is used to compute $`G_i^{n+1},1iM`$ in terms of field values $`G_i^{n1}`$ and $`G_i^n`$. The regularity of $`\mathrm{\Phi }`$ at $`r=0`$ implies that $`G_0^n=0`$ for all $`n`$. The boundary values $`G_{M+1}^n`$ which are required by (34) are supplied by the matching procedure (step 3).

Step 2. Characteristic evolution. The characteristic algorithm is implemented on a uniform grid based on the dimensionless compactified radial coordinate

$$\eta =1\frac{1}{1+r/R_m},\frac{1}{2}\eta 1$$
(35)
so that points at radiative infinity (corresponding to $`\eta =1`$) are included in the grid. In order to include one of the technical problems in matching an interior Cartesian grid to an exterior spherical grid, let there be a small gap between the outer radius $`R_B`$ of the Cauchy grid and the matching radius $`R_m`$ (which is also the inner radius of the characteristic grid). In 3-dimensional Cartesian-spherical matching the outermost Cauchy grid points and the innermost characteristic grid points are necessarily distinct. This is represented here by a gap $`R_mR_B=\kappa \mathrm{\Delta }r`$, where $`\kappa 0`$ is an arbitrary parameter. The characteristic grid consists of the uniformly spaced points $`\eta _\alpha =\frac{1}{2}+\alpha \mathrm{\Delta }\eta `$ ($`0\alpha N_\eta `$), where $`\mathrm{\Delta }\eta =(2N_\eta )^1`$. The retarded time levels $`u=u_n`$ for the characteristic evolution are chosen so that $`u=u_n`$ intersects the time level $`t=t_n`$ of the Cauchy evolution at the matching radius; therefore, $`u_n=t_nR_m`$ and $`\mathrm{\Delta }u=\mathrm{\Delta }t`$. We denote by $`g_\alpha ^n`$ the value of $`g`$ at $`\eta =\eta _\alpha `$, $`u=u_n`$. The initial characteristic data consist of $`g_\alpha ^0`$, $`0\alpha N_\eta `$.

In the $`n`$th iteration of the evolution, we compute the field values at the grid points with $`u=u_n`$ using the values of $`g_\alpha ^{n1}`$, which are known either from initialization or from the previous iteration. As already described, this is done by the characteristic marching algorithm based on the integral identity (17).

At the inner boundary of the characteristic grid ($`\alpha =1`$), the previous scheme must be slightly modified, since $`g_{\alpha 2}^n`$ is not defined. For this initial step, $`PQRS`$ is chosen so that $`\eta _P=\eta _0`$, $`\eta _Q=\eta _1`$, and $`g_R`$, $`g_S`$ are approximated by quadratic interpolation in terms of $`g_0^{n1}`$, $`g_1^{n1}`$, $`g_2^{n1}`$, which have already been computed. Besides these field values, the final evaluation of $`g_1^n`$ still requires the value of $`g_0^n`$, which is supplied by the matching procedure (step 3).

Step 3. Matching. Numerous schemes are possible in the case of spherically symmetric matching. Here I describe one which is stable for a wide range of gap sizes, $`0\kappa 2`$ and works in the more complicated situation with an interior Cartesian grid and an exterior spherical grid.

The required boundary values $`G_{M+1}^n`$ and $`g_0^n`$ are computed by radial interpolations at constant $`t`$, using the field values at points $`A`$, $`B`$, $`E`$, and $`F`$ in Fig. 4. The first two of these field values are already known at the $`n`$th step, while the last two can be obtained by cubic radial interpolations along the previously evolved characteristics $`u=u_{n1}`$ and $`u=u_{n2}`$, respectively. At the initial step, point $`F`$ lies on the characteristic $`u=u_1=\mathrm{\Delta }tR_m`$, which is not evolved by the algorithm; this field value is supplied along with the initial data. Once the Cauchy and characteristic boundary values are computed, a new iteration may be performed starting from Step 1 above.

Since all the interpolations employed in the matching step have fourth order error, the matching algorithm has the same second order global accuracy exhibited by the separate Cauchy and characteristic algorithms, as confirmed by numerical tests.

In summary, CCM has been implemented and tested for nonlinear waves in 3-dimensional space. No special assumption is made about the waves crossing the computational interface $`r=R_m`$ and nonlinear effects in the exterior characteristic domain are automatically taken into account. In numerical experiments CCM converged to the exact solution (with the matching boundary fixed at an arbitrary position) in highly nonlinear problems. For comparison, nonlocal ABC’s yielded convergent results only in linear problems. In terms of both computational cost and accuracy, CCM is a very effective way to solve the 3-dimensional wave equations, with or without a nonlinear term. It is possible to achieve convergence with an ABC by refining the grid and simultaneously enlarging the radius of the outer boundary. However, this is very expensive computationally, especially for small target error in the determination of the radiated signal in a 3-dimensional problems . Because CCM is convergent under grid-refinement alone, for small target error its performance is significantly better than any available alternative. In strongly nonlinear problems. CCM appears to be the only available method which is able to produce numerical solutions which converge to the exact solution with a computational interface located at an arbitrary fixed position.

## V Application to black holes

By definition, a black hole traces out a world tube in space-time which is the boundary of what is visible to an outside observer. As a result, this worldtube is necessarily a characteristic hypersurface traced out in space-time by light rays (the characteristics of general relativity).

A mechanical analogue of a black hole can in principle be constructed by letting a reservoir of water empty through a funnel into a lower reservoir. By arranging the flow velocity in the funnel to exceed the velocity of sound in the upper reservoir, an acoustic version of a black hole results . Sound waves can travel from the upper reservoir to the lower but not in the reverse direction.

The simplest example of a black hole arises in the spherically symmetric collapse of a star whose energy has been depleted to the extent that internal pressure cannot withstand gravity. A point sized black hole first forms at the center of the collapsing star, which then expands into a sphere growing with the velocity of light, thus tracing out a light cone in space-time. Normally a light cone keeps expanding forever. What makes a black hole light cone unique is that gravity halts this expansion and eventually the black hole just hovers in equilibrium at a fixed size. Each point on the spherical black hole still moves along a light ray, but the sphere of light rays is in a delicate gravitational balance between growing and shrinking.

This spherically symmetric black hole was discovered analytically by Schwarzschild as a simple solution to Einstein’s equations. Unfortunately, Schwarzschild’s black hole does not emit gravitational waves (just as a spherically symmetric charge distribution does not emit electromagnetic waves).

The inspiral and merger of a binary system of black holes is a powerful source of gravitational waves, which lie in the frequency band detectable by the new gravity wave observatories. The lack of symmetry necessitates a computational treatment to determine the waveform of the radiated signal. The binary black holes trace out a characteristic hypersurface in space-time so that their simulation by characteristic evolution is a natural approach. Using a characteristic evolution algorithm similar to that which I have described here, except now applied to Einstein’s equation, we have found fascinating results. The binary black holes, rather than initially forming as a point caustic in the Schwarzschild case, form as a cross-over surface where light rays collide. This cross-over surface is itself bounded by a ring of caustics which in a sense mark the merger of the two individual black holes into a single black hole.

The individual black holes form as spheres but, as illustrated in Fig. 5, as the holes approach, their mutual gravitational tidal distortion produces sharp pincers just prior to merger. At merger, the pincers join to form a single temporarily toroidal black hole, as illustrated in Fig. 6. The inner hole of the torus subsequently closes up to produce first a single peanut shaped black hole and finally a spherical black hole. Details of this merger can be viewed at the web site http://artemis.phyast.pitt.edu/animations. We are now using characteristic evolution to compute the gravitational wave signal emitted in such black hole collisions in the anticipation that they will be detected by the new gravity wave observatories.

## VI Acknowledgements

I thank the Pittsburgh Supercomputing Center and NPACI for making computing time available for this research and the National Science Foundation for research support under grant NSF PHY 9510895.


format: text

text:
# Dynamic Cosmic Strings II: Numerical evolution of excited Cosmic Strings

## I Introduction

This is the second paper in a series devoted to the study of dynamic cosmic strings. In the previous paper , henceforth referred to as paper I, we derived the equations of motion for a time dependent cylindrically symmetric cosmic string coupled to a gravitational field with two degrees of freedom. The treatment involved using a Geroch decomposition to reformulate the problem in terms of matter fields on a $`2+1`$ dimensional spacetime and two geometrical variables $`\nu `$ and $`\tau `$ which describe the gravitational degrees of freedom. Unlike the original cylindrical metric the reduced $`2+1`$ spacetime is asymptotically flat and this allows us to implement a numerical treatment of the field equations which includes null infinity as part of the grid. In paper I we presented a Cauchy-characteristic matching (CCM) code that reproduced the results of several vacuum solutions with both one and two degrees of freedom but did not perform entirely satisfactorily in the presence of matter. We have therefore developed a second implicit, purely characteristic code. The numerical details, the testing and the convergence analysis of this implicit code will be described in this paper. We also give a detailed analysis of the interaction between the string and a pulse of gravitational radiation.

After briefly describing the variables used for the coupled system in section II the details of the numerical scheme are presented in section III. We begin by describing the relaxation method we use for the much simpler problem of a static cosmic string in Minkowski spacetime and demonstrate how this approach leads naturally to the implicit scheme used to solve the general case of a dynamical cosmic string coupled to gravity. The testing of the code is described in section IV. This involves comparing it with exact solutions, checking that static cosmic string initial data remain static when they are evolved, and undertaking a time dependent convergence analysis of the code. We are able to show that the code is accurate, stable and shows clear second order convergence. In section V we analyse the interaction between an initially static cosmic string and a Weber–Wheeler type pulse of gravitational radiation. We show that this interaction causes the string fields $`S`$ and $`P`$ to oscillate and examine how this oscillation depends upon both the width and amplitude of the pulse and also upon the constants $`e`$, $`\lambda `$ and $`\eta `$. The key result is that the frequencies are essentially independent of the nature of the Weber–Wheeler wave but are proportional to the masses $`m_S`$ and $`m_P`$ associated with the scalar field $`S`$ and the vector field $`P`$. The string continues to ring after the gravitational pulse has largely radiated away, but the oscillations slowly decay and the variables eventually return to their static values. It is only possible to observe this effect because of the long term stability of the code. Finally in section VI we discuss our results and outline future work.

## II The variables and field equations

We begin with a brief summary of the variables used to describe the spacetime and the string. In cylindrical polar coordinates $`(t,\rho ,\varphi ,z)`$ we may write the line element for a cylindrically symmetric spacetime with two degrees of freedom as a modified version of that given by Jordan, Ehlers, Kundt and Kompaneets ,

$$ds^2=e^{2(\gamma \psi )}(dt^2d\rho ^2)\rho ^2e^{2\psi }d\varphi ^2e^{2(\psi +\mu )}(\omega d\varphi +dz)^2,$$
(1)
where $`\psi `$, $`\omega `$, $`\mu `$ and $`\gamma `$ are functions of $`t`$, $`\rho `$. As shown in paper I, in order to work with an asymptotically flat spacetime it is convenient to make a Geroch decomposition  in which the 4-dimensional metric is replaced by a 3-dimensional one and two auxiliary scalar fields. These scalar fields are the norm of the axial Killing vector $`\nu `$ and the Geroch potential $`\tau `$ and are related to $`\psi `$ and $`\omega `$ by equations (15) and (16) of paper I, while the conformal 3-dimensional line element is given by

$$d\stackrel{~}{\sigma }^2=e^{2(\gamma +\mu )}(dt^2d\rho ^2)\rho ^2e^{2\mu }d\varphi ^2.$$
(2)
The cosmic string is described by a complex scalar field $`\mathrm{\Phi }`$ coupled to a U(1) gauge field $`A_\mu `$, which in cylindrical symmetry can be written as

$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}S(t,\rho )e^{i\varphi },`$ (3)
$`A_\mu `$ $`=`$ $`{\displaystyle \frac{1}{e}}[P(t,\rho )1]_\mu \varphi .`$ (4)
It is also helpful to introduce rescaled coupling constants and variables given by

$`X`$ $`=`$ $`{\displaystyle \frac{S}{\eta }},`$ (5)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{e^2}{\lambda }},`$ (6)
$`r`$ $`=`$ $`\sqrt{\lambda }\eta \rho ,`$ (7)
$`\stackrel{~}{t}`$ $`=`$ $`\sqrt{\lambda }\eta t.`$ (8)
Here $`\eta `$ is the vacuum expectation value of the scalar field while $`\alpha `$ represents the relative strength of the coupling between the scalar and vector field given by $`e`$, compared to the self-coupling of the scalar field given by $`\lambda `$. Critical coupling, for which the masses of the scalar and vector fields are equal, is given by $`\alpha =8`$ .

The field equations for a cosmic string coupled to gravity were derived in paper I in terms of Cauchy coordinates \[Paper I (26)–(33)\] and also in terms of compactified characteristic coordinates \[Paper I (43)–(49)\]. However the fully characteristic numerical scheme also requires characteristic equations in the inner region. In terms of the retarded time $`u=\stackrel{~}{t}r`$ and the radius $`r`$ the field equations are given by

$`\mathrm{}\nu `$ $`=`$ $`\nu _{,r}\mu _{,r}+{\displaystyle \frac{\tau _{,r}^2\nu _{,r}^2}{\nu }}\nu _{,u}\mu _{,r}\nu _{,r}\mu _{,u}+2{\displaystyle \frac{\nu _{,u}\nu _{,r}\tau _{,u}\tau _{,r}}{\nu }}+8\pi \eta ^2\left[2e^{2(\gamma +\mu )}(X^21)^2+e^{2\mu }\nu ^2{\displaystyle \frac{2P_{,u}P_{,r}P_{,r}^2}{\alpha r^2}}\right],`$ (9)
$`\mathrm{}\tau `$ $`=`$ $`\tau _{,r}\mu _{,r}2{\displaystyle \frac{\tau _{,r}\nu _{,r}}{\nu }}\tau _{,u}\mu _{,r}\tau _{,r}\mu _{,u}+2{\displaystyle \frac{\tau _{,r}\nu _{,u}+\tau _{,u}\nu _{,r}}{\nu }},`$ (10)
$`\mathrm{}\mu `$ $`=`$ $`\mu _{,r}^2+{\displaystyle \frac{\mu _{,r}}{r}}{\displaystyle \frac{\mu _{,u}}{r}}2\mu _{,u}\mu _{,r}+8\pi \eta ^2\left[2{\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}(X^21)^2+e^{2\gamma }{\displaystyle \frac{X^2P^2}{r^2}}\right],`$ (11)
$`0`$ $`=`$ $`2\gamma _{,r}+2r\gamma _{,r}\mu _{,r}r\mu _{,rr}+r\mu _{,r}^2{\displaystyle \frac{r}{2\nu ^2}}(\tau _{,r}^2+\nu _{,r}^2)8\pi \eta ^2\left[rX_{,r}^2+{\displaystyle \frac{1}{\alpha }}e^{2\mu }\nu {\displaystyle \frac{P_{,r}^2}{r}}\right],`$ (12)
$`\mathrm{}P`$ $`=`$ $`2{\displaystyle \frac{P_{,r}}{r}}+2{\displaystyle \frac{P_{,u}}{r}}P_{,r}\mu _{,r}+P_{,r}{\displaystyle \frac{\nu _{,r}}{\nu }}+P_{,r}\mu _{,u}+P_{,u}\mu _{,r}{\displaystyle \frac{P_{,r}\nu _{,u}+P_{,u}\nu _{,r}}{\nu }}\alpha {\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}PX^2,`$ (13)
$`\mathrm{}X`$ $`=`$ $`X_{,r}\mu _{,r}X_{,u}\mu _{,r}X_{,r}\mu _{,u}4{\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}X(X^21)e^{2\gamma }{\displaystyle \frac{XP^2}{r^2}},`$ (14)
where $`\mathrm{}`$ represents the flat-space d’Alembert operator

$$\mathrm{}=2\frac{^2}{ur}\frac{^2}{r^2}\frac{1}{r}\left(\frac{}{r}\frac{}{u}\right).$$
(15)
Note that there are two further Einstein equations which are not used in the numerical scheme as they are a consequence of the above equations and their derivatives. These equations are only used to provide a check on the numerical accuracy of the code.

## III Numerical methods

In order to solve the above field equations we have developed two independent codes. The first is based on the Cauchy-characteristic matching code of Dubal et al.  and d’Inverno et al. . This code performs well in the absence of matter and has been used in paper I to study several cylindrically symmetric vacuum solutions. In paper I we have also given details of the convergence analysis and described the modifications with respect to  that lead to long term stability with both gravitational degrees of freedom present. However the CCM code performed less satisfactorily in the evolution of the cosmic string. This is due to the existence of unphysical solutions to the evolution equations (9)–(14) which diverge exponentially as $`r\mathrm{}`$. Controlling the time evolution near null infinity by means of a bump function enabled us to select the physical solutions with regular behavior at $`I^+`$, but the bump function itself introduced noise which eventually gave rise to instabilities. We therefore implemented a second implicit, purely characteristic, code which allowed us to directly apply boundary conditions at the origin as well as null infinity and thus suppress diverging solutions. It is interesting that this problem is already present in the calculation of the static cosmic string in Minkowski spacetime. We will, therefore, first describe the numerical scheme used in the static Minkowskian case where the equations are fairly simple. We then present the modifications necessary for the static and dynamic case coupled to the gravitational field.

### A The static cosmic string in Minkowski spacetime

In (9)–(14) we set the metric variables to their Minkowskian values and all time derivatives to zero to obtain the equations for the static cosmic string in Minkowski spacetime (cf. )

$`r{\displaystyle \frac{d}{dr}}\left(r^1{\displaystyle \frac{dP}{dr}}\right)`$ $`=`$ $`\alpha X^2P,`$ (16)
$`r{\displaystyle \frac{d}{dr}}\left(r{\displaystyle \frac{dX}{dr}}\right)`$ $`=`$ $`X\left[P^2+4r^2(X^21)\right].`$ (17)
The boundary conditions are

$`P(0)=1,`$ $`\underset{r\mathrm{}}{lim}P(r)=0,`$ (18)
$`X(0)=0,`$ $`\underset{r\mathrm{}}{lim}X(r)=1.`$ (19)
In order to cover the whole spacetime with a finite coordinate range, we divide the computational domain into two regions. In the inner region ($`0r1`$) we use the coordinate $`r`$, while in the outer region we introduce the compactified radius

$`y`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{r}}},`$ (20)
which covers the range $`1y0`$ corresponding to the region $`r1`$. It is also useful to combine $`r`$ and $`y`$ into a single radial variable $`w`$ defined by

$`w=\{r\mathrm{for}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}r132/\sqrt{r}\mathrm{for}r>1\text{,}`$ (21)
so that the region $`0w<3`$ corresponds to $`0r<\mathrm{}`$ and infinity is mapped to $`w=3`$ .

In terms of the variable $`y`$, (16) and (17) take the form

$`y{\displaystyle \frac{d}{dy}}\left(y^5{\displaystyle \frac{dP}{dy}}\right)`$ $`=`$ $`4\alpha X^2P,`$ (22)
$`y{\displaystyle \frac{d}{dy}}\left(y{\displaystyle \frac{dX}{dy}}\right)`$ $`=`$ $`4X\left[P^2+4{\displaystyle \frac{(X^21)}{y^4}}\right].`$ (23)
The number of grid points in each region may differ, but each half-grid is uniform. Thus we use a total of $`N:=n_1+n_2`$ grid points where the points labelled $`n_1`$ and $`n_1+1`$ both correspond to the position $`r=1=y`$. The points $`n_1`$, $`n_1+1`$ form the interface between the two regions (see Figure 1). One point will contain the variables in terms of $`r`$, the other in terms of $`y`$. With the computational grid covering the whole spacetime, we now face a two point boundary value problem. Due to the existence of unphysical solutions diverging at $`y=0`$ shooting methods turned out to be unsuitable for solving this problem. On the other hand numerical relaxation, as described in  for example, allows us to directly control the behavior of $`P`$ and $`X`$ at infinity. The form of equations (16), (17) suggests that in order to write them as a first order system we should introduce the auxiliary variables $`Q=r^1P_{,r}`$ and $`R=rX_{,r}`$. The equations may then be written in the form

$`P_{,r}`$ $`=`$ $`rQ,`$ (24)
$`X_{,r}`$ $`=`$ $`{\displaystyle \frac{R}{r}},`$ (25)
$`Q_{,r}`$ $`=`$ $`\alpha {\displaystyle \frac{PX^2}{r}},`$ (26)
$`R_{,r}`$ $`=`$ $`X\left[{\displaystyle \frac{P^2}{r}}+4r(X^21)\right].`$ (27)
The corresponding equations in the outer region are given by

$`P_{,y}`$ $`=`$ $`2{\displaystyle \frac{Q}{y^5}},`$ (28)
$`X_{,y}`$ $`=`$ $`2{\displaystyle \frac{R}{y}},`$ (29)
$`Q_{,y}`$ $`=`$ $`2\alpha {\displaystyle \frac{X^2P}{y}},`$ (30)
$`R_{,y}`$ $`=`$ $`2X\left({\displaystyle \frac{P^2}{y}}+4{\displaystyle \frac{X^21}{y^5}}\right).`$ (31)
Standard second order centered finite differencing results in $`4(N2)`$ non-linear algebraic equations which are supplemented by the 4 boundary conditions (19) and 4 interface relations

$`P_{n_1+1}`$ $`=`$ $`P_{n_1},`$ (32)
$`X_{n_1+1}`$ $`=`$ $`X_{n_1},`$ (33)
$`Q_{n_1+1}`$ $`=`$ $`Q_{n_1},`$ (34)
$`R_{n_1+1}`$ $`=`$ $`R_{n_1}.`$ (35)
We then start with piecewise linear initial guesses for $`P`$ and $`X`$ (and the corresponding derivatives $`Q`$ and $`R`$) and solve the $`4N`$ algebraic equations iteratively with a Newton–Raphson method. Results for various choices of the string parameter $`\alpha `$ are shown in Figure 6 of paper I.

### B The static cosmic string coupled to gravity

From the numerical point of view, the problem of solving for a static cosmic string coupled to gravity through the Einstein equations is virtually identical to that of a static string in Minkowski spacetime. The only difference is the much higher degree of complexity of the equations due to the appearance of the functions $`\nu `$, $`\tau `$, $`\mu `$ and $`\gamma `$ as extra variables. We do not present the equations here, since they may be derived from the fully dynamic case by setting all time derivatives to zero in the relevant equations. The solution is again obtained using the relaxation method described in the previous section. As our initial guess for the metric variables we use Minkowskian values, and for the string variables $`X`$ and $`P`$ we use the previously calculated values for a Minkowskian string with the same string parameters. The results are shown in Figure 7 of paper I.

### C The dynamic code

In the dynamic case all variables $`\nu `$, $`\tau `$, $`\mu `$, $`\gamma `$, $`P`$ and $`X`$ are functions of $`u,r`$ and we have to solve the system (9)–(14) of partial differential equations (PDEs). In order to control the behavior of the solution at infinity, we need a generalisation for PDEs of the relaxation scheme applied to ordinary differential equations (ODEs). In view of the characteristic feature of the relaxation scheme, namely the simultaneous calculation of new function values at all grid points, this generalisation leads directly to implicit evolution schemes as used for hyperbolic or parabolic PDEs. Therefore, the dynamic code is based on the implicit, second order in space and time Crank–Nicholson scheme (see  for example). For each solution step we consider two spatial slices of the grid-type of Figure 1, labelled $`n`$ and $`n+1`$. We then apply centered finite differencing according to

$`f`$ $`=`$ $`{\displaystyle \frac{f_{k+1}^{n+1}+f_k^{n+1}+f_{k+1}^n+f_k^n}{4}},`$ (36)
$`f_{,r}`$ $`=`$ $`{\displaystyle \frac{f_{k+1}^{n+1}f_k^{n+1}+f_{k+1}^nf_k^n}{2\mathrm{\Delta }r}},`$ (37)
$`f_{,u}`$ $`=`$ $`{\displaystyle \frac{f_{k+1}^{n+1}+f_k^{n+1}f_{k+1}^nf_k^n}{2\mathrm{\Delta }u}},`$ (38)
where $`f`$ represents any of our variables. Assuming that all functions are known on slice $`n`$ we arrive at a large set of algebraic equations for the $`f_i^{n+1}`$, similar to the static case, which needs to be supplemented by boundary conditions and interface relations analogous to (32)–(35). Again the system of algebraic equations is solved iteratively with the Newton–Raphson method. The initial guess for the data on the new slice $`n+1`$ is taken from the previous slice and convergence is typically achieved within three iterations. For this purpose we rewrite the dynamic equations (9)–(14) as a first order system. The equations for the variables $`\nu `$, $`\tau `$ and $`X`$ involve radial derivatives which may be written in terms of the second order operator $`\frac{}{r}(r\frac{}{r})`$ so we introduce the corresponding variables $`N=r\nu _{,r}`$, $`T=r\tau _{,r}`$ and $`R=rX_{,r}`$ \[cf. paper I equations (110)–(115)\]. The equation for $`\mu `$ on the other hand involves $`\frac{}{r}(r^2\frac{}{r})`$, while that for $`P`$ involves $`\frac{}{r}(r^1\frac{}{r})`$. We therefore introduce the corresponding variables $`M=r^2\mu _{,r}`$ and $`Q=r^1P_{,r}`$. Finally the equation for $`\gamma `$ involves only one $`r`$ derivative and may be written in first order form without having to introduce any further quantities. In terms of these variables (9)–(14) become:

$`\nu _{,r}`$ $`=`$ $`{\displaystyle \frac{N}{r}},`$ (39)
$`\tau _{,r}`$ $`=`$ $`{\displaystyle \frac{T}{r}},`$ (40)
$`\mu _{,r}`$ $`=`$ $`{\displaystyle \frac{M}{r^2}},`$ (41)
$`P_{,r}`$ $`=`$ $`rQ,`$ (42)
$`X_{,r}`$ $`=`$ $`{\displaystyle \frac{R}{r}},`$ (43)
$`2N_{,u}`$ $`=`$ $`N_{,r}+{\displaystyle \frac{T^2N^2}{r\nu }}+{\displaystyle \frac{NM}{r^2}}+2{\displaystyle \frac{\nu _{,u}N\tau _{,u}T}{\nu }}\nu _{,u}{\displaystyle \frac{\nu _{,u}M}{r}}N\mu _{,u}`$ (45)
$`+8\pi \eta ^2\left[2e^{2(\gamma +\mu )}r(X^21)^2+{\displaystyle \frac{1}{\alpha }}e^{2\mu }\nu ^2(2P_{,u}QrQ^2)\right],`$
$`2T_{,u}`$ $`=`$ $`T_{,r}2{\displaystyle \frac{TN}{r\nu }}+2{\displaystyle \frac{\tau _{,u}N+\nu _{,u}T}{\nu }}+{\displaystyle \frac{TM}{r^2}}\tau _{,u}{\displaystyle \frac{\tau _{,u}M}{r}}T\mu _{,u},`$ (46)
$`2M_{,u}`$ $`=`$ $`M_{,r}+{\displaystyle \frac{M^2}{r^2}}2\mu _{,u}M2r\mu _{,u}+8\pi \eta ^2\left[e^{2\gamma }X^2P^2+2{\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}r^2(X^21)^2\right],`$ (47)
$`2(r+M)\gamma _{,r}`$ $`=`$ $`M_{,r}2{\displaystyle \frac{M}{r}}{\displaystyle \frac{M^2}{r^2}}+{\displaystyle \frac{T^2+N^2}{2\nu ^2}}+8\pi \eta ^2\left[R^2+{\displaystyle \frac{1}{\alpha }}e^{2\mu }\nu r^2Q^2\right],`$ (48)
$`2Q_{,u}`$ $`=`$ $`Q_{,r}{\displaystyle \frac{QM}{r^2}}+Q\mu _{,u}{\displaystyle \frac{Q\nu _{,u}}{\nu }}{\displaystyle \frac{P_{,u}N}{r^2\nu }}+{\displaystyle \frac{QN}{r\nu }}+{\displaystyle \frac{P_{,u}}{r^2}}+{\displaystyle \frac{P_{,u}M}{r^3}}\alpha {\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}{\displaystyle \frac{PX^2}{r}},`$ (49)
$`2R_{,u}`$ $`=`$ $`R_{,r}X_{,u}{\displaystyle \frac{X_{,u}M}{r}}+{\displaystyle \frac{RM}{r^2}}R\mu _{,u}4{\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}rX(X^21)e^{2\gamma }{\displaystyle \frac{XP^2}{r}}.`$ (50)
The corresponding first order system in the outer region is given by

$`\nu _{,y}`$ $`=`$ $`2{\displaystyle \frac{N}{y}},`$ (51)
$`\tau _{,y}`$ $`=`$ $`2{\displaystyle \frac{T}{y}},`$ (52)
$`\mu _{,y}`$ $`=`$ $`2yM,`$ (53)
$`P_{,y}`$ $`=`$ $`2{\displaystyle \frac{Q}{y^5}},`$ (54)
$`X_{,y}`$ $`=`$ $`2{\displaystyle \frac{R}{y}},`$ (55)
$`2N_{,u}`$ $`=`$ $`{\displaystyle \frac{1}{2}}y^3\left(N_{,y}2yNM2{\displaystyle \frac{T^2N^2}{y\nu }}\right)y^2\nu _{,u}MN\mu _{,u}+2{\displaystyle \frac{\nu _{,u}N\tau _{,u}T}{\nu }}\nu _{,u}`$ (57)
$`+8\pi \eta ^2\left[2e^{2(\gamma +\mu )}{\displaystyle \frac{(X^21)^2}{y^2}}+{\displaystyle \frac{1}{\alpha }}e^{2\mu }\nu ^2\left(2P_{,u}Q{\displaystyle \frac{Q^2}{y^2}}\right)\right],`$
$`2T_{,u}`$ $`=`$ $`{\displaystyle \frac{1}{2}}y^3\left(T_{,y}2yTM+4{\displaystyle \frac{TN}{y\nu }}\right)T\mu _{,u}y^2\tau _{,u}M+2{\displaystyle \frac{\tau _{,u}N+\nu _{,u}T}{\nu }}\tau _{,u},`$ (58)
$`2M_{,u}`$ $`=`$ $`{\displaystyle \frac{1}{2}}y^3(M_{,y}2yM^2)2{\displaystyle \frac{\mu _{,u}}{y^2}}2\mu _{,u}M+8\pi \eta ^2\left[e^{2\gamma }X^2P^2+2{\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}{\displaystyle \frac{(X^21)^2}{y^4}}\right],`$ (59)
$`2(y^2M+1)\gamma _{,y}`$ $`=`$ $`y^2M_{,y}+4yM+2y^3M^2{\displaystyle \frac{N^2+T^2}{y\nu ^2}}16\pi \eta ^2\left({\displaystyle \frac{R^2}{y}}+{\displaystyle \frac{1}{\alpha }}e^{2\mu }\nu {\displaystyle \frac{Q^2}{y^5}}\right),`$ (60)
$`2Q_{,u}`$ $`=`$ $`{\displaystyle \frac{1}{2}}y^3\left(Q_{,y}+2yQM2{\displaystyle \frac{QN}{y\nu }}\right)+y^4P_{,u}y^4{\displaystyle \frac{P_{,u}N}{\nu }}{\displaystyle \frac{\nu _{,u}Q}{\nu }}`$ (62)
$`+y^6P_{,u}M+Q\mu _{,u}\alpha {\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}y^2PX^2,`$
$`2R_{,u}`$ $`=`$ $`{\displaystyle \frac{1}{2}}y^3(R_{,y}2yRM)X_{,u}R\mu _{,u}y^2X_{,u}Me^{2\gamma }y^2XP^24{\displaystyle \frac{e^{2(\gamma +\mu )}}{\nu }}{\displaystyle \frac{X(X^21)}{y^2}}.`$ (63)
In order to solve these equations we must supplement them by appropriate initial and boundary conditions. We start by considering boundary conditions on the axis. In general we find the code is more stable if one imposes boundary conditions on the radial derivatives rather than the variables themselves. For the variables $`\nu `$, $`\tau `$ and $`X`$ we therefore impose the required boundary conditions on the initial data, but in the subsequent evolution we impose the weaker condition that their radial derivatives are finite on the axis. This ensures that the evolution equations propagate the axial conditions given on the initial data. For the variable $`\mu `$ we impose the condition that $`M`$ is zero on the axis which is equivalent to the rather weak condition that $`r^2\mu _r`$ vanishes there. The inverse power of $`r`$ in the definition of $`Q`$ makes it unsuitable to specify the value of this quantity at $`r=0`$ so in this case we work with the variable directly and require that $`P=1`$ on the axis. Finally the variable $`\gamma `$ is given by a purely radial equation, so in this case we must specify the value on the axis which is chosen to be zero to ensure elementary flatness. Therefore at $`r=0`$ we require

$`N`$ $`=`$ $`0,`$ (64)
$`T`$ $`=`$ $`0,`$ (65)
$`M`$ $`=`$ $`0,`$ (66)
$`\gamma `$ $`=`$ $`0,`$ (67)
$`P`$ $`=`$ $`1,`$ (68)
$`R`$ $`=`$ $`0.`$ (69)
For the boundary conditions at null infinity we know that regular solutions of the cylindrical wave equation have radial derivatives that decay faster than $`1/r`$ so that we may take the variables $`N`$, $`T`$ and $`R`$, which satisfy a wave type equation, to vanish at $`y=0`$. The asymptotics of $`\mu `$ are slightly different due to the additional power of $`r`$ in the radial derivative (similar to the spherically symmetric wave equation) but for a regular solution $`\mu _{,y}`$ vanishes at null infinity. The $`P`$ equation does not satisfy a wave type equation due to the inverse power of $`r`$ but has asymptotic behavior given by a modified Bessel function. The physically relevant finite solution has exponential decay so in this case one may impose the condition that $`Q=0`$ at $`y=0`$. Hence we require the solution to satisfy the following boundary conditions at $`y=0`$

$`N`$ $`=`$ $`0,`$ (70)
$`T`$ $`=`$ $`0,`$ (71)
$`\mu _{,y}`$ $`=`$ $`0,`$ (72)
$`Q`$ $`=`$ $`0,`$ (73)
$`R`$ $`=`$ $`0.`$ (74)
These boundary conditions are sufficient to determine the solution of the first order system (39)–(63) while suppressing the unphysical solutions which are singular on the axis or null infinity. Note that $`\gamma `$ is determined by the constraint equation (12), which is a first order ODE, and thus only needs one boundary condition.

## IV Testing the code

In paper I we used the CCM code to reproduce several exact vacuum solutions to an accuracy of about $`10^5`$ with second order convergence. We have also checked the properties of the codes for the static cosmic string in Minkowski spacetime and coupled to gravity. Both clearly showed second order convergence. Here we will focus on testing the implicit dynamic code. We have carried out four different independent tests, namely

* Reproducing the non-rotating vacuum solution of Weber and Wheeler ,
* Reproducing the rotating vacuum solution of Xanthopoulos ,
* Using the results for the static cosmic string (paper I) as initial data and checking that the system stays in its static configuration,
* Convergence analysis for the string hit by a Weber–Wheeler wave.

Two additional tests arise in a natural way from the field equations and the numerical scheme. As described above there are two additional field equations which are algebraic consequences of the other field equations. We have verified that these equations are satisfied to second order accuracy ($`\mathrm{\Delta }r^2`$). Furthermore the numerical scheme calculates the residuals of the algebraic equations to be solved, which have thus been monitored in test runs. They are satisfied to a much higher accuracy (double precision machine accuracy), so the total error is dominated by the truncation error of the second order differencing scheme. Another independent test is the comparison with the explicit CCM code which yields good agreement for as long as the latter remains stable. The four main tests are now described in more detail.

### A The Weber–Wheeler wave

In the first test we evolve the analytic solution given by Weber and Wheeler , which describes a gravitational pulse of the + polarisation mode. This solution has two free parameters, $`a`$ and $`b`$, which can be interpreted as the width and amplitude of the pulse. The equations together with a more detailed discussion have been given in paper I.

We prescribe $`\nu `$ as initial data according to the analytic expressions obtained for $`a=2`$ and $`b=0.5`$ and set the other free variables to zero, while $`\gamma `$ is calculated via quadrature from the constraint equation (12). In Figure 2 we show the deviation of the numerical results from the analytical one for $`N=1920`$ grid points (320 points in the inner, 1600 points in the outer region) and a Courant factor of 0.5 with respect to the inner region. The convergence analysis (see below) shows that this number of points provides sufficient resolution in the outer region while still keeping computation times at a tolerable level. All computations presented in this work have been obtained with these grid parameters, unless stated otherwise. The code stays stable for much longer time intervals than shown in the figure, but does not reveal any further interesting features as the analytic solution approaches its Minkowskian values and the error goes to zero.

### B The rotating solution of Xanthopoulos

Xanthopoulos  derived an analytic vacuum solution for Einstein’s field equations in cylindrical symmetry containing both the + and $`\times `$ polarisation mode. Its analytic form in terms of our metric variables and a more detailed discussion has been given in paper I. The solution has one free parameter $`a`$ which is set to one in this calculation.

The error of our numerical results is displayed in Figure 3, where we have used the same grid parameters as in the Weber–Wheeler case. Again we have run the code for longer times and found that the error approaches zero. We conclude that the code reproduces both analytic vacuum solutions with excellent accuracy comparable to that of the CCM code and exhibits long term stability.

### C Evolution of the static cosmic string

The tests described above only involve vacuum solutions, so the matter part of the code and the interaction between matter and geometry has not been tested. An obvious test involving matter and geometry is to use the result for the static cosmic string in curved spacetime as initial data and evolve this scenario. All variables should, of course, remain at their initial values. We have evolved the static string data for our standard grid and the parameter set, $`\alpha =1`$ and $`\eta =0.2`$, which corresponds to a strong back-reaction of the string on the metric. The results are shown in Figure 4. The system stays in its static configuration with high accuracy over a long time interval.

### D Convergence analysis

Our investigation of the interaction between the cosmic string and gravitational waves will focus on the string being hit by a wave of the Weber–Wheeler type. In order to check this scenario for convergence we have run the code for the parameter set $`\eta =0.2`$, $`\alpha =1`$, $`a=2`$, $`b=0.5`$ for different grid resolutions, where $`a`$ and $`b`$ are again the width and amplitude of the Weber–Wheeler wave. In our case it is of particular interest to investigate the time dependence of the convergence to see what resolution we need in order to obtain reliable results for long runs. We calculate the convergence rate in the same way as in the static case (cf. paper I), but this time the $`\mathrm{}_2`$-norm is a function of time. So for each variable we have

$`\mathrm{}_2[\mathrm{\Delta }\mathrm{\Psi }^N](u)`$ $`=`$ $`\sqrt{{\displaystyle \frac{[\mathrm{\Psi }_k^N(u)\mathrm{\Psi }_k^{4320}(u)]^2}{N}}},`$ (75)
where the upper label “4320” indicates that the high resolution reference solution has been calculated for $`N=4320`$ grid points. In Figure 5 we show the convergence factor $`\mathrm{}_2[\mathrm{\Psi }^{1920}]/\mathrm{}_2[\mathrm{\Psi }^{2880}]`$ as a function of $`u`$ for $`\nu `$, $`\mu `$, $`\gamma `$, $`P`$ and $`X`$. The initial data for $`\tau `$ is identically zero for this scenario and stays zero during the evolution. The number of grid points is increased by a factor of 1.5 here (instead of the more commonly used 2) to reduce the computation time. Only points common to all grids have been used in the sum in equation (75). For second order convergence we would expect a convergence factor of $`1.5^2`$. Although the results in Figure 5 show weak variations with $`u`$, second order convergence is clearly maintained for long runs. In each case the outer region contains 5 times as many grid points as the inner region (e.g. $`n_1=320`$, $`n_2=1600`$ for the $`N=1920`$ case). The reason for this is that in the dynamical evolutions $`X`$ and especially $`P`$ exhibit significant spatial variations out to large radii. Due to the compactification, the spatial resolution of our grid decreases as we move towards null infinity and to resolve the spatial variations of the string variables out to sufficiently large radii we therefore have to introduce a large number of grid points in the outer region. No such problems occur in the inner region. If significantly fewer grid points are used in the outer region for this analysis, the convergence properties of the string variables can deteriorate to roughly first order level.

## V Time dependence of the string variables

### A Static cosmic strings excited by gravitational waves

The scenario we are going to investigate in this section is an initially static cosmic string hit by a gravitational wave of Weber–Wheeler type. For this purpose we use the static results with two modifications as initial data. Firstly the static metric function $`\nu _0`$ is multiplied by the exact Weber–Wheeler solution to simulate the gravitational wave pulse. Thus we guarantee that initially the cosmic string is indeed in an equilibrium configuration provided the wave pulse is located sufficiently far away from the origin and its interaction with the string is negligible. Ideally the numerical calculation would start with the incoming wave located at past null infinity. In order to approximate this scenario, we found it was sufficient to use the large negative initial time $`u_0=20`$. The second modification is to calculate $`\gamma `$ from the constraint equation (12) to preserve consistency with the Einstein field equations. In Figure 6 the corresponding initial data for $`\nu `$, $`P`$ and $`X`$ are shown for the parameter set $`\eta =10^3`$, $`\alpha =1`$, $`a=2`$ and $`b=0.5`$. From now on we will refer to these values as “standard parameters” and only specify parameters if they take on non-standard values. Note that $`\tau `$ vanishes on the initial slice in this case and stays identically zero throughout the evolution. The case of rotating gravitational waves hitting a cosmic string will be analysed in a future publication. The time evolution of the “standard configuration” is shown in Figure 7 where we plot $`\nu `$, $`\mu `$, $`\gamma `$, $`P`$ and $`X`$ as functions of $`w`$ at times $`20`$, $`10`$, $`0`$, $`2`$ and $`10`$. While the wave pulse is still far away from the origin, its interaction with the cosmic string is negligible (dotted lines). When it reaches the core region, however, it excites the cosmic string and the scalar and vector field start oscillating (dashed curves). After being reflected at the origin, the wave pulse travels along the outgoing characteristics and the metric variables $`\nu `$, $`\mu `$ and $`\gamma `$ quickly settle down into their static configuration which is close to Minkowskian values for $`\eta =10^3`$. The vector and scalar field of the cosmic string, on the other hand, continue ringing albeit with a different character. Whereas the oscillations of the scalar field $`X`$ are dominant in the range $`r2`$ and have significantly decayed at $`u=10`$ as shown in the figure, the vector field oscillations propagate to large radii and fall off very slowly (solid curves). This behavior is also illustrated in the right panel of Figure 8 which shows a contour plot of $`P`$ as a function of $`(u,r)`$ out to $`r=50`$. We shall see below, that the oscillations of $`P`$ will also decay eventually and the cosmic string will asymptotically settle back into its equilibrium configuration.

### B Frequency analysis

We will now quantitatively analyse the oscillations of the scalar and vector field of the cosmic string. Since we are working in rescaled coordinates, time and distance are measured in units of $`1/\sqrt{\lambda }\eta `$ and frequency in its inverse. To avoid complicated notation, however, we will omit the units from now on unless the meaning is unclear. In order to measure frequencies, we Fourier analyse the time evolution of the scalar and vector field for fixed radius $`r`$. Figure 9 shows $`P`$ and $`X`$ for standard parameters as functions of time at $`r=1`$ together with the corresponding power spectra. The plots for $`X`$ show a characteristic frequency $`f_X=0.43`$ whereas for $`P`$ we find a strong peak at $`f_P=0.16`$. The spectrum for $`P`$, however, also shows a strong mode at $`f=0.43`$. We have calculated similar power spectra for a large class of parameter sets and discovered this effect on numerous occasions — in addition to a strong peak at the characteristic frequency of $`P`$ or $`X`$ there is a second maximum at the characteristic frequency of the other field. In general the characteristic mode of $`X`$ resulted in stronger peaks at smaller radius, that of $`P`$ was stronger at larger radii. We attribute this feature to the interaction between the scalar and vector component of the cosmic string. The variation of the relative strength of the oscillations with radius confirms the corresponding observation in Figure 7. The accuracy of the measurements of the frequencies is limited by the resolution of the Fourier spectra which again is limited by the time interval covered in the evolution and, thus, by computation time. For tolerable computation times we get an accuracy $`\mathrm{\Delta }f0.01`$ which corresponds approximately to one bin in the frequency spectra.

In order to investigate the dependency of the oscillations on $`\alpha `$, $`\eta `$, $`a`$, $`b`$ and the radial position $`r`$, we have varied each parameter over at least two orders of magnitude while keeping the other parameters at standard values. We have found the following dependencies:

* The frequencies of both $`X`$ and $`P`$ did not show any variations with $`\eta `$ for $`\eta <0.1`$. (Note that $`\eta `$ does, however, appear in the units). For larger values of $`\eta `$, the non-linear interaction between string and geometry becomes dominant and we did not find a simple relation between frequency maxima and parameters.
* The variation of the parameters $`a`$ and $`b`$, the width and amplitude of the Weber–Wheeler pulse, has no measurable effect on the frequencies of $`X`$ and $`P`$, but only determined the amplitude of the oscillations. A narrow, strong pulse leads to larger amplitudes.
* For small $`r`$ the oscillations in $`X`$ are stronger, whereas those for $`P`$ dominate at large $`r`$. The frequency values, however, do not depend on the radius. For radii greater than 10 the oscillations in $`X`$ had decayed so strongly that we could no longer measure its frequency.
* To first order approximation the frequency of the scalar field, $`f_X`$, is independent of $`\alpha `$ over the observed range $`0.1\alpha 10`$. $`f_P`$, however, shows a strong dependence on $`\alpha `$ which is also illustrated in Figure 8, where we compare contour plots of $`P`$ obtained for $`\alpha =0.2`$ and $`1`$. The frequency is significantly larger for $`\alpha =1`$.

In Figure 10 we show $`f_P`$ and $`f_X`$ as functions of $`\alpha `$ together with the power laws obtained from a linear regression of the corresponding double logarithmic data. In the range $`5\alpha 8`$ (where $`\alpha =8`$ corresponds to the critical coupling), the expected values of $`f_P`$ and $`f_X`$ become similar and we observed only one maximum in the corresponding power spectra. Therefore a classification with respect to the vector or scalar origin of the frequencies is not obvious. These values (shown by filled lozenges in Figure 10) have not been used in the regression analysis. We obtain power law indices $`\sigma _X=0.00`$ and $`\sigma _P=0.50`$, so that

$`f_X`$ $``$ $`\mathrm{const}.`$ (76)
$`f_P`$ $``$ $`\sqrt{\alpha }.`$ (77)
Whereas the fitted curve for $`f_P`$ coincides with the observed values to high accuracy, Figure 10 shows a small but significant deviation of the measured $`f_X`$ from the fitted constant line. Indeed the sizable range over $`\alpha `$ for which we observe only one frequency indicates significant non-linear interaction similar to the effect of phase locking in non-linear systems of ODEs . If we consider $`f_X`$ and $`f_P`$ to be given by (76) and (77) in terms of the rescaled unphysical coordinates and we transform this back into physical units using $`\alpha =e^2/\lambda `$, we arrive at the following relations for the physical variables

$`f_X`$ $``$ $`\sqrt{\lambda }\eta ,`$ (78)
$`f_P`$ $``$ $`e\eta .`$ (79)
As shown in  up to constant factors $`\sqrt{\lambda }\eta `$ and $`e\eta `$ are the masses of the scalar and the vector field, $`m_X`$ and $`m_P`$, and we conclude that $`X`$ and $`P`$ have characteristic frequencies

$`f_X`$ $``$ $`m_X,`$ (80)
$`f_P`$ $``$ $`m_P.`$ (81)
Since the frequencies for $`X`$ and $`P`$ seem only to depend upon the respective masses we have attempted to confirm these results by considering the oscillations of a cosmic string in two further scenarios. Firstly since the frequencies do not depend upon the Weber–Wheeler pulse we take as initial data the static values for the metric variables but excite the string by adding a Gaussian perturbation to either the $`X`$ or $`P`$ static initial values. The evolution is then computed using the fully coupled system. Secondly since the frequencies do not seem to depend upon the strength of the coupling to the gravitational field we have completely decoupled the gravitational field and considered the evolution of a cosmic string in Minkowski spacetime. The initial data is taken to be that for a static string in Minkowski spacetime with a Gaussian perturbation to either the $`X`$ or $`P`$ values. The evolution is then computed using the equations for a dynamical string in a Minkowskian background \[equations (85) and (86) of paper I\]. In both cases we find the same frequencies, to within an amount $`\mathrm{\Delta }f=0.01`$, that we found in the original case of the fully coupled system excited by a Weber–Wheeler pulse. Furthermore the frequencies did not depend on the location or shape of the field perturbation nor upon the choice of $`X`$ or $`P`$ as the perturbed field.

### C The long term behavior of the dynamic string

The time evolution shown in Figures 7 and 8 indicate that the oscillations of the cosmic string excited by gravitational waves gradually decay and metric and string settle down into an equilibrium state. We have calculated a very long run ($`20u410`$) to investigate the long term behavior in detail. The unphysically large value of $`\eta =0.1`$ is chosen for this calculation in order to guarantee a strong interaction between spacetime geometry and the cosmic string. In Figure 11 we show the difference $`\mathrm{\Delta }f:=f_{\mathrm{evol}}f_{\mathrm{stat}}`$ between the time-dependent $`\nu `$, $`\mu `$, $`\gamma `$ and $`X`$ and their corresponding static results obtained for the same parameters. For the vector field $`P`$ a similar 3-dimensional plot would require an extreme resolution to properly display the oscillations of the vector field (cf. Figure 7). For this reason we proceed differently and calculate the $`\mathrm{}_2`$-norm and the maximum of $`\mathrm{\Delta }P`$ for each slice $`u=\mathrm{const}`$. Both functions are plotted versus time in Figure 11. The incoming wave pulse can clearly be seen as a strong deviation of $`\nu `$ from the static function. The pulse excites the cosmic string and is reflected at the origin at $`u=0`$. The metric variables and the scalar field $`X`$ then quickly reach their equilibrium values. The oscillations in $`P`$ decay on a significantly longer time scale which is also evident in Figures 7 and 8 and the $`\mathrm{}_2`$-norm of $`\mathrm{\Delta }P`$ slowly approaches zero. Significantly longer runs than shown here are prohibited by the required computational time, but the results indicate that $`P`$ will also eventually reach its equilibrium configuration.

## VI Conclusion

In this paper we have described the details of the implicit, fully characteristic, numerical scheme which is used to solve the field equations for a cosmic string coupled to gravity. A feature of the cosmic string equations is that they admit exponentially diverging unphysical solutions. By using a Geroch decomposition it is possible to reformulate the problem in terms of fields which describe the string on an asymptotically flat $`2+1`$-dimensional metric as well as two auxiliary fields $`\nu `$ and $`\tau `$ which describe the gravitational degrees of freedom. We can then introduce a conformally compactified radial coordinate $`y`$ which allows us to include null infinity as part of the numerical grid. As well as avoiding the need to introduce artificial outgoing boundary conditions at the edge of the grid this approach has the advantage that we can also enforce boundary conditions for the string variables at null infinity which rule out the unphysical solutions. The use of the geometrically defined Geroch variables also improves the long term stability of the code compared to the use of metric variables.

The code has been shown to reproduce the results of two exact vacuum solutions, the Weber–Wheeler solution which describes a pulse of gravitational radiation with just the $`+`$ polarisation state, and a solution due to Xanthopoulos which describes a gravitational wave with both the $`+`$ and $`\times `$ polarisation states. The code has also been shown to reproduce the results of the static cosmic string code in that initial data corresponding to a static solution do not change when evolved in time using the dynamical code. For both the exact vacuum solutions and the static initial data the code shows excellent long term stability. Finally a time dependent convergence analysis demonstrates clear second order convergence of the code.

After demonstrating the reliability of the code we use it to analyse the interaction between an initially static cosmic string and a Weber–Wheeler type pulse of gravitational radiation. We find that the gravitational wave excites the string and causes the string variables $`X`$ and $`P`$ to oscillate before the configuration slowly settles back into its equilibrium state. In terms of the unphysical rescaled variables we find that the frequencies of the oscillations are essentially independent of the value of the coupling constant $`\eta `$ and of the width and amplitude of the Weber–Wheeler pulse. We also find that the frequency of $`X`$ is independent of $`\alpha `$ while that of $`P`$ is proportional to $`\sqrt{\alpha }`$. When this result is translated back into the physical units we find that the frequency of the scalar field is proportional to the mass of the scalar field and the frequency of the vector field is proportional to the mass of the vector field. This result is confirmed by investigating two further scenarios. Firstly we consider the evolution of static initial data for the string coupled to the gravitational field, but with a Gaussian perturbation to one of the string variables, and secondly we consider the same thing but in a Minkowskian background with the gravitational field decoupled. In both cases we obtain the same relationship between the frequencies and the mass.

Having investigated the interaction between a Weber–Wheeler type pulse of gravitational radiation and the cosmic string, the next obvious step is to consider the interaction between the string and a pulse of gravitational radiation with both polarisation states present. The code will also be a valuable tool in comparing the numerical results with those that one gets from using perturbation theory or the thin string limit to approximate the behavior of a cosmic string coupled to gravity.

###### Acknowledgements.

We would like to thank Ray d’Inverno for helpful discussions and Denis Pollney for help with GRTensor II.