Strong cosmic censorship: taking the rough with the smooth

It has been argued that the strong cosmic censorship conjecture is violated by Reissner-Nordstr\"om-de Sitter black holes: for near-extremal black holes, generic scalar field perturbations arising from smooth initial data have finite energy at the Cauchy horizon even though they are not continuously differentiable there. In this paper, we consider the analogous problem for coupled gravitational and electromagnetic perturbations. We find that such perturbations exhibit a much worse violation of strong cosmic censorship: for a sufficiently large near-extremal black hole, perturbations arising from smooth initial data can be extended through the Cauchy horizon in an arbitrarily smooth way. This is in apparent contradiction with an old argument in favour of strong cosmic censorship. We resolve this contradiction by showing that this old argument is valid only for initial data that is not smooth. This is in agreement with the recent proposal that, to recover strong cosmic censorship, one must allow rough initial data.


Introduction
The strong cosmic censorship conjecture [1] asserts that, in some physically relevant class of initial data for Einstein's equation (e.g. smooth, complete, asymptotically flat), the maximal Cauchy development is, generically, inextendible. In other words, classical physics is predictable from the initial data. The Reissner-Nordström and Kerr solutions of the vacuum Einstein equation (with vanishing cosmological constant Λ) admit Cauchy horizons.
Consistency with the conjecture requires that such a Cauchy horizon is non-generic: it is expected that, if the initial data is perturbed, then generically the resulting perturbed spacetime will not admit a Cauchy horizon [2][3][4][5].
Making this conjecture precise is surprisingly subtle. 1 Various arguments indicate that, when the initial data is perturbed, the spacetime metric (and other fields) can be extended continuously across a Cauchy horizon [7][8][9][10]. For the Kerr solution, this has been proved recently [11]. So the "C 0 formulation" of the strong cosmic censorship conjecture (where "inextendible" means "inextendible with continuous metric") is false. However, it has also been argued that, generically, curvature invariants diverge at the Cauchy horizon, so the extended spacetime cannot be C 2 there [5]. Hence the C 2 formulation of strong cosmic censorship appears to be true. 2 This is not the end of the story because the total tidal distortion experienced by an observer crossing the Cauchy horizon can remain finite, so the divergence in curvature might not be strong enough to destroy a macroscopic observer [8]. Therefore demanding a C 2 metric appears to be too strong a requirement.
Ultimately, the question of whether or not an observer can cross the Cauchy horizon depends on the equations of motion for the matter that the observer is made of. And of course the observer will have an effect on the geometry determined by the Einstein equation. This motivates formulating the strong cosmic censorship conjecture as the statement that the maximal Cauchy development should be inextendible as a solution of the equations of motion. Since the equations of motion are second order, one might think this implies that the fields should be C 2 . However, one can still make sense of the equations of motion with lower smoothness than this by considering weak solutions. 3 Weak solutions have important physical applications e.g. they describe shocks in a compressible perfect fluid. For the vacuum Einstein equation, a weak solution must have locally square integrable Christoffel symbols in some chart. This leads to Christodoulou's formulation [13] of the strong cosmic censorship conjecture, that, generically the maximal Cauchy development is intextendible as a spacetime with locally square integrable Christoffel symbols. If this is correct then, generically, one cannot extend beyond the Cauchy horizon consistently with the classical equation of motion.
A popular toy model for studying strong cosmic censorship is a linear massless scalar field. In this case, the analogue of the Christodoulou formulation of strong cosmic censorship is that, for generic smooth initial data, at the Cauchy horizon the scalar field will not belong to the Sobolev space H 1 loc of functions that are locally square integrable with a locally square integrable gradient. 4 More informally, the energy of the scalar field will diverge at the Cauchy horizon. Here "energy" refers to the energy on a spacelike surface intersecting the Cauchy horizon, according to an observer with velocity normal to the surface. This linear version of the Christodoulou formulation of the strong cosmic censorship conjecture has been proved to be true for Reissner-Nordström [14] and Kerr [15] black holes with Λ = 0.
The instability of the Cauchy horizon arises from a blue-shifting of perturbations entering the black hole at late time. It was observed long ago that this effect is weaker for Λ > 0 because there is a competing red-shifting of late time perturbations since such perturbations can disperse by falling across the cosmological horizon. 5 This led to the claim that the C 2 version of strong cosmic censorship is violated for near-extremal Reissner-Nordström-de Sitter (RNdS) [16] or Kerr-de Sitter (Kerr-dS) [17] black holes. However, subsequent work [18] argued that this conclusion is invalid because it neglects backscattering of outgoing radiation just inside the event horizon. It was argued that, in the presence of such outgoing radiation, the Cauchy horizon instability is still strong enough to ensure that the C 2 formulation of strong cosmic censorship is respected. 6 Nevertheless, it has been conjectured that the Christodoulou formulation would be violated for near-extremal RNdS and Kerr-dS black holes [6]. As we have discussed above, this formulation seems more relevant than the C 2 formulation.
Interest in this topic has been revived by recent work of Cardoso et al [19]. This work considered linear massless scalar field perturbations of a RNdS black hole. It was found that, for a near-extremal black hole, such perturbations have finite energy at the Cauchy horizon and therefore violate the toy model of strong cosmic censorship discussed above. Going beyond the toy model, one can consider the backreaction of the scalar field on the geometry using nonlinear results of Refs. [20][21][22]. Cardoso et al argued that, at the nonlinear level, such perturbations would respect the C 2 formulation of strong cosmic censorship but, for a near-extremal black hole, the Christodoulou formulation would be violated, in agreement with the conjecture of Ref. [6].
This raises the question of whether the same worrying behaviour is exhibited in the more physical case of Kerr-dS black holes. Surprisingly, the answer appears to be negative: Ref. [23] argued that the Christodoulou formulation of strong cosmic censorship is respected by gravitational (or massless scalar field) perturbations of such black holes, even close to extremality. Thus the evidence suggests that, for Λ > 0, the Christodoulou formulation of strong cosmic censorship is respected by the vacuum Einstein equation but not by the Einstein-Maxwell-massless scalar field equations! 7 Our discussion so far has concerned only perturbations arising from smooth initial data. Very recently, Dafermos and Shlapentokh-Rothman (DSR) [25] have suggested a way of rescuing strong cosmic censorship with Λ > 0, namely to consider initial perturbations 5 For Λ < 0 one would expect the Cauchy horizon instability to be stronger than for Λ = 0 because perturbations outside a black hole decay very slowly. It has been suggested that the C 0 formulation of strong cosmic censorship might be valid for Λ < 0 [11]. 6 There is a problem with this claim which we will discuss below. 7 For massless scalar field perturbations, it has been argued that a near-extremal Kerr-Newman-dS black hole respects strong cosmic censorship provided that it rotates sufficiently rapidly [24]. The latter condition cannot be relaxed because the zero rotation limit gives RNdS, for which strong cosmic censorship is violated.
which are not smooth. As discussed above, the equations of motion can be formulated even with low differentiability. For linear massless scalar field perturbations of RNdS, DSR proved that, generically, the solution at the Cauchy horizon is less regular (in the sense of Sobolev spaces) than the initial data. Now there will be some minimum level of regularity which is acceptable, either physically or mathematically, e.g. for finiteness of energy or (in the nonlinear context) for local well-posedness of the initial value problem. The DSR result suggests a "rough" (i.e. non-smooth) formulation of the strong cosmic censorship conjecture: if one has an initial perturbation with the minimum acceptable level of regularity then, generically, the perturbation at the Cauchy horizon will not have this minimum acceptable regularity [25].
A lack of smoothness of the initial perturbation was already present, although not noticed, in the earlier work of Ref. [18]. As we will show in section 2, the argument of Ref. [18] overlooks a subtlety which implies that this argument only works for initial data that is not C 1 at the event horizon. Thus the work of Ref. [18] does not establish that the C 2 formulation of strong cosmic censorship is respected, because the initial perturbation does not belong to C 2 . Instead, as we will explain, the argument of Ref. [18] is evidence in favour of the rough version of strong cosmic censorship proposed by DSR.
In this paper, we will hammer a few more nails into the coffin of the smooth versions of strong cosmic censorship for RNdS. We will study linearized electromagnetic and gravitational perturbations of a RNdS black hole. Our results assume that the perturbation arises from smooth initial data. We will show that, near extremality, Christodoulou's formulation of strong cosmic censorship is violated by such perturbations. This is analogous to the massless scalar field results of Ref. [19]. However, in contrast with that case, our results show that, in pure Einstein-Maxwell theory, the C 2 version of strong cosmic censorship is also violated near extremality. In fact, generic perturbations arising from smooth initial data can be arbitrarily smooth at the Cauchy horizon. More precisely, if one desires that every perturbation arising from smooth initial data is C r at the Cauchy horizon then this can be achieved by taking the black hole to be close enough to extremality and large enough. Hence, in pure Einstein-Maxwell theory with Λ > 0, not only are the Christodoulou and C 2 formulations of strong cosmic censorship violated (for smooth initial data), but so is the C r formulation for any r ≥ 2! This paper is organized as follows. In section 2 we review the RNdS solution and discuss the arguments of Refs. [16,18]. We will explain the connection between strong cosmic censorship and quasinormal modes of the RNdS solution. In sections 3-6 we discuss linearized electromagnetic and gravitational perturbations of RNdS. We will study these perturbations using the Kodama-Ishisbashi (KI) formalism [26]. In section 3 we determine the condition for a linearized gravitoelectromagnetic perturbation to be extendible across the Cauchy horizon as a weak solution of the equations of motion. In section 4 we give the KI master equations and boundary conditions that we later solve analytically and numerically. We also show that vector-type and scalar-type perturbations in RNdS are isospectral, i.e. they have the same frequency spectrum. In section 5, we show that RNdS gravitoelectromagnetic quasinormal modes fall into three familes, as in the case of the quasinormal modes of a scalar field discussed in [19]. For all of them, there are regimes in the parameter space where we can derive some analytical approximations. We compare them with the exact numerical data and this proves valuable to identify and classify the quasinormal mode families. Finally, in section 6 we present our main results for the spectral gap of gravitational and electromagnetic perturbations. Section 7 contains further discussion of the implications of our results.
where R is the Ricci scalar of the metric g and F = dA is the Maxwell field strength associated to the potential 1-form A. We define the de Sitter radius L by In static coordinates (t, r, θ, φ), the Reissner-Nordström de Sitter (RNdS) solution with mass and charge parameters M and Q is with dΩ 2 2 being the line element of a unit radius S 2 (parametrized by θ and φ) and For an appropriate range of parameters the function f has 3 positive roots r − ≤ r + ≤ r c corresponding to the Cauchy horizon CH + , event horizon H + R and cosmological horizon H + C respectively. We will denote the (positive) surface gravities associated to each of these three horizons as κ − , κ + and κ c , respectively. For any non-extremal RNdS black hole it can be shown that [18] κ − > κ + .
The extremal configuration occurs when κ + and κ − vanish. This happens when Q = Q ext where When presenting many of our results and associated plots we will parametrize the RNdS solution using the dimensionless parameters Q/Q ext and y + . The causal structure of a non-extremal RNdS black hole is shown in Fig. 1. Region I is the region with r + < r < r c between the event horizon and cosmological horizon, i.e. the black hole exterior. Region II is the black hole interior, where r − < r < r + .
In region I we define the tortoise coordinate r * by  and we fix the constant of integration by imposing r * = 0 at r = (r + + r c )/2. We then define Eddington-Finkelstein coordinates in region I by u = t − r * and v = t + r * . In ingoing Eddington-Finkelstein coordinates (v, r, θ, φ) the metric takes the form This metric can be analytically extended into region II so these coordinates cover regions I and II of Fig. 1. We will also make use of Kruskal coordinates near the event horizon. These are defined in region I by and these coordinates also allow the metric to be analytically extended into region II (where U + > 0, V + > 0) as well as two further regions not shown in Fig. 1. The future event horizon H + R is the surface U + = 0. In region II, we have V + = e κ + v and we define the coordinate u in this region by Note that u → +∞ as we approach H + R in either region I or region II. In region II we define t and r * by u = t − r * and v = t + r * . The coordinate r * ranges from −∞ at the event horizons H + L and H + R to +∞ at the Cauchy horizons CH + L and CH + R (see Fig. 1). In region II, the ingoing Eddington-Finkelstein coordinates are smooth at the "left" component CH + L of the Cauchy horizon. We will be interested in the "right" component of the Cauchy horizon CH + R . To introduce coordinates regular there, we use outgoing Eddington-Finkelstein coordinates (u, r, θ, φ). The metric is (2.10) The Cauchy horizon CH + R is the surface r = r − in these coordinates. In region II, we define Kruskal coordinates near the Cauchy horizon as The Cauchy horizon CH + R is the surface V − = 0 in these coordinates. Finally, in region I, we define Kruskal coordinates at the cosmological horizon by (2.12) The future cosmological horizon H + C is the surface V c = 0.

