# Non-linear electrodynamics non-minimally coupled to gravity:
symmetric-hyperbolicity and causal structure.

Érico Goulart  Federal University of São João d’El Rei, C.A.P. Rod.: MG 443,
KM 7, CEP-36420-000, Ouro Branco, MG, Brazil  Santiago Esteban Perez
Bergliaffa  Departamento de Física Teórica, Instituto de Física, Universidade
do Estado de Rio de Janeiro, 20550-013, Rio de Janeiro, Brazil

###### Abstract

It is shown here that symmetric hyperbolicity, which guarantees well-
posedness, leads to a set of two inequalities for matrices whose elements are
determined by a given theory. As a part of the calculation, carried out in a
mostly-covariant formalism, the general form for the symmetrizer, valid for a
general Lagrangian theory, was obtained. When applied to nonlinear
electromagnetism linearly coupled to curvature, the inequalities lead to
strong constraints on the relevant quantities, which were illustrated with
applications to particular cases. The examples show that non-linearity leads
to constraints on the field intensities, and non-minimal coupling imposes
restrictions on quantities associated to curvature.

##  I Introduction

Well-posedness of the initial value problem stands out as a basic requirement
to be satisfied by any relativistic field theory. Broadly speaking, it implies
that solutions for a given problem exist, are unique, and depend continuously
on the initial data. Hence, well-posedness is at the roots of physics, for it
amounts to the predictability power of a given theory. While it is difficult
to decide whether a given non-linear theory has a well-posed initial value
problem, a necessary and sufficient condition for well-posedness around a
given solution is that all the linearised problems obtained by linearising
near such a solution are well-posed, see for instance Kreiss and Lorenz
(1989); Papallo and Reall (2017).

In the linearized regime, well-possedness means that the equation of motion is
hyperbolic. At least three notions of hyperbolicity can be distinguished
Kreiss and Lorenz (1989); Beig (2006): (a) weak hyperbolicity, in which all
the roots of the characteristic equation are real, (b) strong hyperbolicity
implies that there is an energy estimate that sets a bound for the energy of a
solution at a given time in terms of the initial energy 111In fact, strong
hyperbolicity is a necessary and sufficient condition for the initial-value
problem to be well-posed., and (c) symmetric hyperbolicity, which is a
sufficient condition for well-posedness Sarbach and Tiglio (2012). It follows
that symmetric hyperbolicity implies strong hyperbolicity, which in turn
implies weak hyperbolicity.

There are many examples in which the requirement of some kind of hyperbolicity
imposes severe restrictions on the Lagrangian of a given theory. For instance,
as shown in Papallo and Reall (2017), the equations of motion of Lovelock’ss
theory are always weakly hyperbolic for weak fields but not strongly
hyperbolic in a generic weak-field background. The well-posedness of Horneski
theory has been analyzed in Papallo and Reall (2017); Papallo (2017). As shown
in the latter reference, the most general Horndeski theory which is strongly
hyperbolic for weak fields in a generalized harmonic gauge is simply _k-
essence_ coupled to Einstein’s gravity (see also Bernard  _et al._ (2019)).
The well-posedness of scalar-tensor effective field theory was studied in
Kovács and Reall (2020), where it was shown that the equations of motion are
strongly hyperbolic at weak coupling.

We would like to explore here the constraints imposed by the requirement of
well-posedness in the case of nonlinear electromagnetic theories coupled to
gravity. Nonlinear electromagnetism has been widely studied in several
contexts. A non-exhaustive list of applications and references includes black
holes Moreno and Sarbach (2003); Bretón and Perez Bergliaffa (2015); Toshmatov
_et al._ (2018); Magos and Bretón (2020); Falciano  _et al._ (2021),
astrophysics Heyl and Hernquist (2004); Harding and Lai (2006); Turolla  _et
al._ (2015) and cosmology Novello  _et al._ (2004, 2007); Campanelli  _et al._
(2008); Vollick (2008); Kruglov (2015). There are also several articles
devoted to different aspects of the propagation of perturbations in nonlinear
electromagnetic theories, such as Gutierrez  _et al._ (1981); Deser  _et al._
(1999); Goulart de Oliveira Costa and Esteban Perez Bergliaffa (2009); Perlick
_et al._ (2018). The matter of (symmetric) hyperbolicity for nonlinear
electromagnetism minimally coupled to gravity was analyzed in a flat spacetime
background in Abalos  _et al._ (2015), while the hyperbolicity of Maxwell’s
equations with a local (and possibly nonlinear) constitutive law in flat
spacetime was considered in Perlick (2011). Our aim here is to determine the
restrictions that follow from the requirement of symmetric hyperbolicity on
nonlinear electromagnetism non-minimally coupled to gravity, several aspects
of which have been studied in detail in Balakin and Lemos (2005); Balakin  _et
al._ (2008, 2010); Dereli and Sert (2011). We shall restrict to couplings
linear in the curvature, since higher-order couplings produce higher-than-
second order equations for the gravitational field Balakin and Lemos (2005).

As mentioned above, a sufficient condition for well-posedness to hold is that
the system under study admits a symmetric-hyperbolic representation. The
theory of first-order symmetric hyperbolic systems, originally due to
Friedrichs Friedrichs (1954), has been extensively developed (see for instance
Beig (2006) and references therein). To study the evolution of a system we
shall adopt here the modern geometric approach to the subject outlined by
Geroch in Hall  _et al._ (1996), in which covariance is kept during most of
the calculation, instead of using a 3+1 decomposition of spacetime (as for
instance in Perlick (2011)) 222For a list of other approaches to the analysis
of the evolution of minimally coupled electromagnetism, see Abalos  _et al._
(2015).. We shall see that such approach leads to a general form for the
symmetrizer, and to conditions for symmetric hyperbolicity that are easier to
evaluate than those for other types of hyperbolicity.

The structure of the paper is as follows. In Section II the basic notation
used in the equations of motion is presented. Symmetric hyperbolicity, the
related concept of symmetrizer, and the role of the constraints are analyzed
in Section III. The equations defining the characteristic cones will be
deduced in Section IV. Our main result, namely the explicit form of the
matrices that are needed to investigate the symmetric hyperbolicity of any
nonlinear electromagnetic theory non-minimally coupled to gravity can be found
in Section V. In Section VI, some examples of the restrictions imposed by
symmetric hyperbolicity are presented for different theories. Our closing
remarks are presented in Section VII.

##  II Lagrangians and equations of motion

To begin with, let $\mathcal{M}$ denote a smooth four-dimensional spacetime
with a Lorentzian metric $\boldsymbol{g}$ of signature $(+,-,-,-)$. For the
sake of concreteness we assume $\mathcal{M}$ to be also oriented and globally
hyperbolic i.e. $M\cong\mathbb{R}\times\Sigma$, with $\Sigma$ a
codimension-$1$ hypersurface. Sticking to the conventions $R_{ab}\equiv
R^{c}_{\phantom{a}acb}$ and $[ab]\equiv ab-ba$, the canonical decomposition of
the Riemman tensor into its irreducible parts reads

$R^{ab}_{\phantom{a}\phantom{a}cd}=W^{ab}_{\phantom{a}\phantom{a}cd}+\frac{1}{2}\delta^{[a}_{\phantom{a}[c}S^{b]}_{\phantom{a}d]}+\frac{R}{12}g^{ab}_{\phantom{a}\phantom{a}cd},$
(1)

where $W_{abcd}$ is the Weyl conformal tensor, $S_{ab}$ is the traceless part
of the Ricci tensor, $R$ is the scalar curvature and $g_{abcd}=g_{a[c}g_{d]b}$
is the Kulkarni-Nomizu product of the metric with itself. Clearly, each factor
in the decomposition has the same algebraic symmetries as the full Riemann
tensor, that is

$R_{abcd}=-R_{bacd}=-R_{abdc},\quad\quad\quad R_{abcd}=R_{cdab},$ (2)
$R_{abcd}+R_{acdb}+R_{adbc}=0.$ (3)

An arbitrary rank four covariant tensor satisfying Eqs. (2) has 21 independent
components and is often referred to as a double symmetric $(2,2)$-form (the
skew pairs can be interchanged). If, in addition, the tensor satisfies Eq.
(3), it is called an algebraic curvature tensor (see for instance Besse
(1987)).

If we are given a double symmetric $(2,2)$-form on $\mathcal{M}$, say $\chi$,
we may obtain its left-dual $\star\chi$ and right-dual $\chi\star$ by

$\star\chi_{abcd}=\frac{1}{2}\varepsilon_{ab}^{\phantom{a}\phantom{a}pq}\chi_{pqcd},\quad\quad\quad\chi_{abcd}\star=\frac{1}{2}\varepsilon^{pq}_{\phantom{a}\phantom{a}cd}\chi_{abpq},$
(4)

where $\varepsilon_{abcd}$ is the Levi-Civita tensor with
$\varepsilon_{0123}=\sqrt{-g}$. Notice that the place of the $\star$ indicates
the pair of skew indices which are Hodge-dualized. A direct consequence is
that

$\star\star\chi_{abcd}=-\chi_{abcd},\quad\quad\quad\chi_{abcd}\star\star=-\chi_{abcd}.$
(5)

We aim at constructing Lagrangians describing the coupling of gravity with the
electromagnetic field using invariants that are at most linear in the
curvature 333As discussed in Balakin and Lemos (2005), higher-order couplings
lead to equations of motion with derivatives of the metric higher than two.,
and respecting the $U(1)$ gauge invariance of electromagnetism. Recalling that
every algebraic curvature tensor satisfy the Ruse-Lanczos identity, we may
construct the following independent rank four tensors (see the Appendix for
details)

