9512.net
甜梦文库
当前位置:首页 >> >>

A strongly hyperbolic and regular reduction of Einstein's equations for axisymmetric spacet


DAMTP-2005-17

A strongly hyperbolic and regular reduction of Einstein’s equations for axisymmetric spacetimes
Oliver Rinne? and John M Stewart?
Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, United Kingdom (Dated: 9 February 2005) This paper is concerned exclusively with axisymmetric spacetimes. We want to develop reductions of Einstein’s equations which are suitable for numerical evolutions. We ?rst make a Kaluza-Klein type dimensional reduction followed by an ADM reduction on the Lorentzian 3-space, the (2+1)+1 formalism. We include also the Z4 extension of Einstein’s equations adapted to this formalism. Our gauge choice is based on a generalized harmonic gauge condition. We consider vacuum and perfect ?uid sources. We use these ingredients to construct a strongly hyperbolic ?rst-order evolution system and exhibit its characteristic structure. This enables us to construct constraint-preserving stable outer boundary conditions. We use cylindrical polar coordinates and so we provide a careful discussion of the coordinate singularity on axis. By choosing our dependent variables appropriately we are able to produce an evolution system in which each and every term is manifestly regular on axis.
PACS numbers: 04.25.Dm, 04.20.Ex, 04.20.Cv Keywords: Numerical relativity, axisymmetric spacetimes, regularity conditions, (2+1)+1 formalism, Z4 system, constraint-preserving boundary conditions, initial boundary value problem

arXiv:gr-qc/0502037v1 9 Feb 2005

I.

INTRODUCTION

In the early days of numerical relativity attention focussed on spherically symmetric spacetimes. The next simplest subject would have been axisymmetric spacetimes, but very little progress was made in this area because, in polar coordinates adapted to the symmetry, there is a coordinate singularity on axis. Many attempts to deal with this proved unsuccessful, most evolutions became unstable and attention quickly turned, with the rapidly increasing capacity and speed of hardware, to spacetimes without symmetries. For a comprehensive review see e.g., [13]. However axisymmetric situations arise frequently in astrophysics and deserve further study. The ?rst question is, obviously, why were the evolutions unstable? Suppose we assume that in a neighbourhood of the axis the spacetime is regular so that we can introduce locally cartesian coordinates (t, z, x, y ) with the axis given by x = y = 0. Then the Killing vector is ξ = ?y?/?x + x?/?y , and we say that M is axisymmetric if Lξ M = 0. We say that a tensor ?eld M is regular on axis if its components with respect to this chart can be developed as a Taylor series in x and y in a neighbourhood of the axis1 . We may introduce polar coordinates r = x2 + y 2 and ? = arctan(y/x), so that the Killing vector is ξ = (?/??). Because the transformation between cartesian coordinates and polar coordinates (t, z, r, ?) is singular on axis2 , axisymmetric tensor ?elds which are regular on axis may take strange forms when expressed in polar coordinates. For example if Mαβ is symmetric and has these properties then ? ? A B rD r2 F ? ? B C rE r2 G ?, (1) Mαβ = ? 2 ? ? rD rE H + r J r3 K 2 2 3 2 2 r F r G r K r H ?r J where A, B, . . . , K are functions of t, z and r2 . (See appendix A for this and related results.) In particular both the metric and its ?rst time-derivative must take this form. If the evolution scheme fails to preserve precisely the indicated r-dependencies then M becomes irregular on axis and instability is inevitable.

? Electronic ? Electronic 1 2

address: o.rinne@damtp.cam.ac.uk address: j.m.stewart@damtp.cam.ac.uk For numerical purposes analytic and regular are synonymous. Note that the cartesian expression for the Killing vector ξ vanishes on axis, whereas the polar form does not.

2 Clearly the problems lie mainly in the last row and column of (1). In section II we outline a reduction of Einstein’s equations in spacetime M to a 3-dimensional Lorentzian manifold N with chart (t, z, r), e?ectively removing the last row/column. (The relation between Mrr and M?? still needs careful treatment.) Given a solution of the reduced equations in N one can recover the solution in M. However there is no need to do this. All axisymmetric quantities of physical interest in M have their counterparts in N , and so one might as well work in N . Because that manifold is Lorentzian we can perform an ADM decomposition. The details of this (2 + 1) + 1 scheme, with an arbitrary axisymmetric matter source, are given in section II. In the (2 + 1) + 1 system we have a spatial foliation with lapse function α, shift vector η A , and spatial 2-metric HAB . (Here A, B , . . . range over r and z .) Since every 2-metric is conformally ?at it is tempting to set Hrz = 0, Hrr = Hzz = ψ 4 say. The constraint equations imply elliptic equations for ψ and the shift vector components, and a quasi-maximal slicing condition supplies an elliptic equation for α. Although elliptic equations can become computationally expensive, multigrid techniques, if applicable, can reduce signi?cantly the cost. The same choice was made by Choptuik et al [5]. Like those authors we encountered several di?culties, the most severe being that in near-critical collapse situations the discretized version of the energy constraint equation for ψ failed to remain diagonally dominant so that the multigrid algorithm failed to converge. Gar?nkle et al [10] used a conjugate gradient method which gets round this particular problem. If instead we used an evolution equation for ψ there appeared to be a serious violation of the energy constraint. We decided therefore to look at the other extreme, to minimize the number of elliptic equations solved. There are many, apparently di?erent, ways of doing this, but we chose to implement the Z4 extension of Einstein’s equation suggested by Bona et al [2], [3], but adapted to the (2+1)+1 scheme. The Z(2+1)+1 scheme is outlined in section III. There the energy and momentum constraints as well as the Geroch constraint (which arises in the (2+1)+1 decomposition) are converted into evolution equations for the new variables. Pursuing this philosophy we would like to impose dynamical gauge conditions, which are discussed in section IV. These are generalizations of the harmonic gauge condition. However cylindrical polar coordinates are not harmonic, and this means that we can only do this for the lapse function. For simplicity we assume that the shift vector vanishes. We now have a system of pure evolution equations and it is highly desirable to cast them into ?rst order form. This is carried out in section V where we express them in conservation form with sources. We show in section VI that this system is strongly hyperbolic and we write down the conversions between √ conserved and characteristic variables. The characteristic speeds (divided by the lapse function α) are 0, ±1 and ± f where f is a constant parameter in the lapse evolution equation. (The choice f = 1 gives a harmonic time coordinate and a symmetric hyperbolic system.) Knowing the characteristic structure enables us to formulate constraint-preserving boundary conditions for the initial-boundary value problem in section VII. In section VIII we return to our starting point, regularity on axis. By a judicious choice of new dependent variables (which does not a?ect the hyperbolicity) we can write our ?rst order strongly hyperbolic system in a form where each and every term is manifestly regular on axis. Finally we discuss some of the implementation possibilities and applications of our scheme in section IX.
II. THE (2+1)+1 FORMALISM

Spacetime is assumed to be a four-dimensional manifold (M, gαβ ) with signature (? + ++) and a preferred polar coordinate chart (t, r, z, ?). Axisymmetry means that there is an everywhere spacelike Killing vector ?eld ξ = ?/?? with closed orbits. The idea is to perform ?rst a Kaluza-Klein-like reduction from M to the Lorentzian threedimensional manifold N formed by the trajectories of the Killing vector. (This was carried out for vacuum spacetimes by Geroch [11]. Nakamura et al [14, 15] extended the reduction to include general matter sources and then considered the case of a perfect ?uid. Choptuik et al [5] added a massless scalar ?eld source, albeit without rotation. The current analysis includes arbitrary sources and rotation.) After this an ADM-like reduction is applied to N . We shall follow the notation of Nakamura et al. [14, 15], with some, clearly stated, changes. In the following, Greek indices range over t, r, z, ?, lower-case Latin indices over t, r, z , and upper-case Latin indices over r, z . Tensor ?elds M α... β... in M which are both orthogonal to ξ and are axisymmetric, i.e., M α... β... ξ β = M α... β... ξα = . . . = 0, Lξ M α... β... = 0 are said to be tensor ?elds in N . Some basic tensor ?elds in N are the norm of the Killing vector λ2 = gαβ ξ α ξ β > 0 , the (Lorentzian) metric in N , hαβ = gαβ ? λ?2 ξα ξβ ,

3 the Levi-Civita tensor ?αβγ = λ?1 ?αβγδ ξ δ , and the “twist vector” ωα = ?αβγδ ξ β ?γ ξ δ , (2)

where ? is the covariant derivative of gαβ . Note that neither ξ α = δ α ? nor ξα = gα? is in N . However if a solution of the equations in N is given then the solution in M can be reconstructed [11]. We shall argue that a reconstruction is not necessary. All physically relevant quantities in M have their counterparts in N . We de?ne the projections of the energy-momentum tensor as follows, τ = λ?2 ξ γ ξ δ Tγδ , τα = λ?2 hα γ ξ δ Tγδ , ταβ = hα hβ Tγδ . Note that near the axis λ = O(r), and the factors of λ are included in (3) to ensure that the projections are either O(1) or O(r) near the axis, so that regularity conditions (cf. section VIII) are easily imposed. Following the standard ADM procedure, the three-dimensional manifold N is foliated into two-dimensional spacelike hypersurfaces Σ(t) with unit normal na = ?α ?a t and intrinsic metric Hab = hab + na nb . The line element in N takes the form ds2 = ?α2 dt2 + HAB dxA + η A dt dxB + η B dt ,
γ δ