The work of Moss et al.
Strong cosmic censorship for RNdS black holes was first studied by Moss and collaborators in a series of papers. In this section we will review the arguments of Moss et al presented in Refs. [16,27,18]. The analysis of [16] concluded that strong cosmic censorship is violated by some RNdS black holes. However, this conclusion was modified in Refs. [27,18], resulting in the revised conclusion of Ref. [18] that in fact (the C 2 version of) strong cosmic censorship is never violated by RNdS black holes. We will explain why this latter conclusion is valid only if one allows non-smooth initial data. We will consider perturbations by a scalar field Φ although the results of Moss et al apply also to the case of coupled electromagnetic and gravitational perturbations, which we will study later. One can prescribe initial data for the scalar field on the surface Σ of Fig. 1 since this is a Cauchy surface for regions I and II. Equivalently, one can prescribe (characteristic) initial data for the scalar field on the null surface H + L ∪ H − ∪ H − c . We will follow the latter approach. Given generic initial data for the scalar field, we want to know how the field behaves at the Cauchy horizon CH + R . This problem was first investigated by Mellor and Moss (MM) [16]. Their results (rederived below) indicate that the scalar field will fail to be C 1 at the Cauchy horizon if there exists a sufficiently slowly decaying quasinormal mode. More precisely, let α be the spectral gap, i.e. the distance from the real axis in frequency space to the lowest (slowest decaying) quasinormal frequency. Define MM showed that if β < 1 then the scalar field fails to be C 1 at CH + R . When gravitational backreaction is included, the blow-up of the derivatives of Φ at CH + R is expected to cause a blow up of curvature. Thus if β < 1 for all black holes then the C 2 version of strong cosmic censorship is expected to hold. However, by studying quasinormal modes, MM argued that RNdS black holes with |Q| ≈ M have β > 1, so the scalar field is C 1 at CH + R , which is evidence for a violation of the C 2 version of strong cosmic censorship.
MM modified this claim in Ref. [27]. They were motivated by earlier work on a toy model (null dust) [28] which suggested that the analysis of Ref. [16] missed an important effect arising from late-time ingoing radiation propagating along the cosmological horizon H + c . MM argued that, in the presence of such radiation, the scalar field will fail to be Subsequently, Brady, Moss and Myers (BMM) [18] argued that one must also include the effect of scattering of outgoing radiation propagating near H + R and in this case the scalar field will fail to be C 1 at CH + R if β < 1 where β = min(α, κ + , κ c )/κ − . In view of (2.4), this gives β < 1 for any non-extremal RNdS black hole and so BMM concluded that the C 2 version of strong cosmic censorship is always respected.
We will show that the arguments of Refs [27,18] are valid only for initial data which is not smooth, in fact not even C 1 , at, respectively, the future cosmological horizon H + C or future event horizon H + R . Hence this work cannot be regarded as evidence in favour of the C 2 version of strong cosmic censorship because the initial data is not in C 2 . However, we will show that these arguments can be reinterpreted as evidence in favour of the rough version of strong cosmic censorship proposed in Ref. [25]. If one insists on smooth initial data then the original conclusion of MM is still valid: it is simply the quasinormal modes which determine whether or not strong cosmic censorship (in either the C 2 or Christodoulou formulation) is violated.
We will consider solutions which can be written as superpositions of mode solutions. A mode solution has the separable form where Y m is a spherical harmonic. Substituting this into the wave equation or Klein Gordon equation (if the field is massive) one finds that the function R satisfies an equation of the form where the potential V (r) is independent of ω and vanishes exponentially fast as a function of r * as r * → ±∞ in either region I or II. We will start by considering solutions in region II. For real ω, by reformulating (2.15) as an integral equation, one can define two linearly independent solutions with the following behaviour as r * → −∞ in region II [4,16,29,25] 8 : R out,+ gives a scalar field solution Φ smooth on H + L and R in,+ gives a solution smooth on H + R (see Fig. 1). Similarly as r * → ∞ we can define two linearly independent solutions by R out,− ∼ e iωr * , R in,− ∼ e −iωr * , (2.17) 8 We use the notation of Ref. [25] although our mode functions differ from theirs by a factor of r. For us, Rout,+ ∼ e iωr * means Rout,+ = e iωr * R out,+ whereRout,+ is a real analytic function of r for r− < r < rc withRout,+(r+) = 1. and these give scalar field solutions that are smooth at CH + R and CH + L , respectively. We can now write where A and B are the transmission and reflection coefficients for fixed frequency scattering of waves propagating out from H + L andÃ,B are the transmission and reflection coefficients for scattering of waves propagating in from H + R . In region II, initial data can be specified on the characteristic hypersurface H + L ∪ H + R . We assume that the data on H + L is a wavepacket with Fourier transform Z(ω): 19) and the data on H + R is a wavepacket with Fourier transformZ(ω): It follows that the solution in region II is These are, respectively, the parts of Φ that are outgoing and ingoing near the Cauchy horizon. The outgoing part is smooth at CH + R and the ingoing part is smooth at CH + L . We are interested in how smooth the ingoing part is at CH + R where r * → ∞ and we have and hence, taking a derivative w.r.t. the Kruskal coordinate V − that is smooth at CH + R , We now want to examine whether To do this we need to determine whether or not the integral decays faster than e −κ − v as v → ∞. To determine the decay of the integral, we can deform the contour of integration into a line of constant Im(ω) in the lower half complex ω plane. How far we can deform the contour depends on the analyticity properties of the quantity F(ω), which we will now investigate, following [4,16,29].
First, to calculate B andÃ we proceed as follows. For functions f (r * ) and g(r * ) the Wronskian is (a prime denotes a derivative w.r.t. r * ) (2.26) and this is constant (in r) if f, g are solutions of (2.15). We now havẽ where the latter expression follows from evaluating the Wronskian in the denominator at r * → ∞. Similarly, The analyticity properties of the solutions of the radial equation have been determined in Refs. [4,29,25]. The result is that R in,+ (ω, r) can be analytically continued to the complex ω plane, except for simple poles at negative integer multiples of iκ + . Similarly, R out,+ has simple poles at positive integer multiples of iκ + and R out,− has simple poles at negative integer multiples of iκ − . Using (2.4), it follows that, in the lower half-plane, the first pole ofÃ is at −iκ + and the first pole of B is at −iκ − . Consider first the case in which the wavepackets on H + L and H + R are compactly supported in u and v respectively. Then Z(ω) andZ(ω) are entire functions. Using (2.4) it then follows that we can deform our contour of integration to a line of constant Im(ω) in the lower half-plane, until we hit a pole inÃ(ω) at ω = −iκ + . Now, if the wavepacket on H + R is generic thenZ(−iκ + ) = 0 and so this pole will also be a pole of F(ω) with residue proportional toZ(−iκ + ). Hence (2.29) Using (2.4), the above quantity diverges at CH + R where v → ∞. Hence, for generic compactly supported smooth initial data prescribed on H + L ∪ H + R the solution will not be C 1 at CH + R , in apparent support of strong cosmic censorship. It turns out that this argument is too quick because, in the problem of interest, we are not free to prescribe the initial data on H + R . Instead, this data is determined by the solution outside the black hole, i.e. in region I. We will now review the argument of Ref. [16] that shows that in factZ(−iκ + ) vanishes, invalidating the above argument. This analysis will reveal instead that the question of strong cosmic censorship depends on quasinormal modes of the black hole.
First we need to define the mode functions in region I. We define two linearly independent solutions R in,+ and R out,+ of equation (2.15) using exactly the same conditions (2.16) as before except that now these conditions are being applied in region I instead of region II. We define a second pair of linearly independent solutions R in,c and R out,c in region I in terms of their behaviour at the cosmological horizon r * → ∞ R out,c ∼ e iωr * , R in,c ∼ e −iωr * . (2.30) We can expand R in,+ in terms of these solutions as Here, T and R are the transmission and reflection coefficients for scattering of waves incident from H − c . Similarly, we can write whereT andR are the transmission and reflection coefficients for waves progating out of H − . In region I, initial data can be specified on the characteristic hypersurface H − ∪ H − c . We assume that the data on H − is a wavepacket with Fourier transform X(ω): and the data on H − c is a wavepacket with Fourier transformX(ω): It follows that the solution in region I is We can now evaluate this on the event horizon H + R , where r * → −∞. The first term vanishes there provided our initial outgoing wavepacket on H − vanishes on the black hole bifurcation sphere, as it must for the Fourier transform to be well-defined. This leaves where in the final step we evaluated the numerator at r * → ∞. Similarly, Consider initial data which is compactly supported on H − and H − c (w.r.t. u, v respectively), so X(ω) andX(ω) are entire functions. Recall that the analytic continuation of R in,+ has simple poles at negative integer multiples of iκ + . It follows that the analytic continuations of T andR have zeroes at these locations. Hence, for this initial data,Z(ω)(−iκ + ) = 0, as first explained by Mellor and Moss [16].
Recall that the behaviour of Φ near CH + R is determined by the analyticity properties of F(ω). From the above we have We start by considering the case in which the initial data on H + L and H − are compactly supported functions of u, and the initial data on H − c is a compactly supported function of v, so Z(ω), X(ω) andX(ω) are entire functions. In the above expression, the mode functions with poles in the lower half plane are R in,+ (at negative integer multiples of iκ + ), R out,− (at negative integer multiples of iκ − ) and R out,c (at negative integer multiples of iκ c ). However, in F(ω) the poles associated to R in,+ and R out,c will cancel out in the ratios of Wronskians. Therefore singularities of F(ω) in the lower half plane can only arise from the poles in R out,− and where W [R in,+ , R out,c ] = 0 .
This is the condition for R in,+ and R out,c to be linearly dependent, the defining condition of a quasinormal mode. The corresponding values of ω are called quasinormal frequencies.
We see that, for compactly supported initial data, F(ω) is analytic in the lower half-plane except for poles at quasinormal frequencies and at negative integer multiples of iκ − . As discussed above, the spectral gap α is defined as the infimum (smallest value) of −Im(ω) over all quasinormal modes. Deform the contour of integration in (2.24) the line Im(ω) = −α + for arbitrarily small > 0. In other words, we push the contour of integration down until just before it hits the "lowest" (i.e. slowest decaying) quasinormal mode(s). In doing this we may pick up contributions from poles at multiples of −iκ − if these lie closer to the real axis than the lowest quasinormal mode. However, the contribution from such poles to the integral of (2.24) will have v-dependence e −nκ − (for positive integer n), and the contribution to (2.24) will behave as e (1−n)κ − v = (−V − ) n−1 , which is smooth at CH + R . The non-smooth part of (2.24) arises from the integral along the new contour of integration. This integral decays as e −αv for large v. Hence the non-smooth part of (2.24) is proportional to where β is defined in (2.13). If β < 1 then the scalar field is not C 1 at the Cauchy horizon CH + R (where V − = 0). If β < 1/2 then it does not even have locally square integrable derivatives, i.e. it does not have locally finite energy. On the other had, if β ≥ r for some positive integer r then the above result is consistent with the scalar field being C r at the Cauchy horizon. So, for compactly supported initial data, the question of strong cosmic censorship reduces to identifying the most slowly decaying quasinormal modes of the black hole [16].
We now investigate what happens when we relax the condition that the initial data on H − c has compact support. For now we continue to assume compact support on H − and H − L . Solutions arising from such initial data were first considered by Mellor and Moss [27]. They argued that late time ingoing radiation propagating along H + c will lead to an additional pole in F(ω) at ω = −iκ c . Their argument goes as follows. Assume that the wavepacket on H − c is smooth at the cosmological bifurcation sphere B c (U c = V c = 0). The wavepacket must vanish there (otherwise it cannot be built as a superposition of modes as assumed above). Demanding that it does so smoothly leads to the condition Φ ∝ V c as V c → 0 on H − c (i.e. on U c = 0), which implies Φ ∝ e −κcv for large v on H − c . This implies that the Fourier transformX(ω) is analytic in the strip −iκ c < Im(ω) ≤ 0 but X(ω) generically has a simple pole at ω = −iκ c . This is the basis of the claim in Ref. [27] that F(ω) has a pole at ω = −iκ c . However, this claim is incorrect because, in (2.40), this pole inX(ω) is cancelled by a corresponding pole in W [R in,+ , R out,c ] arising from the pole in R out,c at ω = −iκ c . In other words, this pole is cancelled by a corresponding zero in the transmission coefficient T (ω). 9 We see that considering this data with non-compact support on H − c does not change our conclusions above: it is still the quasinormal modes which determine whether or not strong cosmic censorship is violated. However, in making this statement we have assumed that our initial data is smooth at the cosmological bifurcation sphere. If we allow nonsmooth data, as advocated in Ref. [25], then late-time ingoing radiation does lead to a new effect. Consider Φ ∝ e −γκcv for large v on H − c , i.e. Φ ∝ V γ c on H − c , with 0 < γ < 1. Clearly such data is not differentiable at V c = 0, but it has locally finite energy if γ > 1/2 since this is the condition for the gradient of Φ to be locally square integrable. For such data, X(ω) has a pole at ω = −iγκ c and, in the expression for F(ω), this is not cancelled by a zero of T (ω). Hence at CH where δ = γκ c /κ − . Locally finite energy at the Cauchy horizon requires δ > 1/2. If κ c < 2κ − then we can choose γ > 1/2 such that δ < 1/2. 10 In other words, for a RNdS black hole with κ c < 2κ − , ingoing wavepackets with locally finite energy on H − c give solutions whose energy is not locally finite at CH + R . However, we emphasize that such wavepackets are not smooth at the cosmological bifurcation sphere.
Next we consider relaxing the condition that the wavepacket on H + L has compact support. This was important in the argument of Ref. [18] asserting that (the C 2 version of) strong cosmic censorship is respected for any RNdS black hole. Once again, we will first 9 More generally, writing the initial data on H − c as Φ = f (Vc) = f (e −κcv ), for smooth f , taking the Fourier transform and repeatedly integrating by parts one can see thatX(ω) can have poles at negative integer multiples of iκc. These are all cancelled by corresponding zeros in T (ω). 10 More mathematically, the initial data is such that the solution initially belongs to H 1 loc but the solution at CH + R does not belong to H 1 loc .
consider the case of smooth initial data. On H + L ∪ H − (i.e. the line V + = 0) we will assume that the data vanishes at the bifurcation sphere B, i.e. at U + = 0, (which is required for the Fourier transforms Z(ω) and X(ω) to exist as functions) but has non-vanishing derivative there, so for small U + we have, for some constant k (2.43) In region II this gives Φ| H + It follows that Z(ω) and X(ω) both have poles at ω = −iκ + , with equal and opposite residues. Hence it appears that F(ω) will have a pole at ω = −iκ + [18]. But we will now show that the poles in Z and X cancel out in the expression for F(ω). First, note that if there is a pole at ω = −iκ + in F then the residue of this pole is proportional to Recall that R in,+ has a simple pole at ω = −iκ + , i.e.
where g(ω, r) is analytic at ω = −iκ + . The solution R in,+ is obtained by converting (2.15) to an integral equation, and solving by iteration [4,29]. Indeed this is how one sees that it has a simple pole at ω = −iκ + . One can also see from this procedure that the residue h(r) can be expressed as a series in e κ + r * , and is proportional to e κ + r * as r * → −∞. Now, h(r) must satisfy (2.15) with ω = −iκ + . But the solution of (2.15) with behaviour e κr * = e iωr * as r * → −∞ is R out,+ (−iκ + , r). Hence we have 11 for some constant of proportionality c. It turns out that c has opposite signs in regions I and II because of the way we defined R out,+ . To see this, note that e −iωt R in,+ is smooth at H + R hence e −κ + t h(r) should be smooth at H + R . But in region I near H + R we have e −κ + t R out,+ (−iκ + , r) ∼ e −κ + u = −U + whereas in region II we have e −κ + t R out,+ (−iκ + , r) ∼ e −κ + u = +U + . Hence smoothness implies that the constant c has equal magnitude but opposite sign in regions I and II. It follows that, since the numerator is evaluated in region II and the denominator in region I, we have and so the residue (2.44) vanishes. Hence F(ω) does not have a pole at ω = −iκ + . Similarly, it does not have a pole at any negative integer multiple of iκ + . 12 Hence, once 11 At the special values ω = −inκ+ (n = 1, 2, 3, . . .), Rout,+ gives mode solutions that can be smoothly extended through H + R , proportional to U n + near H + R , and the second linearly independent solution of (2.15) gives non-smooth mode solutions involving log U+. 12 One can argue as in footnote 9 that X(ω) and Z(ω) can have poles ω = −inκ+ for n = 1, 2, 3, . . ..
Their residues are related by a factor of (−1) n . This is cancelled by a corresponding factor of (−1) n relating the constant c in regions I and II. The residue in F(ω) then vanishes exactly as for the n = 1 case.
again, we find that for smooth initial data, relaxing the condition of compact support does not lead to anything new, in contrast to the claim of Ref. [18]. The reason that the argument of Ref. [18] fails is that the poles in Z and X at ω = −iκ + cancelled out in F(ω). This cancellation arose because we assumed that the first derivative of Φ was continuous at B, i.e. that the inital data is C 1 there. In order to avoid such a cancellation we have to consider initial data that is not C 1 , i.e. we have to consider rough initial data, as proposed in Ref. [25]. For example, consider initial data which vanishes on H − and H − c , i.e. X(ω) =X(ω) = 0. It follows that the resulting solution will vanish throughout region I. On H + L we take initial data Φ| H + L ∝ U γ + as U + → 0+, where 0 < γ ≤ 1 and hence Clearly our initial data is continuous, but not C 1 , at U + = 0. The resulting solution will fail to be C 1 at U + = 0, i.e. along the event horizon H + R . In terms of u, our data behaves as e −γκ + u as u → ∞ on H + L so Z(ω) has a pole at ω = −iγκ + and hence, even for γ = 1, F(ω) has a pole at the same location. It then follows that at CH + R we have , we see that the solution at CH + R is less smooth than the initial data. In particular, the condition for the initial data to have locally square integrable first derivatives (i.e. finite energy) is γ > 1/2 whereas the condition for the solution at CH + R to have locally square integrable first derivatives is δ > 1/2. For any non-extremal RNdS black hole, we can choose our initial data so that γ > 1/2 but δ < 1/2. Hence one can find an initial wavepacket with finite energy that has infinite energy at the Cauchy horizon. So if we allow such rough initial data then the Christodoulou version of strong cosmic censorship is respected, as argued in Ref. [25].
Once one is prepared to contemplate non-smooth initial data, there is no reason to work with wavepackets to show that this version of strong cosmic censorship is respected. One can work just as well with an outgoing mode solution in region II with complex frequency ω = ω 1 − iγκ + (as was done in Ref. [25] for ingoing mode solutions in region I). In region II, R out,+ can be analytically continued to complex ω, as long as ω is not a positive integer multiple of iκ + . These mode solutions behave as e −iωu near H + R . Now hence such modes vanish on H + R (i.e. U + = 0) if γ > 0. We extend the mode into region I simply by taking it to vanish in region I, i.e. we take vanishing initial data on H − and H − c . At the Cauchy horizon the reflected part of the mode behaves as with δ given by (2.50). As before, (2.4) implies that we can choose γ > 1/2 such that δ < 1/2. The initial data then has locally finite energy but the energy diverges at the Cauchy horizon. 13 In summary, we have seen that the argument of Ref. [18] does not support the strong cosmic censorship conjecture for smooth initial data. However, a modification of this argument can be viewed as supporting the strong cosmic censorship conjecture for rough initial data formulated in Ref. [25]: initial data with locally finite energy generically gives a solution whose energy is not locally finite at the Cauchy horizon.

