format: text

text:
# Stability properties of black holes in self-gravitating nonlinear electrodynamics

## I Introduction

In the last few years, interesting exact solutions ABG-E1 ; ABG-E2 ; ABG-E3 ; ABG-M1 ; Br ; BH ; EH  have been obtained in the context of nonlinear electrodynamics (NED) coupled to gravity. These solutions are interesting because for a suitable choice of the electromagnetic Lagrangian they describe static black hole spacetimes where the black hole does not exhibit the usual singularity at its center, but instead all fields are everywhere regular and fall off in the asymptotic region. The existence of these solutions does not contradict the singularity theorems HE  since the corresponding spacetimes are not globally hyperbolic and their stress-energy tensor does not satisfy the strong energy condition. Examples of regular black hole spacetimes have been given by Bardeen Bardeen . While the corresponding Einstein tensor satisfies the weak energy condition these spacetimes do not satisfy the vacuum Einstein equations. As shown recently by Ayón-Beato and García ABG-M1  the non-vanishing Einstein tensor in the Bardeen model can be associated with the stress-energy tensor of a nonlinear electromagnetic Lagrangian and the Bardeen black holes can be obtained as exact solutions of an appropriate NED coupled to gravity.

An important issue in field theories is the study of stability of a given static or stationary field configuration. Unstable solutions have usually less physical significance than stable ones since they are likely to decay to a stable configuration. For example, it is known that static spherically symmetric regular particle-like solutions in Einstein-Yang-Mills theories are linearly unstable and numerical studies suggest that (depending on the details of the initial perturbations) a perturbation either causes the system to explode or to collapse to a Schwarzschild black hole (which is linearly stable) Zhou-Diss .

In the context of regular black hole solutions in NED, we study the behavior of linear (classical) fluctuations of the gravitational and electromagnetic fields in the domain of outer communication. Basically, there are two possible outcomes: 1) The fluctuations grow exponentially in time which means that the background configuration is (at least linearly) unstable. In this case a natural question to ask is to what configuration (if any) the solution will decay to. 2) The fluctuations are bounded in time which means that the background configuration is linearly stable. This suggests that the perturbed configuration eventually settles down to the background configuration after the possible emission of gravitational and electromagnetic radiation.

In this article, we derive simple conditions on the electromagnetic Lagrangian that guarantee the corresponding black hole solutions to be linearly stable. In order to do so we cast the perturbation equations into a system of coupled wave equations with symmetric potential. Linear stability can then be established if this potential is positive definite. We show that the latter condition holds if the electromagnetic Lagrangian and its first and second order derivatives satisfy appropriate inequalities. We check these inequalities for several of the regular black holes that have been found in self-gravitating NED and so prove that these solutions are linearly stable. However, we also argue that this situation is not generic. That is, one can easily obtain an unstable regular black hole solution by slightly changing the electromagnetic Lagrangian belonging to a regular black hole solution.

A brief review of the theory of NED in general relativity and of the regular black hole solutions is presented in section II. In section III we derive the perturbation equations and cast them into the form of a wave equation with potential. We derive sufficient conditions for this potential to be positive and show in section IV that the latter condition implies the linear stability of the black holes. Conclusions are drawn in section V. A demonstration of the fact that the Bardeen black holes are linearly stable is given in an appendix at the end of this article.

## II Regular black holes in self-gravitating nonlinear electrodynamics

In this section, we first review a few facts from NED and then focus on spherically symmetric solutions. In particular, we give a brief description of the regular black hole solutions obtained in ABG-E1 ; ABG-E2 ; ABG-E3 .

### II.1 Equations of motion

The action for NED in general relativity is

$$S[𝒈,𝑨]=\frac{1}{4\pi }\sqrt{g}d^4x\left(\frac{1}{4}R(F)\right),$$

where $`R`$ denotes the Ricci scalar with respect to the spacetime metric $`g_{\mu \nu }`$, $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$ is the electromagnetic field, and $``$ is an arbitrary function in the invariant $`F\frac{1}{4}F^{\mu \nu }F_{\mu \nu }`$. In Einstein-Maxwell theory, $`(F)=F`$, but here we consider more general choices of electromagnetic Lagrangians.111However we do not consider the case where $``$ depends on the invariant $`F_{\mu \nu }F^{\mu \nu }`$ as well. The field equations are

$`G_\nu ^\mu =2\left[_FF^{\mu \alpha }F_{\nu \alpha }\delta _\nu ^\mu \right],`$ (1)
$`_\mu \left(_FF^{\alpha \mu }\right)=0,`$ (2)
$`_\mu F^{\alpha \mu }=0,`$ (3)
where $`_F/F`$ and $`G_\nu ^\mu `$, $``$ and $``$ denote the components of the Einstein tensor, the covariant derivative and the Hodge dual with respect to $`g_{\mu \nu }`$, respectively.222We use the metric signature $`(,+,+,+)`$ throughout, and Greek letters denote spacetime indices taking values $`(0,1,2,3)`$. The Hodge dual of $`F_{\mu \nu }`$ is given by $`F^{\mu \nu }=\frac{1}{2}\eta ^{\mu \nu \alpha \beta }F_{\alpha \beta }`$, where $`\eta _{\mu \nu \alpha \beta }`$ is completely skew-symmetric with $`\eta _{0123}=\sqrt{g}`$. We have included the Bianchi identities (3) for completeness and later use.

It is also convenient to express the field equations in terms of the field $`P_{\mu \nu }_FF_{\mu \nu }`$ and the Hamiltonian $`2_FF`$ which we consider to be a function of $`P\frac{1}{4}P^{\mu \nu }P_{\mu \nu }`$ SGP :

$`G_\nu ^\mu =2\left[_PP^{\mu \alpha }P_{\nu \alpha }\delta _\nu ^\mu \left(2_PP\right)\right],`$ (4)
$`_\mu P^{\alpha \mu }=0,`$ (5)
$`_\mu (_PP^{\alpha \mu })=0.`$ (6)
However, one should notice that the set of field equations (1,2,3) (called the $`F`$ framework) are equivalent to the set of equations (4,5,6) (called the $`P`$ framework) only if the map $`FP`$ is invertible. Since $`P=_F^2F`$, this requires that

$`0{\displaystyle \frac{P}{F}}`$ $`=`$ $`_F\left(_F+2_{FF}F\right),\text{and}`$
$`0{\displaystyle \frac{F}{P}}`$ $`=`$ $`_P\left(_P+2_{PP}P\right).`$ (7)
Finally, we mention that there is a duality between the $`F`$ framework and the $`P`$ framework: Using the fact that

$$_F_P=1,=2_PP,P^{\mu \alpha }P_{\nu \alpha }=P^{\mu \alpha }P_{\nu \alpha }2\delta _\nu ^\mu P,$$

it is not difficult to see that the equations (4,5,6) can be obtained from the equations (1,2,3) by performing the following transformation Br

$`g_{\mu \nu }`$ $``$ $`\overline{g}_{\mu \nu }g_{\mu \nu }`$
$`F_{\mu \nu }`$ $``$ $`\overline{F}_{\mu \nu }P_{\mu \nu },`$ (8)
$`(F)`$ $``$ $`\overline{}(\overline{F})(P)=(\overline{F}).`$
Notice however that the $`FP`$ duality transformation (8) is only a symmetry of the system if $`\overline{}(F)=(F)`$, i.e. if $`(F)+2F_F(F)=(F)`$. The only Lagrangian that satisfies this identity is $`(F)=\alpha F+\mathrm{\Lambda }/2`$, with $`\alpha `$ and $`\mathrm{\Lambda }`$ two constants. If $`\alpha =1`$ this is the Maxwell Lagrangian with a cosmological term. In general, the $`FP`$ duality connects theories belonging to different Lagrangians. (For a more detailed analysis of duality invariant nonlinear electrodynamics, see Refs. SGP ; GR .)

