Testing the nonlinear stability of Kerr-Newman black holes

# Testing the nonlinear stability of Kerr-Newman black holes

Miguel Zilhão Center for Computational Relativity and Gravitation and School of Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623, USA Departament de Física Fonamental & Institut de Ciències del Cosmos, Universitat de Barcelona, Martí i Franquès 1, E-08028 Barcelona, Spain    Vitor Cardoso CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal. Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada    Carlos Herdeiro Departamento de Física da Universidade de Aveiro and I3N, Campus de Santiago, 3810-183 Aveiro, Portugal    Luis Lehner Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada    Ulrich Sperhake Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Cambridge CB3 0WA, UK Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125, USA Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA
###### Abstract

The nonlinear stability of Kerr-Newman black holes (KNBHs) is investigated by performing numerical simulations within the full Einstein-Maxwell theory. We take as initial data a KNBH with mass , angular momentum to mass ratio and charge . Evolutions are performed to scan this parameter space within the intervals and , corresponding to an extremality parameter () ranging from to . These KNBHs are evolved, together with a small bar-mode perturbation, up to a time of order . Our results suggest that for small , the quadrupolar oscillation modes depend solely on , a universality also apparent in previous perturbative studies in the regime of small rotation. Using as a stability criterion the absence of significant relative variations in the horizon areal radius and BH spin, we find no evidence for any developing instability.

## I Introduction

According to the celebrated uniqueness theorems (see Robinsoon:2004zz (); Chrusciel:2012jk () for reviews), the Kerr-Newman (KN) solution Newman:1965my () describes the most general stationary, regular (on and outside a horizon) single black hole (BH) configuration of Einstein-Maxwell theory. The solution is a 4-parameter family, described by mass , angular momentum , electric charge and magnetic charge . The magnetic charge, however, besides being absent in standard electrodynamics, can be removed by using the electromagnetic duality of the electrovacuum Einstein-Maxwell theory Deser:1976iy (). As such, it is often neglected—as will be the case here.

Even though it is unlikely that the KNBH plays a relevant role in astrophysics Gibbons:1975kk (); Blandford:1977ds (); Barausse:2014tra (), this solution has raised considerable interest since its discovery, as an arena for theoretical investigations. In particular, it provides an ideal testing ground for studying the interplay between gravity and electrodynamics at a nonlinear level and the extent to which fundamental properties of the Kerr space-time are modified by the electromagnetic field.

As for similarities, the KN line element is of course remarkably similar to that of the Kerr solution. In particular, special properties of Kerr also apply to the more general KN spacetime. For instance, the Liouville integrability of the geodesic equations observed in Kerr, is still present in KN Carter:1968rr (). Indeed, KN possesses a hidden constant of motion, which permits the separability of test particle equations. Geometrically, this conserved quantity can be understood from the existence of an irreducible Killing tensor Carter:1977pq (). Yet another consequence of this hidden symmetry is that scalar perturbations, obtained by solving the scalar wave equation in the KN background, are separable Carter:1968ks ().

A different behaviour, on the other hand, is found when considering electromagnetic and gravitational perturbations. Electromagnetic fluctuations are decoupled from gravitational fluctuations in the Kerr geometry, and both separate in an elegant way when using the Newman Penrose formalism Teukolsky:1972my (). These properties allowed for a number of significant results to be achieved for the Kerr geometry, most notably its mode stability Whiting:1988vc (). In contrast, electromagnetic and gravitational perturbations do not decouple in the KN background and need to be studied jointly: a small gravitational fluctuation in such background induces a perturbation in the electromagnetic field which is of the same order of magnitude. The separability of the relevant equations in the KNBH background is a formidable open problem MTB (). This difficulty has prevented the analysis of various physical properties of the KNBH, most notably its mode stability and oscillation properties.

Understanding the stability of a solution to Einstein’s equations plays a central role in assessing the solution’s physical relevance. As such, considerable effort has been devoted towards establishing a proof of the stability of the Kerr solution beyond mode analysis Dafermos:2010hd (). Also, mounting (but certainly partial) evidence for stability has been furnished by a large body of numerical simulations performed over the last decades. These include binary mergers of BHs and/or neutron stars as well as rotating stars undergoing collapse. These efforts have accumulated considerable support for this solution being stable at the nonlinear level as well, at least within the time-scales and regimes probed by these simulations (we refer the reader to Lehner:2014asa () for a further discussion of this point and numerous representative references of relevant examples, and to Cardoso:2014uka () for similar efforts in a broader context).