Recent work on strong cosmic censorship with Λ > 0
For smooth initial data, we have explained why the conclusion of Ref. [16] remains valid and so whether or not strong cosmic censorship is respected can decided by looking at quasinormal modes. However, one deficiency of the above analysis is the assumption that the initial data vanishes at the bifurcation spheres B and B c . This assumption is required so that the Fourier transforms Z(ω), X(ω) andX(ω) are functions, rather than distributions. This assumption has been eliminated by more recent work in the mathematics literature [30], which proves that, for any smooth initial data, if β > 1 then the scalar field is C 1 at the Cauchy horizon and if β > 1/2 then the scalar field has finite local energy at CH + R . The recent numerical study of Ref. [19] showed that massless scalar field perturbations of RNdS black holes always have β < 1 so generic scalar field perturbations are not C 1 at the Cauchy horizon, which supports the C 2 formulation of strong cosmic censorship for the Einstein-Maxwell-massless scalar field theory. However, it was also found that nearextremal RNdS black holes have β > 1/2 and so, for smooth initial data, the Christodoulou version of strong cosmic censorship is violated in this theory.
Surprisingly, this conclusion does not hold for Kerr-dS black holes. Indeed, Ref. [23] showed that Kerr-dS black holes always have β < 1/2 and so, for smooth initial data, the Christodoulou version of strong cosmic censorship is respected for such black holes in Einstein gravity coupled to a massless scalar field. In fact, it was shown that the same result holds for linearized gravitational perturbations so it was argued that the Christodoulou version of strong cosmic censorship, with smooth initial data, is satisfied by the vacuum Einstein equations.
Finally, we should mention the work of Ref. [22]. This studies spherically symmetric perturbations of RNdS in the nonlinear Einstein-Maxwell-scalar field system. For this system, it is proved that the smoothness at CH + R is determined by how fast perturbations decay at late time along the event horizon H + R . Since linear theory should be reliable for determining the latter, this work provides justification for believing that nonlinear effects will not invalidate the conclusions of a linear analysis of the behaviour near CH + R .

Bound for weak solutions of linearized Einstein-Maxwell equations
As discussed previously, the spectral gap α is defined as the infimum (smallest value) of −Im(ω) over all quasinormal frequencies ω. Defining β = α/κ − as in (2.13), we showed above that if β < 1/2 then generic scalar field perturbations arising from smooth initial data do not have locally square integrable derivatives (i.e. locally finite energy) at the Cauchy horizon. What about gravitoelectromagnetic modes? What condition yields a linearized gravitoelectromagnetic perturbation that constitutes a weak solution of the equations of motion at the Cauchy horizon? Is the critical value still β = 1/2? In this section we will show that the answer to the latter question is positive. The analysis is rather technical so the reader may wish to skip to the summary in subsection 3.3.
Coupled linear gravitational and electromagnetic perturbations of RNdS can be studied using the Kodama-Ishisbashi (KI) formalism [26]. This formalism divides linearized gravitoelectromagnetic perturbations into perturbations arising from vector spherical harmonics and those arising from scalar spherical harmonics (there are no tensor spherical harmonics in 4d). We will consider the vector sector first (subsection 3.1) and then the scalar sector (subsection 3.2). The main conclusions are summarized in subsection 3.3.

