The origin of the Reissner-Nordstr\"om de Sitter instability

We give strong numerical evidence for the existence of an instability afflicting six-dimensional Reissner-Nordstr\"om de Sitter (RNdS) black holes. This instability is akin of the Konoplya-Zhidenko instability present in RNdS black holes in seven spacetime dimensions and above. Moreover, we perform a detailed analysis of the near-horizon limit of extremal RNdS black holes, and find that unstable gravitational modes effectively behave as a massive scalar field whose mass violates the AdS$_2$ Breitenl\"ohner-Freedman bound (if and only if $d\geq 6$), thus providing a physical argument for the existence of the instability. Finally, we show that the frequency spectrum of perturbations of RNdS has a remarkable intricate structure with several bifurcations/mergers that appears unique to RNdS black holes


Introduction
Kodama and Ishibashi proved that Reissner-Nordström de Sitter (RNdS) black holes in d = 4 and d = 5 spacetime dimensions are linearly-mode stable [1]. Therefore, the finding by Konoplya and Zhidenko [2,3], that RNdS black holes can be unstable to gravitational perturbations if they live in d ≥ 7 spacetime dimensions (n = d − 2 ≥ 5) came with some surprise. The existence of this instability was further confirmed in [4] where it was also noted that there is a simple (necessary but not sufficient) criterion − originally due to [5] − that predicts the existence of the instability. Essentially, [5] proposes that a system should be unstable whenever an integral (between the event and cosmological horizons) of the Schrödinger potential of the perturbation is negative. The instability of [2] is in these conditions [4]. More recently, this RNdS instability was also studied in the framework of the large d limit of general relativity [6].
There are however some fundamental questions that are left unanswered by [2,4,3,6] and that we address in this manuscript. Firstly, we would like to understand the physical origin of this instability, and in particular why it only appears in higher-dimensions. This is particularly important, since for d = 4 and d = 5 Kodama and Ishibashi proved linearmode stability [1]. It is also important in the context of recent studies on strong cosmic censorship violation in RNdS black holes [7][8][9][10][11][12][13][14][15][16][17][18]. Secondly, [2,4] find that the system is unstable for d ≥ 7 but their analysis leaves the stability properties of the d = 6 case undetermined. Again, stability is proven only for d = 4, 5 so could it be that there is an instability also for d = 6?
In this paper we address these two questions. In section 2 we start by reviewing the RNdS black holes and its relevant sector of perturbations that can be unstable. Then, in section 3 we point out that there is a criterion for instability that predicts and, at the same time justifies, an instability in RNdS black holes. This instability criterion was conjectured by Durkee and Reall [19] and later proved by Hollands and Wald [20]. Strictly speaking it was proven only for non-positive cosmological constant backgrounds, Λ ≤ 0, but, as argued in [20], it should also hold for de Sitter backgrounds (Λ > 0). Applied to our system, this Durkee-Reall instability criterion essentially states that that one should take the nearhorizon limit of the gravitational perturbation master equation about the extremal RNdS black hole. After the limit is taken, the gravitational master equation effectively reduces to a Klein-Gordon equation for a massive scalar field in an AdS 2 background. If this nearhorizon effective mass is smaller than the AdS 2 Breitenlöhner-Freedman (BF) mass bound then the full RNdS extremal geometry should be unstable. By continuity this instability should extend away from extremality. As shown in section 3, this criterion predicts the existence of an instability for RNdS black holes in d ≥ 6. In section 4, we numerically solve the perturbation master equation and we confirm that the instability is indeed present for RNdS black holes in d ≥ 6 dimensions. Therefore, the two main outcomes of our study are: 1) we provide a physical origin of the Konoplya-Zhidenko instability (violation of the AdS 2 BF bound), and 2) we establish the presence of the instability in d = 6.
Adding to these two main results, we will also take the opportunity to explore in detail the properties of the instability. For example, we will establish that the instability is present in a broader region of the parameter space than originally reported in [2,4,3]. Indeed, in one of our studies we will look for the instability directly in the extremal configuration and will find that all extremal RNdS black holes are unstable for d ≥ 6 1 . That is, [2,4,3] identified the presence of the instability only for large values of r + /r c (where r + and r c are the event and cosmological horizons of RNdS) but we find that the instability is actually present in the whole range r + /r c ∈]0, 1[ for extremal RNdS. Our findings are not in conflict with the Durkee-Reall instability criterion, since the latter is known to be a sufficient, but not necessary condition for the existence of an instability. Additionally, we will also analyse the properties of the instability away from extremality and we will also look directly for the onset of the instability. Finally, we will also study in some detail the quasinormal mode structure of the perturbations as we span the two-dimensional parameter space of RNdS black holes. We will find an intricate network of quasinormal mode branches with interesting bifurcations/mergers that, to the best of our knowledge, seem to be unusual in black hole perturbations (at least in Λ ≤ 0 backgrounds).

Reissner-Nordström de Sitter black hole and its perturbations
We work with the Einstein-Maxwell theory, in d = n + 2 spacetime dimensions (n ≥ 2), with a positive cosmological constant Λ described by the action and where R is the Ricci scalar of the metric g, L is the de Sitter length scale, and F = dA is the Maxwell field strength associated with the Maxwell potential A.
A known solution of this theory is the Reissner-Nordström de Sitter (RNdS) black hole. In static coordinates, the gravitational and electric fields of this solution with mass M and charge q parameters are with dΩ 2 n being the line element of a unit radius S n and For an appropriate range of parameters, specified below, f has three real positive roots r − ≤ r + ≤ r c corresponding to the Cauchy horizon CH, event horizon H + and cosmological horizon H C respectively. We can express M and L in terms of r + , r c and q. The temperature of the event and cosmological horizons are, respectively, given by and T c = − f (rc) 4π . When T + = 0 we have an extremal RNdS black hole. This happens for q = q ext with The Einstein-Maxwell equations of motion are invariant under the scaling g → λ 2 g, A → λA and L → λL, with λ ∈ R, which we can use to construct dimensionless quantities in units of r c . Therefore, we choose to parametrize the RNdS solution using the dimensionless parameters q/q ext and y + . We are interested on gravitoelectromagnetic perturbations of RNdS. These were studied in detail by Kodama-Ishibashi (KI) in [1]. Perturbations of (2.2) can be analysed according to how they transform under diffeomorphisms of the S n sphere. There are a total of three families of perturbations that decouple from each other, namely the scalar, vector and tensor perturbations. These perturbations are built from scalar, vector and tensor harmonics on S n , respectively. We are primarily interested in scalar perturbations, which are built from spherical harmonics S ( x), where x collectively parametrise coordinates on the n− sphere. These harmonics are such that with λ S = ( + n − 1) and ≥ 0 being an integer. Modes with = 0, 1 were shown to be pure gauge in [1]. In particular, modes with = 0 describe changes in the mass of the background RNdS black hole while = 1 modes represent translations. Onwards, we shall take ≥ 2.
In the Schwarzschild limit, q = 0, Φ S − (t, r) and Φ S + (t, r) describe, respectively, purely gravitational and purely electromagnetic perturbations [1]. However, when q = 0 the gravito-electromagnetic perturbations are coupled. As described in [1], one can introduce a separation anstaz of the form which introduces the frequency ω. One can then manipulate the Einstein and Maxwell equation to find a decoupled ordinary differential equation for the radial master field Φ − ω (r) (see Eq. (5.59) of [1]), namely: where the potential V − (r; n, r + , r c , q, ) can be found in equations (5.61)-(5.63) of [1].

Near-horizon criterion for instability
As reviewed below, the near-horizon limit of an extremal RNdS black hole is described by the direct product spacetime AdS 2 × S n . One of the main observations of the current manuscript is that the existence of some instabilities of the full RNdS black hole can be inferred from studying the behaviour of the perturbation equation in the near horizon limit. More concretely, in the near-horizon limit the Kodama-Ishibashi master equation (2.7) reduces to an equation that has the form of a Klein-Gordon equation for a massive scalar field in AdS 2 . We will confirm that, according to the Durkee-Reall criterion [19,20], when this near-horizon effective mass violates the AdS 2 Breitenlöhner-Freedman mass bound then the full RNdS geometry is unstable. 2 This happens for d ≥ 6 (n ≥ 4). We will find that the AdS 2 BF bound violation gives a sufficient (but not necessary) criterion for the presence of an instability and also justifies its origin.
To get the near-horizon geometry of the extremal RNdS black hole, one first takes (2.2) with q = q ext and zooms in around the event horizon region by making the coordinate transformations: with L 2 2 = r 2 + 2ny n+1 + − (n + 1)y 2n + − n + 1 (n − 1) −4ny n+1 + + (n + 1) y 2n + + ny 2 + − (n − 1) 2 . (3.1) The near-horizon solution is then obtained by taking ε → 0 yielding (after a U (1) gauge transformation) This geometry is the direct product of AdS 2 ×S n and has a Maxwell potential that is linear in the radial direction. This limiting solution is still a solution of the (n + 2)-dimensional Einstein-Maxwell-dS theory. On the other hand, the AdS 2 metric solves the 2-dimensional Einstein-AdS equations with AdS 2 radius L 2 , R µν = −L −2 2 g µν .
One can now take the near-horizon limit (3.1), together with ω = εL −2 2ω directly on the gravitational master equation (2.7). After taking ε → 0 this yields where 2 2 is the d'Alembertian in AdS 2 and µ 2 = µ 2 (n, y + , ) is the near-horizon effective mass of the system. Its particular expression is long and not enlightening so we do not reproduce it here. It depends on the harmonic number , on the dimension d = n + 2 and on y + . The keypoint is that this near-horizon analysis is expected to provide a criterion for instability [19,20]: whenever this mass is smaller than the AdS 2 Breitenlöhner-Freedman (BF), µ 2 BF L 2 2 = −1/4, one should have an instability of the full RNdS black hole. We further expect this instability to be the one found in [2,4,6] (for d ≥ 7). If so, the violation of the AdS 2 BF bound effectively explains the origin of the instability found in [2,4,6].
In Fig. 1 we set = 2 and plot µ 2 − µ 2 BF L 2 2 as a function of y + = r + /r c for n = 4, 5, 6, 7, 8, 9 (see plot legends). We conclude that for n ≥ 4 there are always values of y + where µ 2 − µ 2 BF L 2 2 < 0 and thus for which the AdS 2 BF bound is violated and an instability should be present. Note that this includes the d = 6 (n = 4) case, which was also analysed in [2,4], but for which no instability was found. For n = 2, 3, i.e. d = 4, 5, one always has µ 2 > µ 2 BF and we do not display these curves in Fig. 1. This result is consistent with the fact that d = 4, 5 RNdS are stable to all linear-mode gravitational perturbations [1].
In the next section we will confirm that the AdS 2 BF bound violation observed in Fig.  1 (for n ≥ 4) indeed gives a sufficient criterion for the presence of an instability and also justifies its origin. But we will also show that this condition is not a necessary condition.  Figure 1. Difference between the effective near-horizon AdS 2 mass µ 2 L 2 2 and the AdS 2 BF bound µ 2 BF L 2 AdS2 = −1/4 as a function of r + /r c . When this quantity is negative one should expect an instability [19,20]. So far we discussed only harmonics with = 2. But we can also ask what happens if we consider modes with ≥ 3. Typically, we find that as the integer increases it becomes harder to get negative µ 2 − µ 2 BF L 2 2 . For example, in n = 4 or n = 5 dimensions, for ≥ 3 the AdS 2 BF bound is no longer violated for any y + . As another example, for n = 6, 7 the harmonics = 2, 3 generate a violation of the AdS 2 BF bound but this is no longer the case for ≥ 4. As a final example, for n = 8, 9 the harmonics = 2, 3, 4 generate a violation of the AdS 2 BF bound but this is no longer the case for ≥ 5.

Setup of the problem
Our aim is to solve the KI linear ODE (2.7), subject to relevant boundary conditions, to find the properties of linear mode perturbations in RNdS. More concretely, we want to fix the spacetime dimension n and the RNdS black hole − described by the dimensionless quantities {y + , q/q ext } − as well as the harmonic quantum number and look for modes that are unstable.
To achieve our aims we have followed a strategy that is anchored on three main studies: I. Generically, we have a non-extremal RNdS black hole with q/q ext < 1. We can solve (2.7) as a quadratic eigenvalue problem for ω to find the eigenvalue/eigenfunction pair ωr c , Φ − ω for a given set of {n, y + , q/q ext , }. For a fixed {n, } we can then repeat the analysis and scan the full 2-dimensional parameter space of RNdS black holes.
II. As argued in section 3 the instability has a near-horizon/extremal origin. So, when present, it should make its appearance at extremality and then extend away from extremality until it eventually shuts down. Therefore, we first consider the RNdS with q = q ext and solve (2.7) as a quadratic eigenvalue problem for the frequency to find the eigenpair ωr c , Φ − ω , for a given {n, y + , } directly at extremality.
III. Alternatively, instead of looking for the instability timescale, we can search directly for the onset of the instability whereby the frequency vanishes, ω = 0. To find this instability threshold, we solve (2.7) as a nonlinear eigenvalue problem for the black hole charge to find the onset charge q = q onset above which RNdS is unstable.
The above three strategies are complementary. In particular, the 2-dimensional surface in the plot {y + , q/q ext , Im(ωr c )} generated using the first study must: i) have the extremal curve with q/q ext = 1 obtained in the second study as a boundary line, and ii) intersect the onset curve ω = 0 of the third study when Im(ωr c ) = 0. So the three complementary studies also provide non-trivial independent checks of the numerical results we obtain.
To proceed, one needs to use a numerical scheme to solve our boundary value problems. For that it is good to first introduce a new compact radial coordinate, namely such that y = 0 describes the black hole horizon, r = r + , and y = 1 marks the location of the cosmological horizon, r = r c . Next, one needs to discuss the issue of the boundary conditions. These boundary conditions are different for the three eigenvalue problems I−III listed above and we discuss them separately.
Consider first the eigenvalue problem I. We have a non-extremal black hole and deformations propagate between the event and cosmological horizons where we impose as boundary conditions that the perturbations are regular in Eddington-Finkelstein coordinates. In particular, we take the modes to be ingoing at the event horizon and outgoing at the cosmological horizon. To find these boundary conditions we first do a Frobenius analysis about y = 0 to find that where T + is the event horizon temperature (its surface gravity divided by 2π) as defined below (2.3) and the tilde (here and in other expressions) is used to state that the quantities are measured in units of r c , e.g., Ingoing boundary conditions at the black hole horizon requires choosing the lower sign in (4.2). A similar analysis around the cosmological horizon, y = 1, yields whereT c is the (dimensionless) cosmological horizon temperature defined below (2.3). Imposing outgoing boundary conditions at the cosmological horizon demands choosing the lower sign in (4.4). We will use a Chebyshev collocation scheme to numerically solve for ωr c , Φ − ω (see [21] for a review on the subject). Hence, we want to perform a field redefinition that automatically enforces the above boundary conditions. This motivates the wavefunction redefinition where, for our choice of boundary conditions, Q ω (y) is a smooth function of y with a regular Taylor series expansion at both y = 0 and y = 1. The boundary conditions for Q ω (y) are then of the Neumann type, i.e. ∂ y Q ω (y) y=0,1 = 0. Consider now the eigenvalue problem II listed above. This time we have an extremal black hole, meaning that at the event horizon f vanishes quadratically. It follows that, instead of (4.2), this time a Frobenius analysis about y = 0 yields β ≡ 2ny + (n + 1)y 2n + − 2ny n+1 Requiring ingoing boundary conditions at the event horizon of the extremal RNdS black hole amounts to choose the upper sign in (4.6). For the extremal RNdS, a Taylor expansion about the cosmological horizon still yields (4.4) and choosing the lower sign in (4.4) still amounts to impose outgoing boundary conditions at the cosmological horizon, alike in the non-extremal case. Much like in the non-extremal case, since we use a pseudospectral collocation scheme, these boundary conditions are best implemented if we introduce the field redefinition such that the above boundary conditions simply translate into Neumann conditions, namely ∂ y Q ω (y) y=0,1 = 0 for the smooth function Q ω (y). Finally, let us consider the nonlinear eigenvalue problem III listed above. In this case we want to find the eigenpair q/q ext , Φ − ω that describes the onset of the instability in the non-extremal RNdS as we scan the y + parameter. The boundary conditions for this problem are straighforward: a Frobenius analysis about y = 0 (y = 1) indicates that we have a term that diverges as A log y (a log(1 − y)). We impose boundary conditions that eliminate these divergent terms: A ≡ 0, a ≡ 0 by taking pure Neumann boundary conditions for Φ − 0 (y).