Much less is known about the stability of the KN solution. Indeed, due to the difficulties mentioned above, progress has only been made recently in studying electromagnetic and gravitational perturbations in either the slow or extreme rotation limits. For instance, thorough analysis of the behavior of perturbations in the slow rotation regime indicate that KN is linearly stable for all values of the charge Pani:2013ija (); Pani:2013wsa (). The methods used for this purpose, however, are not able to probe fast rotating KN solutions and thus require a different strategy. In recent years, interesting perturbative approaches have been developed that exploit the particular structure arising in the near-extreme limit either directly Yang:2013uba (); Yang:2014tla () or through the Kerr/CFT correspondence and the expanded set of isometries arising in such scenarios Porfyriadis:2014fja (); Hadar:2014dpa (). These works are providing incipient evidence for linear stability in near-extremal BHs. In this work we shall explore the nonlinear stability of the KN solution using tools from numerical relativity, which allow us to probe the fast rotating limit (see also East:2013mfa () for studies in the non-charged case).

The formalism we employ here has been described in Refs. Zilhao:2012gp (); Zilhao:2013nda (), which was previously employed to study collisions of charged BHs with equal and with opposite charges in Einstein-Maxwell theory. With this formalism we are able, after choosing appropriate initial data, to analyse the behavior of perturbed KNBHs. Anticipating some of the discussions, we will show that our evolutions reveal no evidence for instabilities.

This paper is organized as follows. In Sec. II we briefly review the formalism used in Zilhao:2012gp (); Zilhao:2013nda () for evolving the Einstein-Maxwell system. Section III addresses the construction of appropriate initial data to describe a (perturbed) KNBH. Section IV describes the diagnostic tools used to monitor the evolution and decide on whether instabilities are present. The numerical results are reported in Sec. V and our conclusions and final remarks are made in Sec. VI.

## Ii Formalism

Following our previous work on collisions of charged BHs in Refs. Zilhao:2012gp (); Zilhao:2013nda (), we consider the enlarged electrovacuum Einstein-Maxwell equations

 Rμν−R2gμν =8πTμν , (1) ∇μ(Fμν+gμνΨ) =−κnνΨ , ∇μ(⋆Fμν+gμνΦ) =−κnνΦ ,

where is the Maxwell tensor and its Hodge dual, is a constant and is the 4-velocity of Eulerian observers. We recover the standard Einstein-Maxwell system when and merely introduce these fields as a means to damp and control violations of the magnetic and electric constraints during the numerical evolution Komissarov:2007wk (); Palenzuela:2008sf (). The electromagnetic stress-energy tensor takes the usual form

 Tμν=14π[FμλFνλ−14gμνFλσFλσ] . (2)

We perform a Cauchy (3+1) decomposition by introducing a 3-metric , and decompose the Maxwell tensor and its dual into the electric and magnetic 4-vectors as

 Fμν =nμEν−nνEμ+ϵμναβBαnβ , (3) ⋆Fμν =nμBν−nνBμ−ϵμναβEαnβ ,

where we use the convention , , .

## Iii Initial data

As already mentioned in the Introduction, the KN solution is defined by three (physical) parameters: mass , spin and electric charge . In Boyer-Lindquist coordinates , the metric and vector potential take the form (see, e.g., Ref. wald1984general ())

 ds2 =−(Δ−a2sin2θρ2)dt2+ρ2Δdr2BL+ρ2dθ2 (4) +(r2BL+a2)2−Δa2sin2θρ2sin2θdϕ2 −2asin2θ(r2BL+a2−Δ)ρ2dtdϕ, A =−QrBLρ2(dt−asin2θdϕ),

where

 ρ2 ≡r2BL+a2cos2θ, Δ ≡r2BL−2MrBL+a2+Q2.

In order to obtain initial data suitable for numerical evolutions using the “moving punctures” technique Campanelli:2005dd (); Baker:2005vv (), we express the solution in terms of a quasi-isotropic radial coordinate . Following Refs. Brandt:1996si (); Krivan:1998td (); Cook:2000vr () we perform the coordinate transformation

 rBL=R(1+M+√a2+Q22R)(1+M−√a2+Q22R),

and the metric then takes the form [ are spatial indices]

 ds2=(−α2+βϕβϕ)dt2+2βϕdϕdt+γ¯i¯jdx¯idx¯j, (5)