Vector-type gravitoelectromagnetic perturbations of RNdS
Vector perturbations of the background (2.2) are described by [26] where f a , H T and A are functions of {x a } = {t, r}. Additionally, D j is the covariant derivative with respect to the unit S 2 metric γ ij and V i is a vector spherical harmonic, i.e. a regular solution of Here, 2 ≡ γ ij D i D j and regularity requires that the eigenvalues k 2 V are quantized as 3) The case v = 1 (k 2 V = 1) is special case since in this case V i is a Killing vector on the S 2 and thus D (i V j) = 0. Consequently, from (3.1) it follows that the metric components δg ij on S 2 are not perturbed.
For v > 1, all the information about the perturbations can be encoded in two gauge invariant variables Ω and A. The latter was introduced in (3.1) while the former is defined in terms of f a , H T via where ab denotes the anti-symmetric tensor on the 2-dimensional orbit spacetime. These two gauge invariant variables obey a coupled system of two master equations [26] 14 Once we have solved (3.5), we can reconstruct the original metric perturbations (3.1a) using the map [26] This determines δg µν up to a gauge transformation (infinitesimal diffeomorphism) corresponding to H T . A convenient choice of gauge is H T = 0. Note that the Maxwell perturbation is gauge invariant in the vector sector [26]. We will be interested in quasinormal modes, for which we have where ≡ v and the frequency ω is determined in terms of v and a radial "overtone" number n = 0, 1, 2, . . .. The quantized spectrum of frequencies is determined requiring that the perturbations are ingoing at the future event horizon H + R and outgoing at the future cosmological horizon H + c (see Fig. 1). In general, quasinormal frequencies are complex, ω = ω R + iω I , with ω I < 0 so that quasinormal modes decay exponentially with time outside the black hole.
For the regularity analysis at H + R it is convenient to work in ingoing coordinates (v, r, θ, φ) since they are regular both in regions I and II of Fig. 1. Then, a quasinormal mode is an analytic function of these coordinates in region I and can be analytically continued into region II. In these ingoing coordinates, a quasinormal mode has time dependence e −iωv , and thus it diverges as v → −∞, i.e. along the red line on Fig. 1. We will determine the frequency spectrum of vector quasinormal modes in Section 4.
As reviewed above, the behaviour at the Cauchy horizon CH + R of a generic perturbation arising from smooth initial data is determined by the lowest quasinormal mode [16]. Therefore we need to determine the smoothness at CH + R of the metric and Maxwell perturbations of our quasinormal modes. To do this, it is convenient to use outgoing coordinates in the black hole interior. Converting (3.7) to these outgoing coordinates in region II yields for some functionsΩ ω andÃ ω . A Frobenius analysis of (3.5) about the right Cauchy horizon CH + R , dictates that there is a pair {Ω (1) , Ω (2) } of linearly independent solutions for 14 In our conventions the parameter κ of [26] is equal to √ 2.
Ω and another pair {A (1) , A (2) } for A. These two pairs of linearly independent solutions behave as , where Ω (1,2) and A (1,2) denote non-vanishing smooth functions at r = r − . The solutions labelled (1) are outgoing at CH + R . These are smooth at CH + R . The solutions labelled (2) are ingoing at CH + R . These are not smooth at CH + R . Our quasinormal mode will be a superposition of the ingoing and outgoing solutions at CH + R . Given the behaviours (3.9) for the master variables, what is the corresponding behaviour of the metric and Maxwell perturbations at the Cauchy horizon? Again, we work in outgoing coordinates {x µ } = {u, r, θ, φ} and write the metric perturbation in these coordinates as δg µν . The KI formalism maintains covariance w.r.t. diffeomorphisms on S 2 and on the transverse 2d orbit space. Hence δg µν takes the same form as in (3.1a) with f a replaced by the quantityf a obtained from f a by the 2d coordinate transformation (t, r) → (u, r), andH T = H T . Choosing the gaugeH T = 0, we find that the two linearly independent solutions forf a have the following behaviour near Cauchy horizoñ . The behaviour at the Cauchy horizon of the Maxwell perturbation δF µν follow straightforwardly from (3.1b) and (3.9b).
Note that the outgoing solutionsf a are not. This holds in the gaugeH T = 0. We will now determine how much smoother we can make the solution using a gauge transformation.
In the vector sector, an infinitesimal gauge vector ξ has a harmonic decomposition Under such gauge transformation the metric perturbation transforms according to 12) and the Maxwell perturbation δF is invariant: see (3.1b) and recall that A is, by construction, a gauge invariant variable. We now assume and we want to choose the coefficients L (k) to make (3.10) as smooth as possible at r = r − . We find that L (0) can be chosen to setf (2;0) r = 0 in (3.10). We can then choose L (1) to setf (2;1) r = 0. But this choice then dictates that f (2) u and H (2) T behave as (r − r − ) iω/κ − because the gauge parameters L (k) with k ≥ 2 do not appear at this order. Altogether, we can find a gauge where the two linearly independent gravitoelectromagnetic solutions at the Cauchy horizon have the leading behaviour where α a = {0, 1} for a = {u, r}, respectively and f a , H T and δ F ai , δ F ij are functions that are smooth at r = r − (recall that δ F ab components are not excited in the vector sector; see (3.1)). At CH + R , our gravitoelectromagnetic quasinormal mode is some linear combination of the smooth outgoing solution (1) and the non-smooth ingoing solution (2). There is no reason for the coefficients in this linear combination to vanish. Therefore, the regularity of the quasinormal mode is determined by the ingoing solutions.
For the vacuum Einstein equation, the regularity of the metric required for a weak solution is that the Christoffel symbols should be square integrable in some chart [13]. By linearizing this condition, or by considering second order perturbation theory [23], the corresponding condition for a linearized metric perturbation to constitute a weak solution is that, in some gauge, the perturbation, and its first derivatives, should be locally square integrable, i.e. the perturbation should belong to the Sobolev space H 1 loc . In Einstein-Maxwell theory, the corresponding statement is that, in some gauge, the metric perturbation should belong to H 1 loc and the Maxwell field strength perturbation should be locally square integrable (i.e. belong to L 2 loc ). From (3.14) we see that we can reach a gauge for which the least smooth components of the metric perturbation behave as δ g (2) Similarly, (3.14) shows that the least smooth components of the Maxwell field strength perturbation behave as δ F (2) ∼ (r − r − ) p−1 (again with p = iω/κ − ). Once again this is locally square integrable if, and only if, 2(γ − 1) > −1 (again with γ = Re(p)). Hence, the condition for a vector-type gravitoelectromagnetic quasinormal mode to constitute a weak solution at the Cauchy horizon is γ > 1/2, i.e.
The above analysis shows that this condition is sufficient for the mode to constitute a weak solution at the Cauchy horizon. We believe it is also a necessary condition, and this can probably be proved along similiar lines to the argument in Ref. [23], exploiting gauge invariance of the KI variables. However, since we are mainly interested in violation of strong cosmic censorship, we will not perform such an analysis here. The above analysis was for the case v > 1. For the special case v = 1, the field H T is not defined since D (i V j) = 0. It follows that the two quantities defined by the RHS of (3.4) are no longer gauge invariant (and thus, neither is Ω). There is a single gauge invariant quantity (denoted by F 1 = ab rD a (f b /r) in (4.8) of [26]) and the map that reconstructs δg µν and δF µν from the gauge invariant quantity is (necessarily) different from the one described above for the > 1 case: in the end of the day A is the only dynamical field although it still obeys the wave equation (3.5b) (with k 2 V = 1) [26]. We have done this analysis and gravitoelectromagnetic field reconstruction 15 and we find that the condition for a v = 1 vector-type gravitoelectromagnetic quasinormal mode to constitute a weak solution at the Cauchy horizon is still given by (3.15).

