1 Introduction

In recent years, advancement in technology brought a revolution in observational astrophysics that make it possible to test some of most intriguing predictions of general relativity [1, 2]. While, due to the uncertainty in the data, the possibility of a modified theory of gravity is not discarded right away, the stature of Einstein’s general relativity as the most successful theory of gravity remains unaltered [3,4,5]. This situation gives confidence to the scientists to propose interesting research work to test the theory in more and more extreme conditions. A possible alternative to this path is to find paradoxes within the theory of general relativity i.e., to check mathematical consistency of the theory. The existence of the Cauchy horizon in Kerr and Reissner–Nordström solutions is one of such paradoxes since the theory loses its predictive power beyond that region. However, soon it was realized that the Cauchy horizons are subject to blue shift instability that can turn them into curvature singularities under the influence of even small perturbations [6,7,8,9]. This phenomena led Penrose to propose the strong cosmic censorship conjecture which can be stated as follows: “for a generic initial data, the maximal Cauchy development (the largest manifold that is uniquely determined by Einstein’s field equations from a given set of initial data ) is inextendible as a \(C^{0}\) metric” [9,10,11,12,13]. This conjecture ensures that the observers who dare to cross the Cauchy horizon are torn apart by the infinite tidal forces. An another way to look at the problem is to consider the effect of a linear perturbations on the spacetime metric [14]. Here, the question of determinism is answered by considering nonlinear effects where the linear perturbations under consideration acts as the source of second order metric perturbations. Adopting the approach mentioned above, a version of strong cosmic censorship conjecture for a Einstein–Maxwell–scalar field system can be stated as follows : “the maximal Cauchy development of the stationary, axisymmetric solutions of a Einstein–Maxwell–scalar field system can not be extended across the Cauchy horizon (provided they exist!) with square integrable Christoffel symbols and scalar fields that live on Sobolov space i.e., \(\Phi \in H^{1}_{loc}\)” [14,15,16,17]. The same conclusion is obtained for the static and spherically symmetric solutions of a Einstein–Maxwell–Dirac field system with the Dirac field \(\varPsi \) replacing the scalar field \(\Phi \) [18, 19].

Although, any contradiction to this version of strong cosmic censorship conjecture is yet to be found for asymptotically flat black holes [9], the same conclusion may not be drawn in the presence of a positive cosmological constant [20]. In this scenario, the effect of blue shift amplification at the Cauchy horizon may be compromised by the exponential decay of the perturbations at event horizon. In fact, recently Cardoso et al. have found a finite parameter space where strong cosmic censorship conjecture gets violated in Reissner–Nordström-de Sitter black holes in presence of a massless, neutral scalar perturbations [21]. Since then, significant amount of works have been done on this topics, further confirming the same for several fundamental fields in a Reissner–Nordström–de Sitter background [18, 19, 22,23,24,25,26,27]. Interestingly enough, strong cosmic censorship conjecture seems to be always respected in rotating black holes in de Sitter background in the presence of scalar fields [14, 24, 27]. Here, the rotation of the black hole plays a crucial role to uphold the conjecture. So it is legitimate to check whether the conjecture still holds true in presence of other fundamental fields in rotating black hole background in order to better understand the effect of the rotation. In this paper, we study the effect of Dirac fields on strong cosmic censorship in Kerr–de Sitter background. From an astrophysical perspective, this situation is extremely important since Kerr–de Sitter solution describes the most realistic black holes in an expanding universe.

The paper is organized as follows: in Sect. 2, we constitute weak solutions of Einstein’s equations in presence of Dirac field and discuss the effect of these solutions on strong cosmic censorship conjecture which we use to find the the criteria for strong cosmic censorship violation in Sect. 3. In Sect. 4, we present the main results of the paper. Finally, with some relevant discussions in last section, we conclude our paper. Throughout the paper, we use units in which \(G = c = \hslash = 1\).

2 Weak solution of Einstein equation in presence of massless Dirac fields