(3)

where α is the lapse function and η A is the shift vector. The extrinsic curvature is given by χab = ?Ha c Hb d n(c|d) = ? 1 2 Ln Hab , where | denotes the covariant derivative compatible with hab . The trace of the extrinsic curvature is abbreviated as χ = χA A . Further variables de?ned in each Σ(t) are the alternating symbol ?AB = nc ?cAB , the ??-component of the extrinsic curvature, K? ? = ?λ?1 na λ,a , and the rotational degrees of freedom, E A = λ?3 ?Ab ωb , B ? = λ?3 na ω a . Again, the last two de?nitions di?er from those in [14, 15] by factors of λ, and we write B ? with an upstairs index. The various projections of the energy-momentum tensor are J? SA ρH JA SAB and, of course τ de?ned in (3). For convenience we also introduce AA ≡ α?1 α,A , LA ≡ λ?1 λ,A . (6) = = = = = ? na τ a , Ha A τ a , na nb τab , ?HA a nb τab , HA a HB b τab ,

(4) (5)

4 We can now write down the projections of Einstein’s equations Gαβ = κTαβ , which split into sets of elliptic constraint equations and hyperbolic evolution equations. The constraint equations are C = CA = C? = 0, where the energy constraint is
1 2 2 AB C≡ 1 + (2) R) ? LA ||A ? LA LA + χK? ? ? 4 λ EA E A + B ? 2 ? κρH , 2 (χ ? χAB χ

(7)

with || denoting the covariant derivative of HAB , The linear momentum constraints are

(2)

RAB the Ricci tensor of HAB and (2) R = (2) RA A the Ricci scalar. (8)

2 ? B CA ≡ χA B ||B ? (χ + K? ? ),A + LB χAB ? LA K? ? ? 1 2 λ B ?AB E ? κJA ,

and the angular momentum or Geroch constraint is
1 A A ? E ||A + 3 C? ≡ 2 2 LA E ? κJ .

(9)

We shall express all the evolution equations in terms of the Lie derivative along the normal direction, Ln = α?1 (?t ? Lη ). The time evolution of the metric variables is given by Ln HAB = ?2 χAB , Ln λ = ?λK? ? . The evolution equations for the extrinsic curvature variables are Ln χAB =
(2) 1 2 ?2 λ ?AC ?BD E C E D ? HAB EC E C ? B ? 2 1 HAB ρH ? SC C ? τ ?κ SAB + 2 ? ? A ||A

(10) (11)

RAB + (χ + K? ? )χAB ? 2χA C χCB ? AA||B ? AA AB ? LA||B ? LA LB , (12) EA E A ? B
?2

Ln K?

?

= K? (χ + K? ) ? L
1 κ ?2

ρH ? S A A + τ .

? LA (L + A ) ?

A

A

1 2 2λ

(13)

The rotation variables evolve according to Ln E A = (χ + 3K? ? )E A + ?AB B ? ,B + ?AB (AB + 3LB )B ? ? 2κS A , Ln B ? = χB ? + ?AB EA AB + ?AB EA||B . (14) (15)

(The principal parts of these equations resemble the axisymmetric Maxwell equations, which justi?es the notation. The Geroch constraint (9) is also a “Maxwell equation”.) Energy-momentum conservation ?α T αβ = 0 implies the following evolution equations for the matter variables: the energy conservation equation Ln ρH = ?J A ||A ? J A (2AA + LA ) + χρH + χAB S AB + K? ? τ + λ2 E A SA , the Euler equations Ln JA = ?SAB ||B + JA (χ + K? ? ) ? SAB (AB + LB ) ? AA ρH + LA τ +λ2 (EA J ? + ?AB S B B ? ) , and angular momentum conservation Ln J ? = ?S A ||A ? S A (AA + 3LA ) + J ? (χ + 3K? ? ) . (18) (17) (16)

In many situations there may be a conserved particle number density N α satisfying Lξ N α = 0. If we introduce the projections ν α = hα β N β , ΣA = H A b ν b , σ = ? ν a na ,

then particle number conservation ?α N α = 0 implies the following evolution equation for σ , Ln σ = ?ΣA ||A ? ΣA (AA + LA ) + σ (χ + K? ? ) . (19)

5
III. THE Z4 EXTENSION

Bona et al. [2] suggested the following covariant extension of the Einstein ?eld equations:
1 T gαβ , Rαβ + 2?(α Zβ ) = κ Tαβ ? 2

(20)

which reduces to the Einstein equations if and only if Zα = 0. We shall assume that Zα shares the axisymmetry, i.e., Lξ Zα = 0. To apply the (2+1)+1 formalism to equations (20), it is convenient to rewrite them as Einstein’s equations Gαβ = κTαβ with a modi?ed energy-momentum tensor Tαβ = Tαβ ? 2 γ ?(α Zβ ) ? 1 . 2 gαβ ?γ Z κ (21)

The Z covector is decomposed into parts parallel and orthogonal to the Killing vector ξ α , Z ? = λ?2 ξ α Zα , Z a = ha α Z α ,

and further parallel and orthogonal to the timelike normal na , θ = ? na Z a , Z A = HA a Z a .

We now compute the modi?ed matter variables corresponding to Tαβ , ?nding τ = τ + κ?1 Ln θ + Z A ||A + (AA ? LA )Z A + (K? ? ? χ)θ ,

SA = SA + κ?1 ?Z ? ,A + B ? ?AB Z B + EA θ , SAB = SAB + κ?1 ? 2Z(A||B ) + 2χAB θ +HAB Ln θ + Z
C ||C

J ? = J ? + κ?1 Ln Z ? + E A ZA ,

(22) + (AC + LC )Z ? (χ + K? )θ
C ?

,

ρH = ρH + κ?1 Ln θ ? Z A ||A + (AA ? LA )Z A + (χ + K? ? )θ . Inserting the above into the standard (2+1)+1 equations, we obtain the Z(2+1)+1 equations. Because of the Lie derivatives in (22), the constraints (7–9) become evolution equations for the Z covector, Ln θ = C + (LA ? AA )Z A + Z A ||A ? (χ + K? ? )θ ,
B ?

JA = JA + κ?1 Ln ZA ? θ,A + 2χAB Z B + AA θ ,

(23) (24) (25)

Ln ZA = CA ? 2χAB Z ? AA θ + θ,A , Ln Z = C? ? E ZA .
A

The main evolution equations are modi?ed as follows, Ln χAB = . . . + 2Z(A||B ) ? 2χAB θ , Ln E A = . . . + 2Z ?,A ? 2E A θ ? 2B ? ?AB ZB ,

Ln K? ? = . . . + 2LA Z A ? 2K? ? θ ,

where . . . denote the right-hand-sides of (12), (13) and (14), respectively. The remaining evolution equations are unchanged.

6
IV. DYNAMICAL GAUGE CONDITIONS

To complete our evolution formalism, we need to prescribe the gauge variables α and η A . It is well-known [6, 8] that the Einstein equations can be written in symmetric hyperbolic form g γδ gαβ,γδ ? 0 (? denoting equality to principal parts) by adopting harmonic gauge, which for the Z4 extension of the ?eld equations becomes g γδ xα ;γδ = 2Z α , where we are treating xα as a scalar ?eld. It would therefore be desirable to choose this or a related gauge for the system at hand. However, these gauge conditions are incompatible with cylindrical polar coordinates; the wave operator contains a term r?1 ?r , and the right-hand-side of the evolution equation for the shift vector becomes singular on the axis. However, we can keep the time component of the harmonic gauge equation, which in (2+1)+1 language reads ?t α = ?α2 (χ + K? ? ? 2θ) . This condition has been generalized to [3] ?t α = ?f α2 (χ + K? ? ? mθ) , (26)

where f and m are free constant parameters. We choose the shift vector to vanish: η A = 0. Apart from simplifying the evolution equations, this choice also makes it much easier to set up outer boundary conditions because a variable shift would mean that the characteristic speeds could change sign at the boundary (see section VI). Our choice for the shift vector means that henceforth Ln = α?1 ?t .
V. FIRST-ORDER REDUCTION

The Z(2+1)+1 equations are ?rst-order in time but second-order in space. In order to analyze their hyperbolicity and to be able to treat them with standard numerical methods for hyperbolic conservation laws, we perform a ?rst-order reduction. We focus on the geometry evolution equations here; for given geometry, the matter evolution equations form a separate strongly hyperbolic system (see appendix D). To eliminate the second spatial derivatives, we introduce new variables for the ?rst spatial derivatives of the twometric, DABC ≡ 1 2 ?A HBC , and we regard AA and LA de?ned by (6) as independent variables. Evolution equations for these new ?rst-order variables can be obtained from (10), (11) and (26) by commuting space and time derivatives: Ln DABC = ?χBC,A , Ln LA = ?K? ? ,A , The two independent traces of DABC are denoted by DI A ≡ DAB B , DII A ≡ DB BA , (27) (28) (29)

Ln AA = ?f (χ,A + K? ? ,A ? mθ,A ) .

where indices are (formally, for D is not a tensor) raised and lowered with the 2-metric HAB . The De Donder-Fock decomposition [6, 8] is used for the Ricci tensor,
(2)

RAB = ?DC AB,C + 2DII (A,B ) ? DI (A,B ) ? 2DCAB DII C

?ΓCAB (2DII C ? DI C ) + 4DCDA DCD B ? ΓACD ΓB CD ,

7 where of course the Christo?el symbols are given by ΓABC = DCAB + DBCA ? DABC . It is then straightforward to write the Z(2+1)+1 equations in conservation form ?t u + αF D (u) Here, the set of conserved variables is u = (HAB , λ, α, DABC , LA , AA , χAB , K? ? , E A , B ? , θ, ZA , Z ? )T , and the components of the ?ux vector are given by F D HAB F Dλ F Dα D F DABC F D LA = = = = = 0, 0, 0, D δA χBC , D δA K ? ? ,
,D

(30)

= αS (u) .

(31)

F D K? ? = L D , F D EA F D B? F Dθ F D ZA = = = =

D II I F D χAB = DD AB ? δ( A 2D B ) + 2ZB ) ? D B ) ? LB ) ? AB ) ,