Numerical results
The near horizon analysis of the extremal RNdS black hole and associated Durkee-Reall criterion [19,20] discussed in section 3 suggests that near-extremal black holes should be unstable for n ≥ 4. To confirm this is indeed the case, we first search for unstable modes directly in the extremal RNdS black hole using the numerical scheme II outlined in the previous section 4.1. We find three main results: 1. The gravitational instability in the RNdS black hole is present when the spacetime dimension satisfies n ≥ 4 (d ≥ 6) and = 2. In particular, it is present for d = 6 (n = 4), as predicted by the near horizon criterion, which is a result that was not established in previous literature. This is explicitly shown for n ≥ 5 in the right panel of Fig. 2 and for n = 4 in left panel of Fig. 2. These plots display the imaginary part of the dimensionless frequency, Im(ωr c ) as a function of y + = r + /r c (the real part of the frequency of the unstable modes vanishes).
2. At extremality, q/q ext = 1, the instability is present for any value of y + = r + /r c i.e. for y + ∈ [0, 1] (when n ≥ 4 and = 2) with Im(ωr c ) → 0 as y → 0 or y → 1 and attaining a maximum in between. This is explicitly shown for n ≥ 5 in the right panel of Fig. 2. For n = 4 it is computationally much harder to generate the associated numerical data but we believe this case should behave much like the n ≥ 5 cases. Note that in previous literature [2,4], the presence of the instability was established only for large values of r + /r c (and n ≥ 5) and later it was (incorrectly) claimed that the instability is present only above a critical value of y + [3].
3. The near horizon Durkee-Reall criterion [19,20] discussed in section 3 provides a sufficient but not necessary condition for the instability. That is to say, our numerical Figure 2. Instability timescale at extremality, q = q ext as a function of r + /r c (the real part of the frequency of these unstable modes vanishes). Left panel: The d = 6 (n = 4) case: the Durkee-Reall criterion is satisfied in between the two black dashed vertical lines. Right Panel: The distinct curves describe different spacetime dimensions, d = n + 2 ≥ 5 (see legend). The instability is stronger for higher n but, at extremality, it is always present for any value of r + /r c in the allowed range ]0, 1[. results show that the system is indeed always unstable whenever the AdS 2 BF bound is violated. But the instability also extends to values of the parameter space where the near-horizon criterion does not signal an instability. This is best seen comparing the analytical near-horizon predictions of Fig. 1 with the actual numerical results of Fig. 2. Recall that both these plots are for = 2 and q = q ext . Typically, the near-horizon analysis of Fig. 1 predicts instability for a finite range of y + . However, we find that, at extremality, the system is unstable in the full range y + ∈ [0, 1]. This is certainly the case for n ≥ 5 (see right panel of Fig. 2) and should also be true for n = 4.
Having established that all, 0 ≤ y + ≤ 1, extremal RNdS black holes are unstable for n ≥ 4 we might then ask how far away from extremality does the instability extend into. To address this question we first search directly for the onset of the instability using the numerical scheme III outlined in the previous section 4.1. This critical charge q onset above which the RNdS solution is unstable is shown in Fig. 3 for n = 5, 6, 7, 8, 9 and = 2. The left panel shows 1 − q onset /q ext as a function of y + while the right panel shows the logarithmic plot of the same quantity to zoom the details of the small y + region. We see that for a given dimension n, q onset /q ext decreases as y + grows from 0 into 1: the instability extends further away from extremality for high y + . On the other hand, for a given y + we see that increasing the dimension n favours the instability in the sense that q onset /q ext becomes smaller as n grows.
To have a broader perspective of the properties of the instability, next we search directly for the instability timescale in the non-extremal RNdS black hole. For this discussion we fix the dimension to be n = 5 and the associated data is shown in Fig. 4 (for other dimensions the plot is qualitatively similar). Recall that non-extremal RNdS black holes are parametrized by q/q ext and r + /r c and these are the two horizontal axes of Fig. 4. On the other hand the vertical axis is the imaginary part of the dimensionless frequency, Im(ωr c ) (the real part of the frequency of the unstable modes vanishes). In this 3-dimensional . Instability timescale for n = 5: imaginary part of the frequency Im(ωr c ) as a function of r + /r c and q/q ext (only shown the parameter space region where the instability is present, i.e. that has Im(ωr c ) ≥ 0). The brown dots in the plane q/q ext = 1 represent instability data collected independently using a numerical code at extremality. The green dots with Im(ωr c ) = 0 represent data collected independently using a numerical code for the onset of the instability.
plot we also show the extremal (brown) curve already displayed in Fig. 2 and the onset (green) curve already shown in Fig. 3 (for n = 5). The fact that the 2-dimensional surface describing the regime where RNdS is unstable ends on the extremal and onset curves obtained using independent numerical codes provides a non-trivial check of our numerical results. This plot reveals the following properties (some were already discussed above): • At extremality, q/q ext = 1, the instability is present in the whole range 0 ≤ y + ≤ 1.
• As we move away from extremality, we see that for a given (large) charge, namely in the range 0.9128 q/q ext ≤ 1 (when n = 5) the system is unstable only for y + above a (non-vanishing) critical value. For smaller charges q the system is stable. In equivalent words, for a given y + , RNdS black holes are unstable if their charge is above q onset , with q onset approaching q ext as y + → 0.
• The maximum strength of the instability is attained for black holes that are close to extremality, but not at extremality. This is better seen in Fig. 5 where we show the instability timescale for three families of RNdS solutions at constant y + (namely y + = 0.70, 0.95, 0.99) as a function of the dimensionless charge ratio q/q ext . We see that typically Im(ωr c ) grows as q/q ext increases but its maximum occurs slight before one reaches the extremal configuration q/q ext = 1. Up until now we have discussed only modes with = 2. This is because for each dimension n ≥ 4 this is either: i) the only harmonic for which the instability is present or ii) the harmonic mode where the instability is stronger. But we can also discuss briefly other harmonics. As discussed in section 3, the near-horizon analysis shows that as increases it becomes harder to get negative µ 2 − µ 2 BF L 2 2 . For example, in n = 4 or n = 5 dimensions, for ≥ 3 the AdS 2 BF bound is no longer violated for any y + . As another example, for n = 6, 7 the harmonics = 2, 3 generate a violation of the AdS 2 BF bound but this is no longer the case for ≥ 4. As a final example, for n = 8, 9 the harmonics = 2, 3, 4 generate a violation of the AdS 2 BF bound but this is no longer the case for ≥ 5. The full numerical analysis confirms these analytical predictions to be correct. In particular, the numerical analysis also concludes that when the instability is present for more than the = 2 harmonic, modes with lower are more unstable. To illustrate this, in Fig. 6 we take an extremal RNdS BH in n = 9 and compare the instability timescale of = 2 and = 3 modes. We see that the = 3 instability timescale is typically two orders of magnitude smaller than the = 2 timescale. We have also explicitly checked that for n = 9 the = 4 harmonic (but not higher 's) is also unstable but with a strength that is ∼ 5 orders of magnitude smaller than the = 2 instability. For example, for y + = 1/2 the timescale of three harmonics are:  Figure 6. Instability timescale at extremality, q = q ext , for n = 9 as a function of r + /r c for the = 2 (blue ) and = 3 (red ) harmonics. The = 2 curve was already shown (with same colour/shape code) in Fig. 2.
So far we have focused our attention on describing key properties of the instability present in RNdS black holes and confirming that we can understand its origin as being due to the violation of a relevant AdS 2 BF bound discussed in section 3. The relevant modes become unstable close to extremality.
Next, we give a broader view of the properties of these perturbations and follow the unstable modes as we move away from extremality all the way down to the Schwarzschild-dS black hole with q = 0. This analysis turns out to reveal an interesting quasinormal mode structure with properties that, to the best of our knowledge, no other known black hole quasinormal mode spectrum exhibits. In our scan of solutions we have found three main regimes that are distinct from each other and thus deserve a separate discussion. Illustrative examples of each of these three cases are the RNdS black holes with y + = 0.70, y + = 0.95 and y + = 0.99 (and n = 5, = 2). The main instability properties near extremality of these three cases were already discussed in Fig. 5. Next we extend the analysis of these modes all the way towards q = 0, i.e. well away from the region where the solutions are unstable and we also discuss "secondary" (less unstable or stable) modes that are nevertheless relevant to have a good overview of the quasinormal mode structure.
The key properties of these three families can be summarized as follows: 1. The case y + = 0.99 represents what typically happens for RNdS black holes that have very large y + (i.e. close to unit). The relevant frequency spectrum for this case is displayed in Fig. 7 (imaginary part) and Fig. 8 (real part of the frequency). In Fig. 7 we identify the dark red disk (•) curve that describes the (most) unstable mode that was already identified in Fig. 5 with red ⊗ (and in Fig. 4) and that has Re(ωr c ) = 0. In particular, the brown disk with q/q ext = 1 in this curve describes the extremal solution also identified in Fig. 2. We see that this curve has a maximum close to extremality (as discussed previously) and, decreasing q/q ext we find that the mode first becomes stable at q/q ext ∼ 0.92 (consistent with Fig. 3) and then reaches a bifurcation point A at q/q ext ∼ 0.861 where it joins the light blue lozenge (♦) and blue diamond ( ) curves. As seen in the inset plot of the left panel, the light blue ♦ branch extends all the way from A down to q = 0. As shown in the companion Fig.  8, the real part of the frequency of this branch is always non-vanishing (Re(ωr c ) = 0) and it decreases monotonically as q increases from zero towards the bifurcation point A with q/q ext ∼ 0.861. At this point A, the light blue ♦ branch bifurcates into the dark red • curve (upper branch) already discussed above and into the blue curve (lower branch). These two branches (• and ) both have Re(ωr c ) = 0. That is to say, a complex eigenvalue 3 bifurcates into two branches that have purely imaginary eigenvalues. Now let us follow the blue diamond branch. We see that it extends from point A with q/q ext ∼ 0.861 up to point B with q/q ext ∼ 0.902. At this point B it merges with the brown inverted triangle curve (that also approaches point B but coming from point D) and the solution now extends for higher q/q ext along the light blue ♦ branch that again has Re(ωr c ) = 0 (see the half-circle BC section of Fig. 8) up to point C with q/q ext ∼ 0.934. At this point C we have again a new bifurcation into two branches: the upper branch with orange triangles ( ) and the lower branch with blue diamonds ( ). Both these branches then extend independently towards extremality (q/q ext = 1) in such a way that Im(ωr c ) grows monotonically from negative into positive values (moreover Re(ωr c ) = 0 for both branches). This is best seen in the zoomed right panel of Fig. 7. For completeness we should also Figure 7. Frequency spectrum (imaginary part) for y + = 0.99 and n = 5, = 2 with zooms in the relevant regions where instability is present.  Fig. 7 where the brown inverted triangle curve continues as the green square ( ). Both curves have Re(ωr c ) = 0. The green square curve extends away from D to higher q/q ext where it has a new turning point but we do not analyse/discuss this further.
Although we have not done an exhaustive scan of the parameter space, we have collected enough evidence that the frequency spectrum for RNdS black holes with y + close to 1 has a bifurcation/merger structure qualitatively similar to the one shown in Fig. 7 and Fig. 8 (some of this evidence comes from Figs. 12-13, 14-15 and 16-17 to be discussed later). It is also important to point out that a bifurcation structure similar to the one seen around point A of Figs. 7-8 also emerges from a large d analysis of the instability: see Fig. 1 of [6]. However, at least for finite values of the dimension, the system seems to have a more complicated bifurcation structure than the one reported in [6]: there are further bifurcations/mergers like B and C in Figs. 7-8.
2. Next we consider the case y + = 0.95 which represents RNdS that have moderately large y + that is however not large enough to be close to unit. The relevant frequency spectrum for this case is displayed in Fig. 9 (imaginary part) and Fig. 10 (real part of the frequency). This case distinguishes from the previous case mainly because the bifurcations/mergers (like A and C of Fig. 7) cease to exist. Instead they are replaced by crossovers of quasinormal mode branches. Take the left panel of Figs. 9 and 10. As before, the light blue lozenge (♦) curve starts at q = 0 and has Re(ωr c ) = 0. This curve extends all the way to extremality (although it becomes very difficult to compute the properties of this curve for q/q ext above 0.98) with an intricate structure best seen in the zoom plots of the right panel of Figs. 9 and 10. Indeed, in the Im(ωr c ) plot the light blue lozenge (♦) curve has some zig-zagged regions that correspond to minima cusp points in the Re(ωr c ) plot: these are regions A, C, E, · · · in Figs. 9 and 10). These regions coincide with the points where the dark red disk (•) and orange triangle ( ) curves (that have both Re(ωr c ) = 0) crossover the light blue lozenge (♦) curve (note that the red disk curve describing the most unstable mode is precisely the one already displayed in Fig. 5 as a dark-red •). It is important to emphasize that A and C describe crossover not intersection points (since the eigenfunctions of the branches that intersect are distinct and points A, C and E in the light blue ♦ curve have Re(ωr c ) = 0 unlike in the • and curves). Comparing with the situation of Fig. 7 and Fig. 8 we can say that the bifurcation/merger points A and C of Fig. Figure 9. Frequency spectrum (imaginary part) for y + = 0.95 and n = 5, = 2 with zooms in the relevant regions where instability is present. 7 become crossover points in Fig. 9. Further note that, although not shown in Fig.  9, there should exist at least another family of modes to the right of the orange curve and "parallel" to it that should crossover the ♦ curve in region E. Finally, we should pointed out that the orange empty ( ) triangle modes of Fig. 9 are not the extension of the orange filled triangle ( ) modes of Fig. 7: this will become clear later when we discuss the curves and of Fig. 12.
3. Finally, we consider the case with y + = 0.70 which is a representative case of what happens with RNdS black holes that have intermediate and small values of y + . The relevant frequency spectrum for this case is displayed in Fig. 11 (imaginary part in left panel and real part in the right panel). This case distinguishes from the previous case because the structure of the spectrum is now very simple: the ziz-zag regions on the Im(ωr c ) plots of Fig. 9, i.e. the cusps in the Re(ωr c ) plots of Fig. 10, have now flatten out completely and we simply have very simple crossovers A and C between the light blue lozenge (♦) (which has Re(ωr c ) = 0) and, respectively, the red disk (•) curve that is unstable (already shown in Fig. 5 as a dark-red •) and the orange triangle ( ) curve.
As observed before the RNdS family of black holes spans a 2-dimensional parameter space that we are taking to be the dimensionless ratios y + = r + /r c and q/q ext . In the sequence of Figs. 7-8, 9-10 and 11 we have always fixed y + and analysed how the frequency spectrum changes as the charge of the RNdS black hole changes. Next, to have a 3-dimensional perspective of the system, we complement this analysis: we choose some relevant RNdS solutions with a fixed charge q/q ext and discuss how their frequency changes as y + changes (we take again n = 5 and = 2). Particularly illustrative cases that further help revealing the properties of the system are q/q ext = 0.95, q/q ext = 0.88 and q/q ext = 0.85. The analysis of these three cases unveils the following properties: 1) We start with a RNdS black hole with q/q ext = 0.95. In Fig. 12 (imaginary part) and Fig. 13 (real part) we show how the frequency ωr c changes as we vary y + . The curve on the top of Fig. 12 with dark red disks (•) is the branch that is unstable for large y + . It confirms, as already seen in Fig. 3, that RNdS with q/q ext = 0.95 become Figure 12. Left and right panels: Frequency spectrum (imaginary part) for q/q ext = 0.95 and n = 5, = 2 with zooms in relevant regions. We use the same shape/colour code (•, ♦, , , ) used in [9][10]. unstable for y + 0.881 (see left panel of Fig. 12). Additionally, out of an infinite family of modes with Im(ωr c ) < 0 for any y + , we further show the two families of modes that have the smallest |Im(ωr c )|. These are: i) the light blue ♦ branch (that has Re(ωr c ) = 0; see Fig. 13) that bifurcates at q/q ext ∼ 0.977 into two branches with Re(ωr c ) = 0, namely the orange triangle ( ) branch and the solid blue diamond ( ) branch, and ii) the orange empty triangle ( ) branch. A zoom in of the region close to y + = 1 is displayed in the right panel of Fig. 12. This plot at constant q/q ext = 0.95 complements well those plots at constant y + shown previously (namely Figs. 7,9,11). For example, in the right panel of Fig. 12, we identify the upper three points •, , with y + = 0.99. These three points are the three points with the same shape/colour code •, , with q/q ext = 0.95 of Fig. 7 (which describes RNdS solutions with y + = 0.99). As another example, in the left panel of Fig. 12, we identify the three points •, , with y + = 0.95 that are the three points •, , with q/q ext = 0.95 of Fig. 9 (which describes RNdS solutions with y + = 0.95). This is a good moment to pause and note, as observed in the end of the discussion of Fig. 9, that the orange (first introduced in Fig. 7) and modes (first introduced in Fig. 9) are indeed distinct. As a final example, in the left panel of Fig. 12, we further identify the three points •, , with y + = 0.70 that are the three points •, , with q/q ext = 0.95 of Fig. 11 (which describes RNdS solutions with y + = 0.70).
2) Next consider RNdS black holes with q/q ext = 0.88. In Fig. 14 (imaginary part) and Fig. 15 (real part) we show the variation of the frequency as we vary y + . From Fig. 3 (green curve) such black holes are stable for all values of y + . It follows that the curve on the top of Fig. 14 with dark red disks (•) has Im(ωr c ) < 0 for any y + : this is the family of modes that become unstable but only for q/q ext 0.881. Additionally, out of an infinite family of modes with Im(ωr c ) < 0 for any y + , we further show the families of modes that have the smallest |Im(ωr c )|. This is the light blue ♦ branch (with Re(ωr c ) = 0; see Fig. 15) that bifurcates at q/q ext ∼ 0.971 into two branches with Re(ωr c ) = 0, namely the blue diamond ( ) branch and the brown inverted triangle ( ) branch. A zoom in of the region close to y + = 1 is displayed in the right panel of Fig. 14. This plot complements previous plots at constant y + (of Figs. 7,9,11). For example, in the right panel of Fig. 14, we identify the upper three points •, , with y + = 0.99 which are the three points •, , with q/q ext = 0.88 in the left panel of Fig. 7 (which describes RNdS solutions with y + = 0.99; we do not show the green square family in Fig. 14). As another example, still in the right panel of Fig. 14, we identify the two points •, ♦ with y + = 0.95 that are the two points •, ♦ with q/q ext = 0.88 of Fig. 9 (which describes RNdS solutions with y + = 0.95; we do not show the orange triangle family in Fig. 14). As a final example, in the left panel of Fig. 14, we also identify the two points •, ♦, this time with y + = 0.70, that are the two points •, ♦ with q/q ext = 0.95 of Fig. 11 (which describes RNdS solutions with y + = 0.70). Figure 14. Left and right panels: Frequency spectrum (imaginary part) for q/q ext = 0.88 and n = 5, = 2 with zooms in relevant regions. We use the same shape/colour code (•, ♦, , ) used in [9][10] 3) Finally, we considered RNdS black holes with q/q ext = 0.85. In Figs. 16 and 17 we display the three family of modes with the lowest |Im(ωr c )| for this RNdS family. Such black holes are stable for all values of y + (see green curve in Fig. 3). In the main plot of the right panel, we identify the blue ♦ point with y + = 0.99 which is also the blue ♦ point with q/q ext = 0.85 in the left panel of Fig. 7 (which describes solutions with y + = 0.99). Also in the main plot of the right panel, we further identify the two points (blue ♦ and dark red •) with y + = 0.95 which are also the two points (blue ♦ and dark red •) with q/q ext = 0.85 in the right panel of Fig. 9 (which describes solutions with y + = 0.95). Finally, in the left panel, we identify the two points (dark red • and blue ♦) with y + = 0.7 which are the two points (dark red • and blue ♦) with q/q ext = 0.85 of Fig. 11 (which describes solutions with y + = 0.70). Note that one of the two branches (blue ♦ and dark red •) has the lowest timescale depending on the window of y + .