where

 γ¯i¯jdx¯idx¯j =ψ4[dR2+R2(dθ2+sin2θdϕ2) +a2hR4sin4θdϕ2], α =(R+RH)(R−RH)R√r2BL+a2(1+σsin2θ), βϕ =−aσsin2θ,βϕ=βϕ/γϕϕ, ψ4 ≡ρ2/R2,h≡(1+σ)/(ρ2R2), σ ≡(2MrBL−Q2)/ρ2.

Here, is the location of the event horizon in the quasi-isotropic coordinate . The nonzero components of the extrinsic curvature take the form

 KRϕ=ψ−2HEsin2θR2,Kθϕ=ψ−2HFsinθR, (6)

where

 HE ≡aMGρ3√r2BL+a2(1+σsin2θ), (7) HF ≡−ασa3cosθsin2θρ, G ≡(r2BL−a2)ρ2+2r2BL(r2BL+a2) −Q2MrBL(2ρ2+a2sin2θ).

The electric and magnetic field can be computed from (3). Its nonzero components are

 ER =QR(2r2BL−ρ2)(r2BL+a2)ρ6√r2BL+a2(1+σsin2θ), (8) Eθ =−2a2QαRcosθsinθρ6, BR =2aQRrBL(r2BL+a2)cosθρ6√r2BL+a2(1+σsin2θ), Bθ =aQα(2r2BL−ρ2)sinθρ6.

We finally transform to Cartesian coordinates . Our initial data then reads

 γijdxidxj=ψ4[dx2+dy2+dz2 +a2h(y2dx2−2xydxdy+x2dy2)], (9) Kij=Λ¯iiΛ¯jjK¯i¯j,Ei=Λi¯iE¯i,Bi=Λi¯iB¯i, (10)

where , in a form analogous to that of Refs. Hannam:2006zt (); Shibata:2009ad (); Shibata:2010wz ().

In order to study the stability of this solution, we follow Shibata:2009ad (); Shibata:2010wz () and introduce a small bar-mode perturbation to the 3-metric and specify the initial conditions for the 3-metric elements as

 ^γij=γij⎡⎢⎣1+Ax2−y2M2e−(R−R0)22R2H⎤⎥⎦−1, (11)

where , is the unperturbed solution given by (9), and is a tunable parameter that localizes the perturbation.

This perturbation is constraint violating111Constraint violations are an inherent consequence of the numerical modelling of spacetimes in general relativity at the level of the numerical discretization error. Following a common approach (see e.g. Shibata:2009ad (); Shibata:2010wz ()), we here add a small perturbation to the initial data in order to trigger an instability more rapidly (if one exists). This perturbation introduces an additional constraint violation at a level well below that due to the discretization but we significantly mitigate this effect by localizing the perturbation well within the horizon.; confining the fluctuation within the horizon (by choosing ) will however produce only a weak gravitational wave signal. Thus, we here choose to monitor quantities that describe the horizon deformation, as explained in the next section. We find that our results, described below, can also be used to understand the gravitational wave signal at large distances. Note also that for our choices of perturbation amplitudes , when looking at, for instance, the Hamiltonian constraint violation, we see no noticeable differences when comparing with non-perturbed cases.

## Iv Diagnostics

We analyze the result of our numerical investigations using the following quantities. (i) The (coordinate invariant) horizon “areal” radius

 rareal=√r2H+a2, (12)

where . (ii) The ratio between the polar and equatorial horizon circumferences

 Cp/Ce=r2Hπ(r2H+a2)∫π0√1+a2r2Hcos2θ dθ, (13)

which, for known and , allows one to determine . Finally, we (iii) quantify the “strength” of the bar-mode perturbation in terms of the following distortion parameters Saijo:2000qt (); Franci:2013mma ()

 η+≡Ixx−IyyIxx+Iyy,η×≡2IxyIxx+Iyy, (14)

where

 Iij=∫Hd3x√γxixj (15)

is the quadrupole moment of the apparent horizon. We also compute the radiation from the system by computing the Newman-Penrose scalar at distances far from the BH. We have found, however, that the behavior of is better suited to analyze the response of the near BH region to the perturbations—which, as discussed in the previous section, are initially concentrated in that region. This is a natural observation as the BH potential barrier essentially traps the induced perturbations in the BH’s vicinity.

To monitor the evolution, we compute the relative difference of both the areal radius of the apparent horizon and the measured BH spin to the known analytic value

 δ[f(t)]≡maxt>5M|f(t)−f0|f0, (16)

