Apples with Apples

Welcome
welcome

Meetings

Methods
fulcrum

Results
bhspin

Legal Notice

Give us your
Feedback

Polarized Gowdy wave testbed

All of the tests described so far considered initial data which were perturbations of a flat background. Here we use a genuinely curved exact solution - a polarized Gowdy spacetime - to test codes in a strong field context. The polarized Gowdy spacetimes are solutions of the vacuum Einstein equations on the 3-torus, and describe an expanding universe containing plane polarized gravitational waves [Gowdy(1971)]. Gowdy spacetimes have previously been used for testing numerical relativity codes by a number of authors [van Putten(1997),New et al.(1998)New, Watt, Misner, and Centrella,Garfinkle(2002)]. They have also extensively been studied in mathematical cosmology; see e.g. [Ringstrom(2003)] for latest results. An extensive analytical and numerical study of Gowdy spacetimes has been carried out by Berger [Berger(2002)].

The polarized Gowdy metric is usually written as

$\displaystyle {ds}^2 = t^{-1/2} e^{\lambda/2} (-{dt}^2 + {dz}^2) + t\,{dw}^2,$

(1)

where

$\displaystyle {dw}^2=e^P {dx}^2 + e^{-P}{dy}^2.$

(2)

Here the time coordinate

is chosen such that time increases as the universe expands. For code testing, it is quite interesting to compare collapsing and expanding situations. We will thus carry out our tests in both time directions. The quantities $\lambda$ and

are functions of

and

only and are periodic in

. The metric is singular at

which corresponds to the cosmological singularity.

With the metric (1), the Einstein evolution equations can be reduced to a single linear equation for [Gowdy(1971)]:

$\displaystyle P_{,tt} + t^{-1} \, P_{,t} - P_{,zz} = 0.$

(3)

The constraint equations become

$\displaystyle \lambda_{,t} = t\, (P_{,t}^2 + P_{,z}^2)$

(4)

and

$\displaystyle \lambda_{,z} = 2 \, t \, P_{,z} \, P_{,t},$

(5)

and correspond to the Hamiltonian (4) and momentum (5) constraints. The general solution to eq. (3) is a sum of terms of the form $\alpha \log t + \beta$ , where $\alpha$ and $\beta$ are real constants, and terms of the form $Z_0(2\pi nt)\cos(2\pi nz)$ and $Z_0(2\pi nt)\sin(2\pi nz)$ , where

is an integer (assuming periodicity of 1 in

) and

is a linear combination of the Bessel functions

and

. We follow [New et al.(1998)New, Watt, Misner, and Centrella] in the choice of the particular solution and set

$\displaystyle P = J_0(2\pi t)\cos(2\pi z),$

(6)

which yields

$\displaystyle g_{xx}=te^P,\ g_{yy}=te^{-P},\ g_{zz}=t^{-1/2}e^{\lambda/2},$

(7)

$\displaystyle K_{xx}$	$\displaystyle =$	$\displaystyle - \frac{1}{2} t^{1/4} e^{-\lambda/4}e^P(1+t P_{,t}),$
$\displaystyle K_{yy}$	$\displaystyle =$	$\displaystyle - \frac{1}{2} t^{1/4} e^{-\lambda/4}e^{-P}(1 - t P_{,t}),$	(8)
$\displaystyle K_{zz}$	$\displaystyle =$	$\displaystyle \frac{1}{4} t^{-1/4} e^{\lambda/4}(t^{-1}-\lambda_{,t}).$

The shift vanishes, and the lapse is given as

$\displaystyle \alpha=\sqrt{g_{zz}} = t^{-1/4}e^{\lambda/4}.$

(9)

Using our choice for

(6), the constraint eqs. (4,5) yield an expression for $\lambda$ :

$\begin{displaymath}\begin{array}[c]{r c l} \lambda &=& -2\pi t J_{0}(2\pi t) J_{... ...\pi )\bigr]- 2\pi J_{0}(2\pi ) J_{1}(2\pi )\bigr\}. \end{array}\end{displaymath}$

(10)

Note that $A = \int_z t P_{,t} \, dz$ is a constant of motion and can be used to monitor the accuracy of a code. In our case

is set to zero by the choice of initial data.

**Figure 1:** The quantities and $\lambda$ appearing in the Gowdy metric (1) are plotted versus and .
$\includegraphics[width=20pc]{GowdyP.eps}$ $\includegraphics[width=20pc]{GowdyLambda.eps}$

**Figure 2:** The extrinsic curvature components $K_{xx}$ and $K_{zz}$ , as given in Eqs. (8), and at the coordinate origin are plotted versus time .
$\includegraphics[width=30pc]{GowdyK.eps}$

Figures showing , $\lambda$ and extrinsic curvature components, constructed from the analytic formulas, are given in Figures 1 and 2. While slowly decays to zero, $\lambda$ shows a secular linear growth due to the cosmological expansion, and both and $\lambda$ exhibit gravitational wave oscillations. Note that although the individual extrinsic curvature components do not exhibit a fixed sign, is negative and decays in absolute value, consistent with the cosmological expansion. The linear growth of $\lambda$ leads to exponential growth in the metric component $g_{zz}$ . This makes an evolution with standard ADM variables much harder than evolving the Gowdy quantity .