Scalar-type gravitoelectromagnetic perturbations of RNdS
Scalar perturbations of the background (2.2) take the form [26] δg with f ab , f a , H T , H L , E and E b being functions of {x a } = {t, r} and ab is the anti-symmetric unit tensor. Moreover, E 0 = Q/r 2 was introduced in (2.3) and we have defined The scalar spherical harmonics S, and the associated scalar-type vector harmonic S i and traceless scalar-type tensor harmonic S ij are defined by (note that S i i = 0) The eigenvalues are quantized as Harmonics with s = 0 are non dynamical -they correspond to variations of the black hole parameters M, Q. Harmonics with s = 1 are special because S ij vanishes for these harmonics. For now we assume s > 1 and comment on the case s = 1 at the end of this section. Gauge invariant variables for the scalar perturbations are E, E a − already introduced in (3.16) − and, for s > 1, F and F ab defined as [26] The Bianchi identity requires that F b a is traceless, The reader can find the full details in the discussions (4.8)-(4.15) and (4.31)-(4.33) of [26].
The equations of motion imply that the gauge invariant quantities E and E a can be expressed in terms of a single KI master variable A as On the other hand, introducing where here and henceforward, we assume that all perturbed quantities Q(t, r) have the Fourier decomposition Q(t, r) = e −iωt Q(r) with ω being the associated frequency. The KI master variables Φ and A obey the following coupled system of equations [26] f where f is defined in (2.3). The potential V S and source term S Φ (Φ, A) are lengthy expressions given in equations (5.42)-(5.44) of [26]. The auxiliary quantity P Z is given in (C.8) of [26]. Given a solution of the above equations we will need to reconstruct the metric and Maxwell field perturbations in terms of the master variables Φ and A. For that, we first write the variables X, Y and Z in terms of Φ and A and their derivatives as [26] where the coefficients P X0 , P X1 , P XA , P Y 0 , P Y 1 , P Y A and P Z are functions of r that can be found in equations (C.4)-(C.10) and (C.11)-(C.16) of [26]. It follows from the equations of motion, including the Bianchi identity (3.21), that f ab and H L can be written as a function of X, Y, Z (i.e. of Φ, A, their radial first derivative and ω) and of f a , H T and their radial derivatives. To simplify our task (and without prejudice since we will consider gauge transformations later) we can fix the gauge as Then, the metric functions f ab and H L depend only on X, Y, Z. That is to say, via (3.26) and the Bianchi identity (3.21) they can be written solely in terms of the master variables Φ, A and their radial derivative as (3.28) We will find the frequency spectrum of scalar quasinormal modes in Section 4. But first we must discuss the behaviour of the scalar-type perturbations at the Cauchy horizon CH + R . The discussion of regularity at this null hypersurface proceeds very similarly to the vector sector case. Namely, the master variables Ω and A for the scalar quasinormal modes also admit the Fourier decomposition (3.7) and, when analytically continued into region II and converted to outgoing coordinates, these master variables also behave as (3.8). Moreover, a Frobenius analysis of (3.25) around CH + R , dictates that there is a pair {Ω (1) , Ω (2) } of linearly independent solutions for Ω and another pair {A (1) , A (2) } for A. These two pairs of linearly independent solutions still behave as described in (3.9) (with the identification ≡ s ).
Given these behaviours for the KI scalar master variables Ω and A, we can now find the behaviour of the metric and Maxwell perturbations for the outgoing and ingoing modes near CH + R . Just as we did for vector-type perturbations, in region II we transform to outgoing coordinates {x µ } = {u, r, θ, φ} in which the metric perturbation δg ab takes the same form as in (3.16), with f a replaced by the quantityf a obtained from f a via the coordinate transformation from (t, r) to (u, r), f ab is similarly replaced byf ab , butH L = H L and H T = H T are unchanged. Similarly, the Maxwell perturbation δF µν is written in terms of E a andẼ = E.
Choosing the gauge (3.27) (which translates intof a = 0,H T = 0), we find that the outgoing (smooth) and ingoing (non-smooth) solutions forf ab ,H L ,Ẽ andẼ a have, respectively, the following expansions about the Cauchy horizon: rr ). We will now show that we can make the outgoing solution smooth, and the ingoing solutions smoother at the Cauchy horizon with a gauge transformation. In the scalar sector, an infinitesimal gauge vector ξ has the harmonic decomposition Under such gauge transformation the metric and Maxwell perturbations transform according to We assume the following expansions for the functions appearing in the gauge transformation and we now try to choose the constants {N a , L (k) } to eliminate as much as we can the leading terms in (3.29) that are responsible for the lack of smoothness at the Cauchy horizon. Consider first the ingoing solution (labelled by superscript (2) ). We find that a choice of P and also to eliminate the leading term, proportional to (r a becomes non-zero as a result of the gauge transformation; the term (r − r − ) iω/κ − −1 in f (2) r has a contribution due to P (0) r and another due to L (0) ). We can now choose P and to eliminate the term proportional to (r r . But this choice then dictates that the leading term of f (2) uu , f (2) u and H (2) L,T is (r − r − ) iω/κ − because these terms in these quantities do not depend on the higher order gauge parameters and we have no more gauge freedom to avoid such powers.
Consider now the outgoing solution (labelled by superscript (1) ) in (3.29). With a choice of gauge parameters {N and also to eliminate all the terms (r − r − ) 0 log(r − r − ) that typically appear in the fields f L,T and δ F T (these fields become non-zero as a result of the gauge transformation). But with this choice it follows that the leading term of f (1) uu and f (1) u is (r − r − ) 0 since these terms do not depend on the higher order gauge parameters, i.e. we have no further gauge freedom to eliminate such terms. After these gauge transformations, the electromagnetic fields δ F (1) ab , δ F (1) ai also behave as (r − r − ) 0 .
Altogether, our analysis shows that we can find a gauge where the two linearly independent gravitoelectromagnetic solutions at the Cauchy horizon have the leading behaviour where α a = {0, 1} for a = {u, r} (respectively) and α ab = {0, 1, 1} for ab = {uu, ur, rr} (respectively), and f a , f ab , H L,T and δ F ab , δ F ai are smooth functions that depend on ω and s (recall that δ F ij is not excited in the scalar sector; see (3.16)). Note that the outgoing solution is manifestly smooth at the Cauchy horizon. As explained above, for a weak solution we need the metric perturbation and its first derivative to be locally square integrable, and the Maxwell field strength perturbation to be locally square integrable. Using the above results, we can repeat the argument we used for vector-type perturbations to see that the condition for a scalar-type quasinormal mode to be extendible as a weak solution across the Cauchy horizon is exactly the same condition (3.15) that we obtained for vector-type perturbations.
Finally, in this section we have so far assumed s > 1. Harmonics with s = 1 are special because S ij vanishes for these harmonics; as a consequence, the field H T is not defined. It follows that, for s = 1, the fields F and F ab defined in (3.20) are no longer gauge invariant [26]. Additionally, the Bianchi identity no longer implies (3.21) and it turns out that only the electromagnetic field is dynamical [26]. For our purposes, a pragmatic way to deal with this s = 1 case, as suggested in [26], is to impose (3.21) as a gauge condition and then fix a residual gauge freedoom at our convenience. 16 We can then reconstruct the gravitoelectromagnetic fields δg µν and δF µν in this particular gauge following steps similar to those described above for the s > 1. Finally, we add again gauge transformations to make our solutions smoother. In the end of the day, we find that the condition for a s = 1 scalar-type gravitoelectromagnetic quasinormal mode to constitute a weak solution at the Cauchy horizon is still given by (3.15).

Conclusions
We have shown that the condition for a linearized gravitoelectromagnetic mode solution to be extendible as a weak solution across the Cauchy horizon is (3.15). We define β in terms of the spectral gap α as in (2.13). If β < 1/2 then there exists a quasinormal mode which violates (3.15). One can add an arbitrary multiple of this quasinormal mode to any other linear perturbation. Hence if β < 1/2 then a generic linear perturbation cannot be extended as a weak solution across the Cauchy horizon. So if β < 1/2 then the Christodoulou formulation of strong cosmic censorship is respected.
Conversely, if β > 1/2 then all quasinormal modes respect (3.15). Since the behaviour at the Cauchy horizon is determined by the slowest decaying quasinormal mode, in this case, any linearized gravitoelectromagnetic perturbation arising from smooth initial data can be extended across CH + R as a weak solution of the equation of motion, so the Christodoulou version of strong cosmic censorship is violated for smooth initial data.
Finally, we can consider extendibility in C r . By this we mean that there exists a gauge so that, at CH + R , the metric is C r and the Maxwell field strength is C r−1 (so the Maxwell potential is C r in some gauge). It is easy to see from the above analysis that a quasinormal mode is extendible in C r across CH + R if −Im(ω)/κ − ≥ r. Thus, in Einstein-Maxwell theory, the C r version of strong cosmic censorship is respected if β < r and violated if β > r.

Computing the gravitoelectromagnetic quasinormal modes
In this section, we first discuss (subsection 4.1) the Kodama-Ishisbashi (KI) master equations [26] and boundary conditions of the quasinormal mode problem that we later solve analytically and numerically. We will also prove that vector-type and scalar-type modes of RNdS have the same frequency spectrum, i.e. they are isospectral (subsection 4.2).

Vector-type modes
The vector equations (3.5) describe a pair of coupled ODEs for the gauge invariant variables Ω and A. They can be rewritten as a pair of two decoupled ODEs for a pair of master variables Φ ± . These are linear combinations of the original gauge invariant variables, namely Φ ± = a ± r −1 Ω + b ± A (4.1) where a ± and b ± are functions of M, Q, given in equations (4.35)-(4.36) of [26]. Under (4.1), (3.5) tranform into the KI vector master equations where the potentials are given by When Q = 0, Φ − and Φ + are simply proportional to Ω and A, respectively. Thus, in the neutral limit, Φ − and Φ + represent, respectively, the gravitational and electromagnetic modes of the Schwarzschild black hole. Note that Φ + modes have V = 1, 2, 3 . . . whereas Φ − modes have V = 2, 3, 4, . . .. Vector quasinormal modes are solutions of (4.2) that obey ingoing boundary conditions at the black hole horizon and outgoing boundary conditions at the cosmological horizon. More concretely, at the black hole horizon r = r + a Frobenius analysis yields the expansion where Φ is either Φ + or Φ − . Regularity at the event horizon, which follows from demanding a smooth expansion in ingoing coordinates (v, r, θ, φ) around H + R , requires that we discard the solution with the positive sign. Similarly, a Frobenius expansion at the cosmological horizon r = r c yields the two possible solutions and imposing outgoing boundary conditions at the cosmological horizon H c R requires that we discard the irregular solution with plus sign. We are thus lead to introduce the field redefinition: whereΦ ± (r) is a smooth function at r = r + and at r = r c . This effectively imposes the desired boundary conditions since our numerical method can only search for smooth functionsΦ ± (r). Inserting (4.6) into (4.2) we get a pair of decoupled ODEs forΦ ± . Each of these ODEs is quadratic in the frequency ω. That is to say, for each we have to solve a quadratic eigenvalue problem to find the eigenvalue ω and the associated eigenfunctionΦ − (or ω and Φ + ). The boundary conditions forΦ ± (r) follow directly from doing a Taylor expansion of the master equation about the black hole and cosmological horizons. These reveals that at both horizons we have a Robin boundary condition, i.e. of the type (4.7) where Q +,1 , Q +,0 , Q c,1 and Q c,0 are known functions which are at most second order polynomials in ω.
It is also convenient to use a radial coordinate whose range is independent of the black hole parameters. We define such that y ∈ [0, 1] with y = 0 (y = 1) corresponding to the event (cosmological) horizon. 17 The resulting equation forΦ − (orΦ + ) can now be solved using a pseudospectral grid discretization (with the methods reviewed in [31]) as a standard quadratic eigenvalue problem or employing a Newton-Raphson algorithm. In the former method one writes the equation as a quadratic eigenvalue problem for the frequency ω, which is then solved using Mathematica's built-in routine Eigensystem. More details of this method and the discretization scheme can be found e.g. in [32]. The second method is based on an application of the Newton-Raphson root-finding algorithm, and is detailed in [33,31]. The advantage of the first method is that it gives all modes simultaneously. The second method computes a single mode at a time, and only when a seed is known that is sufficiently close to the true answer. However, this method is much quicker as both the size of the grid and numerical precision increases, and can be used to push the numerics to extreme regions of the parameter space.

Scalar-type modes
The pair of coupled ODEs (3.25) for the scalar gauge invariant variables Φ and A can be rewritten as a pair of two decoupled ODEs for a pair of scalar master variables Φ ± . The latter are given by the linear combinations where a ± and b ± are functions of M, Q, given in equations (5.57)-(5.58) of [26]. Inserting (4.9) into (3.25) yields the KI scalar master equations where the potentials V s± are given by equations (5.60)-(5.63) of [26]. When Q = 0, Φ − is proportional to Φ and Φ + is proportional to A. Hence, in the neutral limit, Φ − and Φ + represent, respectively, the gravitational and electromagnetic scalar modes of the Schwarzschild black hole. Note that Φ + modes have S = 1, 2, 3 . . . whereas Φ − modes have S = 2, 3, 4, . . .. Scalar quasinormal modes are solutions of (4.10) that obey ingoing boundary conditions at the black hole horizon and outgoing boundary conditions at the cosmological horizon. The analysis of these boundary conditions is very much similar to the one done for the KI vector sector. In fact equations (4.4) to (4.7) and the subsequent discussion apply without change to the scalar sector of perturbations.

Isospectrality
As discussed in previous sections, gravitoelectromagnetic perturbations of RNdS black holes come in two classes: vector-type and scalar-type. Although they obey two seemingly distinct equations of motion, it turns out they have the same quasinormal mode spectra. For this reason, the spectrum of quasinormal modes of RNdS black holes is said to be isospectral. This is a classical result in the context of asymptotically flat RN black holes, which was first uncovered by Chandrasekhar in [34]. It turns out the same result applies in the context of RNdS black holes, but with more involved algebra.
Just as in [34], we start by noting that the scalar potential V s± (r) − introduced in (4.10) − can be written in the following compact manner where , (4.12c) and f (r) is given in (2.3). Rather remarkably, the vector potential (4.3) takes a similar form with the same quantities defined in (4.12).

Because of this simple relation between the scalar and vector potentials, one can relate solutions of the vector equation to solutions of the scalar equation (and vice versa), via the map
where, momentarily, we added the subscripts s and v to distinguish between scalar and vector perturbations. Maps between solutions might not take physical solutions into physical solutions since one has to check that the maps preserve the relevant boundary conditions. This is the case (i.e. the map (4.14) preserves the boundary conditions) for asymptotically flat RN black holes and RNdS black holes, but it is not the case for RN black holes with anti de-Sitter boundary conditions [35]. For this reason isospectrality occurs in the former two cases, but not in the latter. Note that the differential map (4.14) alone is not enough to guarantee that the critical β bound (3.15) found for vector-type modes also holds for scalar-type perturbations, since the two types of metric perturbations are orthogonal to each other. For this reason, in Section 3 we had to do the analysis that finds the bound (3.15) for the vector and scalar-type of perturbations independently. We concluded that it turns out that (3.15) holds for both sectors.

Classifying the families of quasinormal modes and analytical results
Cardoso et al found that massless scalar field quasinormal modes of RNdS can be classified into three families [19]. We find that the same is true for gravitoelectromagnetic quasinormal modes. The three families are 1) "photon sphere" modes, 2) "de Sitter" modes and 3) "near-extremal" modes. The "photon sphere" modes are identified in the geometric optics limit, 1, and are related to the properties of the unstable circular photon orbits in the equatorial plane of the black hole background (subsection 5.1).
The de Sitter modes reduce, when M and Q vanish, to the gravitational and electromagnetic quasinormal modes of de Sitter spacetime (subsection 5.2). Finally, the "nearextremal" modes have their wavefunction peaked near the horizon and an approximate expression for these modes (strictly valid in the extremal limit) can be obtained analysing the perturbations in the near-horizon geometry of a near-extremal RNdS black hole (subsection 5.3).
In the previous Section 4.2 we found that the spectra of vector-type and scalar-type of quasinormal modes is isospectral. It follows that for each family of modes we just have to consider two sectors (not four) of perturbations corresponding to perturbations for each of the gauge invariant variables Φ − and Φ + . As a test of our numerical code, we did several checks (i.e. for different black holes) that the frequency eingenvalues of the vector-type equation of motion are indeed the same as those that solve the scalar-type equation of motion.
In this section we will obtain approximate analytical expressions for the three families of modes (that are valid at least in a certain region of the RNdS parameter space). Then we compare these analytical results with the exact data that results from our numerical search of the frequency spectra in the full RNdS parameter space 0 ≤ y + ≤ 1 and 0 < Q/Q ext ≤ 1.