where is the analytic value. We choose to evaluate the maximum from onward to remove possible large fluctuations due to our initial perturbation (11). Finally we monitor convergence of the solution via standard numerical analysis and the behavior of the constraints to ensure truncation errors remain small throughout the simulation’s time span.

## V Numerical results

We numerically integrate the Einstein-Maxwell system using fourth-order spatial discretization with the Lean code Sperhake:2006cy (). This code is based on the Cactus Computational toolkit cactus (), the Carpet mesh refinement package Schnetter:2003rb (); carpet () and uses AHFinderDirect for tracking apparent horizons Thornburg:2003sf (); Thornburg:1995cp (). Lean uses the BSSN formulation of the Einstein equations Shibata:1995we (); Baumgarte:1998te () with the moving puncture method Campanelli:2005dd (); Baker:2005vv (). We refer the interested reader to Ref. Sperhake:2006cy () for further details on the numerical methods, and to Zilhao:2013nda () for the tests performed with the Einstein-Maxwell implementation.

We evolve the Einstein-Maxwell system of equations (1) for several different charge and spin values until and monitor both the areal radius and the polar to equatorial horizon circumferences ratio. As a practical measure, we consider a configuration to be stable if: (i) during the course of the numerical evolution, the BH areal radius and spin [the latter inferred through Eq. (13)] vary by less than a few percent—consistent with the perturbation—with respect to the analytic value and (ii), their time-dependence show an attenuating behavior. For visual guidance of the parameter space explored in this work, we display in Fig. 1 the extremality curve together with the distribution of spin and charge values of the simulations performed. The simulations performed include several configurations close to extremality plus additional ones far from this regime for comparison purposes.

In Table 1 we list the simulations performed with the corresponding physical parameters used. Note that, except for two instances (runs a0.990_q0.1_A0.0005 and a0.975_q0.2_A0.0005, both of these being cases where ), the relative variation in is always smaller than (and for most cases even smaller than ), which gives us confidence in the accuracy of our numerical evolution since these are consistent with corresponding results for the Schwarzschild case (). The larger variation observed in the two mentioned cases (and, to a lesser degree, also in the a0.992_q0.0_A0.0005 and a0.994_q0.0_A0.0005 runs) is due to a small but steady growth in observed from onward. We saw similar behaviour in other simulations accompanied by a steady increase of the Hamiltonian constraint violations with time. In all such circumstances, this behaviour was successfully cured with an increase in the numerical resolution used. We believe that this is happening in all aforementioned cases: the growth in the measured horizon area is merely telling us that more resolution is needed should we want to accurately evolve such near-extremal configurations () for longer times. The already very high resolution used in such cases effectively limits our ability to do so, however.

### v.1 Nonlinear stability of Kerr-Newman spacetimes

Figures 24 summarize our results. In Fig. 2 we plot the time evolution of the deformation parameters (14) for a “typical” case corresponding to .

The behavior of consists of a sum of damped sinusoids and decays away on timescales of order , consistent with linearized predictions for the ringdown timescale Berti:2009kk (). For neutral or static BHs, the ringing frequency and damping times of the fluctuations match well linearized calculations of quasinormal frequencies Berti:2009kk (). All our simulations display this same behavior: initial fluctuations are damped away. This is one of the main messages of our work: for the parameters we studied, the KN geometry appears to be nonlinearly stable against such perturbations on the timescales examined herein—thus indicating any possible instability should have a secular growth associated to it.

In Fig. 3, we show the Hamiltonian constraint violations for a example. Note that after a brief transient early on, the constraint violation is not significantly growing in time, and that it bears the overall pattern observed for typical numerical BH evolutions.

### v.2 Universality of oscillation modes

Our results indicate a surprising universal relation between the oscillation frequency and damping times of the fluctuations, namely that for large , spacetimes with the same behave in a similar way. This is summarized in Fig. 4 where we show the evolution of for three different values of which share the same , and have . The lines corresponding to the different cases overlap almost perfectly. For comparison, another value of with , but with is exhibited, for which the curve is slightly displaced from the previous ones. We note that if this agreement holds throughout the entire range of charge and mass, this would imply that the characteristic or quasinormal frequencies of these BHs satisfy

 ω=ω(a/√M2−Q2), (17)