Discussion
We believe that our study reveals key aspects of the gravitational instability of Reissner-Nordström de Sitter black holes originally found in [2] and further studied in [4,6]. The fundamental novel results that we establish can be summarized as follows: • RNdS black holes are unstable when d ≥ 6 (see e.g. Fig. 2). For d ≥ 7 this instability was first established in [2,4,6]. In addition, we find it is also present for d = 6 (although with a much lower timescale).
• We have established the physical origin of the instability: it is present because in the near-horizon limit of extremal RNdS, the unstable gravitational modes effectively behave as a massive scalar field whose mass violates the AdS 2 Breitenlöhner-Freedman bound (if and only if d ≥ 6; see Fig. 1). By continuity, the instability then extends away from extremality.
• The instability criterion of the previous item is known as the Durkee-Reall conjecture [19] later proved by Hollands and Wald [20] when Λ ≤ 0. Our findings provide good numerical evidence for the Durkee-Reall conjecture with de Sitter asymptotics (Λ > 0), as already argued in [20]. It would be interesting to formally prove this result for Λ > 0 using the methods of [20].
• Our results also confirm that the Durkee-Reall instability criterion [19,20] provides a sufficient but not necessary condition for instability. For example, for extremal black holes the instability criterion typically predicts instability only for certain windows of the dimensionless horizon radii ratio r + /r c (see Fig. 1) but we find that actually the instability is present for any value of this ratio (at extremality and for d ≥ 6): see e.g. Fig. 2. In particular, this also means that there is no critical minimum value of y + for the existence of the instability unlike it was claimed in [3].
• RNdS black holes are parametrized by the dimensionless ratios y + ≡ r + /r c and q/q ext . For a given dimension n and y + , we have found the onset charge q onset /q ext above which RNdS becomes unstable: see Fig. 3.
• In addition to finding the instability timescale at extremality (Fig. 2) and the onset charge (Fig. 3) we have spanned the 2-dimensional parameter space of RNdS to produce the 3-dimensional plot of Fig. 4 that displays the instability timescale as a function of y + ≡ r + /r c and q/q ext . In particular, as best seen in Fig. 5, after its onset the instability first increases as we approach extremality until it reaches a maximum near but before extremality. Then its strength decreases slightly as we approach further and reach extremality.
• The instability is present for d ≥ 6 for the harmonic mode = 2. Typically, the instability then tends to get weaker and even disappear as the harmonic number grows (see Fig. 6). For example: in d = 6, 7 only the = 2 mode is unstable; in d = 8, 9 = 3 is also unstable; and in d = 10, 11 = 2, 3, 4 modes are unstable but not ≥ 5.
• If we follow the unstable modes to regions of the parameter space where the instability shuts down and if we include the second and/or third families of quasinormal modes with lowest |Im(ωr c )| we find an intriguing network of mode bifurcations/mergers that seems to be absent in black holes with Λ ≤ 0. This spectrum seems so unique and intriguing that we dedicated a special study to it. In the sequence of Figs There are quite a few natural extensions of our work. In particular, it would be interesting to frame the quasinormal mode structure we found in this manuscript in light of the three families of quasinormal modes used to study the strong cosmic censorship conjecture [7] for initial data close to RNdS black holes [8][9][10][11][12][13][14][15][16][17][18]. We should note, however, that once the instability sets in, a novel black hole geometry is likely to form [22]. Strong cosmic censorship should then be studied by investigating how slowly generic perturbations decay around the new hypothetical background geometry. Away from extremality, we should be able to settle this issue by studying the quasinormal mode spectrum of the relevant RNdS black hole. grant ST/J005673/1, STFC capital grant ST/H008586/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National e-Infrastructure.