The fate of strong cosmic censorship conjecture relies on the possibility of finding a solution of Einstein equation at the Cauchy horizon. Even if the metric is not differentiable (but continuous!) at the Cauchy horizon, one can still make sense of Einstein equation there by constituting a weak solution of the equation [14]. This can be understood by considering the effect of linear perturbation on the spacetime. Consider a massless Dirac field which satisfies the equation \(\widehat{\mathscr {D}}\varPsi =0\), triggers a perturbation in the spacetime. Here \(\widehat{\mathscr {D}}\) is the Dirac operator that acts on the spinor \(\varPsi \). Let the Dirac field act as a first order perturbation which induces a second order perturbation of the metric, denoted by \(h_{\mu \nu }^{(2)}\) which satisfies the following equation

$$\begin{aligned} \widehat{\mathscr {O}}h_{\mu \nu }^{(2)}=8\pi ~T_{\mu \nu }^{\varPsi } \end{aligned}$$
(1)

where, \(\widehat{\mathscr {O}}\) is a second order differential operator and \(T_{\mu \nu }^{\varPsi }\) is the stress-energy tensor for the Dirac field which can be expressed in the following form

$$\begin{aligned} \begin{aligned} T_{\mu \nu }^{\varPsi }&=\frac{i}{2}~\left[ \overline{\varPsi } \gamma _{(\mu }\nabla _{\nu )}\varPsi -\nabla _{(\mu }\overline{\varPsi }\gamma _{\nu )}\varPsi \right] \\&\quad -\frac{i}{2}~g_{\mu \nu }\left[ \overline{\varPsi }\gamma ^{\lambda } \nabla _{\lambda }\varPsi -\nabla _{\lambda }\overline{\varPsi }\gamma ^{\lambda }\varPsi \right] . \end{aligned} \end{aligned}$$

Even when \(h_{\mu \nu }^{(2)}\) is not differentiable at Cauchy horizon, we can have a solution of Einstein equation by multiplying Eq. (1) with a smooth, symmetric tensor \(K^{\mu \nu }\). By performing integration by parts, we obtain the following equation

$$\begin{aligned} \int _{\mathscr {M}}~d^{4}x~\sqrt{-g}\left( h_{\mu \nu }^{(2)}~\mathscr {L}^{\dagger }~K^{\mu \nu }\right) =\int _{\mathscr {M}}~d^{4}x~\sqrt{-g}\left( K^{\mu \nu }~T_{\mu \nu }^{\varPsi }\right) \end{aligned}$$
(2)

where, \(\mathscr {M}\) is the maximal Cauchy development i.e. the region of spacetime uniquely determined by a set of generic initial data [10, 15] and \(\mathscr {L}^{\dagger }\) is the adjoint of the operator \(\mathscr {L}\). If this equation is satisfied for any smooth, symmetric function \(K^{\mu \nu }\), there exists a weak solution of Einstein equation provided that both sides of the equation remains finite. The requirement is fulfilled when \(\varPsi \) belongs to \(H^{1}_{loc}\). When such solutions exists at the Cauchy horizon, we can extend the metric across it which leads to the breakdown of strong cosmic censorship conjecture.

3 Dirac equation in Kerr–de Sitter spacetime and the criteria of violation of strong cosmic censorship conjecture

In the previous section, we have found that the presence of Dirac fields in Kerr–de Sitter spacetimes can lead to the violation of strong cosmic censorship, if the spinor field \(\varPsi \) belongs to \(H^{1}_{loc}\) at Cauchy horizon. In this section, we closely inspect this condition. We solve the Dirac equation near the Cauchy horizon of Kerr–de Sitter black holes and then rewrite the condition in terms of black hole parameters. We start with a Kerr–de Sitter black hole spacetime in Boyer–Lindquist coordinate \((t,r,\theta ,\phi )\) whose line element can be expressed as follows [28]

$$\begin{aligned} ds^{2}= & {} -\frac{\varDelta _{r}}{(1+\alpha )^{2}\Sigma }\left[ dt-a \sin ^{2}\theta ~d\phi \right] ^{2}+\Sigma \left[ \frac{dr^{2}}{\varDelta _{r}}+\frac{d\theta ^{2}}{\varDelta _{\theta }}\right] \nonumber \\&+\frac{\varDelta _{\theta }\sin ^{2}\theta }{(1+\alpha )^{2}\Sigma }\left[ a dt-(r^{2}+a^{2})~d\phi \right] ^{2} \end{aligned}$$
(3)