which for small charge can also be written as , where we defined following Refs. Pani:2013ija (); Pani:2013wsa (). This prediction was tested against linearized calculations in the slowly-rotating regime from Pani:2013ija (); Pani:2013wsa (), where frequencies are expressed as . Translated into this notation, universality as described by Eq. (17) would imply that for both the real and imaginary components, which is to very good precision the result presented in Table I of Pani:2013ija (); Pani:2013wsa () for modes.

While such universality seems to hold only for quadrupolar modes (and again, the linearized calculations of Refs. Pani:2013ija (); Pani:2013wsa () are also consistent with universality for only), the mere existence of such property is intriguing and adds to the isospectrality found in linearized studies Pani:2013ija (); Pani:2013wsa ().

Such universality is not an artifact of horizon-deformation measures. Our results indicate that the gravitational-wave signal at large distances (in particular the components of the scalar ) shares the same characteristics.

Recently, an analytical formalism to compute the quasinormal mode spectra of (weakly) charged black holes has been introduced Mark:2014aja (). With it, the extent of this seemingly universal behavior can be scrutinized. This has confirmed such behavior for large spin values (see also Hod:2014uqa ()), but it degrades considerably at low ones markhuan ().

## Vi Conclusions

In this paper we have used the techniques developed in Zilhao:2012gp (); Zilhao:2013nda () for performing BH evolutions in Einstein-Maxwell theory to study the nonlinear stability of the Kerr-Newman BH for a variety of parameters and, in particular, for rapidly spinning BHs. On the timescales explored here (a few hundred ), we have seen no evidence for instabilities in any of the simulations performed. We are able to measure the spin of the BH with high accuracy (see Table 1), which varies only within the expected margin for numerical error, indicating that the solution is stable.

In order to trigger potential instabilities, we have considered an initial perturbation of a particular type: a bar-mode perturbation in the metric coefficients. We do not expect, however, that other types of qualitatively different initial perturbations—like Brill or Teukolsky waves (see e.g. Hilditch:2013cba (), for a recent study using this type of initial data in moving puncture gauge) will give different results; otherwise, an instability would appear to require very specific perturbations neither contained in our bar mode nor in the numerical noise of the initial data. As such, our nonlinear analysis reinforces previous linear results Pani:2013ija (); Pani:2013wsa (); Civin:2014bha () on the stability of the nonextremal KNBH. This contrasts with the instability found for extremal KNBHs Aretakis:2012ei (); Reiris:2013efa (). Thus, the latter, albeit continuously connected to nonextremal KNBHs in parameter space, seem qualitatively disconnected in terms of physical properties.

Our results have also uncovered, in the large rotation regime, a new class of universality for the quadrupolar quasinormal modes of these BHs: they depend solely on the combination , a feature which had been observed previously in the perturbative regime of slow-rotation. The significance of such results is unclear, but together with the isospectrality—observed also in the slow-rotation regime—hints at deeper relations at work also in rotating and charged geometries.

###### Acknowledgements.
We would like to thank Paolo Pani, Zachary Mark and Huan Yang for useful discussions. M.Z. is supported by NSF grants OCI-0832606, PHY-0969855, AST-1028087, and PHY-1229173. V.C. acknowledges financial support provided under the European Union’s FP7 ERC Starting Grant “The dynamics of black holes: testing the limits of Einstein’s theory” grant agreement no. DyBHo–256667. L.L. acknowledges support by NSERC through a Discovery Grant and CIFAR. U.S. acknowledges support by the FP7-PEOPLE-2011-CIG CBHEO Grant No. 293412, the STFC Grant No. ST/I002006/1, the XSEDE Grant No. PHY-090003 by the National Science Foundation, the COSMOS Shared Memory system at DAMTP, University of Cambridge, operated on behalf of the DiRAC HPC Facility and funded by BIS National E-infrastructure capital Grant Nos. ST/J005673/1 and ST/J001341/1 and STFC Grant Nos. ST/H008586/1 and ST/K00333X/1, and the Centro de Supercomputacion de Galicia (CESGA) under Grant No. ICTS-2013-249. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development Innovation. This work was supported by the NRHEP 295189 FP7-PEOPLE-2011-IRSES Grant, and by FCT-Portugal through projects PTDC/FIS/116625/2010, CERN/FP/123593/2011 and the IF program. Computations were performed on the “Baltasar Sete-Sois” cluster at IST, the “Blafis” cluster at Universidade de Aveiro, the NICS Kraken Cluster, the SDSC Trestles Cluster, Cambridge’s COSMOS, on the “venus” cluster at YITP, and CESGA’s Finis Terrae.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters