When considering the evolution of magnetohydrodynamics in a general relativistic setting, there are generally two sets of variables, the primitive variables
P={ρb,P,vi,Bi}and the conservative variables
C={˜ρ,˜τ,˜Si,˜Bi}.The details of the individual variables in each set are discussed more in Step 2. While the primitive variables are the quantities we associate with the usual description of physical properties, the conservative variables have far better numerical stability. As such, IllinoisGRMHD evolves C, and so the following evolution equations are in terms of these variables. The derivation in this document is based on the work by Duez et al (citation below), and I have rederived their results. They label their primitive density as ρ0 (0 for rest mass density), while we label ours as ρb (b for baryonic density). These are the same quantity despite the different labels.
This module is organized as follows
In order to find evolution equations for the conserved quantities for the system, we start by finding a convenient form of the electromagnetic energy-momentum tensor TμνEM, which is defined as
TμνEM=14π(FμλFννλ−14gμνFαβFαβ)where Fμν is the Faraday tensor. We want to examine this with respect to an arbitrary observer with normalized four-velocity ξμ. If we consider the electromagnetic fields in the rest frame of this observer,
ξμEμ=0=ξμBμsince E and B are purely spatial. The Faraday tensor can then be decomposed in terms of E and B:
Fμν=ξμEν−ξνEμ+ξγϵγμνδBδ.Now, we simply want to find an expression for of TμνEM in terms of ξ, E, and B. Let's first consider the relatively simpler term
FαβFαβ=(ξαEβ−ξβEα+ξγϵγαβδγαβBδ)(ξαEβ−ξβEα+ξσϵσαβμBμ)=ξαξαEβEβ−ξαEαEβξβ−EαξαξβEβ+ξβξβEαEα+ξαEβξσϵσαβμBμ−ξβEαξσϵσαβμBμ +ξαEβξγϵγαβδγαβBδ−ξβEαξγϵγαβδγαβBδ+ξσξγϵγαβδγαβϵσαβμBδBμ=−2EβEβ+4ξαξσϵσαβμEβBμ+ξσξγϵγαβδϵσαβμBδBμwhere we have used ξμξμ=−1 and the permutation rules of the Levi-Civita tensor. Further simplification requires abusing the properties of the Levi-Civita tensor, which is given in Appendix A:
FαβFαβ=−2EβEβ+4ξαξσϵσαβμEβBμ+2ξσξγ(δσδδμγ−δσγδμδ)BδBμ=−2EβEβ+4ξαξσϵσαβμEβBμ+2ξσξμBσBμ−2ξσξσBμBμ=−2EβEβ+4ξαξσϵσαβμEβBμ+2BμBμNow, since we are in the rest frame of the observer, the spatial components of the observer's velocity are zero. Then, ξαξσϵσαβμ=ξ0ξσϵσ0βμ=ξ20ϵ00βμ=0
since any components of the Levi-Civita tensor with repeated indices are zero. Therefore,
FαβFαβ=−2EβEβ+2BμBμIn this section, we tackle the term FμλFννλ in the energy-momentum tensor.
FμλFννλ=(ξμEλ−ξλEμ+ξγϵγμλδBδ)(ξνEλ−ξλEν+ξαϵανλβανλBβ)=ξμEλξνEλ−ξμEλξλEν−ξλEμξνEλ+ξλEμξλEν +ξμEλξαϵανλβBβ−ξλEμξαϵανλβBβ+ξγϵγμλδBδξνEλ −ξγϵγμλδBδξλEν+ξγϵγμλδBδξαϵανλβανλBβ=ξμξνEλEλ−EμEν+ξαEλBβ(ξμϵανλβ+ξνϵαμλβ) −ξλξαBβ(ϵανλβEμ+ϵαμλβEν)+ξαξγϵγμλδϵανλβανλBδBβwhere we have again used the relations ξμξμ=−1 and ξμEμ=0=ξμBμ. We've also taken the liberty to relabel dummy indices to simplify the resulting expression. To simplify the final term, we must (unfortunately) once again rely on the generalized Kronecker delta. This time we use the second (far more complicated) expression in Appendix A:
ξαξγϵγμλδϵανλβανλBδBβ=−gνϕξαξγBδBβ(δγα(δμϕδδβ−δμβδδϕ)−δγϕ(δμαδδβ−δμβδδα)+δγβ(δμαδδϕ−δμϕδδα))=−gνϕξαξαBδBβ(δμϕδδβ−δμβδδϕ)+gνγξαξγBδBβ(δμαδδβ−δμβδδα)−gνϕξαξβBδBβ(δμαδδϕ−δμϕδδα)=−gμνξαξαBβBβ+gνδξαξαBδBμ+gνγξμξγBβBβ−gνγξαξγBαBμ−gνδξμξβBδBβ+gμνξαξβBαBβ=gμνBβBβ−BμBν+ξμξνBβBβwhich leaves us with
FμλFννλ=ξμξνEλEλ−EμEν+ξαEλBβ(ξμϵανλβ+ξνϵαμλβ) −ξλξαBβ(ϵανλβEμ+ϵαμλβEν)+gμνBβBβ−BμBν+ξμξνBβBβ=ξμξνEλEλ+(gμν+ξμξν)BβBβ−EμEν−BμBν+ξαEλBβ(ξμϵανλβ+ξνϵαμλβ)where we have again used the fact that ξαξσϵσαβμ=0.
Now we just need to put all these pieces together. The final expression for the electromagnetic stress-energy tensor for an arbitrary observer is
TμνEM|ξ→n=14π(ξμξνEλEλ+(gμν+ξμξν)BβBβ−EμEν−BμBν+ξαEλBβ(ξμϵανλβ+ξνϵαμλβ)−14gμν(−2EβEβ+2BμBμ))=18π(gμν+2ξμξν)(EαEα+BαBα)−14π(EμEν+BμBν)+14π+ξαEλBβ(ξμϵανλβ+ξνϵαμλβ)Now, we can consider two special cases. First, let's consider the normal observer, in which case the normalized 4-velocity becomes the unit normal vector ξ→n
nν=(1α,−βiα)nν=(−α,0).If we also define
ξαϵανλβ≡ϵνλβthen the stress-energy tensor becomes
TμνEM|ξ→n=18π(gμν+2nμnν)(EαEα+BαBα)−14π(EμEν+BμBν)+14πEαBβn(μϵν)αβWe can also consider the observer co-moving with the fluid. We can call this co-moving velocity uμ. In the case of ideal conditions (that is to say, a fluid with perfect conductivity), we have (from Ohm's law) that
uμFμν=0which implies that
Eμ=uνFμν=0Since Eμ=0, the stress-energy tensor simplifies significantly in this frame:
TμνEM|ξ→u=18π(gμν+2uμuν)BβBβ−14πBμBνHowever, variable b is often used instead of Bμ(u), given by
bμ=Bμ(u)√4πIn terms of b, the electromagnetic energy-momentum tensor is given by
TμνEM|ξ→u=b2(12gμν+uμuν)−bμbνBefore discussing the evolution variables, recall the 3+1 decomposition which is used for GRMHD. The metric is given by /ds2=−α2dt2+γij(dxi+βjdt)(dxj+βidt)
In GRMHD simulations, there are two sets of variables which can be used for evolution. The first set, called primitive variables, are
P={ρb,P,vi,Bi}where ρb is the fluid rest-mass density, P is the pressure, vi is the rescaled fluid three-velocity ui/u0, and Bi is the magnetic field. An important sidenote here is that different GRMHD codes define different vi. As an example, the widely used GRHydro thorn in the Einstein Toolkit uses the Valencia formalism
vi=uiW+βiαwhere W=√1−vivi is the Lorentz factor. The quantity vi used by IllinoisGRMHD is not the Valencia three-velocity, so those interested in seeing more about the Valencia formulation of the evolution equations can look to the GRHydro paper and its references.
The second set, the conservative variables, are
C={˜ρ,˜τ,˜Si,˜Bi}where
˜ρ=√−gρbu0=α√γρbu0˜τ=√γnμnνTμν−˜ρ=α2√γT00−˜ρ˜Si=√γSi=α√γT0i˜Bi=√γBiOne can design numerical codes to evolve either set of variables, and we choose to evolve C. Conservative variables are used to evolve these systems because the evolution of primitive variables leads to highly unstable simulations. This change of variables for the sake of stability has some similarity to the change from ADM variables to BSSN variables when evolving space-time, where one set of variables has substantially better numerical stability.
Another question that may arise is why we do not simply evolve the energy variable α2√γT00 and instead subtract ˜ρ. It turns out that evolving T00 directly is inaccurate because the magnetic and internal energy density can be orders of magnitude smaller than the rest mass density (this is mentioned in e.g. the HARM
software paper). Because of this, we subtract the baryonic density evolution equation from the energy evolution equation to find a more accurate and stable evolution equation.
All of our previous derivations, of course, only considered the electromagnetic contribution. We do not derive the hydrodynamic part of the stress-energy tensor and instead merely state that in the ideal MHD limit the full stress-energy tensor takes the form
Tμν=ρbhuμuν+Pgμν+TμνEMOur next step is to derive the actual evolution equations for our designated evolution variables C. To this end, we will use the identity
Γνγν=1√−g∂γ√−gto simplify the equations. We will be frequently using this to condense two terms into one in the following way:
∂νTνμ+ΓνγνTγμ=∂νTνμ+1√−gTγμ∂γ√−g=∂νTνμ+1√−gTγμ∂γ√−g=1√−g(√−g∂νTνμ+Tγμ∂γ√−g)=1√−g∂γ(√−gTγμ)The next two sections derive their equations from the energy-momentum conservation equation
∇νTνμ=0=∂νTνμ+ΓνγνTγμ−ΓγμνTνγRearranging and using the identity for Γνγν, we get
1√−g∂ν(√−gTνμ)=ΓγμνTνγFor the momentum density variable ˜Si, we consider the indices μ∈{1,2,3}. Replacing μ with i to find the evolution equation,
1√−g∂ν(√−gTνi)=ΓγiνTνγ1√−g∂t(√−gT0i)+1√−g∂j(√−gTji)=12Tνγgγβ(gβν,i+giβ,ν−giν,β)∂t˜Si+∂j(α√γTji)=α√γ2Tβνgβν,i+α√γ2Tβν(giβ,ν−giν,β)∂t˜Si+∂j(α√γTji)=α√γ2Tβνgβν,iThe energy evolution equation comes from the same conservation equation as the previous section, but with μ=0. However, we choose to start with μ raised instead of lowered. Then,
0=∇νTμν=∂νTμν+ΓμσνTσν+ΓνσνTμσ=1√−g∂ν(√−gTμν)+ΓμσνTσνSetting μ=0,
1√−g∂ν(√−gT0ν)=−Γ0σνTσν∂t(α√γT00)+∂i(α√γT0i)=−α√γΓ0σνTσν=−α√γ2Tσνg0β(gβσ,ν+gβν,σ−gσν,β)=−α√γ2Tσνg0β(2gβσ,ν−gσν,β)where the final step follows from the fact that Tσν is symmetric. Now, we will consider the right-hand side for various components of Tσν. For these derivations, we will use derivations relating to the extrinsic curvature in Appendix B. For T00,
−α√γ2T00g0β(2gβ0,0−g00,β)=−α√γ2T00[2(g00g00,0+g0ig0i,0)−g00g00,0−g0ig00,i]=−α√γ2T00(g00g00,0+2g0ig0i,0−g0ig00,i)=−α√γ2T00[−α−2∂t(β2−α2)+2βiα2∂tβi−βiα2∂i(β2−α2)]=−√γ2αT00[−∂t(β2−α2)+2βi∂tβi−βi∂i(β2−α2)]=−√γαT00(α∂tα+βiα∂iα−βiβj∂iβj)=√γT00(βiβjKij−∂tα−βi∂iα)where we have used the identity in Appendix B. Next, we look at the mixed term T0i+Ti0:
−α√γ2T0ig0β(2gβ0,i−g0i,β+2gβi,0−gi0,β)=−α√γT0i[g00(g00,i+g0i,0)+g0j(g0j,i+gij,0)−(g00g0i,0+g0jg0i,j)]=−α√γT0i[g00g00,i+g0j(g0j,i+gij,0−g0i,j)]=−α√γT0i[−α−2∂i(β2−α2)+βjα2(∂iβj+∂tγij−∂jβi)]=−√γαT0i[2α∂iα+βj(+∂tγij−∂iβj−∂jβi)]=−√γαT0i[2α∂iα+βj(−2αKij−2Γkijβk)]=2√γT0i(βjKij−∂iα)+2√γαT0iβjΓkijβk=2√γT0i(βjKij−∂iα)Finally, for Tij we have
−α√γ2Tijg0β(2gβi,j−gij,β)=−α√γ2Tij[2(g00g0i,j+g0kgki,j)−g0βgij,β]=−α√γ2Tij[2(g00∂jβi+g0k∂jγki)−g0β∂βγij]=−α√γ2Tij[g00(2∂jβi−∂tγij)+g0k(2∂jγki−∂kγij)]=−α√γ2Tij[−α−2(2∂jβi−∂tγij)+βkα2(2∂jγki−∂kγij)]We can see that, thanks to the symmetry of Tij,
Tijβk(2∂jγki−∂kγij)=Tijβngkn(∂jγki+∂iγkj−∂kγij)=2TijΓkijβkUsing this and using the definition of Kij to replace ∂tγij,
−α√γ2Tijg0β(2gβi,j−gij,β)=√γ2αTij(2∂jβi+2αKij−∂iβj−∂jβi+2Γkijβk−2Γkijβk)=√γ2αTij(2αKij+∂jβi−∂iβj)=√γTijKijwhere the final step again follows from the symmetry of Tij. Putting all this together and multiplying through by α,
α∂t(α√γT00)+α∂i(α√γT0i)=α√γT00(βiβjKij−∂tα−βi∂iα) +2α√γT0i(βjKij−∂iα)+α√γTijKijα√γT00∂tα+2α√γT0i∂iα+α∂t(α√γT00)+α∂i(α√γT0i)=α√γT00(βiβjKij−βi∂iα)+2α√γT0iβjKij −α√γT0i∂iα+α√γTijKij∂t(α2√γT00)+∂i(α2√γT0i)=α√γ[(T00βiβj+2T0iβj+Tij)Kij−(T00βi+T0i)∂iα]We have finally arrived at the evolution equation. However, we still need to get it in terms of the evolution variable ˜τ=α2√γT00−˜ρ. To do so, we can add the fluid density evolution equation to the left-hand side. We also define the source term
s=α√γ[(T00βiβj+2T0iβj+Tij)Kij−(T00βi+T0i)∂iα]Then, the energy evolution equation becomes
∂t(α2√γT00)+∂i(α2√γT0i)−∂t˜ρ−∂i(˜ρvi)=s∂t˜τ+∂i(α2√γT0i−˜ρvi)=sTo find a convenient start for the derivation, I find the dual of Maxwell's equation. Recall that the dual of the Faraday tensor is
F∗μν=12ϵμναβFαβAlso, since the covariant derivative of the metric is 0,
ϵαλμνFμν;λ=∇λ(ϵαλμνFμν)=2∇λF∗αλTherefore, Maxwell's equation becomes
0=ϵαλμνF[μν;λ]=ϵαλμν(Fμν;λ−Fλν;μ+Fνλ;μ−Fμλ;ν+Fλμ;ν−Fνμ;λ)=ϵαλμνFμν;λ+ϵαλμν∇μ(Fνλ−Fλν)+ϵαλμν∇ν(Fλμ−Fμλ)−ϵαλμν∇λFνμ=2∇λF∗αλ+2ϵαλμν∇μFνλ+2ϵαλμν∇νFλμ−ϵαλμν∇λFνμBy exploiting the symmetries of the Levi-Civita tensor,
0=2∇λF∗αλ+2ϵαμνλ∇μFνλ+2ϵανλμ∇νFλμ+ϵαλνμ∇λFνμ=2∇λF∗αλ+4∇μF∗αμ+4∇νF∗αν+2∇λF∗αλ=12∇λF∗αλ0=∇λF∗αλ=∂λF∗αλ+ΓαβλF∗βλ+ΓλβλF∗αβ0=1√−g∂λ(√−gF∗αλ)where the first Γ term dissapears because it is summing over the multiplication of symmetric and anti-symmetric objects.
We can now find the magnetic equations from this. Taking the time component simply gives us the no-monopole constraint. To see this, we need to first consider the individual components of the dual. By the definition of B,
Bμ=12ϵμνβαnνFαβ=nνF∗νμSince B for the normal observer is purely spatial (Bμnμ=0), this implies that F∗00=0. Then, the time component of the magnetic equations is
0=1√−g∂λ(√−gF∗0λ)=∂i(α√γF∗0i)=∂i(α√γBiα)0=∂i(˜Bi)Before examining the spatial components of Maxwell's equations, we will do some preliminary work to make our lives easier later. In the co-moving frame,
Fμν=uμEν−uνEμ+uγϵγμνδBδ=uγϵγμνδBδThen, the dual is
F∗μν=12ϵμναβFαβ=12ϵμναβuγϵγαβδBδ=12ϵαβμνϵαβγδuγBδ=(δμδδνγ−δμγδνδ)uγBδ=uνBμ−uμBνwhere we have again used the Levi-Civita identity from Appendix A. Next, we need to find a relationship between the magnetic field of the normal observer and the co-moving observer. For clarity, let the co-moving magnetic field be Bμ(u) and the normal magnetic field be Bμ. Then, we can define a projection operator
Pμν=gμν+uμuνNaturally, the projection of the co-moving magnetic field should simply project back into the same field:
PμνBν(u)=(δμν+uμuν)Bν(u)=Bμ(u)where we have used the orthogonality relation uνBν(u)=0. Projecting the normal observer's magnetic field,
PμνBν=PμνnαF∗αν=Pμνnα(uνBα(u)−uαBν(u))=nα(δμν+uμuν)(uνBα(u)−uαBν(u))=nα(uμBα(u)+uμuνuνBα(u)−uαBμ(u)−uμuνuαBν(u))=nαuμBα(u)(1+uνuν)−nαuα(Bμ(u)+uμuνBν(u))=−αu0Bμ(u)⇒Bμ(u)=Bμ+uμuνBναu0Since we will only need the spatial component for our purposes,
Bi(u)=Bi+uiuνBναu0=Bi+uiujBjαu0=Biαu0+viujBjαwhere vi is the conservative variable ui/u0. Finally, the spatial components of the dual of Maxwell's equations gives the induction equation for the magnetic field:
0=1√−g∂ν(√−gF∗iν)0=∂t(α√γF∗i0)+∂j(α√γF∗ij)0=∂t(√γBi)+∂j(α√γ[ujBi(u)−uiBj(u)])0=∂t˜Bi+∂j(α√γ[uj(Biαu0+viukBkα)−ui(Bjαu0+vjukBkα)])0=∂t˜Bi+∂j(vj˜Bi−vi˜Bj+√γu0[ujuiukBk−uiujukBk])0=∂t˜Bi+∂j(vj˜Bi−vi˜Bj)In the previous sections, we have derived the evolution equations for the conservative variables \mathbf{C}. To summarize, these are
∂t˜ρ+∂i(˜ρvi)=0∂t˜τ+∂i(α2√γT0i−˜ρvi)=s∂t˜Si+∂j(α√γTji)=α√γ2Tβνgβν,i∂t˜Bi+∂j(vj˜Bi−vi˜Bj)=0where
s=α√γ[(T00βiβj+2T0iβj+Tij)Kij−(T00βi+T0i)∂iα]In order to simplify the various combinations of Fμν, contractions of the Levi-Civita tensor are required. First, summing over the first indices of two such tensors yields
ϵαβσμϵαβγδ=−δβσμβσνβγδ=−δββδσμγδ+δβγδσμβδ−δβδδσμβγ=−4(δσγδμδ−δσδδμγ)+δβγ(δσβδμδ−δσδδμβ)−δβδ(δσβδμγ−δσγδμβ)=−4(δσγδμδ−δσδδμγ)+2(δσγδμδ−δσδδμγ)=2(δσδδμγ−δσγδμδ)Second, summing over the third index (with all other components raised) yields
ϵγμλδϵανλβανλ=gατgνϕgβθϵγμλδϵτϕλθ=−gατgνϕgβθδγμδτϕθ=−gατgνϕgβθ(δγτδμδϕθ−δγϕδμδτθ+δγθδμδτϕ)=−gατgνϕgβθ(δγτ(δμϕδδθ−δμθδδϕ)−δγϕ(δμτδδθ−δμθδδτ)+δγθ(δμτδδϕ−δμϕδδτ))The source term of the energy equation involves the extrinsic curvature. Defined in terms of the 3+1 formalism quantities, the extrinsic curvature is Kij=−12α(∂tγij−∇iβj−∇jβi)=−12α(∂tγij−∂iβj−∂jβi+2Γkijβk)
As for how the extrinsic curvature enters into the energy equations, it involves several relations between it and the other quantities which we must show. First, consider the contraction of the Christoffel symbol with βiβj:
βjΓkijβk=βjβk(γki,j+γkj,i−γij,k)=βjβkγkj,i=α2(γjk−gjk)γkj,i∝γjk∂iγkj−gjk∂iγkj∝γjk∂iγkj−γjk∂iγkj∝γjk∂iγkj−γjk∂iγkj=0The change gjk→γjk is allowed because the 3-metric is identical to the 4-metric with only spatial indices. We can use this to find a relationship between the partial derivative of βj and the extrinsic curvature:
βiβj∂iβj=βiβj(2αKij+∂tγij−∂jβi+2Γkijβk)=βiβj(2αKij+∂tγij−∂iβj+2Γkijβk)=βiβj(2αKij+∂tγij−∂iβj)⇒βiβj∂iβj=βiβj(αKij+12∂tγij)⇒βiβj∂iβj=βiβjαKijwhere in the second step we use the symmetry of βiβj. The final step uses the same trick as with βjΓkijβk, just with the derivative being ∂t instead of ∂i.
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Derivation_of_GRMHD_Evolution_Equations")
Created Tutorial-Derivation_of_GRMHD_Evolution_Equations.tex, and compiled LaTeX file to PDF file Tutorial- Derivation_of_GRMHD_Evolution_Equations.pdf