where,

$$\begin{aligned}&{\varDelta _{r}}(r)=(r^{2}+a^{2})(1-\frac{\Lambda r^{2}}{3})-2M r,\quad {\varDelta _{\theta }}(\theta )=1+\frac{\Lambda a^{2}}{3}\cos ^{2}\theta , \\&{\alpha }=\frac{\Lambda a^{2}}{3},\quad {\varrho }=r+i a\cos \theta \ ,\quad {\Sigma }=\varrho \varrho ^{*} \end{aligned}$$

here, M is the mass of the black hole, a is the black hole rotation parameter and \(\Lambda >0\) is the cosmological constant. Throughout the paper, superscript ‘\(*\)’ denotes the complex conjugate of a quantity. Since we want to check the validity of strong cosmic censorship conjecture in presence of positive cosmological constant, we choose the values of the black hole parameters M, a and \(\Lambda \) in such a way that the spacetime possesses three distinct horizons. The position of Cauchy, event and cosmological horizon which we denote by \(r_{-}\), \(r_{+}\) and \(r_{c}\), respectively, can be found by solving the equation \(\varDelta _{r}(r)=0\). The properties of Dirac fields in this spacetime can be best understood in the framework of Newman–Penrose formalism [29, 30]. Here, we choose our null tetrad to be [28, 31, 32]

$$\begin{aligned} \begin{aligned} l^{\mu }&=\Bigg [\dfrac{(1+\alpha )(r^2+a^2)}{\varDelta _{r}},1,0,\frac{(1+\alpha )a}{\varDelta _{r}}\Bigg ]\\ n^{\mu }&=\frac{1}{2\Sigma }[(1+\alpha )(r^2+a^2),-\varDelta _{r},0,(1+\alpha )a]\\ m^{\mu }&=\frac{1}{\sqrt{2\varDelta _{\theta }}\varrho }\Bigg [i a (1+\alpha )\sin \theta ,0,\varDelta _{\theta },\frac{i(1+\alpha )}{\sin \theta }\Bigg ]\\ \bar{m}^{\mu }&=m^{*\mu } \end{aligned} \end{aligned}$$
(4)

in the \((t,r,\theta ,\phi )\) coordinate. We can easily verify that the only non-vanishing inner product combination is given by the normalization condition, \(\mathbf{l}\cdot \mathbf{n}=-1\) and \(\mathbf{m}\cdot \bar{\mathbf{m}}=1\). The advantage of this choice is that the tetrad vectors are regular across the Cauchy horizon. To see that we write the tetrad in outgoing Eddington–Finkelstein coordinate \((u,r,\theta ,\varphi )\) by using the following transformation [33]

$$\begin{aligned}&{du}=dt-\frac{(1+\alpha )(r^{2}+a^{2})}{\varDelta _{r}}dr,\quad {d\varphi }=d\phi -\frac{(1+\alpha )a}{\varDelta _{r}}dr \end{aligned}$$

Under this transformation the tetrad vectors take the following form

$$\begin{aligned} \begin{aligned} l^{\mu }&=[0,1,0,0]\\ n^{\mu }&=\frac{1}{\Sigma }\Bigg [(1+\alpha )(r^2+a^2),\frac{-\varDelta _{r}}{2},0,(1+\alpha )a\Bigg ]\\ m^{\mu }&=\frac{1}{\sqrt{2\varDelta _{\theta }}\varrho }\Bigg [i a (1+\alpha )\sin \theta ,0,\varDelta _{\theta },\frac{i(1+\alpha )}{\sin \theta }\Bigg ]\\ \bar{m}^{\mu }&=m^{*\mu }.\\ \end{aligned} \end{aligned}$$
(5)

Note that, the tetrad vectors are regular at the Cauchy horizon.