### II.2 Spherically symmetric solutions

We now focus on spherically symmetric and purely electric solutions in the $`P`$ framework, equations (4,5,6). From this, one can also obtain purely magnetic solution in the $`F`$ framework (but with a different Lagrangian) by applying the inverse of the $`FP`$ duality transformation (8). In the remainder of this section, we assume that the Hamiltonian $`(P)`$ is negative and satisfies the weak field limit

$$(P)P,\text{as}P0,$$

although this condition will not be a crucial assumption in the stability analysis below.

A spherically symmetric spacetime can locally be written as $`\stackrel{~}{M}\times S^2`$ where $`\stackrel{~}{M}`$ is a two-dimensional pseudo-Riemannian manifold. The metric takes the form

$$𝒈=\stackrel{~}{g}_{ab}dx^adx^b+r^2d\mathrm{\Omega }^2,$$

where $`d\mathrm{\Omega }^2=d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2`$ is the standard metric on the two-sphere $`S^2`$ and $`\stackrel{~}{𝒈}=\stackrel{~}{g}_{ab}dx^adx^b`$ is a metric of signature $`(1,1)`$, where the indices $`a,b`$ label some coordinates on the manifold $`\stackrel{~}{M}`$. It is not difficult to see that

$$P_{ab}=\frac{Q}{r^2}\stackrel{~}{\eta }_{ab},$$
(9)
with the other components of $`P_{\mu \nu }`$ being zero, satisfies the NED equations (5,6). Here, $`\stackrel{~}{\eta }_{ab}`$ is the volume element corresponding to $`\stackrel{~}{𝒈}`$ and $`Q`$ is a constant that represents the electric charge.

Next, Einstein’s field equations (4) reduce to

$`G_{ab}`$ $`=`$ $`2\stackrel{~}{g}_{ab},`$
$`G_{Ab}`$ $`=`$ $`0,`$
$`G_{AB}`$ $`=`$ $`2g_{AB},`$
where capital indices refer to the angular variables $`\vartheta `$ and $`\phi `$. In Schwarzschild coordinates,

$$\stackrel{~}{𝒈}=S^2Ndt^{\mathrm{\hspace{0.17em}2}}+\frac{dr^{\mathrm{\hspace{0.17em}2}}}{N},$$
(10)
these equations yield $`S=S(t)`$ (which can be set to $`1`$ by an appropriate redefinition of $`t`$) and

$$_tm=0,_rm=r^2(P),P=\frac{1}{2}\frac{Q^2}{r^4},$$
(11)
where the mass function $`m`$ is defined by $`N=12m/r`$. In particular, it follows from (11) that spacetime is static (as long as $`N>0`$), i.e. Birkhoff’s theorem HE  holds.

Given the Hamiltonian $`(P)`$, equations (11) can be integrated to yield the mass function $`m(r)`$:

$$m(r)=M+_r^{\mathrm{}}r^2𝑑r,$$
(12)
where $`M`$ is an integration constant representing the ADM mass. As a consequence of the weak field limit, we have $`m(r)MQ^2/2r`$ as $`r\mathrm{}`$ and the metric is asymptotically flat with the same asymptotic as the Reissner-Nordström metric. If the Hamiltonian $`(P)`$ is such that the integral in equation (12) converges as $`r0`$, we can choose the ADM mass such that $`m(r=0)=0`$. In particular if $`\mathrm{\Lambda }/2`$ as $`r0`$ where $`\mathrm{\Lambda }`$ is a positive constant, one obtains $`m(r)\mathrm{\Lambda }r^3/6`$, that is, the metric approaches the one of de Sitter space as $`r`$ approaches zero. Hamiltonians with a similar type of behavior at $`r=0`$ and $`r=\mathrm{}`$ have been constructed in Ref. ABG-E1 ; ABG-E2 ; ABG-E3 . As an explicit example, we consider the one that has been used in Ref. ABG-E1 :

$$(P)=P\frac{13x}{(1+x)^3}\frac{3}{2Q^2s}\left(\frac{x}{1+x}\right)^{5/2},$$
(13)
where $`x=\sqrt{2Q^2P}=Q^2/r^2`$ and $`s=|Q|/2m`$. The corresponding mass function is obtained from (12):

$$m(r)=M\frac{1}{2}\frac{Q^2r^3}{(r^2+Q^2)^2}+m\frac{r^3}{(r^2+Q^2)^{3/2}}m.$$

Thus, if $`M=m`$, the metric is regular at $`r=0`$ and so is the only non-trivial component of the electric field $`E=F_{tr}=_PP_{tr}`$:

$$E=Qr^4\left(\frac{r^25Q^2}{(r^2+Q^2)^4}+\frac{15}{2}\frac{M}{(r^2+Q^2)^{7/2}}\right).$$

Furthermore, if the charge to mass ratio $`(2s)`$ is small enough, the function $`N(r)`$ intersects the $`r`$-axis at some values $`r_c`$, $`r_h`$ with $`r_c<r_h`$, $`r_h`$ and $`r_c`$ representing the radius of an event and Cauchy horizon, respectively. Properties of the Penrose diagrams belonging to these solutions have been discussed in ABG-E1 . As expected, the corresponding spacetimes possess the same global structure as the Reissner-Nordström black holes except that there is no singularity at $`r=0`$; $`r=0`$ is the origin of the spherical coordinates.

We end this section after some comments. First, it is clear that there is a wide class of Hamiltonians that have the qualitative properties described above, since one could deform slightly the Hamiltonian (13) without affecting its asymptotic properties. Next, as observed in Br ; Br-Lett , regular electrically charged black holes are not well-defined when interpreted in the $`F`$ framework. What happens is that since $`F=E^2/2`$ is negative but vanishes as $`P`$ goes to zero or infinity, it must exhibit a minimum at some $`P^{}`$. Therefore, the condition (7) is violated at at least one point and the electromagnetic Lagrangian $`(F)`$ suffers from branching, see Ref. Br . Next, while the stress-energy tensor corresponding to (13) satisfies the weak energy condition (because $`0`$) it violates the strong energy condition since spacetime approaches the one of de Sitter near the horizon. This, together with the fact that the spacetimes possess Cauchy horizons, explains why there is no contradiction to the singularity theorems. Finally, magnetically charged solutions can also be obtained by using the $`FP`$ duality transformation (8), as mentioned above. An example of a magnetically charged solution is provided by the Bardeen model which has been interpreted as a self-gravitating nonlinear magnetic monopole in ABG-M1 .

## III Perturbation equations

In this section we investigate the stability of the solutions described in section II with respect to non-spherically symmetric dynamical linear fluctuations outside the event horizon. In order to do so, we use a recently introduced gauge-invariant perturbation formalism S-Diss ; ST  and generalize the equations obtained by Moncrief Moncrief  to nonlinear electromagnetic fields. These equations (which we shall call the pulsation equations) are a set of wave-like coupled equations with symmetric potential acting on gauge-invariant and unconstrained quantities which therefore capture the purely physical degrees of freedom associated to the linear fluctuations. Our perturbation formalism, which is described in more details in ST  for vacuum perturbations and in S-Diss  for the Einstein-Maxwell case, has the additional advantage to yield the pulsation equations for any spherically symmetric coordinates on the background.

Notice that an alternative approach is given by the Newman-Penrose formalism (the background has a metric of type D). However, we have not succeeded in decoupling the equations along the lines of C-Book . As we will see at the end of this section, one does not necessarily have the equivalence between the spectra belonging to the pulsation operators in the odd- and even-parity sectors when the electromagnetic Lagrangian is nonlinear. This might be an indication that the Newman-Penrose approach is difficult when nonlinear matter fields are coupled to the metric.