Photon sphere family of modes and its geometric optics limit
In this subsection we will find an analytical expression for the photon sphere quasinormal modes in the geometric optics limit, i.e. in the WKB limit → ∞. We find that this analytical expression gives an imaginary part of the frequency that matches very well the numerical results even for = 1 (the real part is not such a good approximation for low ). Our geometric optics results are independent of the spin of the perturbing field and so they should agree with the geometric optics results for massless scalar field photon sphere modes in Ref. [19]. Consider a null geodesic x µ (τ ) of a RNdS black hole. By spherical symmetry there is no loss of generality in assuming that the geodesic is confined to the equatorial plane θ = π/2. There are conserved quantities associated to the Killing fields K = ∂/∂t and χ = ∂/∂φ: e ≡ −K µẋ µ and j ≡ χ µẋ µ , where the dot represents derivative with respect to the affine parameter τ . This giveṡ The radial motion is governed byṙ where we can check that r + ≤ r s ≤ r c .
The orbital angular velocity (Kepler frequency) of our null circular photon orbit can now be computed using (5.2), (5.5) and (5.7) yielding We now have to compute the largest Lyapunov exponent λ L , measured in units of t, associated with perturbations of an unstable circular photon orbit r(τ ) = r s . This is done considering perturbations r(τ ) = r s + δr(τ ) of the radial geodesic equation (5.3). Small deviations obey the linearized equation which has solution being the desired (largest) Lyapunov exponent. Note that C is an integration constant and the unstable photon orbit parameters r s and b s are given in terms of the RNdS parameters {L, M, Q} by (5.7). Finally, one can reconstruct the spectrum of the photon sphere family of quasinormal modes with 1 using [36][37][38][39][40][41][42][43][44] ω WKB ≈ Ω c − i n + 1 2 λ L , (5.11) where n = 0, 1, 2, . . . is the radial overtone. Note that this geometric optics/WKB approximation is universal in the sense that it is blind the particular sector of perturbations we look at. That is, it is expected to be a good approximation to both photon sphere modes Φ ± (or for a massless scalar field [19]). Note that, at this order, Im(ω WKB ) is independent of (assuming 1) while Re(ω WKB ) does depend on . One might wonder whether next-to-leading order corrections to this result might change significantly (5.11), especially near extremality. However, the corrections to Im(ω) are of order 1/ so, for any fixed background, the corrections to Im(ω) can be made arbitrarily small by taking sufficiently large. 18 So the WKB results for Im(ω) should be reliable for sufficiently large .
We can now analyse −Im(ω WKB )/κ − . In the left panel of Fig. 2 we plot this quantity for n = 0 (which yields the smallest value) as a function of the horizon radii ratio y + = r + /r c and charge ratio Q/Q ext . Over most of the RNdS moduli space we have −Im(ω WKB )/κ − < 1/2. Since we expect our result for to be exact as → ∞, we must therefore have β < 1/2 over most of the RNdS moduli space [19]. Thus the Christodoulou version of strong cosmic censorship is respected by most RNdS black holes. However, for any fixed y + , there is always a critical value for Q/Q ext (close to extremality) above which −Im(ω WKB )/κ − > 1/2. So there is the possibility of a violation of strong cosmic censorship by near-extremal RNdS black holes.  We will now compare the WKB prediction with our numerical results for the quasinormal frequencies of photon sphere modes. In the right panel of Fig. 2 we compare the n = 0 WKB result with our numerical results for −Im(ω)/κ − for the Φ − photon sphere quasinormal modes (with n = 0). From the plot we see that, when = 10, the WKB prediction is in excellent agreement with our numerical results. In fact even for = 2 the plot shows that the WKB prediction is in very good agreement with our numerical results. This agreement extends to other values of y + not shown in the plot. Note that, as expected, the agreement is very good for the imaginary part of the frequency but not so good for the real part (not shown in the plot). As a check of our numerical computations we have also confirmed that we reproduce some (the ones we searched for in our tests) of the quasinormal frequencies listed in [16] (note that this reference only computed what we call photon sphere modes).
Recall that to compute β defined in (2.13) we need to determine the spectral gap α. To determine α we need to find the slowest decaying quasinormal mode, i.e. the one with the smallest value of −Im(ω). We will now discuss which of the photon sphere modes has the smallest value of −Im(ω). There are two types of photon sphere modes: one corresponding to Φ − and another to Φ + . Our numerical results indicate that, for each type, the lower and n modes dominate. Therefore the slowest decaying photon sphere mode must be one of the following (with n = 0): (1) Φ − , = 2, or (2) Φ + , = 1.
Which of these two modes decays most slowly? For most of the black hole parameter space we find that the Φ + modes with = 1 decay most slowly. To illustrate this, in the left panel of Fig. 3 we plot −Im(ω)/κ − vs Q/Q ext at fixed y + for Φ + modes with = 1 (and n = 0) and Φ − modes with = 2 (and n = 0). We see that Φ + modes with = 1 typically have lower −Im(ω)/κ − (for fixed background parameters) than Φ − modes with = 2. However, there are small islands in the parameter space where the opposite occurs: see curve y + = 0.1 (red disks/circles) for Q/Q ext 0.9. A similar conclusion is reached from the right panel of Fig. 3. Here we plot the same modes but this time for RNdS with fixed Q and varying y + . We see that typically the Φ + , = 1 modes dominate over the Φ − , = 2 modes. However, for small y + there is a crossover and the = 2 modes become dominant.
These crossovers will not be a problem for our purposes. For each RNdS black hole we will compute numerically the two types (Φ ± ) of photon sphere quasinormal mode and then pick the one with lowest −Im(ω). This can then be compared with the results from the other families (dS and near-extremal) of quasinormal modes in order to calculate the spectral gap.

de Sitter family of modes
In the de Sitter limit, M = 0, Q = 0, the master equations for Φ + and Φ − are the same. To find the spectrum, we just need to take (4.2) or (4.10) and set M = 0, Q = 0. Using the radial coordinate (4.8) this yields the master equation where we have introduced the dimensionless frequencyω = ω r c (with r c = L for the dS solution). Note that = 1, 2, 3, . . . for electromagnetic modes Φ + and = 2, 3, 4, . . . for gravitational modes Φ − . The general solution of (5.12) is for arbitrary amplitudes A and B, with 2 F 1 (a, b, c; z) being the Gaussian Hypergeometric function. At the origin this solution behaves as Φ ± y=0 ≈ A y +1 + B y − and regularity at y = 0 thus requires that we set B = 0. On the other hand, a Taylor expansion about the cosmological horizon y = 1 yields . (5.14) Requiring outgoing boundary conditions demands that we discard the (1 − y) iω 2 solution. This can be done using the property Γ(−n) = ∞, n ∈ N 0 , i.e. requiring that Γ with = 1, 2, 3, · · · for Φ + modes and at = 2, 3, · · · for Φ − modes. So far we have restricted our attention to the dS limit (M = 0 = Q) of the RNdS solution. Naturally, RNdS has quasinormal modes Φ ± that in the dS limit reduce to (5.15). These are what we call the dS family of RNdS quasinormal modes. Numerically we find that these modes have purely imaginary frequencies and their wavefunctions are localized near the cosmological horizon. Fig. 4 shows some numerical results for the dS family of modes. For concreteness we do this illustration for modes Φ − with { , n} = {2, 0}. In the left panel we fix Q/Q ext and we plot the imaginary part of the frequency Im(ω r c ) as a function of the dimensionless ratio y + = r + /r c . By definition, dS quasinormal frequencies must approach (5.15) as y + → 0 and this is indeed the case (see red diamond). Note that the frequency changes substantially with y + . However, if we instead fix y + and vary Q then we find that the frequency does not change that much as Q/Q ext increases from 0 up to 1. This is illustrated in the right panel Imaginary part of the frequency as a function of Q/Q ext for fixed y + = 0.01 (green squares), y + = 0.05 (brown diamonds) and y + = 0.1 (black disks).  Figure 5. Left panel: de Sitter gravitoelectromagnetic mode Φ − with = 2 and n = 0: −Im(ω)/κ − as a function of Q/Q ext for fixed y + = 0.01 (green squares), y + = 0.05 (brown diamonds) and y + = 0.1 (black disks). Right panel: The ratio between the frequency Im(ω dS ) of de Sitter mode of the left panel with y + = 0.01 and the imaginary part of the geometric optics WKB frequency prediction (5.11) for the photon sphere modes of the same black holes. We see that, for a small black hole, −Im(ω dS ) is smaller than −Im(ω WKB ) for the full range of Q/Q ext .
of Fig. 4. This is similar to what was found for massless scalar field quasinormal modes in Ref. [19]. In particular, note that the result (5.15) works well for any small (y + 1) black hole, independently of Q.
Ultimately we will be interested in the ratio −Im(ω)/κ − . In the left panel of Fig. 5, we plot this quantity for the modes displayed in the right panel of Fig. 4. This plot illustrates that for the dS family, −Im(ω)/κ − can attain large values well above 1/2 or 1. The reason we choose to display data with small y + is because this is the region where the slowest decaying quasinormal modes belong to the dS family (as will be clear later, in Fig. 8).
For a small black hole, we can compare our analytical formula (5.15) for the slowest decaying ( = 1, n = 0) de Sitter modes with our WKB prediction (5.11) for the photon sphere modes. The latter is strictly valid for 1 but we found it worked well even for small . We find that in this small black hole limit, the de Sitter modes always decay more slowly than the WKB prediction for the photon sphere modes. This is illustrated in the right panel of Fig. 5 for the black hole family with y + = 0.01 (the same green square solutions shown in the left panel of the same figure). Thus for small black holes the = 1, n = 0 de Sitter mode is the slowest decaying mode belonging to either the de Sitter or photon sphere families.

Near-extremal family of modes and its near-horizon limit
The third family of quasinormal modes for RNdS black holes is called the near-extremal family since these modes are continuously connected to modes that can be identified in the near-horizon limit of the (near-)extremal RNdS solution, i.e. as r − → r + . The analytical analysis of the near-extremal modes of this subsection (and the near-Nariai modes of the next one) is very much inspired by ideas from Appendix A of [45] and [46,47]. This family of near-extremal modes is also present in the case of massless scalar field perturbations of a RNdS black hole [19].
In this subsection we will first perform an approximate analytical calculation of the near-extremal quasinormal modes using the near-horizon limit. We will then compare this to numerical results for these modes.
It is convenient to define the dimensionless quantities where σ ≥ 0 vanishes at extremality. The idea is to use the manifest SL(2, R) symmetry of the AdS 2 × S 2 near horizon geometry of an extremal RNdS black hole to simplify our calculation. The modes we seek, in the near extremal limit, are supported near the black hole horizon. So the limit we want to take has to accomplish two things: approach extremality, and zoom in near the black hole horizon. This can be achieved by sending σ → 0 while keeping z = x/σ fixed. We can anticipate that ω will vanish linearly as σ, so we define ω r c = δω σ and solve for δω in what follows.
We set and expand (4.2) − or (4.10) since the vector-type and scalar-type modes are isospectral − to leading order in σ. The resulting equation takes a simple form where we defined ϕ = y + Ξ δω , Note thatφ depends on δω, butη ± does not. The expression forη ± is easily shown to be real, and it is then manifestly positive. This will play an important role in what follows. Equation (5.18) can be readily solved in terms of Gaussian Hypergometric functions 2 F 1 via the following combination ± + 2 iφ ; 1 + 2 iφ ; z , (5.20) whereĈ (1) ± andĈ (2) ± are integration constants to be fixed via boundary conditions and a (1) We want to impose ingoing boundary conditions at the event horizon, i.e. regularity in ingoing Eddington-Finkelstein coordinates. This is equivalent to settingĈ (2) ± = 0. Next we need to impose a boundary condition at large −z. In principle this should be done by matching to a solution that is outgoing at the cosmological horizon. But we will follow the simpler approach of simply demanding that the solution vanishes at large −z. This can be motivated by the observation that near-extremal modes are highly localized near the event horizon and are therefore very small at large −z. Ultimately the justification for this boundary condition is that it gives quasinormal frequencies that match very well our numerical results.
At large negative values of z, we get The expansion (5.22) diverges at large positive values of (−z) because of the term proportional to (−z) This can be avoided if we set of the Gamma functions in the denominator to have a pole, which occurs for Γ(−n), with n ∈ N 0 = {0, 1, 2, . . .}. In particular, we quantize the frequency by demanding b (2) with n ∈ N 0 . This equation can be readily solved for δω and hence for ω: which simplifies considerably when written in terms of κ − : whereη ± is defined in (5.19). Note that these quasinormal frequencies are purely imaginary and that they all have −Im(ω)/κ − > 1/2. Which of these modes decays most slowly?
The imaginary part of the frequency increases with overtone number n so consider the fundamental (n = 0) modes. For given , we haveη − <η + , so the Φ − modes decay more slowly than the Φ + modes. It can also be checked that, for any y + ,η ± is an increasing function of . It follows that the slowest decaying modes covered by the above analysis are either the Φ − modes with = 2 or the Φ + modes with = 1 (as there are no Φ − modes with = 1). It is easy to show from (5.26) that it is always the Φ − modes with = 2 which decay the most slowly. The above calculation is, at best, valid only in the near-extremal limit, σ 1, and for small frequencies, |ω r c | 1. In the derivation of (5.26) we have only used the properties of the RNdS near-horizon geometry but no use of the full geometry or its far region was made. So we might question the validity of this approximation. To address this question, in Fig. 6 we compare (5.26) with the exact numerical data for the quasinormal mode family (with purely imaginary frequency) that we henceforth call the near-extremal modes. For illustrative purposes, we do this for the Φ − mode with = 2 and radial overtone n = 0. In the left panel of Fig. 6, we fix Q/Q ext = 0.999 and we plot −Im(ω)/κ − as a function of y + . Since we are very close to extremality we expect that (5.26) should be a good approximation. This is indeed what we find. The red dots representing the numerical data agree very well with the green curve corresponding to (5.26). On the other hand, as expected, the analytical approximation (5.26) becomes increasingly poor as we move away from extremality, i.e. as Q/Q ext moves further away from unity. This is illustrated in the right panel of Fig. 6, where we fix y + = 0.5 and see that the prediction (5.25) (green dashed curve) is an excellent approximation when Q ≈ Q ext but quickly becomes a bad approximation as Q decreases.
The validity of the approximation that leads to (5.26) was also tested in the following way. The fact that we just use the near-horizon geometry to get (5.26) suggests that these quasinormal modes have to be localized near the event horizon and very quickly decay away Figure 6. Near-extremal modes for the gravitoelectromagnetic mode Φ − with = 2 and n = 0. In both plots, the dashed green line is the analytical prediction (5.26), or (5.25), and the red dots are our numerical results. Left panel: −Im(ω)/κ − as a function of y + for near-extremal modes at fixed Q/Q ext = 0.999. The dashed blue curve is the WKB prediction (5.11) for the photon sphere modes (also with Q/Q ext = 0.999). This WKB blue curve continues to increase monotonically as y + decreases. Right panel: Im(ω r c ) as a function of Q/Q ext for near-extremal modes at fixed y + = 1/2. The dashed blue curve is again the WKB prediction (5.11) for the photon sphere modes. We see that for a wide range of charge Q the photon sphere modes decay more slowly than the near-extremal modes but, above a critical charge ratio of Q/Q ext ∼ 0.98, the opposite happens.
from it. Our numerical results confirm that this is the case: the numerical near-extremal mode wavefunctions are indeed localized near the event horizon, r = r + , becoming more localized as extremality is approached.
In summary, we find that the analytical prediction (5.25) works very well for nearextremal modes of near-extremal black holes. It is interesting to compare this analytical prediction, for the dominant Φ − , = 2 modes, to the extremal limit of our WKB prediction (5.11) for the photon sphere modes. This comparison is shown in the left panel of Fig. 6 for Q/Q ext = 0.999. If we go even closer to extremality then the blue curve moves to the right, and −Im(ω WKB )/κ − diverges in the extremal limit. Thus we see that, sufficiently close to extremality, the near-extremal modes always decay more slowly than the WKB prediction for the photon sphere modes. Thus, to the extent that the WKB prediction is valid at small (and, as we have seen, it seems to work well), our analytical results predict that, in a neighbourhood of extremality, the Φ − , = 2 near-extremal modes should be the slowest decaying modes belonging to either the near-extremal or photon sphere families. This is further illustrated in the right panel of Fig. 6 where we are at fixed y + = 0.5 and vary Q/Q ext : as we approach extremality, there is a critical value of the charge ratio above which the near-extremal modes indeed become more slowly decaying than the WKB photon sphere modes.
We can also compare the near-extremal family of modes to the de Sitter family. For the slowest decaying de Sitter modes, we see from Fig. 4 (right panel) that Im(ωr c ) does not vary much as we approach extremality. It follows that −Im(ω)/κ − diverges for the de Sitter modes as we approach extremality. This ratio remains finite for the near-extremal modes, hence the near-extremal modes decay more slowly than the de Sitter modes in a neighbourhood of extremality.
In summary, a combination of analytical and numerical calculations indicates that, in a neighbourhood of extremality, the slowest decaying quasinormal mode across all families is the Φ − near-extremal mode with = 2 and n = 0. Furthermore, we have an analytical prediction from (5.25) for the frequency of this mode. Hence (5.25) gives us an analytical prediction for the behaviour of β as we approach extremality. This is the green curve in the left panel of Fig. 6. We will discuss the implications of this below.