In Newman–Penrose formalism, the equation for a massless Dirac field \(\varPsi \) can be written as four coupled differential equations as follows [29]

$$\begin{aligned} \begin{aligned} (D+\epsilon -\rho )F_{1}+(\bar{\delta }+\pi -\alpha )F_{2}&=0\\ (\varDelta +\mu -\gamma )F_{2}+(\delta +\beta -\tau )F_{1}&=0\\ (D+\epsilon ^{*}-\rho ^{*})G_{2}-(\delta +\pi ^{*}-\alpha ^{*})G_{1}&=0\\ (\varDelta +\mu ^{*}-\gamma ^{*})G_{1}-(\bar{\delta }+\beta ^{*}-\tau ^{*})G_{2}&=0.\\ \end{aligned} \end{aligned}$$
(6)

where, \(D=l\cdot \nabla \), \(\varDelta =n\cdot \nabla \), \(\delta =m\cdot \nabla \), \(\bar{\delta }=\bar{m}\cdot \nabla \) are the directional covariant derivative along the tetrad vectors and \(\alpha \), \(\beta \), \(\gamma \), \(\epsilon \), \(\pi \), \(\rho \) and \(\tau \) are the spin coefficients (for details see [29, 30]). Here, \(F_{1}\), \(F_{2}\), \(G_1\), \(G_2\) denote the spinor components such that \(\varPsi =(F_1,F_2,-G_2,G_1)^{T}\). For our choice of tetrad (given by Eq. (4)), the non-vanishing spin coefficients are given by [31, 32]

$$\begin{aligned} {\rho }= & {} -\frac{1}{\varrho ^{*}} ,\qquad {\tau }=-\frac{i~a~\sqrt{\widetilde{\varDelta _{\theta }}}}{\sqrt{2}\varrho ^{2}} ,\qquad {\pi }=\frac{i~a~\sqrt{\widetilde{\varDelta _{\theta }}}}{\sqrt{2}(\varrho ^{*})^{2}},\\ {\mu }= & {} -\frac{\varDelta _{r}}{2\Sigma \varrho ^{*}} ,\qquad {\gamma }=\mu +\frac{1}{4\Sigma }\frac{d\varDelta _{r}}{dr},\\ {\beta }= & {} \frac{1}{2\sqrt{2}\rho \sin \theta }\frac{d\sqrt{\widetilde{\varDelta }_{\theta }}}{d\theta } ,\qquad {\alpha }=\pi -\beta ^{*}. \end{aligned}$$

where, \(\widetilde{\varDelta _{\theta }}=\varDelta _{\theta }\sin ^{2}\theta \). Due to presence of timelike and angular Killing vectors, the Dirac field can be decomposed as \(\varPsi (t,r,\theta ,\phi )=e^{-i\omega t}e^{i m \phi }(f_{1}, f_{2},-g_{2},g_{1})^{T}\), where \(f_{1}\), \(f_{2}\), \(g_1\), \(g_2\) are functions of r and \(\theta \) only. Moreover, if we take the following transformation [29, 32]

$$\begin{aligned} {f_{1}(r,\theta )}= & {} \frac{R_{-}(r)~S_{-}(\theta )}{\varrho ^{*}},\qquad {f_{2}(r,\theta )}=R_{+}(r)~S_{+}(\theta )\\ {g_{1}(r,\theta )}= & {} R_{+}(r)~S_{-}(\theta ),\qquad {g_{2}(r,\theta )}=\frac{R_{-}(r)~S_{+}(\theta )}{\varrho } \end{aligned}$$

Equation (6) can be decomposed into radial and angular parts which can be written as follows [31, 32]

$$\begin{aligned}&{\mathscr {D}_{-}R_{-}}=\lambda \sqrt{\varDelta _{r}}R_{+},\qquad {\mathscr {D}_{+}\sqrt{\varDelta _{r}}R_{+}}=\lambda R_{-} \end{aligned}$$
(7)
$$\begin{aligned}&{\sqrt{\varDelta _{\theta }}\mathscr {L}_{-}S_{-}}=\lambda S_{+},\qquad {\sqrt{\varDelta _{\theta }}\mathscr {L}_{+}S_{+}}=-\lambda S_{+} \end{aligned}$$
(8)

