Hyperbolicity of Divergence Cleaning and Vector Potential Formulations of GRMHD

TL;DR

This paper analyzes the hyperbolicity of divergence cleaning and vector potential formulations in general relativistic magnetohydrodynamics, demonstrating their strong hyperbolicity and suitability for well-posed initial value problems in numerical relativity.

Contribution

It shows that divergence cleaning and vector potential formulations are strongly hyperbolic, making them preferable for numerical simulations in GRMHD.

Findings

Divergence cleaning formulation is strongly hyperbolic.

Vector potential formulations can be made strongly hyperbolic.

Preferred formulations for numerical relativity are divergence cleaning and vector potential.

Abstract

We examine hyperbolicity of general relativistic magnetohydrodynamics with divergence cleaning, a flux-balance law form of the model not covered by our earlier analysis. The calculations rely again on a dual-frame approach, which allows us to effectively exploit the structure present in the principal part. We find, in contrast to the standard flux-balance law form of the equations, that this formulation is strongly hyperbolic, and thus admits a well-posed initial value problem. Formulations involving the vector potential as an evolved quantity are then considered. Carefully reducing to first order, we find that such formulations can also be made strongly hyperbolic. Despite the unwieldy form of the characteristic variables we therefore conclude that of the free-evolution formulations of general relativistic magnetohydrodynamics presently used in numerical relativity, the divergence cleaning and vector potential formulations are preferred.

Tables

Table 1. Table 1: Overview of the uppercase and lowercase quantities.

	Uppercase	Lowercase
Unit normal	$N^a=W\left(n^a+v^a\right)$	$n^a=W\left(N^a+V^a\right)$
Boost vector	$V^a$	$v^a={\hat{v}}^a/W$
Lorentz factor	$W={\left(1-V^a{V}_a\right)}^{-1/2}$	$W={\left(1-v^a{v}_a\right)}^{-1/2}$
Projector	${}_{}{}^{\text{(N)}}\gamma _{}^a{}_b{}^{}=g^a{}_b{}^{}+N^a{N}_b$	${\gamma }^a{}_b{}^{}=g^a{}_b{}^{}+n^a{n}_b$
Boost metric	${}_{}{}^{\text{(N)}}𝕘_{ab}^{}:={}_{}{}^{\text{(N)}}\gamma _{ab}^{}+W^2{V}_a{V}_b$	${𝕘}_{ab}:={\gamma }_{ab}+{\hat{v}}_a{\hat{v}}_b$
Inverse boost	${}_{}{}^{\text{(N)}}\left({𝕘}^{-1}\right)_{}^{ab}={}_{}{}^{\text{(N)}}\gamma _{}^{ab}-V^a{V}^b$	${\left({𝕘}^{-1}\right)}^{ab}={\gamma }^{ab}-v^a{v}^b$

Table 2. Table 2: Summary of the various unit spatial vectors appearing in our  2 + 1 2 1 2+1 decomposed equations, plus their associated projection operators.

	Uppercase	Lowercase
Unit normal	$N^a$	$n^a$
Spatial 1-form	$S_a$	${𝕤}_a$
Spatial vector	$S^a={}_{}{}^{\text{(N)}}\gamma _{}^{ab}{S}_b$	${\widehat{𝕤}}^a={\left({𝕘}^{-1}\right)}^{ab}{𝕤}_b$
Norm	$S_a{S}^a=1$	${𝕤}_a{\left({𝕘}^{-1}\right)}^{ab}{𝕤}_b=1$
Projector	${}_{}{}^{\text{Q}}⟂_a^b={}_{}{}^{\text{(N)}}\gamma _{}^b{}_a{}^{}-S^b{S}_a$	${}_{}{}^{𝕢}⟂_a^b={\gamma }_a^b-{\widehat{𝕤}}^b{𝕤}_a$
Index notation	${}_{}{}^{\text{Q}}⟂_A^B$	${}_{}{}^{𝕢}⟂_{\widehat{𝔸}}^{𝔹}$

Equations

$$n^a=W(N^a+V^a),N^a=W(n^a+v^a),$$

$${\gamma }^b{}_a={\delta }^b{}_a+n^b{n}_a,{}^{\text{(N)}}{\gamma }^b{}_a={\delta }^b{}_a+N^b{N}_a,$$

$${\nabla }_u\mathbf{U}={\mathbf{A}}^p{\nabla }_p\mathbf{U}+\mathcal{S},$$

$${\mathbf{P}}^S\equiv {\mathbf{A}}^p{S}_p.$$

$$T_{\text{fluid}}^{ab}={\rho }_{0}h{u}^a{u}^b+g^{ab}p,$$

$$T_{\text{em}}^{ab}=F^{ac}{F}^b{}_c-\frac{1}{4}{g}^{ab}{F}_{cd}{F}^{cd},$$