D f (χ + K? ? ? mθ) , F D AA = δ A

1 D F D Z? = ? 2 E .

?2H AD Z ? ? ?AD B ? , ??AD EA , DI D ? DII D + LD ? Z D , D ? χA D + δ A (χ + K? ? ? θ) ,

S (u) is a source term containing no derivatives, SHAB Sλ Sα SDABC SLA SAA SχAB = = = = = = = ?2χAB , ?λK? ? , ?f α(χ + K? ? ? mθ) , 0, 0, 0, AC DC AB + A(A ?2DII B ) + DI B ) + LB ) + AB ) ? 2ZB )

?LA LB ? AA AB + ΓC AB (LC + AC ) + χAB (χ + K? ? ) ? 2χA C χCB 2 C D C ?2 ) ?1 2 λ ?AC ?BD E E ? HAB (EC E ? B ?2θχAB ? 2DCAB DII C + 4DCDA DCD B ? ΓACD ΓB CD , 2 A ?2 = K? ? (χ + K? ? ) + LA (2Z A ? LA ? DI A ) ? 1 ) ? 2K? ? θ 2 λ (EA E + B
C ?1 2 κ(ρH ? SC + τ ) , C II C C ?κ SAB + 1 ? DI C ) 2 HAB (ρH ? SC ? τ ) ? ΓCAB (2Z + 2D

SK? ?

SE A = (χ + 3K? ? ? 2θ)E A + ?AB B ? (3LB ? 2ZB + DI B ) ? 2κS A + (4DII A ? 2AA )Z ? , SB ? = χB ? + ?AB EA DI B , 1 1 I Sθ = AA (DI A ? DII A + LA ? Z A ) + DABC DABC ? 2 ΓABC ΓABC ? 2 D A DI A +(LA ? AA + DI A )Z A ? (χ + K? ? )θ ,

2 AB 2 A ?2 ?LA (LA + DI A ) + 1 ) + χK? ? ? 1 ) ? κρH 2 (χ ? χAB χ 4 λ (EA E + B

SZA = ?AB χA B + AA (χ + K? ? ? θ) + χAB (DI B + LB ? 2Z B ) ? ΓC AB χC B 2 ? B ?LA K? ? ? 1 2 λ B ?AB E ? κJA ? AA θ , SZ ? =
1 A I 2 E (D A

+ 3LA ? 2ZA ? AA ) ? κJ ? .

8 Note that HAB , λ and α have vanishing ?uxes and thus trivially propagate along the normal direction. The rotation variables (E A , B ? , Z ? ) form a decoupled subsystem analogous to Maxwell’s equations.
VI. HYPERBOLICITY

To investigate the hyperbolicity of the Z(2+1)+1 system, we pick a unit covector ?A and de?ne an orthogonal covector πA = ?AB ?B , so that ?A ?A = πA π A = 1 , ?A π A = 0 . (32)

Thus (?A , π A ) form an orthonormal basis for T (Σ). Projection along ? and π is denoted as V ≡ V A ?A , It can be veri?ed that the Jacobian matrix B ≡ ?F ?u V ⊥ ≡ V A πA . (33)

corresponding to the ?ux in the ?-direction is real, diagonalizable and has complete sets of left and right eigenvectors for arbitrary ?A , i.e., the system is strongly hyperbolic. The characteristic speeds are3 λ0 = 0 , λ± 1 = ±1 = ± (speed of light) , f (gauge speed) .

λ± f

For f 1, all characteristic speeds are causal, for f = 1 (harmonic slicing) they are all physical. The characteristic variables (left eigenvectors dotted into the conserved variables) are given by ? Normal modes (eigenvalue 0) l0,1 l0,2 l0,3 l0,4 l0,5 l0,6 l0,7 = = = = = = = f m(D ⊥⊥ + L ? D⊥⊥ ? Z ) ? f (D + D ⊥⊥ + L ) + A , f m(D⊥ + L⊥ ? D ⊥ ? Z⊥ ) ? f (D⊥ + D⊥⊥⊥ + L⊥ ) + A⊥ , D⊥ , D⊥⊥ , D⊥⊥⊥ , L⊥ , A⊥ ,

along with the zeroth-order variables HAB , λ and α. ? Light cone modes (eigenvalue ±1)
± l1 ,2 ± l1 ,3 ± l1 ,4 ± l1,5 ± l1 ,6 ± l1 ,1 = θ ± (D

⊥⊥

= χ



= E⊥ ? B? ,

= E ? 2Z ? ,

= K? ± L ,

= χ⊥⊥ ± D
?

±

1 2 (A⊥

+ L ? D⊥⊥ ? Z ) , + D⊥

⊥⊥ ,

? D⊥⊥⊥ + L⊥ ? 2Z⊥ ) ,

3

Note that a factor of α has been taken out of the ?uxes (31). Note also that for a nonzero shift vector ηA , a term ?α?1 η would have to be added to the above characteristic speeds.

9 √ ? Lapse cone modes (eigenvalue ± f )
± lf = A ?

f (m ? 2) (D f ?1

⊥⊥

+ L ? D⊥⊥ ? Z ) f (m ? 2) +2 θ . f ?1 (34)

±

f χ

+ χ⊥⊥ + K? ? ?

The inverse transformation from characteristic to conserved variables is given in Appendix B. In the special case f = 1, we must set m = 2 to enforce strong hyperbolicity, and the singular term (m ? 2)/(f ? 1) in (34) is to be replaced with an arbitrary ?xed constant (e.g. 0 for simplicity). In this case, the system is even symmetric hyperbolic, i.e. B is symmetrizable with a symmetrizer that is independent of the direction ?A . This is not surprising because the choice of parameters f = 1, m = 2 in (26) corresponds to harmonic slicing. An energy for the system is given by E = χAB χAB + λABC λABC + (K? ? + χ ? 2θ)2 + AA AA +VA V A + K? ? 2 + LA LA + EA E A + B ? 2 + 4Z ? 2 , where VA ≡ AA + DI A + LA ? 2DII A ? 2ZA , A λA BC ≡ DA BC + δ( B VC ) . E is positive de?nite and its time derivative is (considering principal parts only) a total divergence.
VII. CONSTRAINT-PRESERVING BOUNDARY CONDITIONS

In the Z(2+1)+1 formalism, the standard (2+1)+1 constraints (7–9) are replaced with the algebraic constraints θ = ZA = Z ? = 0 . (35)

If those hold at all times, the (2+1)+1 constraints are automatically satis?ed by virtue of the evolution equations for the Z vector (23–25). On the initial time level, one imposes (35) and solves the (2+1)+1 constraints (7–9) so that both the Z vector and its time derivative vanish. In addition, the outer boundary conditions have to be chosen so that no constraint-violating modes propagate into the computational domain at any time. ? αβ = 0. If we A propagation system for (35) follows from the contracted Bianchi identities, which imply ?α T impose the standard matter evolution equations ?α T αβ = 0 in addition, (21) yields the following homogeneous wave equation for the Z covector: ?β ?β Zα + Rαβ Z β = 0 . (36)

To obtain the (2+1)+1 reduction of (36), we take an additional time derivative of the evolution equations for the Z vector (23–25). To principal parts, we have
B L2 n θ ? ?B ? θ 2 B Ln ZA ? ?B ? ZA ? L2 nZ

? 0, B IB ? 1 + LB ) ? ? B (AA + DI A + LA ) 2 ?B ?A (A + D ? 0. ?2(?C DBC A ? ? B DII A ) ,

(37)

? ?B ? Z

B

?

Note that the right-hand-sides are zero if the ordering constraints ?D DABC = ?A DDBC , ?B LA = ?A LB , ?B AA = ?A AB

(38)

are satis?ed, which corresponds to commuting of partial derivatives in the de?nition of the ?rst-order variables. Although this holds analytically, it may not obtain numerically, and so the terms have to be included in the subsidiary system (37).

10 We choose an orthonormal basis ?A , π A such that ?A is normal to the outer boundary under consideration. Indices contracted with ?A (π A ) are denoted by (⊥). Equations (37) decompose into
2 2 L2 n Z ? ? Z ? ?⊥ Z 2 2 L2 n θ ? ? θ ? ?⊥ θ ? 0 ,

? ? 2 D⊥⊥

+1 2 ? ?⊥ (A⊥ ? D⊥
1 2 +2 ?⊥ (?A + D 1 2 2 ? (?A⊥ ? D⊥

L2 n Z⊥

? ? Z⊥ ?

2

2 ?⊥ Z⊥

?

1 ? ?⊥ (A + D +2 ⊥

+ D⊥⊥⊥ ? L⊥ ) ?D
⊥⊥

?D

+ D⊥⊥⊥ ? 2D
⊥⊥



+ L⊥ )

? L ),

2 +?⊥ D ? L2 nZ

,

? 2D⊥⊥ + L )

?? Z ?

2

?

2 ?⊥ Z?

? 0.

This second-order system can be written in ?rst-order form by introducing new variables for the temporal and spatial derivatives of the fundamental variables. Absorbing boundary conditions can then be imposed in the direction, which read (after replacing the new variables with their de?nitions as ?rst derivatives of the fundamental variables) . Ln θ = ?? θ , (39) . 1 (40) Ln Z = ?? (Z + D⊥⊥ ) ? 2 ?⊥ (A⊥ ? D⊥ + D⊥⊥⊥ ? 2D ⊥ + L⊥ ) , . 1 Ln Z⊥ = ?? Z⊥ ? 2 (A⊥ + D⊥ ? D⊥⊥⊥ + L⊥ ) Ln Z
? 1 ?⊥ (A + D ?2 . = ?? Z ? ,

?D

⊥⊥

? 2D⊥⊥ + L ) ,

(41) (42)

. where = denotes equality at the boundary. The straightforward way of setting up stable outer boundary conditions for the main evolution system is to set all ingoing modes to zero at the outer boundaries of the computational domain (so-called absorbing boundary conditions). For a symmetric hyperbolic system, this leads to a well-posed initial boundary value problem (IBVP) [17, 18]. Our goal is to investigate under which circumstances such boundary conditions imply the constraint-preserving ?rst-order conditions (39–42). Suppose l? is an incoming mode with characteristic speed ?λ (λ > 0). Setting l? to zero at the entire boundary at all times implies that . 0 = ?⊥ l? ≡ L1 (l? ) (43) and . 0 = Ln l? ? λ ? l? ? ?⊥ F ⊥ l? ≡ L2 (l? ) , (44)

. where = denotes equality at the boundary and we have used the general general evolution equation for l? . The relations (43) and (44) will now be used to manipulate the general evolution equations for the Z vector at the boundary. As mentioned in section VI, the evolution system decouples into a non-rotational and a rotational part at the level of principal parts, and we consider the two subsystems separately. Let us start with the general evolution equations for θ and ZA , decomposed with respect to our boundary-adapted basis: Ln θ ? ?? (D ⊥⊥ ? D⊥⊥ + L ? Z ) ? ?⊥ (D⊥ Ln Z ? ?? (χ⊥⊥ + K? ? ? θ) + ?⊥ χ ⊥ , Ln Z⊥ ? ? χ ⊥ ? ?⊥ (χ + K? ? ? θ) . If we add to these suitable combinations of equations (43) and (44),
? Ln θ + = L2 (?l1 ,1 ) ,

?D



+ L ⊥ ? Z⊥ ) ,

? ? ? ? Ln Z⊥ + = L2 (?l1 ,2 ) + L1 (l1,1 ? l1,3 ? lf ) ,

Ln Z

? ? ? ? + = L2 (?l1 ,1 + l1,3 + l1,4 ) + L1 (?l1,2 ) ,

we obtain precisely the constraint-preserving boundary conditions (39–41), provided that we choose harmonic gauge ?2 f = 1, m = 2, m f ?1 = 0. (For di?erent values of f and m, the derivatives tangential to the boundary in the evolution

11 equation for Z⊥ di?er from the constraint-preserving boundary condition (41).) Hence absorbing boundary conditions for the non-rotational part of the evolution system preserve the constraints. However, this is not true for the rotational subsystem. A little calculation shows that the only dissipative boundary condition consistent with (42) is ? . ? . + ? . l1 l1 ,5 = 0 , ,6 = l1,6 ? B = 0 . Fortunately, there is a much simpler solution to this problem: we delete Z ? entirely from the evolution system. Note that setting Z ? ≡ 0 does not break general covariance because our spacetime has a Killing vector and we simply choose Zα ξ α = 0. The remaining rotational subsystem is still symmetric hyperbolic, and the constraint
A C? ? 1 2 E ,A

has zero speed with respect to the boundary:
1 AB ? ? B ,AB = 0 . Ln C? ? 2

The absorbing boundary condition
? ⊥ ? . l1 ,6 = E + B = 0

(45)

leads to a well-posed IBVP for this system. In addition to the Z constraints, there are di?erential constraints related to the de?nition of the ?rst-order variables, DABC = LA = AA =
1 2 ?A HBC , λ?1 ?A λ , α?1 ?A α ,

and the ordering constraints (38). By virtue of the evolution equations (10, 11, 26) and (27, 28, 29), it is easy to see that the di?erential and ordering constraints have zero speed with respect to the boundary. In summary, we have shown that for harmonic gauge, absorbing boundary conditions for the non-rotational part of the main evolution system are automatically constraint-preserving, and that there are two ways of setting up constraint-preserving stable outer boundary conditions for the rotational subsystem. We thus have a well-posed IBVP with constraint-preserving boundary conditions. The result concerning the non-rotational part carries over immediately to the general Z4 system without spacetime symmetries. Constraint-preserving boundary conditions for the Z4 system have recently also been analysed by Bona et al. [4]: in essence, those authors propose to implement the constraint-preserving boundary conditions (39–41) directly. While it is not clear analytically that such boundary conditions lead to a well-posed IBVP, Bona et al. are able to prove a necessary condition and to demonstrate numerical stability.
VIII. REGULARITY ON AXIS

Extra care is needed in axisymmetric situations because the coordinate system that is adapted to the symmetry, cylindrical polar coordinates, is singular on the axis. This leads to rather strong regularity conditions on tensor ?elds, as explained in appendix A. In addition to the axisymmetry, we would like to impose re?ection symmetry about the z -axis so that we only need to evolve one half of the (r, z ) plane. This implies that tensor components are either odd or even functions of z . Let us ?rst deal with one of the regularity conditions for 2-tensors Mαβ , which follows from (1) M?? = 1 + O(r2 ) r2 Mrr near the r-axis. For the metric gαβ this implies λ2 g?? = 2 = 1 + O(r2 ) . 2 r Hrr r grr We enforce this condition by replacing λ with a new variable s ? de?ned by λ = r er
2

s ?

Hrr .

(46)

12 ? via To satisfy the corresponding condition for the extrinsic curvature, we introduce Y χrr ? K? ? = + r2 Y Hrr (note that K?? = λ2 K? ? ). Similary for the energy-momentum tensor, we set τ= Srr + r2 τ ?. Hrr

(47)

? (47) can be viewed as a generalization of those in [5] and We remark that the de?nitions of the variables s ? (46) and Y [10]. The remaining dependent variables are rede?ned by taking out systematically the leading order of r and z : ? rrr = D ? zrr = D Hrr 1 ? H /r , 2 rr,r 1 ? 2 Hrr,z /z , ? ? rz = Hrz /(rz ) , = Hrr , H 1 ? rrz = H ? ? rzz = D D 2 rz,r /r , ? zrz = 1 H ? rz,z /z , D ? zzz = D
2 ?1

Hzz = Hzz 1 ? 2 Hzz,r /r , 1 ? 2 Hzz,z /z ,

?

,

s ?r = s ?,r /r , s ?z = s ?,z /z , ? 1 ?r = α ?z = α α ? = α, A ? α ? ,r /r , A ? α ?,z /z , χ ?rr = χrr , χ ?rz = χrz /(rz ) , χ ?zz = χzz , r r z z ? ? ? ? E = E /r , E = E /z , B = B ? /(rz ) , σ ? = σ, ?r = Zr /r , Z ?z = Zz /z , Z ?? = Z ? , Z ?r = Jr /r , J ?z = Jz /z , J ?? = J ? , ρ ?K = ρH ? σ , J ? r = Σr /r , Σ ? z = Σz /z , S ?r = Sr /r , S ?z = Sz /z , Σ ?rr = Srr , S ?rz = Srz /(rz ) , S ?zz = Szz . S ?= θ, θ

It can now be veri?ed4 that the Z(2+1)+1 equations can be written in terms of the new variables u ? in the form ?(r ) (u ?t u ? + αF ?)
2

,r 2

?(z ) (u + αF ?)

2

,z 2

?(u = αS ?) ,

(48)

? and the source S ? are even in r and z and manifestly regular on the axes (cf. appendix C). where the ?uxes F The above transformation of variables does not a?ect the hyperbolicity of the system. To see this one should note that ? rrr ? Drrr /r , D ? zrr ? Dzrr /z , D ? rzz ? Drzz /r , D ? zzz ? Dzzz /z , D Drrr Dzrr s ?r ? Lr ? /r3 , s ?z ? Lz ? /(r2 z ) , Hrr Hrr ?r ? Ar /r , A ?z ? Az /z , A ? ? (K? ? ? χrr )/r2 , χ ?rz ? χrz /(rz ) , χ ?zz ? χzz , Y Hrr r r z z ? ? ? ? E ? E /r , E ? E /z , B ? B ? /(rz ) , ?? θ, Z ?r ? Zr /r , Z ?z ? Zz /z , Z ?? ? Z ? , θ ? r ? Σr /r , Σ ? z ? Σz /z , Σ ?rr ? Srr , S ?rz S ?r ? Jr /r , J ?r ? Sr /r , S ? rrz ? Drrz /(r2 z ) , D ? zrz ? Dzrz /(rz 2 ) , D

χ ?rr ? χrr ,

?z ? Jz /z , J ?z ? Sz /z , S ?zz ? Szz , ? Srz /(rz ) , S

where ? denotes equality to leading derivative order (the lower-order terms are absorbed in the source terms)5 . Hence the variables occuring in the ?uxes undergo a linear transformation, which does not a?ect the hyperbolicity and leaves the characteristic variables unchanged.

4 5

The calculations are rather lengthy. We used the computer algebra system REDUCE. Note that Hrr has zero ?ux and is thus to be treated as a background quantity.

13 By writing the equations in terms of derivatives with respect to r2 and z 2 instead of r and z , however, we have ultimately to multiply the ?uxes by 2r and 2z , respectively, because ?r2 = 1 ?r 2r

and similarly for z . This means that in (r2 , z 2 ) coordinates, the characteristic speeds behave like ? r and ? z , respectively. For computational purposes it is highly desirable to avoid the non-uniform characteristic speeds which would occur on a (r2 , z 2 ) grid. We therefore recommend working in the original (r, z ) grid and discretizing the derivatives in (48) with respect to r2 and z 2 . This has another, perhaps more fundamental, advantage: in (r, z ) coordinates, we can enforce Neumann boundary conditions for all our variables on the axes (since they are even functions of r and z ), whereas in (r2 , z 2 ) coordinates the boundary conditions on the axes are not known. The authors of [15] faced a similar problem but decided to work on an (r2 , z 2 ) grid, using extrapolation on the axis, which led to numerical instabilities. Choptuik et al [5] and Gar?nkle et al [10] take a slightly di?erent approach. They include terms which are formally irregular on axis, but which are regular once appropriate boundary conditions are imposed. Following an idea of Evans [7] these terms are di?erenced with respect to r2 or r3 . The transformation from conserved to primitive ?uid variables and vice versa can be made manifestly regular on the axes by rede?ning the velocities v ?r = vr /r , v ?z = vz /z

and leaving the remaining primitive variables unchanged. Note in particular that the modi?ed general matter variable Srr Hrr
? Hrr v ? 2 ? = ρhW 2 e2r s
2

τ ? = r ?2 τ ?

2 v ?r Hrr

(49)

is then automatically regular. To compute the characteristic variables (geometry and matter), one can start from the regularized conserved variables u ?, compute the original conserved variables u and evaluate the characteristic variables (section VI and appendix D 3). While this transformation is perfectly regular on the axes, the inverse transformation involves negative powers of r and z . In order to implement outer boundary conditions at r = rmax (z = zmax ), one has to ensure that the inverse transformation in the r-direction (z -direction) is regular at z = 0 (r = 0). This can be achieved by rede?ning the characteristic variables in a similar fashion as was done for the conserved variables above. We choose to make the following replacements ? for the characteristic variables in the r-direction: ? rz zH l0,4 → ? l0,4 = z ?2 l0,4 + √ l0,5 H ,

and the leading order of z is taken out of the remaining characteristic variables ? for the characteristic variables in the z -direction: l0,4 → ? l0,4 = r?2 ? rz rH l0,4 + √ l0,5 H ,

l0,6 → ? l0,6 = r?3 (l0,6 ? l0,5 ) , ± ± ?2 ± ?± l1 (l1,4 ? l1 ,4 → l1,4 = r ,3 ) , and the leading order of r is taken out of the remaining characteristic variables We have checked with REDUCE that the transformation from characteristic variables to conserved variables and vice versa is then manifestly regular at the entire outer boundary. However, the transformation from characteristic variables in the r-direction (z -direction) to conserved variables is still formally singular at r = 0 (z = 0), and no regularization procedure can cure this problem. To understand this, one can note that the characteristic variables in the r-direction (z -direction) do not have a de?nite parity in r (z ). For this reason, numerical methods operating in the space of characteristic variables (typically ones based on the solution of the Riemann problem) appear to be unusable near the axes.

14
IX. CONCLUSIONS

We started out by trying to understand why so many previous attempts to evolve axisymmetric spacetimes failed because of on-axis instabilities. This led to a detailed study of the behaviour of the components of axisymmetric tensors with respect to cylindrical polar coordinates, given that the components with respect to cartesian coordinates were regular in the neighbourhood of the axis. This suggested, most strongly, that we should make a Kaluza-Klein-like reduction to a 3-dimensional Lorentzian spacetime [11], a step taken earlier by Maeda et al [14], [15]. (A di?erent version of this idea, without rotation, has been pursued by [5].) Choptuik et al’s evolution [5] was fully constrained—elliptic equations were solved for the metric components, making use of multigrid for computational e?ciency. We made similar choices and obtained similar results. However in strong ?eld situations, e.g., near-critical collapse of Brill waves, the diagonal dominance of the matrix describing the discretized elliptic operators failed, and so the multigrid iterations failed to converge. We have not found a satisfactory resolution of this problem, and so we went to the other extreme—our evolution algorithm involves no elliptic operators. To achieve this we have used the Z4 formalism [2], [3] adapted to our axisymmetric reduction, an evolution equation for the lapse function and a zero shift gauge. Our ?rst order hyperbolic system of conservation laws (with sources) is strongly hyperbolic and, for one choice of parameters, symmetric hyperbolic. In addition the dependent variables have been carefully chosen so that each and every term in the system is regular on axis. We have allowed for arbitrary matter sources and have presented explicitly the details for a perfect ?uid. Initial conditions have to be imposed via the constraint conditions, which means solving elliptic equations. By means of conformal rescalings on the initial hypersurface one can avoid the diagonal dominance problem on that hypersurface and so multigrid techniques are applicable. One must choose a high resolution shock capturing method for the evolution of the hyperbolic system. Our system contains one subtle disadvantage if we choose to work with numerical algorithms that require a transformation between conserved and characteristic variables. The conversion from characteristic variables to our “regularized” conserved variables, see section VIII, is formally singular on the axis. However, we have shown that the transformation can be used to impose outer boundary conditions in a regular way. Because the evolution does not involve elliptic operators it is straightforward to introduce adaptive mesh re?nement, where appropriate, to enhance the resolution. We shall be using our formalism to evolve a number of interesting axisymmetric spacetimes, including the e?ects of rotation and perfect ?uid sources.

15
APPENDIX A: CONSEQUENCES OF AXISYMMETRY

We want to use a (t, z, r, ?) chart adapted to the Killing vector ξ = ?/??. Unfortunately this chart is singular at the axis r = 0. We shall assume elementary ?atness : in a neighbourhood of the axis we can introduce local cartesian coordinates xA = (x, y ) such that x = r cos ?, y = r sin ? ?? r= x2 + y 2 , y ? = arctan . x

With respect to cartesian coordinates the Killing vector is ξ = ?y?/?x + x?/?y. Notice that this representation is valid everywhere, while ξ = ?/?? is valid only for r > 0. We say that a scalar function f (xA ) is regular on axis if f has a Taylor expansion with respect to x and y about A x = 0 convergent in some neighbourhood of r = 0. (Throughout this section we are ignoring t, z dependencies which are implicit in all calculations.) If f is axisymmetric Lξ f = 0 = ?f /?? ? f = f (r), r > 0. (A1)

If f had a Taylor expansion in r and was regular on axis then the expansion could contain only even powers of r since r = x2 + y 2 , and r has no such expansion. While the conclusion is correct, the argument is fallacious because equation (A1) is invalid at r = 0. Instead we start from K ≡ ?yf,x + xf,y = 0, (A2)

which is valid everywhere. In particular we may di?erentiate (A2) an arbitrary number of times with respect to x and also with respect to y . Then setting xA = 0 we obtain linear recurrence relations between the Taylor coe?cients on axis. These can be solved to show f = n fn (x2 + y 2 )n . Next consider a vector ?eld uα . For a = (t, z ), Lξ ua = 0 implies ?ua /?? = 0. This reduces to the scalar ?eld case and we may deduce ua = ua (r2 ). For ux and uy we have ?ux /?? + uy = 0, The general solution for r > 0 is ux = a(r) cos ? ? b(r) sin ?, However in the cartesian chart we have ?yux ,x + xux ,y + uy = 0, Clearly uA = 0 on axis. We write a = ra etc so that ux = xa ? yb, uy = ya + xb. (A4) ?yuy ,x + xuy ,y ? ux = 0. (A3) uy = a(r) sin ? + b(r) cos ?. ?uy /?? ? ux = 0.

We now regard a and b as unknown functions of x and y to be determined by substituting (A4) into (A3), di?erentiating the latter an arbitrary number of times, and then solving the recurrence relations for the Taylor coe?cients of a and b. Again we ?nd that a and b are even functions of r. Thus in the (t, z, r, ?) chart an axisymmetric vector ?eld which is regular on axis must take the form uα = (A, B, rC, D), (A5)

where A, B , C and D are functions of t, z and r2 . Next consider an axisymmetric covector ?eld ωα which is regular on axis. For a = (t, z ), Lξ ωa = 0 implies ?ωa /?? = 0. This reduces to the scalar ?eld case and we may deduce ωa = ωa (r2 ). For the other indices we ?nd ?yωx,x + xωx,y + ωy = 0, ?yωy,x + xωy,y ? ωx = 0, (A6)

which is equivalent to (A3), interchanging uA and ωA . We therefore deduce the analogue of (A4) and hence, in polar coordinates, ωα = (E, F, rG, r2 H ), (A7)

16 where E , F , G and H are functions of t, z and r2 . Finally we consider a symmetric valence 2 axisymmetric tensor ?eld Mαβ which is both axisymmetric and regular on axis. For (a, b) = (t, z ) we have Lξ Mab = 0 and so Mab = Mab (r2 ). Letting a = (t, z ) we have ?yMax,x + xMax,y + May = 0, ?yMay,x + xMay,y ? Max = 0. (A8) This is essentially the same as (A6) and we may deduce Mar = rAa (r2 ) and Ma? = r2 Ba (r2 ). The remaining Killing equations are ?yMxx,x + xMxx,y + 2Mxy = 0, ?yMyy,x + xMyy,y ? 2Mxy = 0, ?yMxy,x + xMxy,y + Myy ? Mxx = 0.

1 If we introduce new variables u = 1 2 (Mxx + Myy ), v = 2 (Mxx ? Myy ) and w = Mxy then

?yu,x + xu,y = 0 which implies u = u(r2 ). The remaining equations are ?yv,x + xv,y + 2w = 0, For r > 0 these can be written as v,? + 2w = 0, so that v = a(r) cos 2? ? b(r) sin 2?, w = a(r) sin 2? + b(r) cos 2?, w,? ? 2w = 0, ?yw,x + xw,y ? 2v = 0. (A9)

where a and b are arbitrary functions of r. If we set a = a/r2 and b = b/r2 then v = (x2 ? y 2 )a ? 2xyb, w = 2xya + (x2 ? y 2 )b. (A10)

(Note that (A9) and its ?rst derivatives imply that v , w and their ?rst derivatives vanish on axis, consistent with (A10).) Substituting (A10) into (A9) gives x3 a,y ? x2 y (a,x + 2b,y ) ? xy 2 (a,y ? 2b,x ) + y 3 a,x = 0,
2 2 3

(A11) (A12)

x b,y + x y (?b,x + 2a,y ) + xy (?2a,x ? b,y ) + y b,x = 0.

3

Di?erentiating these many times and, proceeding as in the scalar and vector cases we conclude that a and b are functions of r2 . Thus Mxx = u + (x2 ? y 2 )a ? 2xyb, Myy = u ? (x2 ? y 2 )a + 2xyb, Mr? = r3 b, Mxy = 2xya + (x2 ? y 2 )b.

Re-expressing these as polar components we obtain Mrr = u + r2 a, Finally combining all of the results we have ? Mαβ M?? = r2 (u ? r2 a). (A13)

where A, B, . . . , K are functions of t, z and r2 .

? A B rD r2 F ? ? B C rE r2 G ?, =? 2 3 ? ? rD rE H + r J r K r2 F r2 G r3 K r2 H ? r2 J

(A14)

17
APPENDIX B: TRANSFORMATION FROM CHARACTERISTIC TO CONSERVED VARIABLES

The conserved variables are given in terms of the characteristic variables by f (m ? 2) 1 + ? ? + ? 1 + + 1 (l1 = ? l0,1 + ,1 ? l1,1 ) ? 2 (l1,3 ? l1,3 + l1,4 ? l1,4 ) f 2(f ? 1) 1 ? l + + lf , + 2f f 1 (m ? 2) (f m ? 2) ? 1 + (l1,2 ? l1 = ? l0,2 + (l0,3 + l0,5 + l0,6 ) ? l0,7 + 2 ,2 ) , fm 2m 2f m + ? = 1 2 (l1,3 ? l1,3 ) , = l0,3 , = l0,4 , = l0,5 , ? 1 + = 2 (l1,4 ? l1 ,4 ) , = l0,6 , f (m ? 2) + + ? ? 1 , (l ? l1 = ,1 ) + 2 lf + lf 2(f ? 1) 1,1 = l0,7 , f (m ? 2) ? + ? + ? 1 + = + 1 (l1 ,1 + l1,1 ) ? 2 (l1,3 + l1,3 + l1,4 + l1,4 ) 2(f ? 1) 1 + ? + √ lf ? lf , 2 f ? 1 + = 2 (l1,2 + l1 ,2 ) ,
? 1 + 2 (l1,3 + l1,3 ) , ? 1 + 2 (l1,4 + l1,4 ) , ? 1 + 2 (l1,5 + l1,5 ) , ? 1 + 2 (l1,6 + l1,6 ) , ? 1 + ?2 (l1,6 ? l1 ,6 ) , ? 1 + 2 (l1,1 + l1,1 ) , + ? + ? + ? 1 (?l1 ?l0,4 + 2 ,1 + l1,1 + l1,3 ? l1,3 + l1,4 ? l1,4 ) , ? 1 1 + 2 (l0,3 ? l0,5 + l0,6 + l0,7 ) ? 2 (l1,2 ? l1,2 ) ,

D

D



D ⊥⊥ D⊥ D⊥⊥ D⊥⊥⊥ L L⊥ A A⊥ χ

χ



χ⊥⊥ = K? ? = E = E⊥ = B? = θ = Z = Z⊥ =

1 + ? ? l1 Z ? = ? (l1 ,5 ) . 4 ,5

Tensor components can then be computed as, for instance, χAB = ?A ?B χ + 2?(A πB ) χ


+ πA πB χ⊥⊥ .

APPENDIX C: REGULAR CONSERVATION FORM

In this appendix we provide some details of the regular form of the Z(2+1)+1 equations outlined in section VIII. The equations are written in conservation form (48) ?(r ) (u ?t u ? + αF ?)
2

,r 2

?(z ) (u + αF ?)

2

,z 2

?(u = αS ?) ,

where the modi?ed variables u ? are de?ned in section VIII. ? (r 2 ) , F ?(z2 ) and sources S ? have been obtained with the computer algebra system REDUCE and have The ?uxes F been typeset automatically using the TeX REDUCE Interface (TRI). From the same source we generated C++ code

18 using the REDUCE source code optimization packages SCOPE and GENTRAN. The resulting expressions are rather ?. lengthy. To get some idea of what is typically involved we state here the ?uxes and source for the variable Y Let us introduce the shorthand ? zz ? r2 z 2 H ?2 ? rr H H =H rz for the determinant of the 2-metric. The ?uxes are given by ? (r F ? Y ? (z F ? Y
2

)

?2 ? 2 ?1 ? ? rrr r2 z 2 H ? rr ? rzz ? 4H ?1 D ? zrr z 2 H ? rr ? zrz z 2 = 2 H ?1 D Hrz ? 2H ?1 D Hrz + 4H ?1 D

2

)

? rz s ? rr s = 2H ?1 z 2 (?r2 H ?r + H ?z ) ,

? ?1 (?A ?r + 2Z ?r ) + 2H ?1 r2 z 2 H ? rz (r2 H ? ?1 H ? rz s +2H ?r ? s ?z ) , rr rr

and the source term is
?1 ? ?1 ? 2 ? ?2 ? 2 ? ? 4χ ? rr ? rr ? rz Y ? ? = ? κτ ?rz ) Y + H ?1 χ ?rr z 2 H Hrz (r2 H S ? + 2 H ?1 χ ?2 ?rr H rr z Hrr Hrz + 2χ Y ?1 ? 2 2 2 ? ?3 ? 2 ?1 ? ? +H χ ?zz Hrr Y ? 4H Drrr r z Hrr Hrz ?2 ? 2 ?3 ? 4 ? rrz r2 z 2 H ? rr ? rrr ? rr ? rrr D Hrz ? H ?2 D r4 z 4 H Hrz + 4H ?1 D ? zrr z 2 H ? ?2 H ? rz ? rrz r4 z 4 H ? ?2 H ? 3 + 4 H ?1 D ? rrr D ? rrr D +2H ?2 D rr rr rz 2 2 4 ? 2 3 ? 2 ?2 ? ? zzz z H ? rz ? zrr r z H ? H ? ? 2H D ? rrr D +2H Drrr D rr rz

?1 ? 2 ?1 ? zrr r2 z 4 H ? rr ? zrr z 2 H ? rr ? rrz D ? rzz r2 z 2 H ? rz ? 4H ?1 D ? rrz D ? rrz D Hrz ? 6 H ?2 D ? 2 H ?2 D 2 ? 2 4 ? ?2 ? ?2 ? ? ? +4H Drrz Dzrz r z Hrz + 2H Drrz Dzzz z Hrr ?r ? ?1 H ? rz ? 2r2 z 2 H ? ?1 H ? rz Z ? rrz (2r4 z 2 H ? ?1 H ? rz s ?r H +2H ?1 D ?r + r2 z 2 A rr rr rr

? rz + 4r2 z 2 H ? rr s ? rz s ? ?2 H ? rz (?5r4 z 2 H ? rz s ?r H ? rrr H ?z ? 2r2 z 2 H ? ?r ? r2 z 2 A + H ?1 D rr 2 2 2 2 ? ?z + 2H ? rr ) ? rr ? 4z H ? rr Z ?r + 2z A ?z H +4r z Hrz Z ?2 ? 2 2 ? ?2 ? 3 4 2 ? 2 2 ? 2 2 ? ? rz + 2H ? rr ) +H Drrr r z H H (?r z Hrz s ?r + r z Hrr s ?z + 2r z Hrz s ? + 2z 2 H
rr rz

?z + 2z 2 Z ?z ? 1) ?r 2 z 2 s ?z ? z 2 A ?2 ? 2 2 ?2 4 2 ? ?1 ? ?1 ? ?1 ? ? rr ? rr +2H Drrz r z Hrz (r z Hrr Hrz s ?r ? 2r2 z 2 H Hrz s ? ? r2 z 2 s ?z ? z 2 H Hrz ? 1)
?2 ? ? zrz z 2 H ? rr ? zrr z 2 H ? rz ? 4H ?2 D ? rzz D ?2 H ? + H ?2 D Drzz D rzz rr + 4H ?1 ? 2 ?2 ? 2 2 ? 2 ? ? ? rr s ? rz s +H Drzz (?r s ?r ? Ar ? 2? s) + H Drzz r z Hrz (?r Hrz s ?r + H ?z + 2H ?)

?1 ? ?1 ? ? rr ? zrr z 2 (r2 H ? rr ?r H Hrz ? s ?z ) +2H ?1 D Hrz s ?r ? A ?2 ? 4 ?2 4 ? ?1 ? 2 ? ?1 ? ? ?1 H ? rz ) +H Dzrr z H (r H Hrz s ?r ? 2r H Hrz s ? ? r2 s ?z ? 4H rz rr rr rr

? zrz z 2 (A ?r + 2? ? zrz z 4 H ? rz (?r4 H ? rz s ? rr s ? rz s ? rz ) +2H ?1 D s) + 2H ?2 D ?r + r2 H ?z + 2r2 H ? + 2H 1 2 ? 2 ?z 2 ? zzz z 2 H ? rr (r2 H ? rz s ? rr s ? rz s + H ?2 D ?r ? H ?z ? 2H ?) ? e2r s Hz E 2 2 ?1 2 ? 2 2 ? rr ? rz ? E ?r 2 H ? rr ?r H ? rr H ?z H ?z 2 H ? rz ?? 2 H ? rr ? 2z 2 E ?r E s ?r ) ? r2 A +e2r s r (?z 4 E + z 2B
2 ?2 ?? ? ?1 Z ?r + 4H ? ?1 s ? ?r H ? ? ? ?1 ? + 2 A ? ?1 s +2r2 H rr rr ?Zr ? 2θ Y rr ?r Zr + r Y ? 2Ar Hrr s ?1 6 2 ? ?1 ? 2 2 4 2 ? ?1 ? 2 4 2 ? ?1 ? 2 ? +H (?r z H H s ? ? 4r z H H s ?s ?r + 2r z H H s ?r Zr rr rz r rr rz rr rz

?1 ? 2 ? ? r2 z 2 H ? rr s ? rz Y ? rr ? rz s ?r H ?2 Hrz s ? ? 2r 2 z 2 χ ?rz H +2r4 z 2 H ?r s ?z ? 2r2 z 2 A z ?1 ? 2 2 ?1 ? 2 ? ?1 ? 2 ? rr ? rr ? rr ?4r 2 z 2 H Hrz s ? + 4r 2 z 2 H Hrz s ?Zr ? 3r2 z 2 H Hrz s ?r

2 ?1 ? 2 ?1 ? 2 ? rz ? rr ? rz s ? rr ? rz s (r4 z 2 H Hrz s ?r ? r4 H ?r ? 2r2 z 2 H Hrz s ? ? r2 z 2 H ?z + H ?2 z 2 H 2 ? 2 ? 2 ? ?1 ? 2 +r Hrr s ?z + 2r Hrz s ?? z H H ). rr rz

?r ? r2 H ? rz s ?z ? 2r2 z 2 H ? rz s ? rz s ? rz s ?r ?z Z +4r2 z 2 H ?s ?z ? 2r2 z 2 H ?r Z 2 ? 1 2 2 2 2 ? 2 ? ? ?z ? 6z H ? H ? s ? ? + 2z χ ?rz + 2z Hrr s ?z Z +2z Az Hrz s rr rz ? rz s ?z + 4z 2 H ? rz s ? rr s ? rz s ?4z 2 H ?Z ?z + H ?z + 2H ?)

Note that the expressions are manifestly regular on the axes and even in r and z .

19
APPENDIX D: MATTER EVOLUTION EQUATIONS 1. Conservation form

The matter evolution equations (16–19) can clearly be written in conservation form ?t u + αF D (u)
,D

= αS (u) .

Following [1], we replace ρH with ρK = ρH ? σ (kinetic energy) and regard as the set of conserved variables u = (ρK , JA , J ? , σ )T . The ?uxes are F D ρK = J D ? ΣD , F D J A = SA D , F D J? = SD , F D σ = ΣD , (D1)

and the source terms are SρK = (ΣA ? J A )(DI A + LA ) + K? ? (τ ? σ ) + χAB S AB ? J A AA + χρK SJA SJ ? +λ2 E A SA , = ?SAB (AB + LB ) + JA (χ + K? ? ) ? AA ρH + LA τ

Sσ = ?(DI A + LA )ΣA + σ (χ + K? ? ) .
2.

+λ2 (EA J ? + ?AB S B B ? ) , = ?(DI A + 3LA )S A + J ? (χ + 3K? ? ) ,

Perfect ?uid

To evaluate the characteristic structure, we need to specify the matter model. Here, we consider a perfect ?uid with four-velocity uα , normalized such that uα uα = ?1 , rest mass density ρ, pressure p and internal energy ?. The dependence of the pressure on the density and the internal energy is given by the equation of state p = p(ρ, ?) . With those de?nitions, the number density is N α = ρuα and the energy-momentum tensor is given by T αβ = ρhuα uβ + pg αβ , where h is the speci?c enthalpy, h=1+?+ The Lorentz factor is de?ned as W ≡ ? u α nα . (D4) p . ρ (D3) (D2)

20 Observers who are at rest in a slice Σ(t) (i.e., who have four-velocity nα ) measure a coordinate velocity v A = W ?1 h α A u α , and the angular velocity is v ? = W ?1 λ?2 ξα uα . Hence we obtain the familiar relation W = (1 ? v 2 )?1/2 , where v 2 = vA v A + λ2 v ? 2 . The variables w = (vA , v ? , ρ, p, ?, h, W )T are called primitive variables.6 The conserved variables can be expressed in terms of the primitive variables as ρK JA J? σ and the remaining matter variables are τ = SA = SAB = ΣA =
3.

(D5)

= = = =

ρhW 2 ? p ? ρW , ρhW 2 vA , ρhW 2 v ? , ρW ,

(D6)

ρhW 2 λ2 v ? 2 + p , ρhW 2 v ? vA , ρhW 2 vA vB + pHAB , ρW vA .

(D7)

Characteristic decomposition

The characteristic decomposition for 3+1 general relativistic hydrodynamics was ?rst derived by the Valencia group [1]. The application to our (2+1)+1 system is straightforward. Our method di?ers slightly in that we choose a general orthonormal basis (?A , π A ) in two-space as in section VI and project vectors along ? (index ) and π (index ⊥). Following the notation of [9], we introduce a few abbreviations. From the equation of state (D2), we form χ≡ ?p , ?ρ κ≡ ?p , ?? hc2 s ≡ χ+ p κ, ρ2

where cs is known as the sound speed. Also set7 K ?1 = 1 ? C± = v ? V ± , c2 sρ , κ V± = v ? λ± s , 1 ? v λ± s A± = 1 ? v λ± s 1 ? v2 ,

ξ = 1 ? v2 ,

? = h3 W (1 ? K?1 )(C + ? C ? )ξ .

6 7

Note only ?ve of those are independent because of (D2), (D3) and (D5). Our de?nitions of ξ and ? di?er from those in [9] by a factor of λ2 to ensure regularity (see section VIII). We have de?ned K?1 instead of K to allow for the special case of the ultrarelativistic equation of state, for which K?1 = 0. As a consequence, ? above has been multiplied by K?1 and the characteristic variable l0,1 has been divided by K?1 .

21 The system is found to be strongly hyperbolic. The characteristic speeds in the ?-direction are λ0 = v , λ± s = 1 1 ? v 2 c2 s v (1 ? c2 s ) ± cs
2 2 (1 ? v 2 ) (1 ? v 2 c2 s ) ? v (1 ? cs )

.

The characteristic variables (corresponding to the left eigenvectors) are l0,1 = l0,2 l0,3
? ls

W hσ ? W (σ + ρK ) + W (v J + v⊥ J⊥ + λ2 v ? J ? ) , 1 ? K ?1 1 = ?v⊥ (σ + ρK ) + v v⊥ J + (1 ? v 2 )J⊥ , hξ 1 = ?v ? (σ + ρK ) + v ? v J + (1 ? v 2 )J ? , hξ h2 = K?1 hW V ± ξσ + K?1 ? A± ? (2 ? K?1 )v J ? +(2 ? K?1 )V ± W 2 ξ (v J + v⊥ J⊥ + λ2 v ? J ? )

+ (K?1 ? 1) ?v + V ± (W 2 ξ ? 1) ? W 2 V ± ξ (σ + ρK ) .

The inverse transformation (corresponding to the right eigenvectors) is given by σ = J 1 + ? l0,1 + W (v⊥ l0,2 + λ2 v ? l0,3 ) + ls + ls , hW + ? = K?1 v l0,1 + 2hW 2 v (v⊥ l0,2 + λ2 v ? l0,3 ) + hW (C + ls + C ? ls ),

+ ? J ? = K?1 v ? l0,1 + hl0,3 + 2hW 2 v ? (v⊥ l0,2 + λ2 v ? l0,3 ) + hW v ? (ls + ls ), 1 l0,1 + W (2hW ? 1)(v⊥ l0,2 + λ2 v ? l0,3 ) ρK = K ?1 ? hW + ? + ? +hW (A+ ls + A? ls ) ? ls ? ls .

+ ? J⊥ = K?1 v⊥ l0,1 + hl0,2 + 2hW 2 v⊥ (v⊥ l0,2 + λ2 v ? l0,3 ) + hW v⊥ (ls + ls ),

4.

Transformation from conserved to primitive variables

The conserved matter variables (D1) are the ones that are evolved in a numerical algorithm. To compute the remaining matter variables (D7) and the eigenvectors, the primitive variables have to be calculated from the conserved variables as an intermediate step. This transformation is much more involved than the opposite direction (D6). To make it explicit, we have to specify an equation of state. Here, we consider the ultrarelativistic equation of state, p = (Γ ? 1)ρtot = (Γ ? 1)ρ(? + 1) = Γ?1 ρh , Γ (D8)

where ρtot is the total energy density. Suppose we are given the conserved variables, and also form ρH = ρK + σ . Consider the quantity J 2 ≡ JA J A + λ2 J ? 2 . Using (D6), (D8) and (D5), we can express J 2 and ρH in terms of the primitive variables as J2 = ρH Γ p2 W 2 (W 2 ? 1) , Γ?1 Γ = p W2 ? 1 . Γ?1
2

Eliminating W yields an equation for the pressure in terms of conserved variables: p = ?2βρH +
2 2 4 β 2 ρ2 H + (Γ ? 1)(ρH ? J ) ,

22 where β ≡ (2 ? Γ)/4. Next de?ne χA ≡ (Γ ? 1)JA , Γp χ? ≡ (Γ ? 1)J ? , Γp χ2 ≡ χA χA + λ2 χ? 2 .

We identify χA = W 2 v A and χ? = W 2 v ? and hence with (D5) we obtain W ?2 = This now enables us to calculate the velocities, vA = W ?2 χA , v ? = W ? 2 χ? . 1 2 χ2 1 + 4 χ2 ? 1 . (D9)

The form of W ?2 in (D9) guarantees that |vA |, |v ? | < 1. This is most important since evolved speeds greater than unity (i.e. greater than the speed of light) can easily cause the numerical code to crash [16]. Finally, we can calculate the speci?c enthalpy and rest mass energy density from (D6) and (D8), h= J? , σv ? W ρ= Γp . (Γ ? 1)h

ACKNOWLEDGMENTS

We thank Anita Barnes for several useful discussions. OR thanks the Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Golm, Germany for his inclusion in their Visitor Programme. OR acknowledges gratefully ?nancial support from the Engineering and Physical Sciences Research Council, the Gates Cambridge Trust and Trinity College Cambridge.

[1] Banyuls F, Font J A, Ib? an ?ez J M, Mart? i J M and Miralles J A 1997 Numerical {3+1} General Relativistic Hydrodynamics: A Local Characteristic Approach Astrophys. J. 476 221–231 ˇ aˇ [2] Bona C, Ledvinka T, Palenzuela C and Z? cek M 2003 General-covariant evolution formalism for numerical relativity Phys. Rev. D 67 104005 ˇ aˇ [3] Bona C, Ledvinka T, Palenzuela C and Z? cek M 2004 Symmetry-breaking mechanism for the Z4 general-covariant evolution system Phys. Rev. D 69 064036 ˇ aˇ [4] Bona C, Ledvinka T, Palenzuela-Luque C and Z? cek M 2004 Constraint-preserving boundary conditions in the Z4 Numerical Relativity formalism Preprint arXiv:gr-qc/0411110 [5] Choptuik M W, Hirschmann E W, Liebling S L and Pretorius F 2003 An axisymmetric gravitational collapse code Class. Quantum Grav. 20 1857–1878 [6] De Donder T 1921 La Gravi?que Einsteinienne Gauthier-Villars, Paris [7] Evans C R 1984 PhD. thesis, University of Texas, Austin, Texas, unpublished [8] Fock V A 1959 The Theory of Space Time and Gravitation Pergamon, London [9] Font J A 2003 Numerical Hydrodynamics in General Relativity Living Rev. Rel. 6 4 [10] Gar?nkle D and Duncan GC 2001 Numerical evolution of Brill waves Phys. Rev. D 63 044011 [11] Geroch R 1971 A Method for Generating Solutions of Einstein’s Equations J. Math. Phys. 12(6) 918–924 [12] Gustafsson B, Kreiss H-O and Sundstr¨ om A 1972 Stability Theory of Di?erence Approximations for Mixed Initial Boundary Value Problems. II Math. Comput. 26 649–686 [13] Lehner L 2001 Numerical relativity: a review Class. Quantum Grav. 18 R25–R86 [14] Maeda K, Sasaki M, Nakamura T and Miyama S 1980 A New Formulation of the Einstein Equations for Relativistic Rotating Systems Prog. Theor. Phys. 63 719–721 [15] Nakamura T, Oohara K and Kojima Y 1987 General Relativistic Collapse of Axially Symmetric Stars and 3D Time Evolution of Pure Gravitational Waves Prog. Theor. Phys. Suppl. 90 13–109 [16] Nielsen D W and Choptuik M W 2000 Ultrarelativistic ?uid dynamics Class. Quantum Grav. 17 733–759 [17] Rauch J 1985 Symmetric positive systems with boundary characteristics of constant multiplicity Trans. Am. Math. Soc. 291 167–187 [18] Secchi P 1996 The initial boundary value problem for linear symmetric hyperbolic systems with characteristic boundary of constant multiplicity Di?. Int. Eq. 9 671–700, Well-posedness of characteristic symmetric hyperbolic systems Arch. Rat. Mech. Anal. 134 155–197


赞助商链接

更多相关文章:
更多相关标签:

All rights reserved Powered by 甜梦文库 9512.net

copyright ©right 2010-2021。
甜梦文库内容来自网络,如有侵犯请联系客服。zhit325@126.com|网站地图