The (coordinate) velocity of light is constant in the coordinates chosen in eq. (1), and with a fixed spatial discretization size $\Delta z$ the Courant condition is consistent with a fixed timestep $\Delta t$ . This makes it convenient to choose the gauge (1) for evolving in the expanding direction. We will see below however, that this leads to exponential growth in the metric component $g_{zz}$ . For the collapsing direction, this would lead to a singularity at , so we will evolve this case with a different slicing as discussed below.

For the forward (expanding) evolution, we set initial data from the exact solution at , which yields initial data of order unity, and evolve with any lapse condition which is equivalent in the continuum limit to the exact lapse given by eq. (9). Due to the exponential growth in the metric variables, such evolutions may crash rather soon but will test the accuracy of a code in a rather harsh situation. In order to evolve in the backward time direction, we choose harmonic time slicing, as has previously been done by Garfinkle [Garfinkle(2002)]. Since harmonic slicing is marginally singularity avoiding [Bona et al.(1997)Bona, Massó, Seidel, and Stela,Alcubierre(2003)], such evolutions should only asymptotically reach the singularity at $t=-\infty$ .

It turns out that it is actually quite simple to write down an exact solution for harmonic slicing, which greatly simplifies the task of choosing appropriate gauge source functions for various formalisms. Starting with the Gowdy metric, as given by eq. (1), we look for a coordinate transformation $(t,x^i) \rightarrow (\tau, x^i)$ , with $t = F(\tau)$ . In the new coordinates, the lapse becomes $\hat \alpha = F(\tau)^{-1/4} \partial_\tau F(\tau) e^{\lambda/4}$ . The harmonicity condition $\Box \, t = 0$ then implies

$\displaystyle e^{-\lambda/2} \left[ F \partial_{\tau\tau} F - \partial_\tau F^2 \right] = \sqrt{ F} \partial_\tau F^3$

(11)

with the solution $F(\tau) = k e^{c \tau}$ , where

and

are free constants. The lapse in this new gauge is

$\displaystyle \hat \alpha(\tau) = c k^{3/4} e^{3 c \tau / 4 + \lambda(F(\tau),z) / 4}.$

(12)

In order to start the collapse slowly, and to simplify initial data, we choose the constants

in such a way that $\hat \alpha = 1$ at the initial time

. Picking a value

for which $J_0(2 \pi t_0) = 0$ , eq. (10) implies that $\hat \alpha$ is independent of

. Using

$\displaystyle \tau_0 = \frac{1}{c} \ln \left( \frac{t_0}{k} \right), \quad \lambda(k e^{c \tau_0}, z) = \lambda_0$

we obtain

$\displaystyle \hat \alpha_0 = c\; t_0^{3/4}\; e^{\lambda_0 / 4}.$

(13)

Given our requirement $\hat \alpha_0 = 1$ , and choosing $t_0 = \tau_0$ , i.e. $F(\tau_0) = \tau_0$ , we get

$\displaystyle c = t_0^{-3/4}\; e^{-\lambda_0 / 4}, \quad k = t_0 e^{-c t_0}.$

(14)

We will choose a particular value of

such that the initial slice is far from the cosmological singularity, but not so far that we have to deal with extremely large numbers. We pick the $20^{\text{th}}$ zero of the Bessel function $J_0(2 \pi t_0)$ , which yields $t_0 \sim 9.8753205829098$ , corresponding to

$\displaystyle c \sim 0.0021195119214617, \quad k \sim 9.6707698127638.$

The values of the metric components found from (1) at

are then $g_{xx} = g_{yy} = t_0$ , $g_{zz} \sim 2.283 \times 10^3$ . This choice challenges a numerical code to accurately track a small effect (the dynamics in $g_{xx}$ , $g_{yy}$ ) together with a larger effect (dynamics in $g_{zz}$ ). Other choices are of course possible, and certainly worth exploring. For the purpose of a standard testbed, which should provide tests which are able to discriminate well between different formulations, the current choice seems appropriate.

The geometry of the grid is chosen analogous to the 1D gauge wave test:

Simulation domain: $z \in [-0.5; +0.5]$ ,
Grid: $z^i = -0.5 + (n-\frac{1}{2}) dz, \quad n=1\ldots 50\rho, \quad dz = 1/(50\rho), \quad \rho \in \mathbb{Z}$
Time step: $dt = dz/4 = 0.005 / \rho$
Run time: , i.e. 1000 crossing times or until code crash.

We output the $L_{\infty}$ and -norms, the maxima and minima, and profiles along the -axis through the center of the grid of $g_{zz}$ , $\alpha$ , , the Hamiltonian constraint and all other nontrivial constraints of the formulation, and some typical evolution variables, depending on the evolution system chosen. We output norms every crossing time, and profiles either every 10 crossing times or once per crossing time for some initial time, depending on the behavior of the solution. We also calculate the $L_{\infty}$ -norms of the difference from the exact solution for $g_{xx}$ and $g_{zz}$ for the expanding direction.

As a sample result we present a comparison of an ADM and a BSSN code for the collapsing direction in Figure 3. While the ADM code shows roughly second order convergence for 1000 crossing times (we show the first 500 for better comparison with the BSSN results), only the lowest resolution BSSN run lasts for 1000 crossing times with the higher resolution runs crashing significantly earlier. The loss of convergence is clear in Figure 3. The poor performance of the BSSN code seems to be rooted in its mixing of components. The comparatively good performance of the ADM code supports the usefulness of this test. Alternative choices of initial data can be made to yield tests with different characteristics, but will not be included in this round of tests.