The reader who is not interested in the derivation of the pulsation equations may skip the remaining of this section.

### III.1 Preliminaries

Our task is to linearize the equations (4,5,6) around a spherically symmetric background, where we have to keep in mind that $`P_{\mu \nu }=_P^1(_\mu A_\nu _\nu A_\mu )`$. In the following, we assume that the function $`_P`$ is regular and positive for all $`r>r_h`$, where $`r_h`$ denotes the radius of the event horizon.

In terms of small perturbations of the metric and the vector potential, $`\delta g_{\mu \nu }`$ and $`\delta A_\mu `$, respectively, the linearized equations become

$`\delta G_{\mu \nu }`$ $`=`$ $`2\{_{PP}P_\mu ^\alpha P_{\nu \alpha }\delta P+_P(2P_{(\mu }^\alpha \delta P_{\nu )\alpha }P_\mu ^\sigma P_\nu ^\rho \delta g_{\sigma \rho })`$ (14)
$`g_{\mu \nu }\kappa _P\delta P\delta g_{\mu \nu }(2_PP)\},`$
$`0`$ $`=`$ $`^\mu \left(\delta P_{\alpha \mu }2\delta g_{\sigma [\alpha }P_{\mu ]}^\sigma +{\displaystyle \frac{1}{2}}g^{\sigma \rho }\delta g_{\sigma \rho }P_{\alpha \mu }\right),`$ (15)
where in the following the function

$$\kappa =1+2_P^1_{PP}P$$
(16)
will play an important role. Notice that in view of the inequality (7) and our assumption on $`_P`$ to be regular and positive, $`\kappa `$ is regular and non-zero outside the event horizon if and only if the portion of the background solution that lies outside the event horizon is well-defined in the $`F`$ framework.

The variation of $`P_{\mu \nu }`$ can be expressed as

$$\delta P_{\mu \nu }=\frac{1}{2}_P^2_{FF}\left(P^{\alpha \beta }\delta F_{\alpha \beta }_PP^{\sigma \alpha }P_\alpha ^\rho \delta g_{\sigma \rho }\right)P_{\mu \nu }+_P^1\delta F_{\mu \nu },$$
(17)
where in terms of $`(P)`$, we have

$$_{FF}=\kappa ^1_P^4_{PP}.$$

General expressions for the perturbed Einstein tensor are found, for instance, in S-Diss .

The perturbed quantities $`\delta g_{\mu \nu }`$ and $`\delta A_\mu `$ are subject to the following gauge transformations: With respect to an infinitesimal coordinate transformation $`\delta x^\mu \delta x^\mu +X^\mu `$, where $`X^\mu `$ is a vector field,

$$\delta g_{\mu \nu }\delta g_{\mu \nu }+\mathrm{\pounds }_Xg_{\mu \nu },\delta A_\mu \delta A_\mu +\mathrm{\pounds }_XA_\mu ,$$
(18)
where $`\mathrm{\pounds }_X`$ denotes the Lie derivative with respect to $`X^\mu `$. Furthermore, the potential $`A_\mu `$ is subject to $`U(1)`$ gauge transformations generated by a scalar $`\chi `$:

$$\delta A_\mu \delta A_\mu +_\mu \chi .$$
(19)
Because of the gauge invariance of the theory, the perturbation equations are invariant with respect to the transformations (18,19). Since we want to separate these transformations from the purely physical degrees of freedom, our strategy is to introduce gauge-invariant combinations of the perturbation amplitudes and to recast all equations in terms of these amplitudes.

Since the background is spherically symmetric, it is convenient to expand all perturbations in spherical tensor harmonics. Perturbations belonging to different angular momentum number $`\mathrm{}`$ and $`m`$ and with different parities decouple from each other. We discuss only the sectors with $`\mathrm{}1`$ since the spherically symmetric sector ($`\mathrm{}=0`$) has already been discussed in the previous section.

### III.2 Odd-parity sector

In each sector $`\mathrm{}m`$ linear perturbations with odd parity are parametrized by two scalar fields $`k`$ and $`\nu `$ and a one-form $`h=h_adx^a`$ on the two-dimensional manifold $`\stackrel{~}{M}`$:

$`\delta g_{ab}=0,\delta g_{Ab}=h_bS_A,\delta g_{AB}=2k\widehat{}_{(A}S_{B)},`$
$`\delta A_\mu dx^\mu =\nu S_Bdx^B,`$
where in terms of the standard spherical harmonics $`Y=Y^\mathrm{}m`$ the transverse vector harmonics $`S_B`$ are defined by $`S_B=(\widehat{}dY)_B=\widehat{\eta }_B^A\widehat{}_AY`$. Here and in the following, quantities with a hat refer to the standard metric $`\widehat{𝒈}=d\mathrm{\Omega }^2`$ on $`S^2`$ and quantities with a tilde refer to the metric $`\stackrel{~}{𝒈}`$ on $`\stackrel{~}{M}`$. It is not difficult to show that $`\nu `$ is invariant with respect to both coordinate and $`U(1)`$ gauge transformations. In the gravitational sector, we can construct the gauge-invariant one-form (see S-Diss ; ST )

$$h^{(inv)}=hr^2d\left(\frac{k}{r^2}\right).$$
(20)
Since the perturbation equations are gauge invariant we can choose a gauge in which $`k=0`$ and $`h^{(inv)}=h`$ to simplify the calculations and then replace $`h`$ by $`h^{(inv)}=hr^2d(r^2k)`$ in the final result.

#### III.2.1 Gravito-electromagnetic equations

Since in the odd-parity sector all scalar perturbations vanish, the expressions (14,15,17) simplify considerably in this case. For example we have $`\delta P=0`$, and so $`\delta P_{\mu \nu }=_P^1(_\mu \delta A_\nu _\nu \delta A_\mu )`$. Explicitly,

$$\delta P_{\mu \nu }=_P^1\left[d\nu S_Bdx^B\mathrm{}(\mathrm{}+1)\nu Yd\mathrm{\Omega }\right]_{\mu \nu },$$
(21)
where we have used $`d(S_Bdx^B)=\widehat{\mathrm{\Delta }}Yd\mathrm{\Omega }=\mathrm{}(\mathrm{}+1)Yd\mathrm{\Omega }`$.

Introducing equation (21) into formula (15), using the background expressions (9) and the fact that $`g^{\alpha \beta }\delta g_{\alpha \beta }=0`$, the linearized electromagnetic equations yield

$$\stackrel{~}{d}^{}\left(_P^1d\nu \right)+\frac{\mathrm{}(\mathrm{}+1)}{r^2}_P^1\nu +Q\stackrel{~}{}d\left(\frac{h}{r^2}\right)=0,$$

where $`\stackrel{~}{}`$ and $`\stackrel{~}{d}^{}\stackrel{~}{}d\stackrel{~}{}`$ denote the Hodge dual and the co-differential, respectively, with respect to the metric $`\stackrel{~}{𝒈}`$.333In a coordinate-dependent notation, one has $`\stackrel{~}{d}^{}d\nu =\stackrel{~}{}^a\stackrel{~}{}_a\nu `$, $`\stackrel{~}{}dh=\stackrel{~}{\eta }^{ab}\stackrel{~}{}_ah_b`$, $`(\stackrel{~}{d}^{}dh)_a=2\stackrel{~}{}^b\stackrel{~}{}_{[a}h_{b]}`$, $`\stackrel{~}{}d\nu _a=\stackrel{~}{\eta }_a^b\stackrel{~}{}_b\nu `$, etc.  Next, from expressions (14) and (21) we find

$$\delta G_{Ab}dx^b=4\left\{Q(\stackrel{~}{}d\nu )r^2h\right\}\frac{S_A}{2r^2},\delta G_{ab}=0,\delta G_{AB}=0.$$