${}^{(1)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv
g^{ab}_{\phantom{a}\phantom{a}cd},\quad\quad\quad^{(2)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv\varepsilon^{ab}_{\phantom{a}\phantom{a}cd},$
(6) ${}^{(3)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv
Rg^{ab}_{\phantom{a}\phantom{a}cd},\quad\quad\quad^{(4)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv
R\varepsilon^{ab}_{\phantom{a}\phantom{a}cd},$ (7)
${}^{(5)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv
W^{ab}_{\phantom{a}\phantom{a}cd},\quad\quad\quad^{(6)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv\star
W^{ab}_{\phantom{a}\phantom{a}cd},$ (8)
${}^{(7)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\equiv\delta^{[a}_{\phantom{a}[c}S^{b]}_{\phantom{a}d]}.$
(9)

Notice that ${\phantom{a}}^{(\Gamma)}\chi_{abcd}$ $(\Gamma=1,2,...,7)$ are
double symmetric $(2,2)$-forms by construction and it can be checked that they
exhaust the interesting possibilities: the first pair, Eqs. (6), involves only
algebraic terms in the metric, while the remaining Eqs. (7)-(9) contain at
most linear terms in the curvature.

Let $F_{ab}=\partial_{[a}A_{b]}$ denote a test electromagnetic field
propagating on $\mathcal{M}$. We shall focus on the propagation properties of
such a field on a fixed gravitational background. The EM field will be
described by the gauge-invariant action functional of the form

$S=\frac{1}{4}\int d^{4}x\sqrt{-g}\ \mathcal{L}(I^{1},I^{2},...,I^{7}),$ (10)

where the factor $1/4$ is introduced for future convenience and the Lagrangian
density is taken as an arbitrary smooth function of the following scalars

$I^{\Gamma}\equiv{\phantom{a}}^{(\Gamma)}H^{ab}F_{ab},\quad{\rm
where}\quad\quad{\phantom{a}}^{(\Gamma)}H^{ab}\equiv\frac{1}{2}{\phantom{a}}^{(\Gamma)}\chi^{ab}_{\phantom{a}\phantom{a}cd}F^{cd}.$
(11)

In analogy with electrodynamics in material media, we call
${\phantom{a}}^{(\Gamma)}H^{ab}$ the $\Gamma$-th induction tensor and
${\phantom{a}}^{(\Gamma)}\chi^{ab}_{\phantom{a}\phantom{a}cd}$ the $\Gamma$-th
constitutive tensor. Particular instances described by Eq. (10) include
Maxwell’s theory, minimally-coupled nonlinear electrodynamics, and the three-
parameter non-minimal Einstein–Maxwell model originating from QED vacuum
polarization in a background gravitational field (see for instance Drummond
and Hathrell (1980)), among others. With the conventions presented above, the
variation of the action with respect to the 4-potential yields a coupled
system of first-order quasi-linear PDEs for the fields, given by

$\nabla_{b}H^{ab}=0,\quad\quad\quad\nabla_{b}\star F^{ab}=0.$ (12)

Here, $\nabla_{a}$ is the covariant derivative compatible with the metric and
the full induction tensor is defined by

$H^{ab}\equiv\displaystyle\sum\mathcal{L}_{\Gamma}{\phantom{a}}^{(\Gamma)}H^{ab},$
(13)

with $\mathcal{L}_{\Gamma}\equiv\partial\mathcal{L}/\partial I^{\Gamma}$, for
conciseness. Notice that the system defined by Eqs. (12) is composed of eight
equations for only six unknowns. In relevant physical situations, the
equations must include six dynamical equations (defined with respect to some
time coordinate, to be identified later) and two constraints, which are to be
imposed on initial data.

##  III Symmetric-hyperbolicity

In order to study whether the equations of motion admit a symmetric-hyperbolic
representation, it is convenient to recast the system of equations Eqs. (12)
in a unified manner as

$K_{A\phantom{a}\beta}^{\phantom{a}m}(x,\Phi)\partial_{m}\Phi^{\beta}+J_{A}(x,\Phi)=0,$
(14)

where $x\in{\mathcal{M}}$, $K_{A\phantom{a}\beta}^{\phantom{a}m}$ is the
principal part of the PDE and $J_{A}(x,\Phi)$ stands for semi-linear
contributions (whose explicit form is unnecessary for our discussion). Here
capital Latin indices $(A=1,...,8)$ stand for the space of multi-tensorial
equations, lowercase Latin indices $(m=0,1,2,3)$ stand for space-time indices,
and greek indices $(\beta=1,..,6)$ for tensorial unknowns. To start with, we
introduce an ordering of the antisymmetric indices to obtain the six possible
collective quantities

$1\rightarrow(01)\quad 2\rightarrow(02)\quad 3\rightarrow(03)\quad
4\rightarrow(32)\quad 5\rightarrow(13)\quad 6\rightarrow(21).$ (15)

Making the identification $\Phi^{\beta}\rightarrow F^{bc}$ and performing
simple manipulations in Eqs. (12), the principal part is then written as

$K_{A\phantom{a}\alpha}^{\phantom{a}m}=\frac{1}{2}\left(X_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m},\
\varepsilon_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m}\right),$ (16)

with

$X_{abcd}\equiv\sum_{\Gamma}\mathcal{L}_{\Gamma}{\phantom{a}}^{(\Gamma)}\chi_{abcd}+4\sum_{\Gamma}\sum_{\Lambda}\mathcal{L}_{\Gamma\Lambda}{\phantom{a}}^{(\Gamma)}H_{ab}{\phantom{a}}^{(\Lambda)}H_{cd}.$
(17)

$X_{abcd}$ consists of a main term involving only first partial derivatives of
the Lagrangian density and a nonlinear term including the Hessian matrix of
the latter. Clearly, if the Lagrangian is a linear combination of the
invariants $I^{\Gamma}$, the last term vanishes and linear equations of motion
follow. More importantly, $X_{abcd}$ is always a symmetric double (2,2)-form
independently on the specific form of the Lagrangian. This is a direct
consequence of the symmetries of the constitutive tensors together with the
symmetry of the Hessian matrix, i.e.
$\mathcal{L}_{\Gamma\Lambda}=\mathcal{L}_{\Lambda\Gamma}$.

In what follows we shall use the covariant approach for first-order symmetric-
hyperbolic systems outlined in Hall  _et al._ (1996): a symmetric
hyperbolization of Eqs. (14) means that there exist a smooth symmetrizer
$h^{A}_{\phantom{a}\alpha}$ and a covector field $n_{m}$, such that:

  1. 1.

$\hat{K}_{\alpha\phantom{a}\beta}^{\phantom{a}m}\equiv
h^{A}_{\phantom{a}\alpha}K_{A\phantom{a}\beta}^{\phantom{a}m}$ is symmetric in
the indices $\alpha,\ \beta$.

  2. 2.

The matrix
$\hat{K}_{\alpha\beta}(n)\equiv\hat{K}_{\alpha\phantom{a}\beta}^{\phantom{a}m}n_{m}$
is positive-definite.

Roughly speaking, the first statement means that it should be possible to
construct from Eqs. (14) a new subsystem of first-order quasi-linear PDEs
given by

$\hat{K}_{\alpha\phantom{a}\beta}^{\phantom{a}m}\partial_{m}\Phi^{\beta}+\hat{J}_{\alpha}=0$
(18)

with $\hat{J}_{\alpha}\equiv h^{A}_{\phantom{a}\alpha}J_{A}$. Clearly, the new
system contains only evolution equations and its dimension is equivalent to
the number of unknown fields. The second statement means that the new system
can be solved uniquely for any given set of initial conditions on a
hypersurface $\Sigma$ with normal covector $n_{m}$. What are the remaining
equations the symmetrizer does not capture? For the system to be consistent
they should not be of the evolution type. In other words, they must be
satisfied automatically once they are satisfied initially i.e., they should be
what we normally call the constraints444See references Hall  _et al._ (1996);
Beig (2006); Kato (1975) for additional details..

The geometrical meaning of the covectors introduced in the previous paragraph
is as follows. At a spacetime point $p\in\mathcal{M}$, the collection of all
covectors satisfying condition (2) is denoted by $S_{p}$. This set defines a
nonempty open, convex cone at $p$ and the tangent vectors $p^{a}\in
T_{p}\mathcal{M}$ such that $p^{a}n_{a}>0$ for all $n_{a}\in S_{p}$ determine
the cone of influence of the physical field, i.e., the maximal speed of
propagation in any given direction. It turns out that the latter is also a
nonempty open, convex cone at $p$.

In the next subsection, we start by showing that a family of symmetrizers
parametrized by a vector field always exists for the system of first order
PDEs given by Eqs. (12). It is important to point out that our result is
general in the sense that it does not depend on the specific form of the
Lagrangian density. We then obtain the conditions that a theory must satisfy
for positive-definiteness (see condition (2) above) to hold. This is first
obtained for a particular covector and henceforth generalized by inspecting
the characteristic varieties (dispersion relations), which are necessarily
well-behaved for symmetric hyperbolic systems.

###  III.1 Symmetrizer

In order to find a symmetrizer, it is convenient to work with projections. In
other words, we seek a multi-tensorial field $h^{A}_{\phantom{a}\alpha}$ such
that the quantity
$\delta\phi^{\alpha}(h^{A}_{\phantom{a}\alpha}K_{A\phantom{a}\beta}^{\phantom{a}m})\delta{\psi}^{\beta}$
is symmetric in $\delta\phi$ and $\delta{\psi}$. Making the identifications

$\delta\phi^{\alpha}\rightarrow
A^{ab},\quad\quad\quad\delta{\psi}^{\alpha}\rightarrow B^{ab},$ (19)

where $A$ and $B$ are generic bivectors, we obtain from Eq. (16), the relation

$K_{A\phantom{a}\beta}^{\phantom{a}m}\delta\psi^{\beta}=\left(\frac{1}{2}X_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m}B^{bc},\
\star B_{a}^{\phantom{a}m}\right).$ (20)

Since there is no known practical procedure to obtain a symmetrizer for an
arbitrary system of PDE’s, hyperbolizations are found, for a sufficiently low
number of dimensions, by solving explicitly the algebraic equations inherent
to the system defined by Eq. (14) and, in higher dimensions, by guessing. Let
us show that the symmetrizer is given by the projection

$\delta\phi^{\alpha}h^{A}_{\phantom{a}\alpha}=\left(A^{a}_{\phantom{a}q}\ ,\
\frac{1}{2}\star X^{a}_{\phantom{a}qrs}A^{rs}\right)t^{q}$ (21)

where $t^{q}$ is an auxiliary vector field and $\star X^{a}_{\phantom{a}qrs}$
is the left Hodge dual as defined before. Indeed, multiplying (21) by (20) one
obtains

$\delta\phi^{\alpha}(h^{A}_{\phantom{a}\alpha}K_{A\phantom{a}\beta}^{\phantom{a}m})\delta{\psi}^{\beta}=\frac{1}{2}\left(X_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m}A^{a}_{\phantom{a}q}B^{bc}+\star
X^{a}_{\phantom{a}qrs}A^{rs}\star B_{a}^{\phantom{a}m}\right)t^{q}$ (22)

Now, defining

$Y^{a}_{\phantom{a}q}\equiv X^{a}_{\phantom{a}qrs}A^{rs},\quad\quad\quad\star
Y^{a}_{\phantom{a}q}\equiv\star X^{a}_{\phantom{a}qrs}A^{rs},$ (23)

and using the well known identity valid for anti-symmetric tensors

$(\star Y^{aq})(\star
B_{am})=-\frac{1}{2}(Y_{ln}B^{ln})\delta_{m}^{\phantom{a}q}+Y_{am}B^{aq},$
(24)

it follows that

$\delta\phi^{\alpha}(h^{A}_{\phantom{a}\alpha}K_{A\phantom{a}\beta}^{\phantom{a}m})\delta{\psi}^{\beta}=\frac{1}{2}\left[X_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m}(A^{a}_{\phantom{a}q}B^{bc}+A^{bc}B^{a}_{\phantom{a}q})-\frac{1}{2}(X_{abcd}B^{ab}A^{cd})\delta^{m}_{\phantom{a}q}\right]t^{q},$
(25)

which is obviously symmetric in $A$ and $B$ since $X_{abcd}=X_{cdab}$.
Therefore, the symmetrizer itself is given by the simple expression

$h^{A}_{\phantom{a}\alpha}=\frac{1}{2}(g^{a}_{\phantom{a}qrz},\star
X^{a}_{\phantom{a}qrz})t^{q},$ (26)

It is remarkable that Eq. (26) is valid independently of the specific content
of the tensor field $X_{abcd}$: in particular it does not matter whether the
equations of motion contain quasi-linear terms or not.

The application of Eq. (26) to Eq. (16), yields the object

$\hat{K}_{\alpha\phantom{a}\phantom{a}\beta}^{\phantom{a}m}=\frac{1}{4}(g^{a}_{\phantom{a}qrz}X_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m}+\star
X^{a}_{\phantom{a}qrz}\varepsilon_{a\phantom{a}\phantom{a}bc}^{\phantom{a}m})t^{q},$
(27)

which, after straightforward algebraic manipulations, becomes

$\hat{K}_{\alpha\phantom{a}\phantom{a}\beta}^{\phantom{a}m}=-\frac{1}{4}\left(g_{q[a}X_{b]\phantom{a}\phantom{a}cd}^{\phantom{a}m}+g_{q[c}X_{d]\phantom{a}\phantom{a}ab}^{\phantom{a}m}+\delta^{m}_{\phantom{a}q}X_{abcd}\right)t^{q}.$
(28)

Notice that this equation is indeed symmetric in the exchange of antisymmetric
indices $ab\Leftrightarrow cd$ and that the auxiliary vector field $t^{q}$
remains (up to now) arbitrary, so we can use it at our disposal. This
concludes the first step of our task.

###  III.2 Positive definitness

Let us next investigate whether the above symmetrizer constitutes a true
hyperbolization of the equations of motion. This will be achieved if we manage
to find a covector field $n_{m}$ such that the characteristic matrix

$\hat{K}_{\alpha\beta}(t,n)=-\frac{1}{4}[t_{[a}X_{b]\phantom{a}\phantom{a}cd}^{\phantom{a}m}+t_{[c}X_{d]\phantom{a}\phantom{a}ab}^{\phantom{a}m}+t^{m}X_{abcd}]n_{m}$
(29)

is positive definite, _i_.e.

$\delta\phi^{\alpha}\hat{K}_{\alpha\beta}\delta\phi^{\beta}>0,$ (30)

for all nonzero vectors $\delta\phi^{\alpha}\in\mathbb{R}^{6}$. A direct
calculation gives the equivalent inequality

$-\left[t_{a}X_{b\phantom{a}\phantom{a}cd}^{\phantom{a}m}A^{ab}A^{cd}+\frac{1}{4}t^{m}X_{abcd}A^{ab}A^{cd}\right]n_{m}>0,$
(31)

where $A^{ab}$ is an arbitrary, but nonzero bivector. Notice that the latter
is a fully covariant expression in the sense that it does not depend on any
particular choice of coordinates. It will restrict, however, the choices of
admissible covector and auxiliary vector fields.

In order to proceed, we shall decompose $A^{ab}$ and $X_{abcd}$ with respect
to the auxiliary vector field $t^{q}$. To do so, assuming that $t^{q}$ is
timelike, future-directd and normalized, we write the bivector as

$A_{ab}=\mathfrak{a}_{[a}t_{b]}+\varepsilon_{ab}^{\phantom{a}\phantom{a}cd}t_{c}\mathfrak{b}_{d},$
(32)

with

$\mathfrak{a}_{a}\equiv A_{ab}t^{b},\quad\quad\mathfrak{b}_{a}\equiv\star
A_{ab}t^{b},\quad\quad\mathfrak{a}_{a}t^{a}=0\quad\quad\mathfrak{b}_{a}t^{a}=0.$
(33)

Similarly, for any double $(2,2)$-form, we have555See for instance Senovilla
(2000) for similar decompositions.

$X_{abcd}=\left\\{g_{abpq}\left(g_{cdrs}\mathbb{P}^{pr}+\varepsilon_{cdrs}\mathbb{Q}^{pr}\right)+\varepsilon_{abpq}\left(g_{cdrs}\mathbb{R}^{pr}+\varepsilon_{cdrs}\mathbb{S}^{pr}\right)\right\\}t^{q}t^{s},$
(34)

with the 2-index tensors given by

$\mathbb{P}_{ab}\equiv X_{acbd}t^{c}t^{d},\quad\mathbb{Q}_{ab}\equiv-
X_{acbd}\star t^{c}t^{d},\quad\mathbb{R}_{ab}\equiv-\star
X_{acbd}t^{c}t^{d},\quad\mathbb{S}_{ab}\equiv\star X_{acbd}\star t^{c}t^{d},$
(35)

and orthogonal to $t^{q}$ by construction. In general, each of the latter has
9 independent components since the number of independent components of
$X_{abcd}$ is 36. However, in the particular case in which
$X_{abcd}=X_{cdab}$, the following simplified relations are valid:

$\mathbb{P}_{ab}=\mathbb{P}_{ba},\quad\quad\mathbb{S}_{ab}=\mathbb{S}_{ba},\quad\quad\mathbb{Q}_{ab}=\mathbb{R}_{ba},$
(36)

so that $\mathbb{P}_{ab}$ and $\mathbb{S}_{ab}$ have 6 independent components
each, while $\mathbb{Q}_{ab}$ (or, equivalently, $\mathbb{R}_{ab}$) have the
remaining 9 components. Notice that, in the context of electrodynamics in
material media, these tensors would be related to the permittivity,
permeability and magneto-electric cross terms of the medium Balakin and Zayats
(2018).

Using the above decompositions in Eq. (31), the new covariant inequality

$(n_{m}t^{m})[\mathbb{P}_{pr}\mathfrak{a}^{p}\mathfrak{a}^{r}-\mathbb{S}_{pr}\mathfrak{b}^{p}\mathfrak{b}^{r}]>2(n_{m}t^{n})\varepsilon^{mq}_{\phantom{a}\phantom{a}np}[\mathbb{R}_{qr}\mathfrak{a}^{p}\mathfrak{a}^{r}+\mathbb{S}_{qr}\mathfrak{a}^{p}\mathfrak{b}^{r}]$
(37)

follows, which is to be satisfied for any $3$-vectors $\mathfrak{a}^{p}$ and
$\mathfrak{b}^{q}$, not vanishing simultaneously. Furthermore, since our
considerations are essentially algebraic we may restrict to an arbitrary point
$p\in\mathcal{M}$. In order to complete our symmetric hyperbolization, it
suffices to find a specific covector $n_{a}\in T^{*}_{p}\mathcal{M}$
satisfying inequality (37) since, if this is the case, there will be a unique
connected, open and convex cone, $S_{p}$, containing the initial covector, and
this cone will exhaust all possibilities. The natural choice here is the
covector $n_{a}=g_{ab}t^{b}$, since the right-hand-side of Eq. (37) will
vanish. Now, since $n_{m}t^{m}>0$ for the latter, we obtain the equivalent
inequalities

$\mathbb{P}_{pr}\mathfrak{a}^{p}\mathfrak{a}^{r}>0,\quad\quad\quad\mathbb{S}_{pr}\mathfrak{b}^{p}\mathfrak{b}^{r}<0,$

for all nonzero vectors $\mathfrak{a}^{p},\mathfrak{b}^{q}$. In words, a
symmetric hyperbolization is achieved for a covector $n_{a}$ if the 2-index
tensors $\mathbb{P}_{ab}$ and $\mathbb{S}_{ab}$, obtained from $X_{abcd}$ and
$\star X_{abcd}\star$ via contractions with $n^{a}$, satisfy the above
inequalities.

For the sake of completness, let us a repeat the calculations in an adapted
frame at $p$, such that

$g_{ab}(p)=\eta_{ab},\quad\quad\quad t^{q}=\delta^{q}_{\phantom{a}0}.$ (38)

Using Eq. (34), it can be shown that $X_{abcd}$ has the following block matrix
display:

$X_{\alpha\beta}=\left(\begin{array}[]{cc}\mathbb{P}&\mathbb{Q}\\\ \\\
\mathbb{Q}^{T}&\mathbb{S}\end{array}\right),$ (39)

where $\mathbb{P}$, $\mathbb{Q}$ and $\mathbb{S}$ denote the covariant
$3\times 3$ matrices constructed with the corresponding tensors in the obvious
way. In order to compute Eq. (29) in matrix form it is convenient to define
the auxiliary $3\times 3$ matrices (a similar notation was used in Perlick
(2011))

$\textbf{A}^{1}=\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&1\\\
0&-1&0\end{array}\right),\quad\quad\textbf{A}^{2}=\left(\begin{array}[]{ccc}0&0&-1\\\
0&0&0\\\
1&0&0\end{array}\right),\quad\quad\textbf{A}^{3}=\left(\begin{array}[]{ccc}0&1&0\\\
-1&0&0\\\ 0&0&0\end{array}\right).$