Nariai modes
RNdS black holes have three horizons, r − , r + and r c . In the previous subsection we considered the extremal limit where r − → r + . There is however another interesting limit − the Nariai limit − which occurs when r + → r c . The surface gravity remains non-zero in this limit. It is natural to wonder wether there is a fourth family of RNdS quasinormal modes that reduce to Nariai quasinormal modes in this limit.
For massless scalar field perturbations, the results of Ref. [19] suggest that these "Nariai modes" are a subset of photon sphere modes, rather than constituting a distinct fourth family of modes. In the Appendix, we will show that this is indeed the case for gravitoelectromagnetic modes. Therefore we do not need to consider the Narai modes as a distinct family.

Results
As explained above, for each type of perturbation (Φ + or Φ − ) we expect quasinormal modes to fall into three families (dS, photon sphere and near-extremal). Furthermore, from the discussion above, we expect that the slowest decaying quasinormal modes for each family and each type of perturbation to be given by the modes with the lowest allowed value of for that type of perturbation (this will be illustrated later in Table 1 for a particular black hole). Therefore our numerical calculations of quasinormal modes have focused on the two gravitoelectromagnetic sectors {Φ − , = 2} and {Φ + , = 1} since other sectors are expected to give more rapidly decaying modes.
As an example of how we classify the quasinormal modes emerging from our numerical calculations, we will consider the family of "lukewarm" RNdS black holes [48,49]. This is the 1-parameter subfamily of RNdS black holes that are in thermal equilibrium since the temperature of the event and cosmological horizons are the same i.e. κ + = κ c 19 . It turns out that this is equivalent to M = |Q| [48]. For a lukewarm hole Q Q ext = 1 1 + y + 3y 2 + + 2y + + 1 1 + 2y + (6.1) with Q/Q ext = 1/ √ 2 ∼ 0.707 for y + = 1 and Q/Q ext = 1 for y + = 0. We have discretized the lukewarm RNdS family with a numerical grid of 100 points for 0 ≤ y + ≤ 1, and we searched for the full spectra of frequencies solving each one of the relevant two perturbation equations as a quadratic eigenvalue problem for ω 2 . To evaluate the numerical convergence of our results we then took the frequency spectrum of each lukewarm solution and inserted it as a seed in a Newton-Raphson code, and we progressively increased the number of grid points along the radial direction 0 ≤ y ≤ 1 − see (4.8) − until we got the desired accuracy for the quasinormal frequency.  Figure 7. Results for the Φ − quasinormal modes with = 2 for lukewarm RNdS black holes.
Left panel: The filled marks identify the fundamental (n = 0) modes of the three families, namely photon sphere (black disks), near-extremal (red diamonds), and de Sitter (blue squares). The black circles represent the next 15 photon sphere overtones (n = 1, · · · , 15) and the 16 blue dotted lines represent the WKB approximation Im (ω WKB ) (n), n = 0, · · · , 15, for the photon sphere modes (valid for 1). The red diamond (in the de Sitter curve) represents the n = 0 pure de Sitter frequency Im(ω r c )| dS = −3. The green triangle (in the near-extremal curve) represents the n = 0 analytical approximation Im(ω r c )| NE = −2 in the limit where Q = Q ext , which for lukewarm RNdS occurs when y + → 0. Right panel: The three families of fundamental (n = 0) quasinormal modes. Here we plot −Im(ω)/κ − against Q/Q ext . The colour code is the same as for the left panel.
As an example, in the left panel of Fig. 7 we give our results for the imaginary part of the frequency for the Φ − modes with = 2. The black disks are the fundamental (n = 0) photon sphere modes. This identification emerges from the fact that they match the geometric optics/WKB approximation (5.11) for Im(ω WKB ) (blue dotted line). These modes also have Re(ω r c ) = 0 which distinguishes them from the purely imaginary dS and near-extremal modes. In the same figure, below this n = 0 photon sphere curve, we identify a total of 15 more curves with black circles. From the left/top to the right/bottom these are the photon sphere overtones n = 1, 2, · · · , 15. This identification follows from: 1) the fact that they match the geometric optics/WKB approximation (5.11) (see the associated 15 blue dashed curves 20 ), and 2) the number of radial zeros in the real and imaginary parts of the associated eigenvectors increases by one unit as n increases by one unit. For clarity of our presentation we decided not to plot the photon sphere modes with n ≥ 16. From the figure the reader can however understand that these accumulate on the right side of the plot. 21 Also on the left panel of Fig. 7 we also see a line of red diamonds. This is the fundamental (n = 0) near-extremal mode of the lukewarm RNdS family. 22 This identification emerges from the fact that: 1) these frequencies are purely imaginary, 2) they converge to Im(ω r c )| NE = −2 in the lukewarm extremal limit y + → 0 (see the green triangle), as dictated by the analytical analysis (5.26), and 3) the eigenvectors of these modes (real functions) are very localized near the event horizon.
Also on the left panel of Fig. 7 there is a curve of blue squares. This is the n = 0 de Sitter family of modes because: 1) these modes are purely imaginary, 2) they converge to Im(ω r c )| dS = −3 as y + → 0 (see the red diamond), in agreement with the analytical analysis (5.15), and 3) the eigenvectors of these modes (real functions) are very much localized near the cosmological horizon. 23 To conclude our analysis of the left panel of Fig. 7, the numerical solution of the quadratic eigenvalue problem gives the full spectrum of eigenfrequencies and associated eigenvectors. We have identified each family of modes that appears in the spectrum using the information discussed in section 5. All the numerical data fits in one of the three classes of modes (de Sitter, photon sphere or near-extremal). Still in the lukewarm family of RNdS, we did a similar analysis for the other relevant sector of perturbations, {Φ + , = 1}, with similar results.
Recall, that we are studying the quasinormal spectra of RNdS to find the spectral gap α in order to calculate β defined by (2.13). To calculate α we need to determine the slowest decaying quasinormal mode across the two types of perturbation (i.e. Φ + and Φ − ) in all three families of quasinormal modes. We can illustrate this with the lukewarm family of RNdS black holes. Focus first on the sector of perturbations {Φ − , = 2} already studied in the left panel of Fig. 7. Clearly, for our purposes, it is enough to compare the leading (n = 0) overtone −Im(ω)/κ − for the three families of modes. This is done in the right panel of Fig. 7. We see that for lukewarm black holes and in the {Φ − , = 2} sector, photon sphere modes have the lowest −Im(ω)/κ − for Q/Q ext 0.955. However, for 0.955 Q/Q ext ≤ 1 the slowest decaying modes are the near-extremal ones. The de Sitter modes are irrelevant for the spectral gap discussion of lukewarm black holes. This analysis still does not identify β for the lukewarm family. For that, we have to repeat the analogue of the right panel of Fig. 7 for the other sector {Φ + , = 1} of perturbations and β is then the minimum of −Im(ω)/κ − over the two sectors of quasinormal modes.
Moving away from the lukewarm family, we will now describe our results for the full moduli space of RNdS black holes. We have spanned the full parameter space 0 ≤ y + ≤ 1 however remarkable that the 1 approximation (5.11) is so accurate for = 2. 21 Without much effort, i.e. without increasing the resolution beyond the value required to have the accuracy desired for the leading overtones, we were able to capture the first ∼40 photon sphere overtones. 22 The higher, n ≥ 1, near-extremal overtones have lower Im(ω rc) and are not shown. 23 The higher, n ≥ 1, de Sitter overtones have lower Im(ω rc) and are not shown. and 0 ≤ Q/Q ext ≤ 1 using a numerical grid with 100 points along y + and another 100 points along Q/Q ext . That is to say, we have computed the {Φ + , = 1} and {Φ − , = 2} quasinormal modes for a total of 10 4 RNdS black holes. Where necessary we further zoomed in a particular region of parameter space, e.g. near Q/Q ext ∼ 1 and/or y + ∼ 0 or y + ∼ 1. Again, all the numerical modes were identified as belonging to one the three families of modes (de Sitter, photon sphere or near-extremal). It is in this sense that we are confident that, for each of the 10 4 RNdS black holes that we studied, the frequency spectra of quasinormal modes belongs to one of the three families discussed in section 5 and no fourth family exists.   The discontinuities in the derivatives of these curves occur across the boundaries of the different regions A, B, C. Note that β > 2 sufficiently close to extremality, and large near-extremal black holes can have arbitrarily large β. The black curve is the analytical prediction for near-extremal modes and the black disks correspond to lukewarm black holes.
Our main results for the spectral gap are presented in Fig. 8. In the left panel we show a density plot where we plot β = α/κ − as a function of the horizon ratio y + = r + /r c and charge ratio Q/Q ext . We identify three regions A, B, C separated by three black curves. In region A the spectral gap is dominated by the de Sitter modes. That is, in this region, the slowest decaying quasinormal mode is a de Sitter mode. This region A extends all the way down to Q → 0, i.e. de Sitter modes dominate the region of parameter space described by very small values of y + . On the other hand, in region B it is the photon sphere modes that dominate. Finally, in region C, i.e. in a band of parameter space around extremality Q/Q ext ∼ 1, it is the near-extremal modes that dominate.
The left panel of Fig. 8, also shows a red dashed curve. This curve identifies solutions with β = 1/2 and, above it, we have a region of parameter space near extremality where the solutions have β > 1/2 (see also the density plot legend). It follows from the discussion of section 3.3 that, in this region, the Christodoulou version of strong cosmic censorship is violated (for smooth initial data) by gravitoelectromagnetic perturbations. These results are similar to the results for massless scalar field perturbations presented in Ref. [19]. However, there is an important qualitative difference between our results and the results for massless scalar field perturbations. In the massless scalar case one always has β < 1 [19]. But in our case we can have β > 1. This is apparent in the right panel of Fig. 8, which plots β against y + for different values of Q/Q ext . The black curve corresponds to the analytical prediction (5.26) for near-extremal modes Φ − with n = 0 and = 2. From the discussion at the end of section 5.3 we expect this analytical prediction to be reliable as we approach extremality. Our plot shows that this analytical result does indeed give an accurate prediction for the value of β close to extremality. From the plot we see that, not only that do near-extremal black holes have β > 1/2, but in fact they have β > 2, which implies (section 3.3) that the C 2 version of strong cosmic censorship is violated for smooth initial data. In fact, for any r, by taking y + large enough we can find a near-extremal black hole for which β > r (the appropriate value of y + can be determined from (5.26)). Hence, for any r, the C r version of strong cosmic censorship is violated for smooth initial data. 1 2 Finally, in Table 1 we present detailed numerical results for a particular near-extremal black hole which violates the C 2 version of strong cosmic censorship. For this particular example we have computed not just the {Φ + , = 1} and the {Φ − , = 2} modes but also the {Φ + , = 2} modes and both types of mode with = 3, 4. From the table we see that the slowest decaying mode for this particular black hole is the {Φ − , = 2} mode, in agreement with the discussion at the end of section 5.3. This black hole has κ − r c = 0.098005 and so β = 0.200374/0.098005 = 2.04.

Taking the rough with the smooth
We have reviewed the reason why quasinormal modes determine the behaviour at the Cauchy horizon of linear perturbations arising from smooth initial data. By calculating the gravitoelectromagnetic quasinormal modes of RNdS black holes we have shown that, the Christodoulou and C 2 formulations of strong cosmic censorship are always violated close to extremality, and, for any r, the C r formulation is violated close to extremality for a sufficiently large black hole. Thus gravitoelectromagnetic perturbations exhibit a much worse violation of strong cosmic censorship than the massless scalar field perturbations considered in Ref. [19].
We emphasize that this violation of strong cosmic censorship in Einstein-Maxwell theory does not occur in pure Einstein gravity. Ref. [23] showed that any non-extremal Kerr-dS black hole has slowly decaying photon sphere gravitational quasinormal modes which ensure that the Christodoulou version of strong cosmic censorship is respected for smooth initial data.
As we have discussed above, Dafermos and Shlapentokh-Rothman (DSR) have shown that one can rescue strong cosmic censorship for RNdS black holes at the expense of considering rough initial data [25]. We have explained how a lack of smoothness of the initial data is also required to make sense of the older argument of Ref. [18] in favour of strong cosmic censorship.
What are we to make of this? Should we allow rough initial data? In physics we often assume that it is sufficient to work with smooth initial data. However, in some theories, smooth initial data can lead to a rough solution. For example, a shock can form in a compressible perfect fluid. Once we accept the existence of shocks, it is natural to weaken the regularity of our initial data to allow for shocks present initially. So for a fluid it is natural to allow rough initial data. However, in Einstein-Maxwell (-scalar field) theory, if we start with smooth initial data then the solution will remain smooth throughout the domain of dependence of this data. Shocks do not form dynamically. So we are not forced to consider rough initial data.
On the other hand, rough initial data can be approximated by a sequence of smooth initial data labelled by an integer n, and all with the same energy as the rough data. The sequence of smooth solutions arising from such data will be close to the rough solution in a region of spacetime that becomes larger as n → ∞, and approaches the Cauchy horizon in this limit (this follows from Cauchy stability of the equations of motion). DSR's rough version of strong cosmic censorship indicates that one can find a sequence such that the energy at the Cauchy horizon diverges as n → ∞. Hence, even for smooth perturbations, the energy at the Cauchy horizon is not bounded by the initial energy. Even if the energy of a smooth perturbation does not diverge at the Cauchy horizon, it can still become arbitrarily large there.
Maybe for some reason one would want the initial data not just to have finite energy but also that the first k derivatives are square integrable, i.e. the initial data has finite H k norm. For example, such a condition might arise from the requirement that the leading higher derivative corrections to the equations of motion are negligible initially. DSR's rough version of strong cosmic censorship implies that there exist smooth initial data whose H k norm on a spacelike surface intersecting the Cauchy horizon is not bounded by the H k norm of the initial data. This suggests that generic smooth initial data for which the leading higher derivative corrections are negligible will give a solution for which the leading higher derivative corrections become large near the Cauchy horizon. This does seem to capture the physics of the strong cosmic censorship hypothesis, namely that there is always a breakdown of effective field theory at a Cauchy horizon.

Comments on quantum effects
The analysis of this paper has been entirely classical. In this section we will discuss the role of Hawking radiation [50] in enforcing strong cosmic censorship. Recall that the behaviour at the Cauchy horizon is determined by the late-time behaviour of the black hole solution. So we need to discuss the effects of Hawking radiation on this late time behaviour. In de Sitter spacetime, we have to account for Hawking radiation both from the black hole horizon and from the cosmological horizon [51].
Consider first pure Einstein-Maxwell theory. In this case there are no charged particles and so Hawking radiation cannot change the charge of the black hole. If the black hole has a higher temperature than the cosmological horizon then it will radiate photons and gravitons and its temperature will decrease. If it has a lower temperature than the cosmological horizon then it will absorb photons and gravitons emitted by the cosmological horizon and the black hole temperature will increase. Thus Hawking radiation will drive the black hole towards a lukewarm solution for which the black hole and the cosmological horizon have equal temperatures, i.e. κ + = κ c [48].
We can approximate the late time solution as a (slightly perturbed) lukewarm solution and the behaviour near the Cauchy horizon will be determined by the behaviour near the Cauchy horizon of a lukewarm black hole. Fig. 8 (right panel) shows that small lukewarm black holes have 1/2 < β < 2 and so (in pure Einstein-Maxwell theory) they violate the Christodoulou formulation of strong cosmic censorship (for smooth initial data) but not the C 2 formulation. Thus it appears that Hawking radiation does not rescue the Christodoulou version of strong cosmic censorship in pure Einstein-Maxwell theory.
However, there is another way in which quantum effects can influence the geometry, namely via vacuum polarization. At late time, one would expect the quantum state of fields outside the black hole to approach the Hartle-Hawking state in the lukewarm black holes background. In this state, the results of calculations in a 2d toy model [52] (with conformally coupled quantum fields) indicate that T µν diverges at the Cauchy horizon. This divergence is proportional to (−V − ) −2 , which is not locally integrable at the Cauchy horizon hence one cannot make sense of the semi-classical Einstein equation G µν = 8π T µν there, even in the sense of weak solutions. This suggests that quantum effects may rescue strong cosmic censorship. It would be interesting to confirm this with a calculation of T µν in the Hartle-Hawking state near the Cauchy horizon of a lukewarm black hole.
Of course, in the real world there exist charged particles e.g. electrons, that an electrically charged RNdS black hole can emit as Hawking radiation, and thereby decrease its charge. If the radiation of charged particles is rapid compared to the radiation of uncharged particles then the black hole will first lose most of its charge, and then evaporate away completely. If the radiation of charged particles is slow compared to the radiation of uncharged particles then the latter would tend to push the black hole onto the lukewarm family of solutions as above. The emission of charged particles would then cause the charge gradually to decrease whilst remaining within the lukewarm family. But ultimately the black hole would evaporate away completely. Note that this conclusion does not depend on the mass of the charged particles. This is because, unlike in flat spacetime, particles of any mass can escape the black hole by tunnelling through the potential barrier separating the event horizon from the cosmological horizon. In other words, the mass of the particle is redshifted away at the cosmological horizon.
It seems that Hawking radiation of charged particles will ensure that strong cosmic censorship is respected. However, one could also imagine a magnetically charged RNdS hole, perhaps formed by pair creation in de Sitter spacetime [48]. By performing an electromagnetic duality rotation, our results on gravitoelectromagnetic perturbations of electrically charged RNdS holes map to identical results for magnetically charged holes. If there are no magnetically charged particles then such black holes will evolve via Hawking radiation to lukewarm holes, which will behave as discussed above.
The procedure described in section 5.3 also applies to the current Nariai analysis as long as we do the identifications x → X and σ → δ in these formulas. We want modes that are regular at X = 0 (which corresponds to have outgoing boundary conditions at the cosmological horizon in the full geometry) and the condition that the solutions should decay at large X quantizes the frequencies. The latter condition is poorly motivated but it gives results that agree well with our numerics.
The frequencies (A.2) of near-Nariai modes have a real and imaginary part. This analytical approximation is strictly valid in the near-Nariai limit, δ 1 (i.e. y + → 1), Q Q ext and for small frequencies, |ω r c | 1. So what are these modes? Do they represent a fourth family of modes in RNdS?
To answer this question we attempted different strategies. In one of them we fix the black hole parameters and the quantum number and we solve the perturbation master equation as an eigenvalue problem to find the frequencies that are allowed in the background. After identifying the frequencies − including the first few overtones n ≥ 0 − that describe the 1) de Sitter, the 2) photon sphere and 3) near-extremal modes we do not find evidence of a new fourth family of modes. In a second approach, we use a Newton-Raphson algorithm whereby we give directly (A.2) as a seed (in a region of parameter space, i.e. y + ∼ 1, where it is a good approximation). Again, such a code does not converge to a new fourth family of modes. Instead, this Newton-Raphson code always converges for the family of modes that we have already identified as being the photon sphere of modes. Moreover, this happens not only when we search for the leading radial overtone, n = 0 in the seed (A.2), but also for the first few other overtones that we attempted (n = 1, 2, 3).
We consider that our experiments give good evidence to support the claim that there is no fourth family of quasinormal modes that can be associated to a Nariai origin. Instead, the Nariai frequencies simply give a good approximate description of photon sphere modes in the limit where y + → 1. These conclusions are best illustrated in Fig. 9. In the left panel, we fix y + = 0.99 and the dashed orange curve describes the analytical Nariai expression (A.2) prediction whereas the dashed blue line is the WKB prediction (5.11) for the photon sphere modes. The black diamonds represent the outcome of our Newton-Raphson search when we give the Nariai frequency (A.2) as a seed. Both the near-Nariai and WKB photon sphere predictions agree very well with the numerical data although, as expected, the near-Nariai result works less well at large Q/Q ext . In the right panel of Fig. 9, we fix Q/Q ext  Figure 9. Photon sphere family of modes Φ − with = 2, n = 0 and the Nariai limit. In both plots, the dashed orange line refers to the Nariai analytical prediction (A.2) for −Im(ω Nariai )/κ − (with = 2, n = 0), while the dotted blue curve is the analytical photon sphere prediction (5.11) for −Im(ω WKB )/κ − (with n = 0). Left panel: −Im(ω)/κ − as a function of Q/Q ext at fixed y + = 0.99, i.e. r + = 0.99 r c . Right panel: −Im(ω)/κ − as a function of y + at fixed Q/Q ext = 0.4995. Note that, as expected, the Nariai analytical prediction is a good approximation only near y + ∼ 1. It seems to describe the y + ∼ 1 limit of the photon sphere modes (black diamonds and WKB dotted blue line). and vary y + . Again, the black diamonds represent the outcome of our Newton-Raphson search when we give the Nariai frequency (A.2) as a seed. As expected, the near-Nariai prediction (A.2) is very good for 1 − y + 1 but quickly gets worst as y + decreases. The black diamonds turn out to be exactly the photon sphere modes that we had already found in an independent analysis. This is confirmed by the agreement with the WKB prediction (5.11). The results presented in this plot are qualitatively the same for any other value of the charge ratio Q/Q ext (and we did a fine-tunned search which spanned the full interval 0 < Q/Q ext < 1).
Analysis similar to the one displayed in Fig. 10 further reinforce our conclusion. In this figure, we take Q/Q ext = 0.0999 and we focus our attention in the interval 0.98 < y + < 1, i.e. very close to the Nariai limit y + → 1. We display the modes we obtain with a Newton-Raphson search when we give analytical Nariai expression (A.2) as a seed. In the left panel we plot the imaginary part of the frequency, while in the right panel we plot the real part of the frequency. The left panel exemplifies again what we already know: as discussed in the previous cases, both the WKB expression (blue dotted curve) and the Nariai expression (orange dashed curve) give good approximations for Im(ω r c ) when y + ∼ 1 and as we move away from the Nariai limit, the analytical expression (A.2) starts being a less good approximation. On the other hand, the right panel of Fig. 10 shows that the analytical Nariai expression (A.2) yields an approximation for the real part of the frequency that is actually even better than the WKB approximation (5.11), as long as 1 − y + 1. However, we would expect that including higher order terms in 1/ would improve the accuracy of  Figure 10. Photon sphere family of modes Φ − with = 2, n = 0 and the Nariai limit for RNdS black holes with Q/Q ext = 0.0999. In both plots, the dashed orange line refers to the Nariai analytical prediction (A.2) for = 2, n = 0, while the dotted blue curve is the analytical geometric optics/WKB photon sphere prediction (5.11) for = 2, n = 0. Left panel: Im(ω r c ) as a function of y + close to the Nariai limit y + ∼ 1. Right panel: Re(ω r c ) as a function of y + close to the Nariai limit y + ∼ 1.
the WKB prediction, which is already remarkably accurate given that we are working with = 2. To conclude, we have shown that the Nariai result (A.2) simply describes the photon sphere family of quasinormal modes in the y + → 1 limit. In the case of a massless scalar field perturbation of RNdS, the analysis of [19] reached the same conclusion.