Using this in connection with the corresponding expressions for the perturbed Einstein tensor in the odd parity sector (see Ref. S-Diss ),

$`\delta G_{Ab}dx^b`$ $`=`$ $`\left\{\stackrel{~}{d}^{}\left[r^4d\left({\displaystyle \frac{h}{r^2}}\right)\right]+\left(\mathrm{}(\mathrm{}+1)2+r^2G_B^B\right)h\right\}{\displaystyle \frac{S_A}{2r^2}},`$
$`\delta G_{AB}`$ $`=`$ $`\stackrel{~}{d}^{}h\widehat{}_{(A}S_{B)},`$ (22)
we end up with the manifestly gauge-invariant equations

$`\stackrel{~}{d}^{}\left[r^4d\left({\displaystyle \frac{h^{(inv)}}{r^2}}\right)\right]+\lambda h^{(inv)}`$ $`=`$ $`4Q\stackrel{~}{}d\nu ,`$ (23)
$`\stackrel{~}{d}^{}\left[_P^1d\nu \right]+{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}_P^1\nu `$ $`=`$ $`Q\stackrel{~}{}d\left({\displaystyle \frac{h^{(inv)}}{r^2}}\right),`$ (24)
$`\stackrel{~}{d}^{}h^{(inv)}`$ $`=`$ $`0,`$ (25)
where here and in the following $`\lambda (\mathrm{}1)(\mathrm{}+2)`$. Equation (25) is void if $`\mathrm{}=1`$ since in this case, $`\widehat{}_{(A}S_{B)}=0`$.

#### III.2.2 Pulsation equations ($`\mathrm{}2`$)

We first consider the case $`\mathrm{}2`$. Then, equation (25) allows to express $`h^{(inv)}`$ in terms of a potential $`\mathrm{\Psi }`$,

$$h^{(inv)}=\frac{1}{\sqrt{\lambda }}\stackrel{~}{}d(r\mathrm{\Psi }),$$

where the factors $`1/\sqrt{\lambda }`$ and $`r`$ turn out to be convenient. Introducing also $`\mathrm{\Phi }=\sqrt{4_P^1}\nu `$, we can integrate equation (23) and eventually obtain the following pulsation equations444When integrating equation (23), an integration constant is set to zero by appropriately choosing the free constant in the potential $`\mathrm{\Psi }`$.:

$`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Psi }+\left[{\displaystyle \frac{\lambda }{r^2}}+r\stackrel{~}{\mathrm{\Delta }}\left({\displaystyle \frac{1}{r}}\right)\right]\mathrm{\Psi }{\displaystyle \frac{\sqrt{4\lambda _P}Q}{r^3}}\mathrm{\Phi }=0,`$
$`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Phi }{\displaystyle \frac{\sqrt{4\lambda _P}Q}{r^3}}\mathrm{\Psi }+\left[{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}+_P^{1/2}\stackrel{~}{\mathrm{\Delta }}_P^{1/2}+{\displaystyle \frac{4Q^2}{r^4}}_P\right]\mathrm{\Phi }=0,`$
where $`\stackrel{~}{\mathrm{\Delta }}\stackrel{~}{d}^{}d=\stackrel{~}{}^a\stackrel{~}{}_a`$ is the Laplace operator with respect to the metric $`\stackrel{~}{𝒈}`$. These equations have the form of a coupled system of wave equations with symmetric potential. We can rewrite them into the form

$$\left(𝑫+𝑺\right)\left(\begin{array}{c}\mathrm{\Psi }\\ \mathrm{\Phi }\end{array}\right)=0,$$
(26)
where

$`𝑫`$ $`=`$ $`\left(\begin{array}{cc}r\stackrel{~}{}^a\frac{1}{r^2}\stackrel{~}{}_ar& 0\\ 0& _P^{1/2}\stackrel{~}{}^a_P^1\stackrel{~}{}_a_P^{1/2}\end{array}\right),`$ (29)
$`𝑺`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}\left(\begin{array}{cc}\lambda & \frac{\sqrt{4\lambda _P}Q}{r}\\ \frac{\sqrt{4\lambda _P}Q}{r}& \mathrm{}(\mathrm{}+1)+\frac{4Q^2}{r^2}_P\end{array}\right).`$ (32)
As we will show in the next section in more detail, stability follows from the fact that the spatial part of the operator $`𝑫`$ and the matrix $`𝑺`$ are manifestly positive.

#### III.2.3 Pulsation equations ($`\mathrm{}=1`$)

In the sector $`\mathrm{}=1`$, $`\lambda `$ vanishes and equation (23) can be integrated to yield

$$r^4\stackrel{~}{}d\left(\frac{h}{r^2}\right)=4Q\nu +6J,$$

where $`J`$ is an integration constant. Using this into equation (24) gives

$$_P\stackrel{~}{d}^{}\left[_P^1d\nu \right]+\left(\frac{2}{r^2}+4_P\frac{Q^2}{r^4}\right)\nu =6_P\frac{QJ}{r^4}.$$
(33)
The general solution to this equation is the sum of a particular solution $`\nu _1`$ and the general homogeneous solution. The particular solution, which is proportional to $`J`$, can be chosen to be a function of $`r`$ only. Therefore, $`\nu _1`$ represents a stationary axial excitation of the background configuration. For $`_P=1`$ an explicit solution is given by $`\nu _1=J/r`$ which corresponds to the Kerr-Newman solution to linear order in the rotation parameter $`J/M`$ S-Diss . For a NED, we have not found an explicit expression for $`\nu _1`$, but we anticipate that $`\nu _1`$ may give rise to slowly rotating generalizations of the solutions found in ABG-E1 ; ABG-E2 ; ABG-E3 . The homogeneous solutions represent electromagnetic radiation which propagates in a stable way since the potential on the left-hand side of equation (33) is positive (see also the next section).

### III.3 Even-parity sector

Gravitational perturbations are more complicated in the even-parity sector than in the odd-parity one. The metric is now parametrized by a symmetric $`2`$ tensor, $`H_{ab}`$, a one-form $`Q_a`$ and two scalar fields $`K`$ and $`G`$ in the following way:

$`\delta g_{ab}`$ $`=`$ $`H_{ab}Y,`$
$`\delta g_{aB}`$ $`=`$ $`Q_a\widehat{}_BY,`$
$`\delta g_{AB}`$ $`=`$ $`Kg_{AB}Y+Gr^2\left(\widehat{}_A\widehat{}_B{\displaystyle \frac{1}{2}}\widehat{g}_{AB}\widehat{\mathrm{\Delta }}\right)Y,`$
where $`G`$ is not present for $`\mathrm{}=1`$ since $`\left(\widehat{}_A\widehat{}_B\frac{1}{2}\widehat{g}_{AB}\widehat{\mathrm{\Delta }}\right)Y=0`$ in this case. For $`\mathrm{}2`$ one can introduce the following gauge invariants GS :

$$H_{ab}^{(inv)}=H_{ab}2\stackrel{~}{}_{(a}p_{b)},K^{(inv)}=K2\frac{p^a\stackrel{~}{}_ar}{r}+\frac{1}{2}\mathrm{}(\mathrm{}+1)G,$$

where $`p_a=Q_a(r^2\stackrel{~}{}_aG)/2`$. With respect to an infinitesimal coordinate transformation parametrized by the vector field $`X=\xi ^aY_a+f\widehat{}^AY_A`$ we have

$$p_ap_a+\xi _a,GG+2f.$$
(34)
Since $`G`$ and $`p_a`$ can both be set to zero by an appropriate coordinate transformation, we can assume that these quantities vanish when deriving the pulsation equations below. In this gauge, $`H_{ab}^{(inv)}=H_{ab}`$, $`K^{(inv)}=K`$ and $`Q_a`$ and $`G`$ vanish. Since the perturbation equations are gauge invariant, it is sufficient to replace $`H_{ab}`$ by $`H_{ab}^{(inv)}`$ and $`K`$ by $`K^{(inv)}`$ in the final expressions and so obtain the equations in an arbitrary gauge.

The electromagnetic potential is parametrized according to

$$\delta A_\mu dx^\mu =\alpha Y+\mu dY,$$

with $`\alpha =\alpha _bdx^b`$ a one-form and $`\mu `$ a scalar on $`\stackrel{~}{M}`$. The corresponding electromagnetic field is

$$\delta F_{\mu \nu }=_\mu \delta A_\nu _\nu \delta A_\mu =\left[d\alpha Y+(d\mu \alpha )dY\right]_{\mu \nu },$$

from which it is clear that the amplitude $`\widehat{\alpha }=\alpha d\mu `$ is $`U(1)`$-invariant. Using the background expression (9) and $`\mathrm{\pounds }_X𝑭=di_X𝑭=d(r^2_PQ\stackrel{~}{\eta }_{ab}X^adx^b)`$ we find that with respect to an infinitesimal change of coordinates,

$$\widehat{\alpha }\widehat{\alpha }+_P\frac{Q}{r^2}\stackrel{~}{}(\xi _adx^a).$$

In view of the transformations (34) we see that we can construct a coordinate invariant one-form which reduces to $`\widehat{\alpha }`$ if $`p_a=0`$.

#### III.3.1 Electromagnetic equations

Having discussed the gauge issues, we now proceed to the derivation of the perturbation equations. First, a small calculation using (17) yields

$$\delta P_{\mu \nu }=\left[\stackrel{~}{}\pi Y+\gamma dY\right]_{\mu \nu },$$

with

$`\pi `$ $`=`$ $`{\displaystyle \frac{_P^1}{\kappa }}\stackrel{~}{}d\widehat{\alpha }+{\displaystyle \frac{\kappa 1}{2\kappa }}{\displaystyle \frac{Q}{r^2}}H,`$ (35)
$`\gamma `$ $`=`$ $`_P^1\widehat{\alpha },`$ (36)
where $`\kappa `$ is defined in (16) and $`H=H_a^a`$. In terms of the scalar $`\pi `$ and the one-form $`\gamma `$ we can write the electromagnetic equations as

$`\stackrel{~}{}d\left[r^2\pi +{\displaystyle \frac{Q}{2}}(H2K)\right]\mathrm{}(\mathrm{}+1)\gamma `$ $`=`$ $`0,`$ (37)
$`\stackrel{~}{d}^{}\gamma `$ $`=`$ $`0.`$ (38)
Equation (38) is the integrability condition for equation (37) and can be used in order to introduce a scalar potential $`\varphi `$ according to

$$\gamma =\stackrel{~}{}d\varphi .$$
(39)
Equation (37) can then be integrated to obtain the algebraic relation555The integration constant is absorbed into the definition of $`\varphi `$.

$$r^2\pi =\mathrm{}(\mathrm{}+1)\varphi +\frac{Q}{2}(H2K).$$

Using the equations (35) and (36) one then obtains the following pulsation equations for the electromagnetic field

$$_P^1\stackrel{~}{d}^{}\left(_Pd\varphi \right)+\kappa \frac{\mathrm{}(\mathrm{}+1)}{r^2}\varphi +\frac{Q}{2r^2}\left(2\kappa KH\right)=0.$$
(40)
If we turn off the gravitational amplitudes (i.e. if $`H=K=0`$), we see that the sign of the potential depends on the function $`\kappa `$ defined in (16). In particular, the equations are stable as long as $`\kappa 0`$. In the next section we will show that this is also a necessary condition for stability, i.e. we will show that if $`\kappa <0`$ in an interval, the equations admit exponentially growing modes for sufficiently large $`\mathrm{}`$.

#### III.3.2 Gravitational equations

Since the electromagnetic field couples back to the gravitational field, we have to consider the linearized Einstein equations (14) as well. The harmonic decomposition of the linearized Einstein tensor can be found in S-Diss . The result is quite complicated, but the structure of the equations becomes somewhat more transparent if one splits the tensor $`H_{ab}`$ into its trace and traceless part and then introduces the one-form

$$C=H_{ab}^{TF}(\stackrel{~}{}^ar)dx^b,$$

where $`H_{ab}^{TF}=H_{ab}\frac{1}{2}\stackrel{~}{g}_{ab}H`$ denotes the traceless part of $`H_{ab}`$. The relevant components of the linearized Einstein tensor may then be expressed in the following way:

$`(\delta G)_{AB}^{TF}`$ $`=`$ $`S\left(\widehat{}_A\widehat{}_B{\displaystyle \frac{1}{2}}\widehat{g}_{AB}\widehat{\mathrm{\Delta }}\right)Y,`$
$`\stackrel{~}{g}^{ab}\delta G_{ab}`$ $`=`$ $`\left(T+G^{ab}H_{ab}^{TF}\right)Y,`$
$`\delta G_{Ab}dx^b`$ $`=`$ $`\left(U+U_{BG}\right){\displaystyle \frac{1}{2}}\widehat{}_AY,`$
$`(\delta G)_{ab}^{TF}r^{|a}dx^b`$ $`=`$ $`\left(V{\displaystyle \frac{1}{2}}G^{ab}H_{ab}^{TF}dr+G_{ab}(C^a{\displaystyle \frac{r}{2}}\stackrel{~}{}^aK)dx^b\right)Y,`$
where

$$U_{BG}=\frac{r}{2N}\left(G^{ab}dr+(\stackrel{~}{\eta }^{ac}G_c^b)\stackrel{~}{}dr\right)H_{ab}^{TF},$$

and

$`S`$ $`=`$ $`{\displaystyle \frac{1}{2}}H,`$ (41)
$`T`$ $`=`$ $`{\displaystyle \frac{2}{r}}\stackrel{~}{d}^{}C{\displaystyle \frac{2}{r^2}}\stackrel{~}{𝒈}(C,dr)+\stackrel{~}{\mathrm{\Delta }}K+{\displaystyle \frac{4}{r}}\stackrel{~}{𝒈}(dK,dr){\displaystyle \frac{\lambda }{r^2}}K{\displaystyle \frac{\lambda +4}{2r^2}}H,`$ (42)
$`U`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left[(\stackrel{~}{d}^{}C)dr+(\stackrel{~}{}dC)\stackrel{~}{}dr\right]d\left(K+{\displaystyle \frac{1}{2}}H\right)+H{\displaystyle \frac{dr}{r}},`$ (43)
$`V`$ $`=`$ $`(\stackrel{~}{d}^{}C){\displaystyle \frac{dr}{r}}+{\displaystyle \frac{1}{r}}d\stackrel{~}{𝒈}(C,dr)+{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{2r^2}}C`$ (44)
$`+`$ $`{\displaystyle \frac{1}{2}}\stackrel{~}{\mathrm{\Delta }}Kdrd\stackrel{~}{𝒈}(dK,dr)+\left(\stackrel{~}{\mathrm{\Delta }}r{\displaystyle \frac{N+1}{2r}}\right)dK+{\displaystyle \frac{N}{2r}}dH.`$
Here, $`N=(\stackrel{~}{}^ar)(\stackrel{~}{}_ar)`$ agrees with the function $`N`$ defined in the previous section, see equation (10). Note that equation (41) is void when $`\mathrm{}=1`$ since then $`\widehat{}_A\widehat{}_BY=3\widehat{g}_{AB}Y`$.

In our case $`G_{ab}=2\stackrel{~}{g}_{ab}`$ is proportional to $`\stackrel{~}{g}_{ab}`$, so we have $`G^{ab}H_{ab}^{TF}=0`$ and $`U_{BG}=0`$. Next, it is not difficult to see from the right-hand side of equation (14) that $`\delta G_{AB}`$ is proportional to $`g_{AB}`$ so that $`S=0`$. Then, it follows from equation (41) that the trace of $`H_{ab}`$ vanishes, $`H=0`$. (When $`\mathrm{}=1`$, $`H`$ can be set to zero by an appropriate coordinate transformation). The remaining components of the right-hand side of (14) yield the expressions

$`T`$ $`=`$ $`4_P{\displaystyle \frac{Q}{r^2}}\pi ,`$ (45)
$`U`$ $`=`$ $`4_P{\displaystyle \frac{Q}{r^2}}d\varphi ,`$ (46)
$`V`$ $`=`$ $`rdK,`$ (47)
where equation (39) has been used in the derivation for $`U`$. Following the vacuum case S-Diss ; ST  we first look at the component $`\stackrel{~}{𝒈}(U,\stackrel{~}{}dr)`$: Using equations (43) and (46), we find

$$\stackrel{~}{}d[CrdK]=4\stackrel{~}{}d\left[_P\varphi \frac{Q}{r^2}dr\right],$$

where we have used the identity $`\stackrel{~}{𝒈}(\alpha ,\stackrel{~}{}dr)=\stackrel{~}{}(\alpha dr)`$ for any one-form $`\alpha `$. This motivates the introduction of a one-form $`Z`$ defined by

$$Z=CrdK+4_P\varphi \frac{Q}{r^2}dr.$$

Since $`dZ=0`$ we can introduce a potential $`\zeta `$ with $`Z=d\zeta `$. Next, we look at the combination $`2VTdr`$. Using equations (42,44) and (45,47), we find

$`d\left[2r\stackrel{~}{𝒈}(Z4_P\varphi {\displaystyle \frac{Q}{r^2}}dr,dr)\right]+\mathrm{}(\mathrm{}+1)Z`$ (48)
$`+`$ $`r\left(\lambda +3(1N)+2r^2\right)dK+\left(\lambda +4_P{\displaystyle \frac{Q^2}{r^2}}\right)Kdr=0.`$
By virtue of the background equations we have

$`_r\left[r\left(\lambda +3(1N)+2r^2\right)\right]`$ $`=`$ $`\lambda +6_rm+6r^2+2r^3_P_rP`$
$`=`$ $`\lambda +4_P{\displaystyle \frac{Q^2}{r^2}},`$
and since $`Z=d\zeta `$, we can integrate equation (48) to obtain666The integration constant is absorbed into the definition of $`\zeta `$.

$$2r\stackrel{~}{𝒈}(d\zeta ,dr)+\mathrm{}(\mathrm{}+1)\zeta +\left(r\lambda +6m+2r^3\right)K8N_P\frac{Q}{r}\varphi =0.$$
(49)
Finally, we compute the combination $`rT+\stackrel{~}{𝒈}(U,dr)`$ which yields, using the equations (42,43) and (45,46),

$$\stackrel{~}{d}^{}Z\frac{2}{r}\stackrel{~}{𝒈}(Z,dr)\left(\frac{\lambda }{r}+4_P\frac{Q^2}{r^3}\right)K+4\left[\stackrel{~}{d}^{}\left(_P\frac{Q}{r^2}dr\right)+2N_P\frac{Q}{r^3}\frac{\mathrm{}(\mathrm{}+1)}{r^2}\frac{Q}{r}_P\right]\varphi =0.$$
(50)
This equation may be used in order to express the scalar $`K`$ in terms of the other perturbation amplitudes.

To summarize, we have the perturbation equations (40) for the electromagnetic perturbation potential $`\varphi `$ and equation (49) for the gravitational perturbation potential $`\zeta `$ where $`H=0`$ and $`K`$ can be eliminated using equation (50). The next task is to bring these equations in form of a wave-like equation with symmetric potential which is more suitable for a discussion of stability.

#### III.3.3 Gravito-electromagnetic equations

The pulsation equations can be brought into the form of a wave-like equation with symmetric potential after the following transformations:

$$\mathrm{\Psi }=\sqrt{\lambda }\frac{\zeta }{a+\lambda },\mathrm{\Phi }=\sqrt{4_P}\left(\varphi \frac{Q}{r}\frac{\zeta }{a+\lambda }\right),$$
(51)
where $`a=3r\stackrel{~}{\mathrm{\Delta }}r4r^2=6m/r+2r^2`$. For the moment, we restrict ourselves to the case $`\mathrm{}2`$. Note that $`a`$ is positive outside the horizon of a regular black hole (see below).

In terms of the new potentials defined in (51), a long calculation yields the final pulsation equations

$$\stackrel{~}{\mathrm{\Delta }}\left(\begin{array}{c}\mathrm{\Psi }\\ \mathrm{\Phi }\end{array}\right)+\left(\begin{array}{cc}V_{11}& \frac{\sqrt{4\lambda _P}Q}{r^3}W\\ \frac{\sqrt{4\lambda _P}Q}{r^3}W& V_{22}\end{array}\right)\left(\begin{array}{c}\mathrm{\Psi }\\ \mathrm{\Phi }\end{array}\right)=0,$$
(52)
where

$`V_{11}`$ $`=`$ $`{\displaystyle \frac{1}{r^2(a+\lambda )}}\left[\mathrm{}(\mathrm{}+1)\lambda 2N\lambda +ar\stackrel{~}{\mathrm{\Delta }}r\right]+{\displaystyle \frac{2N\lambda b}{r^2(a+\lambda )^2}},`$
$`V_{22}`$ $`=`$ $`\kappa {\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}+{\displaystyle \frac{4_PQ^2}{r^4(a+\lambda )}}\left(\lambda +1N2r^2+4N\kappa \right)+_P^{1/2}\stackrel{~}{\mathrm{\Delta }}_P^{1/2}+{\displaystyle \frac{8N_PQ^2b}{r^4(a+\lambda )^2}},`$
$`W`$ $`=`$ $`{\displaystyle \frac{1}{a+\lambda }}\left(\lambda +1N2r^2+2N\kappa \right)+{\displaystyle \frac{2Nb}{(a+\lambda )^2}},`$
and where we have defined

$$b=\lambda +4_P\frac{Q^2}{r^2}.$$

Here, we have also used the background expressions

$`r\stackrel{~}{\mathrm{\Delta }}r=r_rN`$ $`=`$ $`1N+2r^2,`$
$`r_P^1\stackrel{~}{d}^{}(_Pdr)`$ $`=`$ $`1N2r^2+2N\kappa .`$
The equations (52) generalize the pulsation equations obtained by Moncrief Moncrief  to NED coupled to gravity. (For $`=P`$ the equations above reduce to the one obtained in Moncrief .)

#### III.3.4 Pulsation equations ($`\mathrm{}2`$)

We now proceed to show under which conditions the potential in the equation (52) is positive semi-definite.

We first show that the function $`a=6m/r+2r^2`$ is positive outside an event horizon: Let $`r=r_h`$ denote the radius of the event horizon. Then, $`m(r_h)=r_h/2`$ and $`0r_h_rN(r_h)=2r_h^2(r=r_h)+1`$. Therefore, $`a(r=r_h)2`$. Next, using the field equations $`_rm=r^2`$, we have $`_r(ra)=4_PQ^2/r^2>0`$. This shows that $`a`$ cannot become negative or zero for $`r>r_h`$. In order to prove that the potential is positive, we first absorb a term $`r^1\stackrel{~}{\mathrm{\Delta }}r`$ appearing in $`V_{11}`$ and the term $`_P^{1/2}\stackrel{~}{\mathrm{\Delta }}_P^{1/2}`$ appearing in $`V_{22}`$ into the differential operator. We may then rewrite the pulsation equations into the form

$$\left(𝑫+𝑺\right)\left(\begin{array}{c}\mathrm{\Psi }\\ \mathrm{\Phi }\end{array}\right)=0,$$
(53)
where

$`𝑫`$ $`=`$ $`\left(\begin{array}{cc}\frac{1}{r}\stackrel{~}{}^ar^2\stackrel{~}{}_a\frac{1}{r}& 0\\ 0& _P^{1/2}\stackrel{~}{}^a_P\stackrel{~}{}_a_P^{1/2}\end{array}\right),`$ (56)
$`𝑺`$ $`=`$ $`\left(\begin{array}{cc}\frac{\lambda }{r^2(a+\lambda )}\left[c_1+\frac{2Nb}{a+\lambda }\right]& \frac{\sqrt{4\lambda _P}Q}{r^3(a+\lambda )}\left[w+\frac{2Nb}{a+\lambda }\right]\\ \frac{\sqrt{4\lambda _P}Q}{r^3(a+\lambda )}\left[w+\frac{2Nb}{a+\lambda }\right]& \kappa \frac{\mathrm{}(\mathrm{}+1)}{r^2}+\frac{4_PQ^2}{r^4(a+\lambda )}\left[c_2+\frac{2Nb}{a+\lambda }\right]\end{array}\right),`$ (59)
with

$`c_1`$ $`=`$ $`\lambda +1N2r^2,`$
$`c_2`$ $`=`$ $`c_1+4N\kappa ,`$
$`w`$ $`=`$ $`c_1+2N\kappa .`$
Assume that $`0`$ and $`\kappa 0`$ for all $`r>r_h`$. It then follows that $`a`$ ,$`b`$ and $`c_1`$ are positive since $`N1`$. Therefore, the diagonal elements of $`𝑺`$ are positive. In order to show that $`𝑺`$ is positive definite, it remains to check the sign of its determinant. A small calculation yields

$$\frac{r^4(a+\lambda )^2}{\lambda }det𝑺=\kappa \mathrm{}(\mathrm{}+1)\left[(a+\lambda )c_1+2\lambda N\right]+8N\kappa _P\frac{Q^2}{r^2}\left[\mathrm{}(\mathrm{}+1)2N\kappa \right].$$

This is manifestly positive if

$$0<2N\kappa \mathrm{}(\mathrm{}+1)$$
(60)
(if $`\kappa =0`$ $`𝑺`$ is only positive semi-definite.) In particular, this condition is satisfied in linear electrodynamics since then $`\kappa =1`$. Provided that the weak field limit holds and that $`\kappa `$ is regular and positive near the event horizon, we have $`N\kappa 1`$ as $`r\mathrm{}`$ and $`N\kappa 0`$ as $`rr_h`$, so in those limits the condition (60) is satisfied.

#### III.3.5 Pulsation equations ($`\mathrm{}=1`$)

For $`\mathrm{}=1`$ it does not make sense to rescale $`\zeta `$ by $`\sqrt{\lambda }`$ since $`\lambda =0`$ in this case. Considering $`\mathrm{\Psi }=\zeta /a`$ instead of the $`\mathrm{\Psi }`$ defined in (51), we see that the equation for $`\mathrm{\Phi }`$ decouples when $`\mathrm{}=1`$. Repeating the arguments above, it is easy to check that the corresponding potential is positive if $`\kappa >0`$. The equation for $`\mathrm{\Psi }`$ can be rewritten as

$$\frac{1}{r}\stackrel{~}{}^a\left[r^2\stackrel{~}{}_a\left(\frac{\mathrm{\Psi }}{r}\right)\right]=\frac{\sqrt{4_P}Q}{r^3a}\left[w+\frac{2Nb}{a}\right]\mathrm{\Phi }.$$
(61)
On the other hand, the perturbation amplitudes are not necessarily gauge invariant when $`\mathrm{}=1`$. Indeed, as shown in S-Diss , the potentials can be seen to transform according to

$$\zeta \zeta +raf,\varphi \varphi +Qf$$

with respect to an infinitesimal coordinate transformation parametrized by the vector field $`X=r^2\stackrel{~}{}^afY_a+f\widehat{}^AY_A`$ where $`f`$ satisfies $`\stackrel{~}{}^a(r^2\stackrel{~}{}_af)=0`$ (This makes sure that we remain in a gauge with $`H=0`$.) Therefore, $`\mathrm{\Psi }\mathrm{\Psi }+rf`$ while the electromagnetic potential $`\mathrm{\Phi }`$ is gauge invariant. Having solved for $`\mathrm{\Phi }`$, equation (61) can be used in order to obtain the gauge-invariant part of $`\mathrm{\Psi }`$.

## IV Stability analysis

In the previous section, we have shown that linear fluctuations around a spherically symmetric and purely electric solution are governed by a wave-like equation with symmetric potential of the form

$$\left(𝑷\stackrel{~}{}^a𝑷^2\stackrel{~}{}_a𝑷+𝑺\right)u=0,$$
(62)
where $`\stackrel{~}{}`$ denotes the covariant derivative with respect to the metric $`\stackrel{~}{𝒈}`$ and where $`𝑷`$ is a positive definite symmetric matrix and $`𝑺`$ is a symmetric matrix. The vector-valued function $`u`$ on $`\stackrel{~}{M}`$ is a gauge-invariant combination of the perturbed metric and electromagnetic fields. We will show below that linear stability follows if the potential $`𝑺`$ is positive definite. In the previous section, we have shown that this is the case if the following conditions hold outside the event horizon:

$`<0,_P>0,`$ (63)
$`0<N\kappa 3,`$ (64)
where $`\kappa =1+2_P^1_{PP}P`$. The first two conditions (63) imply that the weak energy condition is satisfied outside the horizon, while condition (64) implies that the function $`F(P)`$ is well defined outside the horizon, see (7). We have checked that these conditions are verified for the Hamiltonian (13) used in ABG-E1 , for several values of the parameter $`s`$. However, it does not seem to be difficult to modify the function $`(P)`$ slightly and to violate, for example, the condition $`\kappa >0`$ in a small interval outside the event horizon. As we will show below, the background solution is unstable if $`\kappa <0`$ at some point outside the event horizon.

One way of analyzing the behavior of the solutions to equation (62) is by introducing Schwarzschild coordinates $`(t,r)`$, see equation (10), in which case (62) can be written as

$$\left(_t^{\mathrm{\hspace{0.33em}2}}𝑷_r_{}𝑷^2_r_{}𝑷+N𝑺\right)u=0,$$
(65)
where $`_r_{}N_r`$ denotes the derivative with respect to the tortoise coordinate. (Notice that $`𝑷`$ does not depend on $`t`$.) One then considers perturbations which are smooth and vanish at the boundary points $`r=r_h`$ ($`r_{}=\mathrm{}`$) and $`r=\mathrm{}`$ ($`r_{}=\mathrm{}`$). The spatial part of the wave operator in (65) is symmetric with respect to the $`L^2`$ product

$$(u,v)=_{\mathrm{}}^{\mathrm{}}u^{}v𝑑r_{}$$

for such perturbations. Since that operator is also real, self-adjoint extensions exist (in this case the operator is self-adjoint on some Sobolev space), and the linear stability can be analyzed along the lines of Wald ; VG  using the tools from spectral theory. As a first step one can look at solutions of (65) with time dependency of the form $`ue^{i\omega t}`$, in which case

$$\omega ^2u=\left(𝑷_r_{}𝑷^2_r_{}𝑷+N𝑺\right)u.$$
(66)
If the operator on the right-hand side of (66) is positive, which is the case if $`𝑺`$ is positive definite, it follows that $`\omega `$ is real and positive and there are no unstable modes (if $`\omega =0`$ there might be modes with linear growth in $`t`$). On the other hand, if $`𝑺`$ is negative the operator is not necessarily positive. If one can find an eigenfunction with negative eigenvalue, $`\omega `$ is complex and one has a solution that grows exponentially in time. While it is difficult to solve the eigenvalue problem explicitly in general, there are simple analytical methods to show the existence of at least one unstable mode VG . In fact, if one can find a function $`u`$ with $`(u,u)<\mathrm{}`$ and $`(u,Hu)<0`$, this establishes the existence of at least one negative eigenvalue and the background solution is linearly unstable.

However, the analysis of the eigenvalue problem (66) only gives information on the point spectrum (or on quasi-normal frequencies, but the corresponding modes $`u`$ do not satisfy $`(u,u)<\mathrm{}`$). Generalizing the results obtained in Wald  on should be able to show the following stronger stability result: If the potential $`N𝑺`$ is strictly positive for all $`r>r_h`$ and if $`u(t=0,.)`$, $`_tu(t=0,.)`$ are smooth and have compact support, there is a bound $`C<\mathrm{}`$ such that

$$|u(t,r_{})|C$$

for all $`t0`$ and $`\mathrm{}<r_{}<\mathrm{}`$. (If $`N𝑺`$ is nonnegative but fails to be strictly positive outside the event horizon, the perturbations might exhibit linear growth in time.)

We now use the above arguments to show that if $`\kappa <0`$, the solution is unstable: Suppose that $`\kappa _1\kappa (r_1)<0`$ for some $`r_1>r_h`$. Assuming that $`\kappa `$ is continuous, there is a $`\delta >0`$ small enough such that $`\kappa <\kappa _1/2`$ on the interval $`(r_1\delta ,r_1+\delta )`$ and $`\kappa <0`$ on the interval $`(r_12\delta ,r_1+2\delta )`$. Consider the test function

$$u=(\mathrm{\Psi },\mathrm{\Phi })=(\sqrt{4\lambda ^1_P}Qf(r),rf(r)),$$

for equation (53) where the function $`f(r)`$ is smooth, vanishes identically outside the interval $`(r_12\delta ,r_1+2\delta )`$ and is identically one for $`r(r_1\delta ,r_1+\delta )`$. It is clear that $`(u,u)<\mathrm{}`$, and $`u^{}𝑺u=\kappa \mathrm{}(\mathrm{}+1)f^2`$. Furthermore, a small calculation gives

$$(u,Hu)4\delta \left(1+\frac{4}{\lambda }\right)K^2\delta |\kappa _1|\mathrm{}(\mathrm{}+1),$$

where $`K`$ is a constant which is independent on $`\mathrm{}`$ but large enough, such that

$$\left|Qr_r_{}\left(_P^{1/2}r^1f\right)\right|K,\left|_P^{1/2}_r_{}\left(_P^{1/2}rf\right)\right|K.$$

It is clear that $`(u,Hu)<0`$ for large enough angular momentum number $`\mathrm{}`$. Therefore, the background configuration is linearly unstable if $`\kappa `$ is negative somewhere outside the event horizon.

## V Conclusion

We have derived the pulsation equations governing linear fluctuations on a spherically symmetric and electrically charged black hole solution in NED coupled to gravity. Under a mild technical assumption (namely that $`_P`$ vanishes nowhere outside the event horizon) the pulsation equations in the odd-parity sector cannot exhibit modes that grow exponentially in time outside the event horizon. However, in the even-parity sector, the stability properties depend strongly on the Hamiltonian $`(P)`$, and one does not necessarily have the equivalence between the spectra in the odd- and even-parity sectors via a supersymmetric transformation like for linear perturbations in Einstein-Maxwell theory C-Book . As a consequence, the quasi-normal frequencies might be different in the odd- and even-parity sectors if the electromagnetic Lagrangian is not linear.

In terms of the dimensionless variable $`x=\sqrt{2Q^2P}=Q^2/r^2`$, sufficient conditions for linear stability are

$``$ $`<`$ $`0,`$ (67)
$`_x`$ $`<`$ $`0,`$ (68)
$`_{xx}`$ $`<`$ $`0,`$ (69)
$`3_x`$ $``$ $`xN_{xx},`$ (70)
for all $`0<x<x_h`$ (i.e. $`r>r_h`$), where $`x_h=Q^2/r_h^2`$ corresponds to the value of $`x`$ at the event horizon. While we have shown that a solution is unstable if $`_{xx}>0`$ for some $`x<x_h`$ (the case $`_{xx}0`$ corresponding to a marginal case where growth which is linear in time might appear) we stress that the conditions (67,68,70) are only sufficient. For example, there might exist stable black holes that satisfy $`_{xx}<0`$ but violate the condition (70).

The linear stability of the purely magnetic solutions considered in ABG-M1 ; Br  can be discussed on the same lines as above. In fact, using the $`FP`$ duality transformation (8), one can easily translate the conditions (67,68,69,70) to obtain the sufficient conditions

$``$ $`>`$ $`0,`$ (71)
$`_y`$ $`>`$ $`0,`$ (72)
$`_{yy}`$ $`>`$ $`0,`$ (73)
$`3_y`$ $``$ $`yN_{yy},`$ (74)
where $`y=\sqrt{2g^2F}`$ and $`g`$ is the magnetic charge.

One can show that the conditions (71,72,73,74) are satisfied for the Bardeen black holes of ABG-M1  for any value of the charge to mass ratio (such that one has a black hole). A proof of this statement is given in the appendix. We have also checked that the conditions (67,68,69,70) hold for the black hole solutions found in ABG-E1 ; ABG-E2 ; ABG-E3  for different values of the charge to mass ratio. This establishes the linear stability of these solutions with respect to linear fluctuations in the domain of outer communication. An interesting question is what happens to perturbations inside the black hole. As expected from the mass inflation scenario PI  fluctuations that are falling into the black hole are likely to destroy the Cauchy horizon and to introduce a singularity inside the hole. These questions will be analyzed further in some future work where we also plan to discuss the stability of particle-like solutions ABMS .

## Acknowledgments

The authors wish to thank Eloy Ayón-Beato and Carlos Lousto for many useful discussions. This work was supported by the Mexican Council of Science and Technology, CONACyT, and the Swiss National Science Foundation.

## Appendix A Linear stability of the Bardeen black holes

The Bardeen black holes are obtained in the $`F`$ framework from the Lagrangian ABG-M1

$$(y)=\frac{3}{2sg^2}\left(\frac{y}{1+y}\right)^{5/2},$$

where $`y=\sqrt{2g^2F}`$ and $`s=|g|/2m`$ is half the monopole charge to mass ratio. Since

$`_y`$ $`=`$ $`{\displaystyle \frac{15}{4sg^2}}{\displaystyle \frac{y^{3/2}}{(1+y)^{7/2}}},`$
$`_{yy}`$ $`=`$ $`{\displaystyle \frac{15}{8sg^2}}{\displaystyle \frac{y^{1/2}}{(1+y)^{9/2}}}(34y),`$
$`f(y)y_y^1_{yy}`$ $`=`$ $`{\displaystyle \frac{34y}{2(1+y)}},`$
it is clear that $`_y>0`$ and so it is sufficient to check the condition

$$3Nf(y)>0$$
(75)
for $`0yy_h`$ where $`y_h`$ is the value of $`y`$ at the event horizon.

Since $`f(y)`$ is monotonously decreasing with $`f(0)=3/2`$ the first part of the inequality (75) is verified. Next, consider the component

$$N=g_{tt}=1\frac{1}{s}\frac{y^{1/2}}{(1+y)^{3/2}}$$

of the metric which determines the location of the event horizon. $`N`$ has a single minimum at $`y=1/2`$, so $`y_h1/2`$. Therefore, the second part of the inequality (75) is also verified for all $`0yy_h`$ and we conclude that the Bardeen black holes are stable.