Taking into account that the transposition relation
$(\textbf{A}^{k})^{T}=-\textbf{A}^{k}$ holds, a direct calculation gives:

$\hat{K}_{\alpha\beta}(n)=\left(\begin{array}[]{cc}n_{0}\mathbb{P}-n_{k}(\mathbb{Q}\textbf{A}^{k}-\textbf{A}^{k}\mathbb{Q}^{T})&n_{k}\textbf{A}^{k}\mathbb{S}\\\
\\\ -n_{k}\mathbb{S}\textbf{A}^{k}&-n_{0}\mathbb{S}\end{array}\right).$ (40)

The above relation may be thought as a linear map from $T^{*}_{x}\mathcal{M}$
to the 21-dimensional space of symmetric $6\times 6$ matrices $\mbox{Sym}_{6}$
. The set of all symmetric positive definite matrices forms an open convex
cone in $\mbox{Sym}_{6}$ with apex on the origin. It turns out that the image
of the particular covector $n_{m}=\eta_{mn}t^{n}$ will lie inside this cone
whenever

$\mathbb{P}\succ 0,\quad\quad\quad\quad\mathbb{S}\prec 0.$ (41)

since $n_{m}=(1,0,0,0)$ in our frame. Therefore, if Eqs. (41) are satisfied
pointwisely, symmetric hyperbolicity is guaranteed for the corresponding
covector field. In other words, if the auxiliary vector field $t^{a}$ is
vorticity-free, then initial data given on a hypersurface $\Sigma$ with normal
covectors $t_{a}$ are uniquely evolved away from the hypersurface. As will be
shown in Sect. VI, the requirements set by Eqs.(41) may yield severe
constraints on the Lagrangian density, the intensity of the curvature tensor
or the electromagnetic field.

###  III.3 Constraints

We now discuss the remaining equations, which are left aside by the
symmetrizer. To do so, we recall that, acording to Geroch’s formalism, a
constraint is a tensor $c^{An}$ such that

$c^{An}K_{A\phantom{a}\alpha}^{\phantom{a}m}+c^{Am}K_{A\phantom{a}\alpha}^{\phantom{a}n}=0.$
(42)

It is straightforward to show that, in our case, this tensor is given by

$c^{An}=(xg^{an},yg^{an})$ (43)

where $x,y\in\mathbb{R}$. Indeed, multiplying Eq. (43) by (16) one obtains

$c^{An}K_{A\phantom{a}\alpha}^{\phantom{a}m}=\frac{1}{2}(xX^{nm}_{\phantom{a}\phantom{a}\phantom{a}bc}+y\varepsilon^{nm}_{\phantom{a}\phantom{a}\phantom{a}bc})$
(44)

which is obviously antisymmetric in $n$ and $m$. Here, the emergence of two
real numbers reveals that the vector space of constraints is actually
2-dimensional, as expected.

From the above calculation one concludes that the constraints are complete, in
the sense that the number of constraint equations plus the number of evolution
equations equals the number of initial equations. The constraints are
integrable if the equation

$c^{An}\nabla_{n}(K_{A\phantom{a}\beta}^{\phantom{a}m}\nabla_{m}\Phi^{\beta}+J_{A})=0,$
(45)

is identically satisfied solely due the algebraic structure of the principal
symbol, independently of the equations (14) of the original system. We leave
for the reader to verify that this is indeed the case. This means that the
original system of equations is equivalent to a symmetric-hyperbolic one with
two additional integrable constraints.

##  IV Characteristic cones

Suppose that we manage to find a symmetric hyperbolization with constraints,
as described in the previous sections. We have seen that this choice, however,
is far from unique. Indeed, by the continuity of Eq.(40), any small
deformation of $t_{a}$ will be such that condition (2) of symmetric
hyperbolization is satisfied. It turns out that the set of all admissible
covectors $S_{p}$ in $T_{p}^{*}\mathcal{M}$ is determined by the unique
connected, open, convex, positive cone containing the initial covector Beig
(2006) . Its existence is related to the hyperbolicity of the characteristic
polynomial, defined by

$p(n)\equiv\mbox{det}\left(\hat{K}_{\alpha\beta}(n)\right)$ (46)

which is a homogeneous multivariate polynomial of degree $6$ in our case.
Recall that, at a spacetime point, such a polynomial is called hyperbolic in a
direction $t_{a}$ if $p(t_{a})>0$ and the univariate polynomial
$p(u_{a}+\lambda t_{a})$ only has real roots for all covectors $u_{a}\neq
t_{a}$. Now, the vanishing set of the characteristic polynomial will define an
algebraic variety: the cone of characteristic conormals (or characteristic
cone, for brevity). In general, it will consist of different codimension 1
sheets which may be nested, intersect along lines, or even coincide.
Geometrically, hyperbolicity in the direction of $t_{a}$ is the requirement
that every line parallel to $t_{a}$ intersects this algebraic variety at
exactly $``6"$ points (counting multiplicities). Clearly, this condition
severely constrains the topology of the characteristic cones, thus
guaranteeing well-behaved propagation for small wavy excitations. In
particular, it was shown by Garding Gårding (1959) that the closure of the
connected component of $t_{a}$ in the set $\\{n_{a}|p(n_{a})\neq 0\\}$ is
necessarily convex: the hyperbolicity cone of the polynomial. That symmetric
hyperbolicity implies the hyperbolicity of the characteristic polynomial is
direct. To see this one simply observes that the equation

$\mbox{det}(\hat{K}_{\alpha\beta}(a+\lambda t))=0$ (47)

characterizes the eigenvalues of the quadratic form
$\hat{K}_{\alpha\beta}(a_{m})$ relative to the metric
$\hat{K}_{\alpha\beta}(t_{m})$ \- and these eigenvalues have to be real (see
Beig (2006) for further details).

In order to compute the characteristic polynomial explicitly, we substitute
Eq. (40) into Eq. (46) and use Schur’s determinant identity to obtain the
product of determinants

$p(n)=-\mbox{det}(\mathbb{S})q(n),$ (48)

where

$q(n)\equiv\mbox{det}\left(n_{0}^{2}\mathbb{P}-n_{0}n_{k}(\mathbb{Q}\textbf{A}^{k}-\textbf{A}^{k}\mathbb{Q}^{T})-n_{k}n_{l}\textbf{A}^{k}\mathbb{S}\textbf{A}^{l}\right).$
(49)

Notice that $q(t)$ is necessarily positive definite, since $\mathbb{P}>0$ and
$\mathbb{S}<0$ for this particular covector. Interestingly, the determinant in
Eq. (49) has been calculated several times in literature, see e.g. Perlick
(2011),Rivera and Schuller (2011). In particular, it can be shown that it
factorizes as

$q(n)=(n_{m}t^{m})^{2}P(n),$ (50)

where

$P(n)=\frac{1}{24}\varepsilon_{a_{1}a_{2}a_{3}a_{4}}\varepsilon_{b_{1}b_{2}b_{3}b_{4}}X^{a_{1}a_{2}b_{1}c_{1}}X^{c_{2}a_{3}b_{2}c_{3}}X^{c_{4}a_{4}b_{3}b_{4}}n_{c_{1}}n_{c_{2}}n_{c_{3}}n_{c_{4}}.$
(51)

In other words, the sixth order polynomial always reduces to a product of a
quadratic polynomial and a quartic polynomial. Since the vanishing set of the
quadratic polynomial gives a non-compact variety (plane) inconsistent with the
constraint equations, the characteristic cone is given by the covectors
$k_{a}$ that satisfy the fourth-order Fresnel equation

$P(k)\sim\left(\star
X_{pq}^{\phantom{a}\phantom{a}ar}X^{bpcs}X^{dq}_{\phantom{a}\phantom{a}rs}\star\right)k_{a}k_{b}k_{c}k_{d}=0,$
(52)

The fourth rank tensor defined by the terms between parentheses is called the
Kummer tensor, whereas its totally symmetrized version is usually called the
Tamm-Rubilar tensor. If symmetric hyperbolicity holds, the above polynomial is
necessarily hyperbolic in the direction of $t_{a}$, its vanishing set
determining the causal structure of the theory up to a conformal factor. In
specific situations, where $X_{abcd}$ is sufficiently simple, Eq. (52) will
reduce to the more familiar product of quadratic polynomials 666A general
_Ansatz_ leading to bi-metricity in nonlinear electromagnetism was presented
in Visser  _et al._ (2003).

$P(k)\sim\left(g^{ab}_{(1)}k_{a}k_{b}\right)\left(g^{ab}_{(2)}k_{c}k_{d}\right)=0.$
(53)

In these cases, symmetric hyperbolicity guarantees that the rank-2
contravariant tensors $g^{ab}_{(1)}$ and $g^{ab}_{(2)}$ are nondegenerate,
necessarily of Lorentzian type and with the same signature. Furthermore, if
these tensors coincide 777Covariant conditions on the Fresnel surface for
birefringence to be absent were derived in Itin (2005). , then only one of
them must be considered in the dispersion relation (reduced polynomial) ,
_i_.e.,

$P_{red}(k)\sim g^{ab}_{(1)}k_{a}k_{b}=0.$ (54)

It is important to point out that the abovementioned result stating that
symmetric hyperbolicity leads to a cone strongly constraints the possible
shape of the Fresnel surfaces, defined in a convenient 3-space: they must be
topologically equivalent to those obtained from the 4-d hyperbolicity cone.
Hence, open surfaces are not allowed by symmetric-hyperbolic propagation
888For examples of Fresnel surfaces in Maxwell´s theory non-minimally coupled
to gravity, see Balakin and Zayats (2018). .

##  V General result

We have seen that symmetric hyperbolicity imposes restrictions on admissible
physical theories. Importantly, they are always expressed in terms of matrix
inequalities which are much easier to guarantee than to check the
hyperbolicity of the corresponding characteristic polynomial. In order to
obtain explicit expressions for the matrix inequalities, it is convenient to
introduce the projection tensor $h_{ab}=g_{ab}-t_{a}t_{b}$, which projects
arbitrary tensors onto the “rest spaces” orthogonal to $t^{q}$, i.e,

$h_{ab}=h_{ba},\quad\quad\quad
h_{ac}h^{c}_{\phantom{a}b}=h_{ab},\quad\quad\quad h_{ab}t^{b}=0.$ (55)

Similarly, we decompose the electromagnetic 2-form, the Weyl tensor and the
traceless part of the Ricci tensor as

$F_{ab}=E_{[a}t_{b]}+\varepsilon_{abcd}t^{c}B^{d},$ (56)
$W_{abcd}=\left\\{g_{abpq}\left(g_{cdrs}\mathcal{E}^{pr}-\varepsilon_{cdrs}\mathcal{B}^{pr}\right)-\varepsilon_{abpq}\left(g_{cdrs}\mathcal{B}^{pr}+\varepsilon_{cdrs}\mathcal{E}^{pr}\right)\right\\}t^{q}t^{s},$
(57) $S_{ab}=St_{a}t_{b}+Q_{(a}t_{b)}+N_{ab},$ (58)

with the following definitions

$E_{a}\equiv F_{ab}t^{b},\quad\quad\quad B_{a}\equiv\star F_{ab}t^{b},$ (59)
$\mathcal{E}_{ab}\equiv
W_{acbd}t^{c}t^{d},\quad\quad\quad\mathcal{B}_{ab}\equiv\star
W_{acbd}t^{c}t^{d},$ (60) $S\equiv S_{ab}t^{a}t^{b},\quad\quad Q_{a}\equiv
h_{a}^{\phantom{a}b}S_{bc}t^{c},\quad\quad N_{ab}\equiv
h_{a}^{\phantom{a}c}h_{b}^{\phantom{a}d}S_{cd}.$ (61)

According to the latter, the tensor fields
$\\{E_{a},B_{a},Q_{a},\mathcal{E}_{ab},\mathcal{B}_{ab},N_{ab}\\}$ are
automatically orthogonal to the auxiliary vector field.

In order to compute $\mathbb{P}_{ab}$ and $\mathbb{S}_{ab}$ from Eq. (17), the
following relations involving the constitutive tensors are useful

${\phantom{a}}^{(1)}\chi_{abcd}t^{b}t^{d}=h_{ac},\quad\quad{\phantom{a}}^{(2)}\chi_{abcd}t^{b}t^{d}=0,$
(62)
${\phantom{a}}^{(3)}\chi_{abcd}t^{b}t^{d}=Rh_{ac},\quad\quad{\phantom{a}}^{(4)}\chi_{abcd}t^{b}t^{d}=0,$
(63)
${\phantom{a}}^{(5)}\chi_{abcd}t^{b}t^{d}=\mathcal{E}_{ac},\quad\quad{\phantom{a}}^{(6)}\chi_{abcd}t^{b}t^{d}=\mathcal{B}_{ac},$
(64) ${\phantom{a}}^{(7)}\chi_{abcd}t^{b}t^{d}=N_{ac}.$ (65)

Similarly, for the induction tensors, Eq. (11), one obtains

${\phantom{a}}^{(1)}H_{ab}t^{b}=E_{a},\quad\quad{\phantom{a}}^{(2)}H_{ab}t^{b}=B_{a},$
(66)
${\phantom{a}}^{(3)}H_{ab}t^{b}=RE_{a},\quad\quad{\phantom{a}}^{(4)}H_{ab}t^{b}=RB_{a},$
(67)
${\phantom{a}}^{(5)}H_{ab}t^{b}=\mathcal{E}_{ab}E^{b}-\mathcal{B}_{ab}B^{b},\quad\quad{\phantom{a}}^{(6)}H_{ab}t^{b}=\mathcal{E}_{ab}B^{b}+\mathcal{B}_{ab}E^{b},$
(68)
${\phantom{a}}^{(7)}H_{ab}t^{b}=SE_{a}+N_{ab}E^{b}-\varepsilon_{abcd}t^{b}Q^{c}E^{d},$
(69)

Using the above relations and the Ruse-Lanczos identities (see the Appendix
A), it can be checked that the tensors $\mathbb{P}_{ab}$ and $\mathbb{S}_{ab}$
are given by:

$\displaystyle\mathbb{P}_{ab}$ $\displaystyle=$
$\displaystyle(\mathcal{L}_{1}+\mathcal{L}_{3}R)h_{ab}+\mathcal{L}_{5}\mathcal{E}_{ab}+\mathcal{L}_{6}\mathcal{B}_{ab}+\mathcal{L}_{7}N_{ab}$
$\displaystyle+4E_{a}E_{b}(\mathcal{L}_{11}+2R\mathcal{L}_{13})$
$\displaystyle+4E_{(a}B_{b)}\left[\mathcal{L}_{12}+R(\mathcal{L}_{14}+\mathcal{L}_{23})\right]$
$\displaystyle+4B_{a}B_{b}(\mathcal{L}_{22}+2R\mathcal{L}_{24})$
$\displaystyle+4[\mathcal{L}_{15}E_{(a}+\mathcal{L}_{25}B_{(a}][\mathcal{E}_{b)c}E^{c}-\mathcal{B}_{b)c}B^{c}]$
$\displaystyle+4[\mathcal{L}_{16}E_{(a}+\mathcal{L}_{26}B_{(a}][\mathcal{E}_{b)c}B^{c}+\mathcal{B}_{b)c}E^{c}]$
$\displaystyle+4[\mathcal{L}_{17}E_{(a}+\mathcal{L}_{27}B_{(a}][SE_{b)}+N_{b)c}E^{c}-\varepsilon_{b)pqr}Q^{p}t^{q}E^{r}]+...$

$\displaystyle\mathbb{S}_{ab}$ $\displaystyle=$
$\displaystyle-(\mathcal{L}_{1}+\mathcal{L}_{3}R)h_{ab}-\mathcal{L}_{5}\mathcal{E}_{ab}-\mathcal{L}_{6}\mathcal{B}_{ab}+\mathcal{L}_{7}N_{ab}$
$\displaystyle+4B_{a}B_{b}(\mathcal{L}_{11}+2R\mathcal{L}_{13})$
$\displaystyle-4E_{(a}B_{b)}[\mathcal{L}_{12}+R(\mathcal{L}_{14}+\mathcal{L}_{23})]$
$\displaystyle+4E_{a}E_{b}(\mathcal{L}_{22}+2R\mathcal{L}_{24})$
$\displaystyle+4[\mathcal{L}_{15}B_{(a}-\mathcal{L}_{25}E_{(a}][\mathcal{E}_{b)c}B^{c}+\mathcal{B}_{b)c}E^{c}]$
$\displaystyle-4[\mathcal{L}_{16}B_{(a}-\mathcal{L}_{26}E_{(a}][\mathcal{E}_{b)c}E^{c}-\mathcal{B}_{b)c}B^{c}]$
$\displaystyle-4[\mathcal{L}_{17}B_{(a}-\mathcal{L}_{27}E_{(a}][SB_{b)}+N_{b)c}B^{c}-\varepsilon_{b)pqr}Q^{p}t^{q}B^{r}]+...$

where the dots stand for possible nonlinear terms in the irreducible parts of
the curvature tensor, which we shall not discuss in this work. The tensors
$\mathbb{P}_{ab}$ and $\mathbb{S}_{ab}$ given above are the most general ones
when dealing with symmetric hyperbolicity in models of nonlinear
electromagnetism coupled linearly to curvature. The inequalities given in Eqs.
(41), defined in terms of such tensors, lead to drastic restrictions on
admissible theories: there must be a strong compromise between the
electromagnetic quantities $\\{E^{a},B^{a}\\}$, the spacetime curvature
expressed via $\\{\mathcal{E}_{ab},\mathcal{B}_{ab},Q_{a},N_{ab}\\}$ and the
partial derivatives of the Lagrangian density.

##  VI Applications

Let us present next several examples that illustrate the restrictions obtained
so far. For the sake of completeness, we display, in each case, the tensor
$X_{abcd}$, the relevant inequalities and the corresponding characteristic
polynomial. A generic feature of the latter is non-factorization: in general,
the characteristic varieties are described by vanishing sets of fourth-order
polynomials which do not split into products of second-order polynomials.
Needless to say, this fact makes it considerably difficult to guarantee
hyperbolicity straight from the characteristic polynomial. In constrast with
this fact, symmetric-hyperbolicity automatically implies that the
characteristic varieties are necessarily well-behaved.

###  VI.1 Maxwell electrodynamics $\mathcal{L}=-I^{1}$

Our first example is Maxwell’s linear theory, which is the simplest possible
case in our setting. Using Eq.(17), it follows that

$X_{abcd}=-g_{abcd}.$ (72)

Taking the auxiliary vector field $t^{q}$ as timelike, future-directed and
normalized, Eqs.(V) and (V) give the 2-index tensors as

$\mathbb{P}_{ab}=-h_{ab},\quad\quad\quad\mathbb{S}_{ab}=h_{ab}.$ (73)

In a frame such that $g_{ab}=\eta_{ab}$ and $t^{q}=\delta^{q}_{\phantom{a}0}$,
the corresponding $3\times 3$ matrices are

$\mathbb{P}_{ij}=\delta_{ij},\quad\quad\quad\mathbb{S}_{ij}=-\delta_{ij},\quad\quad\quad
i,j=1,2,3.$ (74)

Since they satisfy Eqs.(41) identically, symmetric-hyperbolicity is guaranteed
for all timelike covector fields $n_{m}=\eta_{mn}t^{n}$. A direct calculation
using Eq. (52) leads to

$\displaystyle P(k)$ $\displaystyle\sim$
$\displaystyle(g^{ab}k_{a}k_{b})^{2}.$

Only one of the quadratic polinomials is needed to obtain the usual dispersion
relation for linear electromagnetic waves in vacuum:

$P_{red}(k)\sim g^{ab}k_{a}k_{b}=0.$ (75)

Clearly, the cone of hyperbolicity at $p$, $S_{p}$, is the connected component
of $t_{a}$ i.e., it coincides with the set of future-directed timelike
covectors, as expected.

###  VI.2 Minimally coupled nonlinear electrodynamics

####  VI.2.1 Single invariant: $\mathcal{L}=\mathcal{L}(I^{1})$

If the Lagrangian density is an arbitrary function of the invariant $I^{1}$, a
direct inspection of Eq. (17) gives

$X_{abcd}=\mathcal{L}_{1}g_{abcd}+4\mathcal{L}_{11}F_{ab}F_{cd}.$ (76)

From Eqs.(V) and (V), we obtain the projections

$\displaystyle\mathbb{P}_{ab}=\mathcal{L}_{1}h_{ab}+4\mathcal{L}_{11}E_{a}E_{b},\quad\quad\quad\mathbb{S}_{ab}=-\mathcal{L}_{1}h_{ab}+4\mathcal{L}_{11}B_{a}B_{b},$
(77)

Simple manipulations using Eq.(41) then yield the matrix inequalities

$\mathcal{L}_{1}\delta_{ij}-4\mathcal{L}_{11}E_{i}E_{j}\prec 0,$ (78)
$\mathcal{L}_{1}\delta_{ij}+4\mathcal{L}_{11}B_{i}B_{j}\prec 0.$ (79)

In particular, they imply the condition $\mathcal{L}_{1}<0$, since the theory
should be well-behaved when either the electric or magnetic field vanish.
However, notice that for field intensities violating the inequalities, the
propagation of small wavy excitations about the background field will be, in
general, ill-posed.

A straightforward calculation using Eq. (52) gives the well-known bi-metric
dispersion relation

$P(k)\sim\left(g^{ab}_{(1)}k_{a}k_{b}\right)\left(g^{ab}_{(2)}k_{c}k_{d}\right),$
(80)

with the effective metrics given by

$g^{ab}_{(1)}\equiv\mathcal{L}_{1}^{2}g^{ab},\quad\quad\quad
g^{ab}_{(2)}\equiv-(\mathcal{L}_{1}g^{ab}+4\mathcal{L}_{11}\tau^{ab}),\quad\quad\quad\tau^{ab}\equiv
F^{ac}F^{b}_{\phantom{a}c},$ (81)

as was shown for instance in Novello  _et al._ (2000). Needless to say, if
Eqs. (78) and (79) hold, the effective metrics are well-behaved, and both
possess the same Lorentzian signature. However, it should be clear from our
discussion that symmetric-hyperbolicity requires more than the simple
Lorentzian nature of the effective metrics.

####  VI.2.2 Two invariants: $\mathcal{L}(I^{1},I^{2})$

The case of two invariants is naturally more involved. The symmetric double
$(2,2)$ form becomes

$X_{abcd}=\mathcal{L}_{1}g_{abcd}+4\\{\mathcal{L}_{11}F_{ab}F_{cd}+\mathcal{L}_{12}(F_{ab}\star
F_{cd}+\star F_{ab}F_{cd})+\mathcal{L}_{22}\star F_{ab}\star F_{cd}\\}.$ (82)

Apart from minor details on definitions, this tensor coincides with the so-
called jump tensor obtained by Obukhov and Rubilar in Obukhov and Rubilar
(2002). Notice also that the term involving $\mathcal{L}_{2}$ in Eq. (17) was
discarded, since it is proportional to the Bianchi identity. We then obtain
the tensors (see Eqs.(V) and (V))

$\mathbb{P}_{ab}=\mathcal{L}_{1}h_{ab}+4\\{\mathcal{L}_{11}E_{a}E_{b}+\mathcal{L}_{12}(E_{a}B_{b}+B_{a}E_{b})+\mathcal{L}_{22}B_{a}B_{b}\\},$
(83)
$\mathbb{S}_{ab}=-\mathcal{L}_{1}h_{ab}+4\\{\mathcal{L}_{11}B_{a}B_{b}-\mathcal{L}_{12}(E_{a}B_{b}+B_{a}E_{b})+\mathcal{L}_{22}E_{a}E_{b}\\},$
(84)

which, in the frame described above, lead to the inequalities

$\mathcal{L}_{1}\delta_{ij}-4\\{\mathcal{L}_{11}E_{i}E_{j}+\mathcal{L}_{12}(E_{i}B_{j}+B_{i}E_{j})+\mathcal{L}_{22}B_{i}B_{j}\\}\prec
0,$ (85)
$\mathcal{L}_{1}\delta_{ij}+4\\{\mathcal{L}_{11}B_{i}B_{j}-\mathcal{L}_{12}(E_{i}B_{j}+B_{i}E_{j})+\mathcal{L}_{22}E_{i}E_{j}\\}\prec
0,$ (86)

A tedious calculation using Eq. (52) yields the effective metrics

$\displaystyle
g^{ab}_{(1)}=\mathcal{X}g^{ab}+(\mathcal{Y}+\sqrt{\mathcal{Y}^{2}-\mathcal{X}\mathcal{Z}})t^{ab},$
(87) $\displaystyle
g^{ab}_{(2)}=\mathcal{X}g^{ab}+(\mathcal{Y}-\sqrt{\mathcal{Y}^{2}-\mathcal{X}\mathcal{Z}})t^{ab},$
(88)

with

$\displaystyle\mathcal{X}=\mathcal{L}_{1}^{2}+2\mathcal{L}_{1}(G-\mathcal{L}_{22}F)+(\mathcal{L}_{12}\mathcal{L}_{12}-\mathcal{L}_{11}\mathcal{L}_{22})G^{2},$
$\displaystyle\mathcal{Y}=2\mathcal{L}_{1}(\mathcal{L}_{11}+\mathcal{L}_{22})+4(\mathcal{L}_{12}\mathcal{L}_{12}-\mathcal{L}_{11}\mathcal{L}_{22})F,$
$\displaystyle\mathcal{Z}=(\mathcal{L}_{11}\mathcal{L}_{22}-\mathcal{L}_{12}\mathcal{L}_{12}),$

and we have used the notation $I^{1}\equiv F$, $I^{2}\equiv G$, $t^{ab}\equiv
F^{ac}F^{b}_{\phantom{a}c}$ for conciseness. This result coincides with
Obukhov and Rubilar (2002); De Lorenci  _et al._ (2000), with minor
modifications of notation.

###  VI.3 Nonminimally coupled nonlinear electrodynamics.

Since the equations of motion that follow from Eq. (10) are very general, it
is convenient to consider some particular cases. Hence, in what follows we
assume that the Lagrangian density takes the form

$\mathcal{L}(I^{1},I^{2})+\alpha I^{3}+\beta I^{5}+\gamma I^{7},$ (89)

where $\mathcal{L}(I^{1},I^{2})$ is an arbitrary function of the usual
electromagnetic invariants, and $\alpha,\beta$ and $\gamma$ are
phenomenological parameters to be determined in principle from experiments.
Such a density is still sufficiently general to include most of the relevant
important models present in literature.

Let us examine next some particular cases that follow from Eq.(89).

#### Linear non-minimal model

The corresponding Lagrangian density is given by

$\mathcal{L}=-I^{1}+\alpha I^{3}+\beta I^{5}+\gamma I^{7},$ (90)

since we expect to recover Maxwell’s theory in the flat spacetime regime. This
model describes several possible non-minimal couplings of linear
electrodynamics with gravity, and encompasses for instance the modifications
of Maxwell´s theory due to one-loop vacuum polarization contributions Drummond
and Hathrell (1980). Clearly, the corresponding equations of motion are linear
and, using Eq. (17), one obtains the double symmetric $(2,2)$-form as

$X_{abcd}=(\alpha R-1)g_{abcd}+\beta
W_{abcd}+\gamma(g_{ac}S_{bd}-g_{ad}S_{bc}+g_{bd}S_{ac}-g_{bc}S_{ad}).$ (91)

Notice that the irreducible parts of the Riemann tensor contribute in
different ways to the equations of motion. Using the above prescription in
Eqs.(V) and (V), we obtain the following $2$-index tensors

$\mathbb{P}_{ab}=(\alpha R-1)h_{ab}+\beta\mathcal{E}_{ab}+\gamma N_{ab}$ (92)
$\mathbb{S}_{ab}=(1-\alpha R)h_{ab}-\beta\mathcal{E}_{ab}+\gamma N_{ab}$ (93)

Let us consider next several sub-cases:

##### Scalar curvature coupling $(\alpha\neq 0,\beta=0,\gamma=0)$.

This is by far the simplest type of non-minimal coupling. Indeed, since the
double symmetric $(2,2)$-form in Eq.(91) reduces to

$X_{abcd}=(\alpha R-1)g_{abcd},$ (94)

the relevant $2$-index tensors read

$\mathbb{P}_{ab}=(\alpha R-1)h_{ab},\quad\quad\quad\mathbb{S}_{ab}=(1-\alpha
R)h_{ab}.$ (95)

Hence, the simple matrix inequality

$(1-\alpha R)\delta_{ij}\succ 0$ (96)

follows. Clearly, symmetric-hyperbolicity requires that $\alpha<1/R$, which
may forbid good propagation for sufficiently high curvature, for a given
$\alpha$. Regarding the characteristic cone, Eq. (52) gives

$P(k)\sim(g^{ab}k_{a}k_{b})^{2},$ (97)

which shows that in this case the dispersion relation is governed by the
background metric i.e. the causal structure is not changed by the coupling.

##### Weyl coupling $(\alpha=0,\beta\neq 0,\gamma=0)$.

It follows from Eqs. (92) and (93) that

$\mathbb{P}_{ab}=-\mathbb{S}_{ab}=-h_{ab}+\beta\mathcal{E}_{ab}.$

Although the ensuing inequalities need to be examined on a case by case basis,
they will lead to limits on the components of the electric part of the Weyl
tensor 999Conversely, for a given geometry, the inequalities may furnish
$\beta$-dependent limits on the region of space-time where the propagation is
symmetric-hyperbolic..

In order to compute the dispersion relation explicitly, we first recall that
the Weyl conformal tensor has only one independent Hodge dual i.e. $\star
W_{abcd}=W_{abcd}\star$. A direct calculation using Eq. (52) gives a quartic
equation of the type

$P(k)\sim G^{abcd}k_{a}k_{b}k_{c}k_{d},$ (98)

with the Kummer tensor given by

$G^{abcd}\equiv
g^{ab}g^{cd}-\frac{\beta^{2}}{3}\left(W^{apbq}W^{c\phantom{a}d}_{\phantom{a}p\phantom{a}q}+\frac{1}{4}g^{ab}W^{cpqr}W^{d}_{\phantom{a}pqr}\right)+\frac{\beta^{3}}{6}\star
W_{pq}^{\phantom{a}\phantom{a}ar}W^{bpcs}W^{dq}_{\phantom{a}\phantom{a}rs}\star.$
(99)

Using the following identites due to Debever and Lanczos

$W_{apcq}W_{b\phantom{a}d}^{\phantom{a}p\phantom{a}q}-\star W_{apcq}\star
W_{b\phantom{a}d}^{\phantom{a}p\phantom{a}q}=Ag_{ab}g_{cd}$ (100)
$W_{apqr}W_{b}^{\phantom{a}pqr}=2Ag_{ab}$ (101)

with $A\equiv\frac{1}{8}W_{pqrs}W^{pqrs}$, we obtain the simplified expression

$G^{abcd}=\left(1-\frac{\beta^{2}}{6}A\right)g^{ab}g^{cd}-\frac{\beta^{2}}{3}W^{apbq}W^{c\phantom{a}d}_{\phantom{a}p\phantom{a}q}+\frac{\beta^{3}}{3}W_{pq}^{\phantom{a}\phantom{a}ar}W^{bpcs}W^{dq}_{\phantom{a}\phantom{a}rs}.$
(102)

##### Traceless Ricci couplings $(\alpha=0,\beta=0,\gamma\neq 0)$.

The relevant matrices take the form

$\mathbb{P}_{ab}=-h_{ab}+\gamma N_{ab},$ (103) $\mathbb{S}_{ab}=h_{ab}+\gamma
N_{ab}.$ (104)

As in the previous case, the corresponding inequalities will relate the
coupling constant $\gamma$ to the curvature quantities described by $N_{ab}$.
The Kummer tensor is given by

$\displaystyle G^{abcd}$ $\displaystyle=$ $\displaystyle g^{ab}g^{cd}-\gamma
g^{ab}\
{}^{(7)}\chi^{cpd}_{\phantom{a}\phantom{a}\phantom{a}p}+\frac{\gamma^{2}}{2}\left({}^{(7)}\chi^{apb}_{\phantom{a}\phantom{a}\phantom{a}p}\
{}^{(7)}\chi^{cqd}_{\phantom{a}\phantom{a}\phantom{a}q}-^{(7)}\chi^{apbq}\
{}^{(7)}\chi^{c\phantom{a}d}_{\phantom{a}p\phantom{a}q}\right)$
$\displaystyle\quad\quad\quad\quad\quad+\frac{\gamma^{3}}{6}\
^{(7)}\chi^{apbq}\ {}^{(7)}\chi^{crds}\ {}^{(7)}\chi_{prqs}.$

##### Non-minimally coupled EM in a cosmological background.

Another relevant example is that of the propagation in a cosmological
background described by the flat Friedman-Lemâitre-Robertson-Walker metric,
for which the Weyl tensor is null. In such a case, only the terms associated
to ${\cal L}_{1}$, ${\cal L}_{3}$, and ${\cal L}_{7}$ survive in Eqs. (V) and
(V). Using Einstein´s equations, the traceless part of the Ricci tensor is
given in terms of the matter by

$S_{ab}=T_{ab}-\frac{T}{4}g_{ab},$

where $T_{ab}$ is the energy-momentum tensor and $T$, its trace. In a
convenient tetrad basis, in which $T_{ab}={\rm diag}\;(\rho,p,p,p)$, it
follows that

$\mathbb{P}_{ij}=-\delta_{ij}\left[1+\beta(\rho-3p)-\frac{1}{2}\gamma(\rho+p)\right],$
(105)
$\mathbb{S}_{ij}=\delta_{ij}\left[1+\beta(\rho-3p)+\frac{1}{2}\gamma(\rho+p)\right],$
(106)

Hence, the inequalities (41) lead to

$1+\beta(\rho-3p)>\frac{1}{2}|\gamma|(\rho+p),$

which is trivially satisfied in the case of linear EM. The propagation will be
symmetric hyperbolic if this inequaly is satisfied at all times for which the
model is valid.

##  VII Conclusions

Well-posedness is a basic requirement for any field theory, and is guaranteed
by symmetric hyperbolicity. We have obtained a general form for the
symmetrizer, given in Eq. (26), valid for a general Lagrangian theory. We have
shown that symmetric hyperbolicity leads to a set of two inequalities for the
matrices $\mathbb{P}$ and $\mathbb{S}$, whose elements are determined by a
given theory. Regarding the constraints, we have verified that they are
integrable.

When applied to nonlinear electromagnetism linearly coupled to curvature, the
matrices $\mathbb{P}$ and $\mathbb{S}$ are expressed in terms the fields, the
Lagrangian, and its derivatives, and also of the different quantities
associated to curvature. They lead to strong constraints on the relevant
quantities, which were illustrated with applications to several particular
cases. The examples show that while in the linear theory, no constraint arises
from symmetric hyperbolic propagation, non-linearity leads to constraints on
the field intensities, and non-minimal coupling imposes restrictions on
quantities associated to curvature. In the general case, symmetric
hyperbolicity relates the electromagnetic quantities $\\{E^{a},B^{a}\\}$, the
spacetime curvature expressed via
$\\{\mathcal{E}_{ab},\mathcal{B}_{ab},Q_{a},N_{ab}\\}$ and the partial
derivatives of the Lagrangian density.

The ideas presented here can be applied in other settings, such as
electromagnetism in material media. We plan to return to this problem in a
future publication.

##  Appendix A Ruse-Lanczos identity

Let $\chi_{abcd}$ denote an arbitrary double symmetric (2,2) form at a point
$p\in\mathcal{M}$. Its double Hodge dual is defined as

$\star\chi_{abcd}\star=\frac{1}{4}\varepsilon_{ab}^{\phantom{a}\phantom{a}pq}\varepsilon_{cd}^{\phantom{a}\phantom{a}rs}\chi_{pqrs}.$
(107)

Writing the traces of $\chi_{abcd}$ as

$\chi_{ab}\equiv\chi^{c}_{\phantom{a}acb},\quad\quad\quad\chi\equiv\chi^{a}_{\phantom{a}a},$
(108)

we may construct the symmetric trace-free tensor

$\psi_{ab}\equiv\chi_{ab}-\frac{1}{4}g_{ab}\chi.$ (109)

Using elementary algebraic manipulations, it can be shown that

$\star\chi_{abcd}\star+\chi_{abcd}=g_{a[c}\psi_{d]b}+g_{b[d}\psi_{c]a},$ (110)

which is called the Ruse-Lanczos identity Hall (2004); de Felice and Clarke
(1992). Aplying this identity to the constitutive tensors results in the
following useful relations

$\displaystyle\star^{(1)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star+^{(1)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0,$
$\displaystyle\star^{(2)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star+^{(2)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0,$
$\displaystyle\star^{(3)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star+^{(3)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0,$
$\displaystyle\star^{(4)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star+^{(4)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0,$
$\displaystyle\star^{(5)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star+^{(5)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0,$
$\displaystyle\star^{(6)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star+^{(6)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0,$
$\displaystyle\star^{(7)}\chi^{ab}_{\phantom{a}\phantom{a}cd}\star-^{(7)}\chi^{ab}_{\phantom{a}\phantom{a}cd}=0.$

## References

  * Kreiss and Lorenz (1989) H. Kreiss and J. Lorenz,  _Initial-boundary Value Problems and the Navier-Stokes Equations_, Initial-boundary value problems and the Navier-Stokes equations No. v. 136 (Academic Press, 1989).
  * Papallo and Reall (2017) G. Papallo and H. S. Reall, On the local well-posedness of Lovelock and Horndeski theories, Phys. Rev. D 96, 044019 (2017), arXiv:1705.04370 [gr-qc] .
  * Beig (2006) R. Beig, Concepts of hyperbolicity and relativistic continuum mechanics, Lect. Notes Phys. 692, 101 (2006), arXiv:gr-qc/0411092 .
  * Sarbach and Tiglio (2012) O. Sarbach and M. Tiglio, Continuum and Discrete Initial-Boundary-Value Problems and Einstein’s Field Equations, Living Rev. Rel. 15, 9 (2012), arXiv:1203.6443 [gr-qc] .
  * Papallo (2017) G. Papallo, On the hyperbolicity of the most general Horndeski theory, Phys. Rev. D 96, 124036 (2017), arXiv:1710.10155 [gr-qc] .
  * Bernard  _et al._ (2019) L. Bernard, L. Lehner, and R. Luna, Challenges to global solutions in Horndeski’s theory, Phys. Rev. D 100, 024011 (2019), arXiv:1904.12866 [gr-qc] .
  * Kovács and Reall (2020) A. D. Kovács and H. S. Reall, Well-Posed Formulation of Scalar-Tensor Effective Field Theory, Phys. Rev. Lett. 124, 221101 (2020), arXiv:2003.04327 [gr-qc] .
  * Moreno and Sarbach (2003) C. Moreno and O. Sarbach, Stability properties of black holes in selfgravitating nonlinear electrodynamics, Phys. Rev. D 67, 024028 (2003), arXiv:gr-qc/0208090 .
  * Bretón and Perez Bergliaffa (2015) N. Bretón and S. E. Perez Bergliaffa, On the thermodynamical stability of black holes in nonlinear electrodynamics, Annals Phys. 354, 440 (2015).
  * Toshmatov  _et al._ (2018) B. Toshmatov, Z. Stuchlík, J. Schee, and B. Ahmedov, Electromagnetic perturbations of black holes in general relativity coupled to nonlinear electrodynamics, Phys. Rev. D 97, 084058 (2018), arXiv:1805.00240 [gr-qc] .
  * Magos and Bretón (2020) D. Magos and N. Bretón, Thermodynamics of the Euler-Heisenberg-AdS black hole, Phys. Rev. D 102, 084011 (2020), arXiv:2009.05904 [gr-qc] .
  * Falciano  _et al._ (2021) F. T. Falciano, M. L. Peñafiel, and J. C. Fabris, Entropy bound in Einstein-Born-Infeld black holes, Phys. Rev. D 103, 084046 (2021), arXiv:2103.14109 [gr-qc] .
  * Heyl and Hernquist (2004) J. S. Heyl and L. Hernquist, A QED model for the origin of bursts from SGRs and AXPs, Astrophys. J. 618, 463 (2004), arXiv:astro-ph/0312608 .
  * Harding and Lai (2006) A. K. Harding and D. Lai, Physics of Strongly Magnetized Neutron Stars, Rept. Prog. Phys. 69, 2631 (2006), arXiv:astro-ph/0606674 .
  * Turolla  _et al._ (2015) R. Turolla, S. Zane, and A. Watts, Magnetars: the physics behind observations. A review, Rept. Prog. Phys. 78, 116901 (2015), arXiv:1507.02924 [astro-ph.HE] .
  * Novello  _et al._ (2004) M. Novello, S. E. Perez Bergliaffa, and J. Salim, Non-linear electrodynamics and the acceleration of the universe, Phys. Rev. D 69, 127301 (2004), arXiv:astro-ph/0312093 .
  * Novello  _et al._ (2007) M. Novello, E. Goulart, J. M. Salim, and S. . Perez Bergliaffa, Cosmological Effects of Nonlinear Electrodynamics, Class. Quant. Grav. 24, 3021 (2007), arXiv:gr-qc/0610043 .
  * Campanelli  _et al._ (2008) L. Campanelli, P. Cea, G. L. Fogli, and L. Tedesco, Inflation-Produced Magnetic Fields in Nonlinear Electrodynamics, Phys. Rev. D 77, 043001 (2008), arXiv:0710.2993 [astro-ph] .
  * Vollick (2008) D. N. Vollick, Homogeneous and isotropic cosmologies with nonlinear electromagnetic radiation, Phys. Rev. D 78, 063524 (2008), arXiv:0807.0448 [gr-qc] .
  * Kruglov (2015) S. I. Kruglov, Universe acceleration and nonlinear electrodynamics, Phys. Rev. D 92, 123523 (2015), arXiv:1601.06309 [gr-qc] .
  * Gutierrez  _et al._ (1981) S. A. Gutierrez, A. L. Dudley, and J. F. Plebanski, Signals and Discontinuities in General Relativistic Nonlinear Electrodynamics, J. Math. Phys. 22, 2835 (1981).
  * Deser  _et al._ (1999) S. Deser, J. G. McCarthy, and O. Sarioglu, ’Good propagation’ constraints on dual invariant actions in electrodynamics and on massless fields, Class. Quant. Grav. 16, 841 (1999), arXiv:hep-th/9809153 .
  * Goulart de Oliveira Costa and Esteban Perez Bergliaffa (2009) E. Goulart de Oliveira Costa and S. Esteban Perez Bergliaffa, A Classification of the effective metric in nonlinear electrodynamics, Class. Quant. Grav. 26, 135015 (2009), arXiv:0905.3673 [gr-qc] .
  * Perlick  _et al._ (2018) V. Perlick, C. Lämmerzahl, and A. Macias, Effects of nonlinear vacuum electrodynamics on the polarization plane of light, Phys. Rev. D 98, 105014 (2018), arXiv:1806.00891 [gr-qc] .
  * Abalos  _et al._ (2015) F. Abalos, F. Carrasco, E. Goulart, and O. Reula, Nonlinear electrodynamics as a symmetric hyperbolic system, Phys. Rev. D 92, 084024 (2015), arXiv:1507.02262 [gr-qc] .
  * Perlick (2011) V. Perlick, On the hyperbolicity of Maxwell’s equations with a local constitutive law, J. Math. Phys. 52, 042903 (2011), arXiv:1011.2536 [gr-qc] .
  * Balakin and Lemos (2005) A. B. Balakin and J. P. S. Lemos, Non-minimal coupling for the gravitational and electromagnetic fields: A General system of equations, Class. Quant. Grav. 22, 1867 (2005), arXiv:gr-qc/0503076 .
  * Balakin  _et al._ (2008) A. B. Balakin, V. V. Bochkarev, and J. P. S. Lemos, Non-minimal coupling for the gravitational and electromagnetic fields: Black hole solutions and solitons, Phys. Rev. D 77, 084013 (2008), arXiv:0712.4066 [gr-qc] .
  * Balakin  _et al._ (2010) A. B. Balakin, J. P. S. Lemos, and A. E. Zayats, Nonminimal coupling for the gravitational and electromagnetic fields: Traversable electric wormholes, Phys. Rev. D 81, 084015 (2010), arXiv:1003.4584 [gr-qc] .
  * Dereli and Sert (2011) T. Dereli and O. Sert, Non-minimally Coupled Gravitational and Electromagnetic Fields: pp-Wave Solutions, Phys. Rev. D 83, 065005 (2011), arXiv:1101.1177 [gr-qc] .
  * Friedrichs (1954) K. O. Friedrichs, Symmetric hyperbolic linear differential equations, Communications on Pure and Applied Mathematics 7, 345 (1954), https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160070206 .
  * Hall  _et al._ (1996) G. Hall, J. Pulham, and A. Scottish Universities Summer School in Physics (46, 1995,  _General Relativity: Proceedings of the Forty Sixth Scottish Universities Summer School in Physics, Aberdeen, July 1995_, A NATO advanced study institute (Scottish Universities Summer School in Physics, 1996).
  * Besse (1987) A. Besse,  _Einstein Manifolds_, Classics in mathematics (Springer, 1987).
  * Drummond and Hathrell (1980) I. T. Drummond and S. J. Hathrell, QED Vacuum Polarization in a Background Gravitational Field and Its Effect on the Velocity of Photons, Phys. Rev. D 22, 343 (1980).
  * Kato (1975) T. Kato, The Cauchy problem for quasi-linear symmetric hyperbolic systems, Archive for Rational Mechanics and Analysis 58, 181 (1975).
  * Senovilla (2000) J. M. M. Senovilla, Superenergy tensors, Class. Quant. Grav. 17, 2799 (2000), arXiv:gr-qc/9906087 .
  * Balakin and Zayats (2018) A. B. Balakin and A. E. Zayats, Non-minimal Einstein–Maxwell theory: the Fresnel equation and the Petrov classification of a trace-free susceptibility tensor, Class. Quant. Grav. 35, 065006 (2018), arXiv:1710.08013 [gr-qc] .
  * Gårding (1959) L. Gårding, An inequality for hyperbolic polynomials, Journal of Mathematics and Mechanics 8, 957 (1959).
  * Rivera and Schuller (2011) S. Rivera and F. P. Schuller, Quantization of general linear electrodynamics, Phys. Rev. D 83, 064036 (2011), arXiv:1101.0491 [hep-th] .
  * Visser  _et al._ (2003) M. Visser, C. Barcelo, and S. Liberati, Birefringence versus bimetricity (2003) arXiv:gr-qc/0204017 .
  * Itin (2005) Y. Itin, No-birefringence conditions for spacetime, Phys. Rev. D 72, 087502 (2005), arXiv:hep-th/0508144 .
  * Novello  _et al._ (2000) M. Novello, V. A. De Lorenci, J. M. Salim, and R. Klippert, Geometrical aspects of light propagation in nonlinear electrodynamics, Phys. Rev. D 61, 045001 (2000), arXiv:gr-qc/9911085 .
  * Obukhov and Rubilar (2002) Y. N. Obukhov and G. F. Rubilar, Fresnel analysis of the wave propagation in nonlinear electrodynamics, Phys. Rev. D 66, 024042 (2002), arXiv:gr-qc/0204028 .
  * De Lorenci  _et al._ (2000) V. A. De Lorenci, R. Klippert, M. Novello, and J. M. Salim, Light propagation in nonlinear electrodynamics, Phys. Lett. B 482, 134 (2000), arXiv:gr-qc/0005049 .
  * Hall (2004) G. Hall,  _Symmetries And Curvature Structure In General Relativity_, World Scientific Lecture Notes In Physics (World Scientific Publishing Company, 2004).
  * de Felice and Clarke (1992) F. de Felice and C. J. S. Clarke,  _Relativity on curved manifolds_ (Cambridge University Press, 1992).