where,

$$\begin{aligned} \begin{aligned} \mathscr {D}_{\pm }&=\left( \sqrt{\varDelta _{r}}\partial _{r}{\mp }i\frac{V(r)}{\sqrt{\varDelta _{r}}}\right) \\ \mathscr {L}_{\pm }&=\partial _{\theta }{\mp }\frac{1+\alpha }{\varDelta _{\theta }} H(\theta )+\frac{1}{2\sqrt{\widetilde{\varDelta _{\theta }}}}\dfrac{d\sqrt{\widetilde{\varDelta _{\theta }}}}{d\theta }\\ \end{aligned} \end{aligned}$$

here, \(\lambda \) is a constant of separation, \(V(r)=(1+\alpha )[\omega (r^{2}+a^{2})-a m]\) and \(H(\theta )=(a \omega \sin \theta -m\csc \theta )\).

Since we are interested in the behavior of the Dirac field near the Cauchy horizon, it is convenient to adopt the outgoing Eddington–Finkelstein coordinate \((u,r,\theta ,\varphi )\) [33]. Moreover, if we set \(\widetilde{R}_{+}(r)=(\varDelta _{r})^{\frac{1}{2}}R_{+}(r)\) and \(\widetilde{R}_{-}(r)=R_{-}(r)\), Eq. (7) can be written in a more symmetric form \(\mathscr {D}_{\pm }\widetilde{R}_{\pm }=\lambda \widetilde{R}_{{\mp }}\). Near the Cauchy horizon, the radial equations can be written as follows

$$\begin{aligned} \dfrac{d^{2}\widehat{R}_{\pm }}{dr_{*}^{2}}+2 i(\omega -m\Omega _{-})\dfrac{d\widehat{R}_{\pm }}{dr_{*}}=0 \end{aligned}$$
(9)

where, \(\widehat{R}_{\pm }=\widetilde{R}_{\pm }\exp [-i(\omega r_{*}-m r_{\phi })]\) and \(\Omega _{-}\) corresponds to the angular velocity of the black hole at the Cauchy horizon. Here, \(dr_{*}=(1+\alpha )(r^{2}+a^2)dr/\varDelta _{r}\) and \(dr_{\phi }=(1+\alpha )a~dr/\varDelta _{r}\). The above equations have two independent solutions such that the spinor field can be written as follows

$$\begin{aligned} \begin{aligned} \widehat{\varPsi }^{(1)}(u,r,\theta ,\varphi )&=e^{-i(\omega u- m\varphi )}~\psi (r,\theta )\\ \widehat{\varPsi }^{(2)}(u,r,\theta ,\varphi )&=e^{-i(\omega u- m\varphi )}~\psi (r,\theta )~(r-r_{-})^{p} \end{aligned} \end{aligned}$$
(10)

where, \(\psi =(\frac{\mathscr {R}_{-}~S_{-}}{\varrho ^{*}},\mathscr {R}_{+}~S_{+},-\frac{\mathscr {R}_{-}~S_{+}}{\varrho },\mathscr {R}_{+}~S_{-})^{T}\). Here, \(\mathscr {R}_{\pm }(r)\) represent some smooth functions of r which are non-vanishing at the Cauchy horizon and \(p=i(\omega -m \Omega _{-})/\kappa _{-}\) where \(\kappa _{-}\) corresponds to the surface gravity of the black hole at the Cauchy horizon. Given the solution of Dirac equation near the Cauchy horizon, we need to check whether the Dirac field belongs to \(H^{1}_{loc}\) or not in order to investigate the possibility of having weak solutions of Einstein equations. In other words, we have to check whether \(\partial _{\mu } \varPsi \) is locally square integrable which boils down to investigate the finiteness of the integral of quantity \(\sim (r-r_{-})^{2(p-1)}\). Thus, the condition for \(\varPsi \) to remain in \(H^{1}_{loc}\) then reduces to the following inequality

$$\begin{aligned} \beta \equiv -\frac{\mathbb {I}{\text {m}}(\omega )}{\kappa _{-}}>\frac{1}{2} \end{aligned}$$
(11)

The existence of weak solutions of Einstein’s equation at Cauchy horizon is guaranteed by this condition which leads to a possible violation of strong cosmic censorship conjecture. In our study, we need to focus on the dominant mode contributions only which are the least damped modes of the quasinormal spectrum.

Fig. 1
figure 1

The variation of imaginary part of quasinormal frequency \(-\mathbb {I}{\text {m}}(\omega )\) as a function of \((a/a_\mathrm{max})\) has been presented for different values of \(\Lambda M^{2}\) and angular eigen mode number l. Here, \(a_\mathrm{max}\) denotes the extremal value of the rotation parameter. As evident from each of these plots, the lowest lying modes always correspond to the mode with \(m=l\) regardless of the the values of \(\Lambda M^{2}\)

Full size image

4 Strong cosmic censorship conjecture for Kerr–de Sitter black holes in presence of Dirac field

In this section, we compute the quasi-normal modes of Kerr–de Sitter black hole in presence of Dirac field using numerical methods. But before that let us put the radial and angular equations given by Eqs. (7) and (8) in standard Teukolsky form which read as follows [28]

$$\begin{aligned}&\varDelta _r^{-s}\partial _{r} \bigg [\varDelta _r^{s+1}\partial _{r}R_{+}(r)\bigg ]+\bigg [4 i (\alpha +1) r s \omega +2 (1-\alpha ) s\nonumber \\&\quad -\frac{2 \alpha r^2 (s+1) (2 s+1)}{a^2}+\frac{V(V-i s \partial _{r}\varDelta _{r})}{\varDelta _r}-\lambda ^2\bigg ]R_{+}(r)=0 \end{aligned}$$
(12)

and

$$\begin{aligned}&\Bigg [\partial _{x}(1+\alpha x^{2})(1-x^{2})\partial _{x}+\lambda ^{2}-s(1-\alpha )+\frac{(1+\alpha )^{2}}{\alpha }\xi ^{2}\nonumber \\&\quad -2\alpha x^{2}+\frac{1+\alpha }{1+\alpha x^{2}}\Bigg \{2s\Big (\alpha m-(1+\alpha )\xi \Big )x-\frac{(1+\alpha )^{2}}{\alpha }\xi ^{2}\nonumber \\&\quad -2m(1+\alpha )\xi +s^{2}\Bigg \}-\frac{(1+\alpha )^{2}m^{2}}{(1+\alpha x^{2})(1-x^{2})}\nonumber \\&\quad -\frac{(1+\alpha )(s^{2}+2s m x)}{1-x^{2}}\Bigg ]S_{+}(x)=0 \end{aligned}$$
(13)

where, \(x=\cos \theta \), \(s=1/2\) and \(\xi =a\omega \). Note that, in the non-rotating limit \((a\rightarrow 0)\), the angular Teukolsky equation (13) gives the separation constant as \(\lambda \rightarrow l(l-1)-s^{2}+s\). This can be used to define the angular eigen mode number l which satisfies the following relation, \(l\ge \mathrm {max}(|m|,|s|)\) [34]. These transformed equations allow us to use the method developed by [28], who showed that Eq. (13) can be transformed into Heun’s equation which, in turn, give us a three-term recurrence relation for the angular equation. This three term recurrence relation can be rewritten in terms of a continued fraction equation which we denote by \(P_{1}(\lambda ,\omega )=0\).

Fig. 2
figure 2

The variation of \(\mathbb {I}{\text {m}}(\omega )/\kappa _{-}\) as a function of \((a/a_\mathrm{max})\) has been presented for different values of \(\Lambda M^{2}\). The parameter \(\beta \) corresponds to the least damped modes of the quasi-normal spectrum. Here, we have considered only the modes with \(m=l\). In each of these plots, red, green and blue curves corresponds to value of \( \mathbb {I}{\text {m}}(\omega )/\kappa _{-}\) for \(l=0.5\), \(l=1.5\) and \(l=10.5\) respectively. The vertical lines in the plots presented in upper-left and upper-right corner corresponds to the value of rotation parameter a for which strong cosmic censorship conjecture gets violated. As evident from the plots, for smaller values of \(\Lambda M^{2}\), we can find a parameter space which does not respect the conjecture

Full size image

The radial Teukolsky equation (12) has five regular singularities at \(r_+\), \(r_-\), \(r_{c}\), \(-(r_{+}+r_{-}+r_{c})\) and spatial infinity. The quasi-normal modes are defined as the eigen values of \(\omega \) with \(\varPsi \) satisfying the following boundary condition: there are only outgoing waves at the cosmological horizon \(r_{c}\) and only ingoing waves at the event horizon \(r_{+}\). In order to satisfy the boundary condition, we write \(R_{+}(r)\) as the multiplication of a function y(z) which is regular at the boundary and a factor which is divergent at \(r_{+}\) and \(r_{c}\), i.e.,

$$\begin{aligned} R_{+}(r)=r^{-(2s+1)}\bigg (\frac{r-r_{-}}{r-r_{+}}\bigg )^{s+2i \frac{V(r_{+})}{\varDelta _{r}'(r_{+})}}e^{iB(r)}~y(z) \end{aligned}$$
(14)

where, \(dB/dr=V(r)/\varDelta _{r}\) and \(z=(r-r_{+})(r_{c}-r_{-})/(r-r_{-})(r_{c}-r_{+})\). By inserting Eq. (14) into Eq. (12) and expressing y(z) as a Frobenious series of the form \(\sum _{n=0}^{\infty }a_{n}z^{n}\), we get a seven-term recurrence relation. Due to convergence problem of seven term recurrence relation [35], it is better to reduce the seven-term recurrence relation into a three term recurrence relation using Gaussian elimination method [36, 37]. Similar to the angular equation, this three term recurrence relation can also be written as an infinite continued fraction equation which we denote by \(P_{2}(\lambda ,\omega )=0\). By solving the angular and radial continued fraction equations simultaneously, one gets the desired quasi-normal modes. However, instead of using the Gaussian elimination procedure to find the three-term recurrence relation for the radial equation, we have employed the Mathematica package developed by Jansen [38] to find the quasinormal modes. Here, we have solved the angular equation \(P_{1}(\lambda ,\omega )=0\) and radial equation for y(z) iteratively, taking the approximate value of \(\lambda \) in Ref. [28] as the initial guess value.

In order to understand the effect of Dirac particles on strong cosmic censorship conjecture, we need to look for modes for which the parameter \(\beta \equiv -\mathbb {I}{\text {m}}(\omega )/\kappa _{-}\) becomes greater than 1/2. As stated earlier, we are interested in the dominant modes (least damped modes) of the quasinormal spectrum for a given value of angular eigen mode number l for the calculation of the parameter \(\beta \). From Fig. 1, we can see that the dominant mode is always corresponds to modes with \(m=l\) for any values of mass scaled cosmological constant \(\Lambda M^{2}\) [14, 39]. Moreover, the imaginary part of the quasinormal frequency decreases with the increase of both the rotational parameter a and the cosmological constant \(\Lambda \). In this regard, our result is fully consistent with the results presented in [39]. Moreover, near to the extremity, a rapid decrease of the value of \(\mathbb {I}{\text {m}}(\omega )\) is observed.

The variation of the quantity \(\mathbb {I}{\text {m}}(\omega )/\kappa _{-}\) with respect to the rotational parameter a for different values of \(\Lambda M^{2}\) is presented in Fig. 2. Here, we consider only the least damped modes for a given value of angular eigen mode number i.e. the modes with \(m=l\). It is interesting to see that for smaller values of mass scaled cosmological constant (\(\Lambda M^{2}\approx \mathscr {O}(10^{-3})\)), modes corresponding to \(l=1/2\) dominate the quasinormal spectrum for certain values of a. However, for larger values of \(\Lambda M^{2}\), the eikonal modes (corresponding to the modes with large l value) become dominant. As evident from Fig. 2, strong cosmic censorship conjecture is respected for larger values of \(\Lambda M^{2}\). But for smaller values of mass scaled cosmological constant (\(\Lambda M^{2}\lessapprox 0.02\)), we are able to find a parameter space for which strong cosmic censorship gets violated. Moreover, in this parameter space, the violation is more severe for smaller values of \(\Lambda M^{2}\) since \(\beta \) becomes greater than 1/2 for smaller values of rotational parameter a. So the assurance that strong cosmic censorship is always respected in astrophysical black holes [14, 24] is no longer valid in the presence of Dirac fields.

5 Conclusion

In recent years, several examples are found which suggest a breakdown of determinism in Reissner–Nordström–de Sitter black holes under the influence of several fundamental fields [18, 19, 21,22,23,24,25,26]. Anyway, astrophysically meaningful Kerr–de Sitter black hole solutions seems to respect the conjecture [14, 24]. However, the previous analyses were done considering only the effect of scalar fields. In this paper, we extended the study by considering Dirac fields in Kerr–de Sitter background.

By considering the effect of linear perturbations on the spacetime metric of our interest, we found that there exist weak solutions of Einstein equation at the Cauchy horizon, if the parameter \(\beta \) defined by Eq. (11) become greater than 1/2 which leads to the violation of strong cosmic censorship conjecture. Comparing our result with Ref. [14], where the effect of scalar fields are considered, we see that the criteria for strong cosmic censorship violation remains the same in the presence of Dirac fields also. Once this criteria is obtained, we performed detailed numerical computation to find lowest lying quasinormal modes to determine the value of \(\beta \) as a function of rotational parameter a for different values of mass scaled cosmological constant \(\Lambda M^{2}\). In Figs. 1 and 2, we presented our main results. From these figures, it is clear that the lowest lying quasinormal modes always correspond to the \(m=l\) mode. Moreover, an increased value of rotational parameter results in modes with smaller decay rate. Near to the extremity, a rapid decrease of decay rate is observed. The decay rate also decreases with the increase of \(\Lambda M^{2}\). From Fig. 2, it is clear that the value of the parameter \(\beta \) always remains smaller than 1/2 for larger values of \(\Lambda M^{2}\). Hence, we can conclude that the strong cosmic censorship conjecture is always respected in “large” rotating black holes (black holes with \(\Lambda M^{2}\gtrapprox 0.02\)). However, for smaller values of \(\Lambda M^{2}\), we found a parameter space where \(\beta \) becomes larger than 1/2. This, in turn, implies that in “smaller” black holes, Dirac fields can be smoothly extended beyond Cauchy horizon. Hence, in presence of Dirac fields, even the rotational parameter can not save the strong cosmic censorship conjecture for a certain parameter range. Moreover, for a fixed value of cosmological constant, black holes with smaller masses are more prone to violate the conjecture.

Note that, in our work, we have considered the effect of linear perturbation only. It may be possible that the violation of strong cosmic censorship can be prevented by considering the non-linear or quantum gravitational effects. However, Cardoso et al. have shown that for Einstein–Maxwell–scalar field system, the non-linear perturbations can not save the conjecture [25]. Moreover, a recent study showed that the violation of this conjecture in \(2+1\) dimensional BTZ black holes even when quantum correction terms are added [40]. However, the studies in \(3+1\) dimensional black holes remain inconclusive. In particular, it has shown that the quantum correction can either amplify or suppress the blue shift instability in Reissner–Nordström black holes [41]. It would be interesting to see the those effects on Einstein–Dirac field system, which we leave for future work. Recently, Dafermos and Shlapentokh-Rothman [15] suggested an interesting proposal that the strong cosmic censorship can still be saved if one starts with rough initial data. This idea is further supported by Ref. [42], where the authors have studied a coupled gravitational and electro-magnetic perturbation and showed that in order to save the strong cosmic censorship, one must consider physically reasonable but slightly less smooth initial data. It will be interesting to see whether this proposal can save the conjecture in presence of Dirac field also. However, this is beyond the scope of this paper.