$$h=1+\epsilon +\frac{p}{{\rho }_0}.$$

$$c_s^2=\frac{1}{h}\left(\chi +\frac{p}{{\rho }_0^2}\kappa \right),\chi ={\left(\frac{\partial p}{\partial {\rho }_0}\right)}_{\epsilon },\kappa ={\left(\frac{\partial p}{\partial \epsilon }\right)}_{{\rho }_0}.$$

$$p=p({\rho }_0,\epsilon ),$$

$$F^{ab}$$

$${}^{\ast }{F}^{ab}$$

$${ϵ}^{abcd}=-\frac{1}{\sqrt{-g}}\left[abcd\right],$$

$$T^{ab}={\rho }_0{h}^{\ast }{u}^a{u}^b+p^{\ast }{g}^{ab}-b^a{b}^b,$$

$${\nabla }_a({\rho }_0{u}^a)=0,$$

$${\nabla }_b{T}^{ab}=0,$$

$${\nabla }_b({{}^{\ast }F}^{ab}-g^{ab}\varphi )=-\frac{1}{\tau }{n}^a\varphi ,$$

$${\nabla }_a({\rho }_0{u}^a)=0,{}^{\text{(u)}}{\gamma }_{ab}{\nabla }_c{T}^{bc}=0,$$

$${}^{\text{(u)}}{\gamma }_{ab}{\nabla }_c({{}^{\ast }F}^{bc}-g^{bc}\varphi )=-\frac{W}{\tau }{V}_a\varphi ,$$

$$u_b{\nabla }_c{T}^{bc}=0,u_b{\nabla }_c({{}^{\ast }F}^{bc}-g^{bc}\varphi )=\frac{W}{\tau }\varphi ,$$

$${\nabla }_{u}p=$$

$${}^{\text{(u)}}{\gamma }_{ab}$$

$$+\left(\frac{2}{{\rho }_0{h}^{\ast }}{}^{\text{(u)}}{\gamma }^{[b}{}_a{b}^{p]}{}^{\text{(u)}}{\gamma }_{bc}+\frac{b_a}{{\rho }_{0}h}{}^{\text{(u)}}{\gamma }^p{}_c\right)(\text{mathbbmss}{g}^{-1}{)}^{ce}{\nabla }_p\perp b_e$$

$$+S_a^{(\hat{\mathbf{v}})},$$

$${}^{\text{(u)}}{\gamma }_{ab}(\text{mathbbmss}{g}^{-1}{)}^{bc}{\nabla }_u\perp b_c=$$

$$-{}^{\text{(u)}}{\gamma }^p{}_a{\nabla }_p\varphi +S_a^{(\perp \mathbf{b})},$$

$${\nabla }_u\epsilon =$$

$${\nabla }_u\varphi =$$

$$S^{(p)}=-c_s^{2}W{\rho }_{0}h{}^{\text{(u)}}{\gamma }^d{}_c(\text{mathbbmss}{g}^{-1}{)}^{ce}{\nabla }_d{n}_e-\frac{\kappa W}{\tau {\rho }_0}(b^a{V}_a)\varphi ,$$

$$S_a^{(\hat{\mathbf{v}})}=-W{}^{\text{(u)}}{\gamma }_{ab}(\text{mathbbmss}{g}^{-1}{)}^{be}{\nabla }_u{n}_e+\frac{2W}{{\rho }_0{h}^{\ast }}{}^{\text{(u)}}{\gamma }^{[b}{}_a{b}^{e]}{V}_b{b}^d{\nabla }_d{n}_e$$

$$+\frac{1}{{\rho }_{0}h}{b}_a\left(W{V}^d{b}^e-W(b^c{V}_c){}^{\text{(u)}}{\gamma }^{de}\right){\nabla }_d{n}_e,$$

$$S_a^{(\perp \mathbf{b})}=2W{}^{\text{(u)}}{\gamma }_{ab}{}^{\text{(u)}}{\gamma }^{[b}{}_c{b}^{d]}(\text{mathbbmss}{g}^{-1}{)}^{ce}{\nabla }_d{n}_e$$

$$+2W{}^{\text{(u)}}{\gamma }_{[a}^e{V}_{b]}{b}^b{\nabla }_u{n}_e+\frac{W}{\tau }{V}_a\varphi ,$$

$$S^{(\epsilon )}=-\frac{Wp}{{\rho }_0}{}^{\text{(u)}}{\gamma }^d{}_c(\text{mathbbmss}{g}^{-1}{)}^{ce}{\nabla }_d{n}_e-\frac{W}{\tau {\rho }_0}(b^a{V}_a)\varphi ,$$

$$S^{(\varphi )}=-\left(W{V}^d{b}^e-W(b^c{V}_c){}^{\text{(u)}}{\gamma }^{de}\right){\nabla }_d{n}_e-\frac{W\varphi }{\tau }.$$

$${}^{\text{(u)}}{\gamma }_{ac}(\text{mathbbmss}{g}^{-1}{)}^{cd}{\nabla }_b\perp b_d:=$$

$$+V_a{b}_d(\text{mathbbmss}{g}^{-1}{)}^{de}{\nabla }_b{\hat{v}}_e.$$

$${\mathbf{B}}^{\text{u}}{\nabla }_u\mathbf{U}={\mathbf{B}}^p{\nabla }_p\mathbf{U}+\mathcal{S}.$$

Full text

I Introduction

It is well appreciated BauSha10x ; Shi16 that the numerical modeling of binary neutron star spacetimes plays, and will continue to play, an important role in the new field of gravitational wave astronomy, particularly in the case of multimessenger events. Such simulations are, however, hampered by relatively poor error behavior as compared with their vacuum, black hole counterparts. This is in part because the equations of motion of these models have a more complicated structure than those of pure general relativity, and are hence less well understood, but also because solutions naturally develop nonsmooth features, not to mention the ever-present complication of the stellar surface.

In a recent contribution SchHilBug17 we employed a new tool, the dual-frame (DF) formalism HilRic13 ; Hil15 ; HilRui16 ; HilHarBug16 , to analyze well-posedness of various fluid models. Well-posedness is the weakest necessary condition to require of a set of evolution partial differential equations (PDE) so that numerical approximation to their solutions may be meaningfully sought. The formalism can be used to exploit structure in the field equations and hence simplifies earlier treatments. This should allow more sophisticated results to be shown in the future.

One of the models treated in Ref. SchHilBug17 was (ideal) general relativistic magnetohydrodynamics (GRMHD), taken in two different guises. In the Valencia flux-balance law form AntZanMir05 we found that the field equations are only weakly hyperbolic, and therefore have an ill-posed initial value problem. Here we attend to two flavors of GRMHD untouched by our earlier study, namely the hyperbolic divergence cleaning (HDC) and vector potential (VP) formulations. Our main result is that both are strongly hyperbolic, provided suitable choices are made in the first-order reduction of the latter.

We work in $3+1$ dimensions in geometric units with $c=G=1$. Our calculations were performed primarily with xTensor for Mathematica xAct_web_aastex ; our notebooks are available online in Ref. HilSch18_WebTar .

II Mathematical background

We start with a short overview of the relevant theory, definitions, and results to the PDE analysis and the DF formalism. These are taken in a highly summarized form from Refs. Hil15 ; HilHarBug16 ; SchHilBug17 .

Index notation.

Latin letters $a$–$e$ are used as abstract indices. We also use $p$ as an abstract index, placing it always on the spatial derivative appearing on the right-hand side of our first-order PDE system. The four-metric $g_{ab}$ is the only object permitted to raise and lower indices. The symbol $\partial_{a}$ stands for a flat covariant derivative. Indices $u$, $S$, $s$, ${\hat{\mathbbmss{s}}}$ and ${\mathbbmss{s}}$ label contraction in that slot with $u^{a}$ or $u_{a}$ and so on, respectively. Capital Latin letters $A$–$C$ are taken as abstract indices and denote appliance of the projection operators ${}^{\textrm{{\tiny{Q}}}}\!\!\!\perp$ and ${}^{{\tiny{q}}}\!\!\!\perp$, to be defined later. Similarly, we use indices ${\mathbbmss{A}}$–${\mathbbmss{C}}$ and ${\hat{\mathbbmss{A}}}$–${\hat{\mathbbmss{C}}}$ to denote the application of the projection operator ${}^{{\tiny{\mathbbmss{q}}}}\!\!\!\perp$ over a vector or dual vector, respectively.

DF formalism.

We describe a region of spacetime in two different frames, namely the lowercase and the uppercase frame. We take the lowercase frame as an Eulerian frame, associated with a coordinate basis as is standard in numerical relativity. We denote the future pointing timelike unit normal vector to spatial slices of constant time, as usual, by $n^{a}$. Additionally, we take any three linearly independent vector fields orthogonal to $n^{a}$ to form a basis of the four-dimensional spacetime. Tensors orthogonal to $n^{a}$ are called lowercase spatial, or just lowercase. The uppercase frame consists of a future pointing timelike unit vector $N^{a}$, which is identified in the application below with the fluid four velocity $u^{a}$, plus any three linearly independent vector fields orthogonal to $N^{a}$. Tensors orthogonal to $N^{a}$ are likewise called uppercase spatial, or just uppercase. The future pointing unit vectors of the lower- and uppercase frames can be mutually $3+1$ decomposed as

[TABLE]

with the Lorentz factor $W=(1-V^{a}V_{a})^{-1/2}=(1-v^{a}v_{a})^{-1/2}$. The vectors $v^{a}=\hat{v}^{a}/W$ and $V^{a}$ are the boost vectors orthogonal to $n^{a}$ and $N^{a}$, respectively. We define projection operators by

[TABLE]

which are also denoted as the lowercase and uppercase spatial metrics, respectively. By definition, the relations $\gamma^{b}{}_{a}n_{b}=0$, ${}^{\textrm{\tiny{(N)}}}\!\gamma^{b}{}_{a}N_{b}=0$ hold. We define furthermore the lowercase and uppercase boost metrics and their inverses, which are presented in Table 1.

PDE analysis.

We consider a quasilinear system of first-order evolution PDEs, in this case GRMHD with HDC, written in the form

[TABLE]

with the covariant derivative along the streamlines of the fluid elements $\nabla_{u}\equiv u^{a}\nabla_{a}$ of the vector of evolved variables, called the state vector $\mathbf{U}$, on the left-hand side. On the right-hand side, the covariant derivative of the state vector is contracted with the principal part $\mathbf{A}^{p}$, $\mathbf{A}^{a}u_{a}=0$. The symbol $\boldsymbol{\mathcal{S}}$ stands for the source term which does not affect the level of hyperbolicity. We need only analyze the system of evolution equations for the matter variables, since they are minimally coupled to the Einstein equations for the components of the metric tensor.

Strong hyperbolicity.

For the hyperbolicity analysis, we have to perform a $2+1$ decomposition against lowercase and/or uppercase spatial vectors and their respective orthogonal spatial projectors. The relevant quantities are defined in Table 2. Taking an arbitrary uppercase unit spatial 1-form $S_{a}$, we define the uppercase principal symbol of the system (3) as

[TABLE]

We call the system (3) weakly hyperbolic, if for each $S_{a}$ the eigenvalues of $\mathbf{P}^{S}$ are real. We call the system (3) strongly hyperbolic, if the system is weakly hyperbolic and for each $S_{a}$ the principal symbol $\mathbf{P}^{S}$ has a complete set of right eigenvectors written as columns in a matrix $\mathbf{T}_{S}$ and there exists a constant $K>0$, independent of $S_{a}$, such that $\mathbf{|}\mathbf{T}_{S}|+|\mathbf{T}_{S}^{-1}|\leq K$. Similar definitions are made if we $3+1$ decompose the system against $n^{a}$ rather than $u^{a}$, and the initial value problem, where data are given at $t=0$, can be well-posed only if it satisfies these lowercase strong hyperbolicity conditions GusKreOli95 ; SarTig12 ; Hil13 .

Frame and variable independence of hyperbolicity.

If the uppercase eigenvalues of the principal symbol fulfill the inequality $|\lambda_{\text{N}}||V|<1$ then strong hyperbolicity is independent of the chosen frame SchHilBug17 . By the form of the energy-momentum tensor of GRMHD, see below, a naturally preferred frame is the fluid rest frame. Therefore, in the PDE analysis in Sec. IV, we work exclusively in the uppercase frame, taken to be the fluid rest frame, $N^{a}\equiv u^{a}$; hence the $3+1$ decomposition in Eq. (3), and in the following, of the equations against the fluid four velocity $u^{a}$ and the orthogonal projector ${}^{\textrm{\tiny{(u)}}}\!\gamma^{a}{}_{b}$. In numerical applications, particular sets of variables, such as the primitive or conservative sets are used. In our analysis, we make a choice of variables which differs slightly from those. Our variables are however related to the code variables by a regular transformation, across which hyperbolicity is unaffected.

III Basics of GRMHD

A brief review of the basic definitions, equations, and assumptions of GRMHD with HDC is now given, following Refs. Ani90a ; AntMirMar10 ; MoeMunFab14 . Presently, the focus lies on the mathematical structure of the equations, and thus we suppress some (important) physical insight and statements. We use Lorentz-Heaviside units for electromagnetic quantities with $\varepsilon_{0}=\mu_{0}=1$, where $\varepsilon_{0}$ and $\mu_{0}$ are the vacuum permittivity (or electric constant) and permeability (or magnetic constant), respectively. Motivated by the arguments given in the previous section, we work exclusively in the uppercase (fluid) frame and thus, present the system of equations in a form so adjusted.

The energy-momentum tensor of GRMHD consists of an ideal fluid part,

[TABLE]

with the four velocity of the fluid elements $u^{a}$, rest mass density $\rho_{0}$, specific enthalpy $h$, and pressure $p$; plus the standard electromagnetic energy-momentum tensor

[TABLE]

with the Faraday electromagnetic tensor field (or for short field strength tensor) $F^{ab}$. The specific enthalpy $h$ can be expressed in terms of $\rho_{0},\ p,$ and the specific internal energy $\varepsilon$ as

[TABLE]

The local speed of sound $c_{s}$ is defined by the relation

[TABLE]

We assume an equation of state (EOS) of the form

[TABLE]

with $p>0$ given, satisfying furthermore that the local speed of sound lies in the range $0<c_{s}\leq 1$.

Using the ideal MHD condition, where the electric conductivity tends to infinity while the electric four-current remains bounded, the field strength tensor and its dual become

[TABLE]

respectively, where we introduced the uppercase magnetic field vector $b^{a}$, satisfying $u_{a}b^{a}=0$, and the Levi-Cività tensor

[TABLE]

where $g$ is the determinant of the spacetime metric $g_{ab}$, $\left[abcd\right]$ is the completely antisymmetric Levi-Cività symbol, and $2{{}^{*}F^{ab}}=-\epsilon^{abcd}F_{cd}$ holds. Note that we use the sign convention of Ref. AlcDegSal09 .

Taking the sum of Eqs. (5) and (6), and substituting the field strength tensor (10), the total energy-momentum tensor of GRMHD may be written as

[TABLE]

with $h^{*}=h+b^{2}/\rho_{0}$, $p^{*}=p+b^{2}/2$, and shorthand $b^{2}=b^{a}b_{a}$.

The covariant system of evolution equations is given by the conservation of the number of particles and the conservation of energy momentum,

[TABLE]

plus the relevant Maxwell equations

[TABLE]

which are already augmented by the terms with the scalar field $\phi$ to drive the Gauss constraint. Since $b^{a}$ has only three free components this equation now gives an evolution equation for $b^{a}$ and $\phi$. Elsewhere the notation $\kappa=\tau^{-1}$ is employed. The constant $\tau$ is the timescale for the exponential driving toward the Gauss constraint for the magnetic field. Typically $\phi$ is set to 0 in the initial and boundary conditions DedKemKro02 .

IV Hyperbolicity analysis of GRMHD with HDC

Projecting Eqs. (14)-(16) along the four velocity of the fluid $u^{a}$ and perpendicular to it by ${}^{\textrm{\tiny{(u)}}}\!\gamma^{a}{}_{b}$, the nine evolution equations which determine the time evolution of the GRMHD system with HDC are

[TABLE]

supplemented with an EOS (9). In the limit of $\phi\rightarrow 0$ we find the uppercase Gauss constraint: ${}^{\textrm{\tiny{(u)}}}\!\gamma^{bc}\nabla_{b}b_{c}=u_{c}\nabla_{b}{{}^{*}F}^{bc}=0$.

Taking Eq. (17) and performing algebraic manipulations similar to the investigation of other formulations of GRMHD in Ref. SchHilBug17 , we derive the evolution equations for the pressure,

[TABLE]

the boost vector,

[TABLE]

the magnetic field,

[TABLE]

the specific internal energy,

[TABLE]

and finally the scalar field variable,

[TABLE]

The sources are given by

[TABLE]

The auxiliary magnetic vector $\bot{b}_{c}$ is defined by the relation

[TABLE]

As usual, square brackets around indices denote antisymmetrization, so that $2\hat{v}^{[a}b^{b]}=\hat{v}^{a}b^{b}-\hat{v}^{b}b^{a}$. We have shown explicitly that the set of equations (18)-(22) is, up to nonprincipal terms, which we have not carefully checked, simply a linear combination of the formulation of GRMHD with HDC used in numerical applications, see, for example, Ref. MoeMunFab14 . This verification can be found in the notebook that accompanies the paper HilSch18_WebTar .

Writing Eqs. (18)-(22) in a vectorial form with state vector $\mathbf{U}=(p,\hat{v}_{a},\bot{b}_{a},\varepsilon,\phi)^{T}$, we obtain, in the notation of Ref. SchHilBug17 , the principal part in the form

[TABLE]

Let $S_{a}$ be an arbitrary unit spatial uppercase $1$-form, $S_{a}S^{a}=1$, and ${{}^{\textrm{{\tiny{Q}}}}\!\!\!\perp^{b}\!{}_{a}}:={}^{\textrm{\tiny{(u)}}}\!\gamma^{b}{}_{a}-S^{b}S_{a}$ be the associated orthogonal projector. Let ${\mathbbmss{s}}_{a}$ and ${{}^{{\tiny{\mathbbmss{q}}}}\!\!\!\perp^{b}\!{}_{a}}$ be their lowercase projected versions, ${\mathbbmss{s}}_{a}=\gamma^{b}{}_{a}S_{b}$, ${{}^{{\tiny{\mathbbmss{q}}}}\!\!\!\perp^{b}\!{}_{a}}:=\gamma^{b}{}_{a}-(\mathbbmss{g}^{-1})^{bc}{\mathbbmss{s}}_{c}{\mathbbmss{s}}_{a}$. Decomposing ${}^{\textrm{\tiny{(u)}}}\!\gamma^{b}{}_{a}$ and $\gamma^{b}{}_{a}$ against $S_{a}$ and ${\mathbbmss{s}}_{a}$, respectively, Eq. (24) can be written as

[TABLE]

where $\simeq$ denotes equality up to nonprincipal terms and uppercase spatial derivatives transverse to $S^{a}$. The uppercase principal symbol is $\mathbf{P}^{S}=\mathbf{B}^{S}=$

[TABLE]

with the state vector ordered as,

[TABLE]

The characteristic polynomial $P_{\lambda}$ for the principal symbol (26) is calculated to

[TABLE]

with the quadratic polynomial for Alfvén waves

[TABLE]

and the quartic polynomial for the magnetosonic waves

[TABLE]

Comparing Eq. (30) with our earlier results for the flux-balance law formulation of GRMHD in Ref. SchHilBug17 , we see that the linear polynomial associated with the Gauss constraint is replaced by the quadratic polynomial $1-\lambda^{2}$. The entropy, Alfvén, and slow and fast magnetosonic uppercase eigenvalues remain the same, as before, and are given by

[TABLE]

respectively, with shorthands

[TABLE]

The remaining two speeds can be associated with the scalar field and the longitudinal magnetic field MoeMunFab14 , and are given by

[TABLE]

Since all uppercase eigenvalues have absolute value smaller than or equal to 1, the relation $|\lambda_{\text{u}}||V|<1$ is satisfied, so we may analyze hyperbolicity independently of the frame SchHilBug17 . Therefore, we analyze the characteristic structure of the principal symbol in the uppercase frame and the result of the analysis applies directly to the numerically used system (in the lowercase).

Continuing the characteristic analysis, we find the left entropy, scalar field, and longitudinal magnetic field, Alfvén, and magnetosonic eigenvectors being

[TABLE]

respectively, where we defined the antisymmetric uppercase two- and three-Levi-Cività tensors as ${{}^{\textrm{{\tiny{(S)}}}}\!\epsilon^{{A}{B}}=S_{d}{}^{\textrm{{\tiny{(u)}}}}\!\epsilon^{d{A}{B}}=u_{c}S_{d}\,{}^{\textrm{{\tiny{Q}}}}\!\!\!\perp^{{A}}\!_{a}{}^{\textrm{{\tiny{Q}}}}\!\!\!\perp^{{B}}\!_{b}\epsilon^{cdab}}$. We employ furthermore the shorthands

[TABLE]

The right eigenvectors can be computed and are presented in the same order,

[TABLE]

We have introduced in the magnetosonic eigenvectors the orthogonal magnetic field vector as $b_{\perp}^{a}={}^{\textrm{{\tiny{Q}}}}\!\!\!\perp_{\ b}^{a}b^{b}$ with $b_{\perp}^{2}=b_{\perp}^{a}b^{\perp}_{a}=b^{A}b_{A}$. As for the prototype algebraic constraint free formulation treated in Ref. SchHilBug17 , rescaled versions of the left eigenvectors (34) and right eigenvectors (36) can be derived. They form complete sets of nine linearly independent eigenvectors under type I, type II, and type II*′* degeneracies AntMirMar10 ; Sch18 . The rescaling can be found in the notebook provided in Ref. HilSch18_WebTar . Thus, as long as $p=p(\rho_{0},\varepsilon)>0$ and $0<c_{s}<1$ hold, the formulation of GRMHD with HDC as given above forms a strongly hyperbolic system of equations.

In the limit of $c_{s}\rightarrow 1$, it can be shown, that the fast magnetosonic waves collide pairwise with the waves associated to the scalar field and longitudinal magnetic field, in the case of which the system is only weakly hyperbolic. The limiting procedure can be found in the provided notebook. This is a consequence of taking the divergence cleaning to happen at the speed of light. By the simple replacement $\phi\rightarrow c_{\phi}^{-2}\phi$, $c_{\phi}>0$ in Eq. (22), the divergence cleaning speed becomes $\lambda_{\pm}=\pm c_{\phi}$. For $c_{\phi}>1$, strong hyperbolicity is also guaranteed in the limiting case $c_{s}=1$. This strategy does however place a nontrivial upper limit on the speed of flows that can be managed with the method, as strong hyperbolicity will break down for sufficiently fast flows. See Ref. SchHilBug17 for details. By modifying the lowercase equations directly it may be possible to avoid this shortcoming, too.

Finally, we want to present the uppercase rescaled characteristic variables for GRMHD with HDC. They are valid for all degeneracies, and are given by

[TABLE]

with $\{\text{m}_{1},\text{m}_{2}\}$ equal to $\{\text{s},\text{f}\}$ or $\{\text{f},\text{s}\}$. The abbreviations in Eq. (37) are given by

[TABLE]

where for type II and even for type II*′* degeneracy we take ${{Q}_{1}}^{a}$ and ${{Q}_{2}}^{a}$ such that in the degenerate limit we have

[TABLE]

For further explanations concerning degeneracies and rescaling, see also Ref. Sch18 .

Using the recovery procedure given in Ref. SchHilBug17 , the lowercase characteristic quantities such as eigenvalues and eigenvectors can be derived. The calculation can be found in the notebook HilSch18_WebTar , but results in rather long expressions which we suppress here. Both the lowercase left magnetosonic eigenvectors and the lowercase right eigenvectors associated with the scalar field and longitudinal magnetic field eigenvalues have a particularly complicated structure, for which a useful simplification seems difficult. In applications it may therefore be appropriate to compute the characteristics numerically.

V Discussion of formulations of GRMHD with VP

The formulations of GRMHD we have thus far considered use the magnetic field as an evolved variable. Another possibility is to introduce the four-vector potential instead GiaRezBai10 ; EtiPasLiu12 ; EtiPasHaa15 . In practice, the potential is then $3+1$ decomposed. Such formulations have the advantage that the Gauss constraint is satisfied by construction, and in this sense can be considered a type of constrained- rather than free-evolution. On the other hand one obtains a system of equations which is a priori not, from the PDE point of view, minimally coupled to the gravitational field equations. The resulting evolution equations for the GRMHD variables are moreover themselves not in first-order form, but rather first order in time and second order in space, and there is an additional gauge degree of freedom. Different choices in this freedom may have different PDE properties as the principal part of the evolution system is altered. We follow Ref. EtiPasLiu12 and focus on the Lorenz gauge, but similar comments hold elsewhere. Strong hyperbolicity of first order in time, second order in space systems can be defined GunGar05 ; HilRic13a by the requirement that there exists a first-order reduction which satisfies the definition given for first order PDEs in Sec. II. Therefore, we must reduce the governing system of equations as in Eq. (3), by introducing reduction variables. There are two natural ways to go about this.

The first, naive, possibility is to introduce reduction constraints $c_{ab}=d_{ab}-\gamma^{c}{}_{a}\gamma^{d}{}_{b}\partial_{c}A_{d}$, which should vanish, for the lowercase spatial derivatives of the lowercase spatial part of the vector potential $A_{a}$, and likewise for the electric potential. The reduction variables $d_{ab}$ should satisfy also the ordering constraint,

[TABLE]

and similarly for the electric potential reduction variables. The reduction constraints must then be added to the equations of motion to remove all second spatial derivatives. Besides that, both the reduction and ordering constraints can be added freely to try and find a hyperbolic reduction. Such a reduction does not use the special structure of the Maxwell equations, does not utilize the fact that the original system satisfies the Gauss constraint by construction, and is not minimally coupled to the evolution equations for the geometric variables. Worse, the resulting principal symbol does not have a clear structure, which makes the analysis very difficult.

The less obvious option is to bring back the magnetic field as a reduction variable for the curl of the spatial vector potential by defining a reduction constraint,

[TABLE]

In this reduction we need not introduce a reduction variable to the electric potential as it appears with at most one spatial derivative. Part of the analog of the ordering constraint in such a reduction turns out to be simply the Gauss constraint,

[TABLE]

A generic PDE system does not allow a reduction of this type, in which new variables that only capture part of the spatial derivatives are introduced. Due to the gauge freedom of the Maxwell equations however the ‘longitudinal’ part of the vector potential does not appear elsewhere in the remaining equations of motion, and so we can close the evolution system using only $B_{a}$. Note that such a restricted reduction does have consequences on the norms in which rigorous estimates would be demonstrated, and also that as usual first derivatives of the metric here are nonprincipal.

Ultimately we end up with evolution equations for the matter variables which are minimally coupled to the Einstein equations. Naively writing out the lowercase principal symbol of the matter variables we can obtain moreover a block-diagonal structure,

[TABLE]

where block $\mathbf{A}$ denotes the principal symbol of the system of evolution equations of the spatial part of the vector potential and the electric potential, whereas $\mathbf{B}$ can be rendered identical to the principal symbol of the prototype algebraic constraint free formulation of GRMHD investigated in Ref. SchHilBug17 . Here, crucially, we rely on the fact that, as it is not to be used in applications, this formal first-order reduction need not be of a flux-balance form, and therefore we can add the ordering constraint $C$ as desired. The upper right block vanishes trivially and the lower left block vanishes by appropriate choice of reduction. We showed already that prototype algebraic constraint free formulation of GRMHD is strongly hyperbolic in the lowercase frame, with an EOS of the form (9) and $0<c_{s}\leq 1$, so all that remains is to show that the block $\mathbf{A}$ satisfies the conditions for strong hyperbolicity. This was done already in Ref. EtiPasLiu12 , but with the use of the reduction variable $B_{a}$ we can give a slightly cleaner treatment. The lowercase principal symbol can be read off from,

[TABLE]

Note that in Eq. (48) the term $D_{a}A_{b}-D_{b}A_{a}$ is written in terms of the reduction variable $B_{a}$ and does not contribute to the principal part. Let $s^{a}$, $s_{a}s^{a}=1$, be a unit spatial lowercase vector and $q^{a}{}_{b}=\gamma^{a}{}_{b}-s^{a}s_{b}$ be the orthogonal projector. The characteristic variables associated with this block are hence

[TABLE]

with speeds $\pm 1$, respectively, and

[TABLE]

with speed [math] for the two orthogonal directions to $s^{a}$. The calculation is provided in a notebook that accompanies the paper HilSch18_WebTar .

VI Conclusion

In previous work SchHilBug17 we examined two formulations of ideal GRMHD, and showed that a formulation similar to that studied in Refs. AniPen87 ; Ani90a , which we call the prototype algebraic constraint free formulation is strongly hyperbolic. Unfortunately, this formulation is not in the flux-balance law form desirable for the application of standard numerical methods. Turning to GRMHD in flux-conservative form, we found the system to be only weakly hyperbolic. This formulation of GRMHD hence has an ill-posed initial value problem. Fortunately, two popular, applicable, alternative formulations of GRMHD were left untreated by that analysis. Presently, we have addressed this shortcoming with the outcome first, that formulations of GRMHD with HDC DedKemKro02 ; AndHirLie06 ; MoeMunFab14 are indeed strongly hyperbolic as long as the sound speed is suitably bounded $0<c_{s}<1$. In fact, it is straightforward to achieve hyperbolicity also in the case $c_{s}=1$ by changing the speed of the cleaning in the formulation. Second, we have shown that by a careful reduction to first-order, formulations of GRMHD with VP EtiPasLiu12 can also be rendered strongly hyperbolic whenever $0<c_{s}\leq 1$. The latter result is a corollary of strong hyperbolicity of the prototype algebraic constraint free formulation. Here we have discussed only the Lorenz gauge choice, but our results carry over trivially to generalized Lorenz gauge, in which there is a modification by source terms, and a natural treatment will be very similar in other cases, too.

Both HDC and the VP formulations were introduced as strategies to control Gauss-constraint violation in applications. Another popular approach, called constrained transport (CT) EvaHaw88 ; BalSpi99 ; Tot00 , uses a carefully constructed discretization so that in a particular approximation the constraint is identically satisfied. There is some subtlety in precisely what continuum PDE should be analyzed given such a constrained evolution, but supposing that the constraints are identically satisfied, they may again be added arbitrarily to the evolution equations, and strong hyperbolicity can again be achieved, in the restricted, constraint-satisfying phase space, as a corollary of hyperbolicity of the prototype algebraic constraint free formulation.

In Ref. SchHilBug17 we discussed two minimally coupled formulations of resistive GRMHD with HDC, one with and one without the evolution of the charge density $q$. Both were found to be only weakly hyperbolic. A natural question is therefore whether the use of the VP approach could cure this problem. Replacing the divergence cleaning variables by $A_{a}$ and $\Phi$, and making a minimally coupled first-order reduction as we did for GRMHD, one arrives with a lower block triangular structure in the principal symbol, with the lower-right block $\mathbf{C}$ being precisely a sub-block of the principal symbol of the original formulation of RGRMHD. Neither of the original two formulations were strongly hyperbolic because $\mathbf{C}$ was not diagonalizable. Consequently, the vector potential formulations are also not strongly hyperbolic. Thus, at least if we insist on taking only minimally coupled first-order reductions, use of a VP reformulation of RGRMHD does nothing to circumvent weak hyperbolicity of RGRMHD.

For numerical applications we therefore have the clear conclusion that, by the fundamental requirement of well-posedness, HDC and VP formulations (and likely also CT schemes) are preferred over their older variant which should henceforth be avoided. From the PDE point of view it is, at this stage, difficult to choose between the favored formulations. One might be tempted to argue in favor of the vector potential formulation, as indeed it is true that there the characteristic structure, inherited from the prototype algebraic constraint free formulation, is simpler, but this is not a principle advantage. In the future it is hoped that the characteristic structure uncovered by our analysis can be put to good use in numerical work in both systems.

Acknowledgements.

We are grateful to Sebastiano Bernuzzi and Bruno Giacomazzo for useful discussions and comments. This work was partially supported by the FCT (Portugal) IF Program No. IF/00577/2015 and the GWverse COST action Grant No. CA16104.