1. Introduction
It is well known that quasibound states are among the most important and fundamental phenomena of physics. There has been much attention devoted to revealing characteristics of quasibound states in various black hole backgrounds in recent years. A large number of methods have been proposed. Such as the WKB method [1–3], the continued-fraction method [1, 4–6], numerical method [1], the matrix method [7–10], the pseudospectral method [11, 12], the asymptotic method [13–15], and the Heun-function method [16–20]. Despite extensive discussions, quasibound state problems in an equivalent medium have not been studied to date. However, this problem is important because analog spacetimes can provide a useful testing tool.
The idea of analog spacetimes can be traced back almost a century ago. In 1920, Eddington [21] pointed out that the gravitational effect on light is equivalent to that produced by a suitable distribution of refractive media. In 1923, Gordon [22] studied the inverse problem and introduced a notion of effective metric to describe the effect of a refractive index on the propagation of light. The metric was to become known as the Gordon metric. The first model on analog spacetime experiment was proposed in 1981 by Unruh [23]. Since then a vast number of research papers on analog gravity have been published, for instance in acoustics [24–29], electromagnetic wave guides [30], Bose–Einstein condensates [31–34], liquid helium [35, 36], and graphene [37–39], etc. The developments of this field before 2011 have been introduced in the excellent review paper [40]. In recent years, especially with the technology breakthrough of transformation optics, the photonic chips were introduced as analogical gravitational systems. From which the simulation experiments have been realized on the gravitational lens [41], Einstein ring [42], topological space of cosmic string [43], and accelerating particles in curved spacetime [44]. These studies show that the analog gravity systems provide a new experimental platform for studying the nature of spacetime. Unlike electromagnetic waves, gravitational waves cannot be controlled with current technology. A natural question is whether researchers can mimic gravitational waves in laboratories. To our best knowledge, no experiments or theories of analog gravitational waves were ever given for strong gravitational background.
Analogy is one of the basic thinking methods in the process of understanding objective things, its physical basis is that different physical systems obey the same dynamical evolution equations. For instance, base on the equivalence between Klein–Gordon equation for the near-horizon Schwarzschild metric and the motion equation of sound waves in a convergent fluid flow, Unruh [23] established the acoustic analog of a black hole. Therefore, similarity analyses provide cross-fertilization of ideas among different branches of science, and are the theoretical basis of analogical experiment research.
Although optical media that mimic Schwarzschild spacetime have been successfully fabricated [41], none of the currently published research results on quasibound states of Schwarzschild black holes can be experimentally verified in terrestrial laboratories using these media. This is because these studies did not use the isotropic coordinate system which is crucial for fabricating the optical medium that mimics Schwarzschild spacetime. It is important to note that using different coordinate systems means examining the quasibound states of particles in Schwarzschild spacetime from the perspective of different observers. In this paper, we investigate quasibound states for massless scalar (spin s = 0), Weyl neutrino (spin s = 1/2), electromagnetic (spin s = 1), massless Rarita–Schwinger (spin s = 3/2), and gravitational (spin s = 2) waves in Schwarzschild spacetime described using the isotropic coordinate system. This spacetime can be considered an equivalent medium analogous to Schwarzschild spacetime, which we refer to as the Schwarzschild equivalent medium. We expect that our study not only solves the quasibound state problem in Schwarzschild equivalent mediums but also gives the fundamental theory for the simulation of gravitational waves in laboratories. Throughout our discussion, the particles associated the above mentioned waves are collectively called the massless spin particles for short.
This paper is organized as follows. In section 2 , we study the nature of Schwarzschild equivalent mediums. In section 3 , we derive the master equation governing massless spin particles in Schwarzschild equivalent mediums, and show that the angular functions of the solutions are the spin-weighted spheroidal harmonics. In section 4 , we study the radial equation in Schwarzschild equivalent mediums, and discuss the form of solutions near and far the event horizon. In section 5 , we obtain remarkably simple formulae for the frequencies of the quasibound States, and show that the quasibound frequencies of the massless spin particles for same statistical properties are identical at the same quantum numbers N and l. Finally, section 6 contains a further discussion and some concluding remarks.
2. Schwarzschild equivalent mediums
In isotropic coordinates, the Schwarzschild metric is given by [45]
$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & \frac{{(1-M/2r)}^{2}}{{(1+M/2r)}^{2}}{\rm{d}}{t}^{2}-{\left(1+\frac{M}{2r}\right)}^{4}\\ & & \times \,({\rm{d}}{r}^{2}+{r}^{2}{\rm{d}}{\theta }^{2}+{r}^{2}{\sin }^{2}\theta {\rm{d}}{\varphi }^{2}),\end{array}\end{eqnarray}$
where M is the mass of the black hole. The event horizon in isotropic coordinates is located at rH = M/2. It is well known that, for the Schwarzschild metric in the standard form, the most obvious characteristic at the event horizon is the reversal there of the roles of t and r as timelike and spacelike coordinates. However, when the isotropic coordinates are adopted, this argument fails. Metric (1) shows that the t direction is timelike and the r direction is spacelike, both in the region r > M/2 and in the region r < M/2. Actually, the metric (1) does not cover the interior region of the Schwarzschild black hole; it covers the exterior region twice [46].In 1923, Gordon [22] introduced the metric describing an equivalent medium, which is
$\begin{eqnarray}{g}_{\mu \nu }={\eta }_{\mu \nu }+\left(\frac{1}{{n}^{2}}-1\right){u}_{\mu }{u}_{\nu },\end{eqnarray}$
where ημν is the flat metric, n is the refractive index, and uμ is the 4-velocity of the medium. This is the famous Gordon metric. In original Gordon metric, n and uμ are constants. To simulate a real spacetime, the following generalization of the Gordon is widely used [47] $\begin{eqnarray}{g}_{\mu \nu }={{\rm{\Omega }}}^{2}\left[{\eta }_{\mu \nu }+\left(\frac{1}{{n}^{2}}-1\right){u}_{\mu }{u}_{\nu }\right].\end{eqnarray}$
Here Ω is called the conformal factor, and the quantities Ω, n, and uμ are allowed to be space and time-dependent.In the rest frame of the medium uμ = (1, 0, 0, 0), equation (3 ) reads
$\begin{eqnarray}{g}_{\mu \nu }={{\rm{\Omega }}}^{2}\left[{\eta }_{\mu \nu }+\left(\frac{1}{{n}^{2}}-1\right){\delta }_{\mu }^{0}{\delta }_{\nu }^{0}\right].\end{eqnarray}$
In spherical polar coordinates r, θ, φ, the line element of a conformal Gordon spacetime can be written as $\begin{eqnarray}{\rm{d}}{s}^{2}={{\rm{\Omega }}}^{2}(r)\left[\frac{1}{{n}^{2}(r)}{\rm{d}}{t}^{2}-({\rm{d}}{r}^{2}+{r}^{2}{\rm{d}}{\theta }^{2}+{r}^{2}{\sin }^{2}\theta {\rm{d}}{\varphi }^{2})\right].\end{eqnarray}$
The comparison of equation (5 ) with equation (1 ) shows that, for the Schwarzschild metric, the conformal factor and the index of refraction take the forms
$\begin{eqnarray}{\rm{\Omega }}={\left(1+\frac{M}{2r}\right)}^{2},\end{eqnarray}$
and $\begin{eqnarray}n={\left(1+\frac{M}{2r}\right)}^{3}{\left(1-\frac{M}{2r}\right)}^{-1}.\end{eqnarray}$
Metric (5 ) with equations (6 ) and (7 ) is also called the Schwarzschild equivalent medium. This metric is important because it has been realized in an optical medium, which portray a good analogical gravitational system. This means that if a spin particle has Maxwell electromagnetic wavelike properties in Schwarzschild spacetime, then propagation of electromagnetic waves in the analogical gravitational system can be used to simulate the properties this spin particle. This is a main motivations for our work.
Starting from the calculation of the spin coefficients [see equations (A1 ) and (A3 )] we can prove that a Schwarzschild equivalent medium is of Petrov type D by Goldberg-Sachs theorem [48]. The values of the non-vanishing Weyl scalar ψ2 and the Ricci scalar R [for definition see equations (A4 ) and (A5 )] can be given by
$\begin{eqnarray}{\psi }_{2}=\frac{1}{6{{\rm{\Omega }}}^{2}}\left[-\frac{{n}^{{\prime\prime} }}{n}+2{\left(\frac{{n}^{{\prime} }}{n}\right)}^{2}+\frac{1}{r}\frac{{n}^{{\prime} }}{n}\right],\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}R & = & \frac{2}{{{\rm{\Omega }}}^{2}}\left[3\left(-\frac{{{\rm{\Omega }}}^{{\prime\prime} }}{{\rm{\Omega }}}+\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}}\frac{{n}^{{\prime} }}{n}-\frac{2}{r}\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}}\right)+\frac{{n}^{{\prime\prime} }}{n}-2{\left(\frac{{n}^{{\prime} }}{n}\right)}^{2}\right.\\ & & \left.+\,\frac{2}{r}\frac{{n}^{{\prime} }}{n}\right],\end{array}\end{eqnarray}$
where the prime denotes the derivative with respect to r.3. Unified wave equation for massless spin particles
Now we will consider the perturbation equation for massless fields of various spin s ≤ 2 on the Schwarzschild equivalent medium background. In general, the field equations of spin 1/2 1, 3/2, and 2 are not accurately separable, but in all the type-D metrics, the massless field equations of spin 1/2, 1, 3/2, and 2 can be decoupled in the case of perturbations [49–52]. The findings indicate that each spin state of a particle corresponds to a equation, resulting in a total of eight equations, as detailed in Appendix B . To enhance the identification of similarities among their solutions, the solutions must be articulated in a unified form. This necessitates the integration of these equations into a single statement (master equation).
In Appendix B , each pair of equations (B3 )–(B6 ) represents the Weyl neutrino, electromagnetic, massless Rarita–Schwinger, and gravitational perturbations, respectively. In each pair, the first equation corresponds to the spin state p = s, while the second corresponds to p = −s. A careful examination of (B3 )–(B6 ) readily allows the equations describing the spin states p = s to be unified into the following formA1 ) of Appendix A . D, Δ, and δ are the directional derivatives defined by the equations A2 ) and (A3 ). Substituting equations (12 ), (A2 ) and (A3 ) into equations (10 ) and (11 ), and applying the transformation 10 ) and (11 ) into a single master equation, which can be written in the elegant form
$\begin{eqnarray}\begin{array}{l}\{[D-(2s-1)\varepsilon +\bar{\varepsilon }-2s\rho -\bar{\rho }]({\rm{\Delta }}-2s\gamma +\mu )\\ -\,[\delta +\bar{\pi }-\bar{\alpha }-(2s-1)\beta -2s\tau ](\bar{\delta }+\pi -2s\alpha )\\ -\,(2s-1)(s-1){\psi }_{2}\}{{\rm{\Phi }}}_{s}=0,\end{array}\end{eqnarray}$
and those for the spin states p = −s into the following form $\begin{eqnarray}\begin{array}{l}\{[{\rm{\Delta }}+(2s-1)\gamma -\bar{\gamma }+2s\mu +\bar{\mu }](D+2s\varepsilon -\rho )\\ -\,[\bar{\delta }-\bar{\tau }+\bar{\beta }+(2s-1)\alpha +2s\pi ](\delta -\tau +2s\beta )\\ -\,(2s-1)(s-1){\psi }_{2}\}{{\rm{\Phi }}}_{-s}=0.\end{array}\end{eqnarray}$
Here, Φs and Φ−s denote the wave functions corresponding to the particle’s spin states p = s and p = −s, respectively. The Greek letters such as ρ represent the spin coefficients, whose definitions are provided in equation ( $\begin{eqnarray}D={l}^{\mu }{\partial }_{\mu },\quad {\rm{\Delta }}={n}^{\mu }{\partial }_{\mu },\quad \delta ={m}^{\mu }{\partial }_{\mu },\quad \bar{\delta }={\bar{m}}^{\mu }{\partial }_{\mu },\end{eqnarray}$
where lμ, nμ, mμ and ${\bar{m}}^{\mu }$ are the null tetrads. For Schwarzschild equivalent mediums, the null tetrads and spin coefficients are given by equations ( $\begin{eqnarray}{{\rm{\Phi }}}_{p}={({\rm{\Omega }}r)}^{(p-s)}{{\rm{\Psi }}}_{p}.\end{eqnarray}$
We can combine equations ( $\begin{eqnarray}\left[({{\rm{\nabla }}}^{\mu }+p{L}^{\mu })({{\rm{\nabla }}}_{\mu }+p{L}_{\mu })-4{p}^{2}{\psi }_{2}+\frac{1}{6}R\right]{{\rm{\Psi }}}_{p}=0,\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{rcl}{L}^{t} & = & \frac{n}{{{\rm{\Omega }}}^{2}}\left(\frac{{n}^{{\prime} }}{n}+\frac{1}{r}\right),\\ {L}^{r} & = & \frac{1}{{{\rm{\Omega }}}^{2}}\left(2\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}}-\frac{{n}^{{\prime} }}{n}-\frac{1}{r}\right),\\ {L}^{\theta } & = & 0,\\ {L}^{\varphi } & = & -\,\frac{1}{{{\rm{\Omega }}}^{2}{r}^{2}}\frac{{\rm{i}}\cos \theta }{{\sin }^{2}\theta }.\end{array}\end{eqnarray}$
Here, ∇μ denotes the covariant derivative in the metric gμν.Equation (14 ) tells us that the wave functions of different massless spin particles $\Psi$p obey the same dynamical equation (also called wave equation). In view of the completely different nature of these particles, it is surprising (and delightful) that their dynamical equations have such a similar structure in Schwarzschild equivalent mediums. Evidently, when p = 0, equation (14 ) is just the (conformally invariant) massless scalar field equation. Therefore, equation (14 ) governs not only the massless fields of spin 1/2, 1, 3/2, and 2, but also the scalar field. Equation (14 ) is the fundamental result, from which we can study the similarity of the wave functions for different massless spin particles in the Schwarzschild equivalent medium.
The fact that the exact solution of equation (14 ) has the form
$\begin{eqnarray}{{\rm{\Psi }}}_{p}={{\rm{e}}}^{-{\rm{i}}\omega t}{S}_{p}(\theta ,\varphi ){R}_{p}(r).\end{eqnarray}$
The equation for Sp is
$\begin{eqnarray}\begin{array}{l}\left[\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial }{\partial \theta }\right)+\frac{1}{{\sin }^{2}\theta }\frac{{\partial }^{2}}{\partial {\varphi }^{2}}+\frac{2{\rm{i}}p\cos \theta }{{\sin }^{2}\theta }\frac{\partial }{\partial \varphi }\right.\\ \left.-\,{p}^{2}{\cot }^{2}\theta +p+(l-p)(l+p+1)\Space{0ex}{3.3ex}{0ex}\right]{S}_{p}(\theta ,\varphi )=0.\end{array}\end{eqnarray}$
Then the angular part of the wave function is a spin-weighted spherical harmonic [53, 54], namely $\begin{eqnarray}\begin{array}{l}{S}_{p}(\theta ,\varphi ){\,\equiv \,}_{p}{Y}_{lm}(\theta ,\varphi )\\ \,=\,{\left[\frac{(2l+1)}{4\pi }\frac{(l+m)!(l-m)!}{(l+p)!(l-p)!}\right]}^{1/2}{{\rm{e}}}^{{\rm{i}}m\varphi }{\left(\sin \frac{\theta }{2}\right)}^{2l}\\ \,\times \,\displaystyle \sum _{k}\left(\begin{array}{l}l-p\\ \quad k\end{array}\right)\left(\begin{array}{l}\quad l+p\\ k+p-m\end{array}\right){(-1)}^{l-k-p}{\left(\cot \frac{\theta }{2}\right)}^{2k+p-m},\end{array}\end{eqnarray}$
where l, and m are integers satisfying the inequalities l ≥ s, −l ≤ m ≤ l. The spin-weighted spherical harmonics satisfy the orthonormality and completeness relations $\begin{eqnarray}\begin{array}{l}{\displaystyle \int }_{0}^{2\pi }{\rm{d}}\varphi {\displaystyle \int }_{0}^{\pi }{\rm{d}}{\theta \,}_{p}{Y}_{lm}{(\theta ,\varphi )\,}_{p}{\bar{Y}}_{l^{\prime} m^{\prime} }(\theta ,\varphi )\sin \theta =\,{\delta }_{ll^{\prime} }{\delta }_{mm^{\prime} },\end{array}\end{eqnarray}$
and $\begin{eqnarray}\displaystyle \sum _{m=-l}^{l}{| \,}_{p}{Y}_{lm}(\theta ,\varphi ){| }^{2}=\frac{2l+1}{4\pi },\end{eqnarray}$
respectively.It is obvious that the actual shape of the metric functions affects only the radial function Rp(r), which has the form
$\begin{eqnarray}{R}_{p}(r)=A\,{{\rm{\Omega }}}^{2p-1}{n}^{(1-2p)/2}{r}^{-(p+1)}g(r),\end{eqnarray}$
where A is some constant, the functions g(r) has to satisfy $\begin{eqnarray}\begin{array}{l}\frac{{{\rm{d}}}^{2}g(r)}{{\rm{d}}{r}^{2}}+\left[{\omega }^{2}{n}^{2}+2{\rm{i}}\omega pn\left(\frac{{n}^{{\prime} }}{n}+\frac{1}{r}\right)-\frac{1}{3}(2p-1)(2p+1)\right.\\ \left.\,\times \,\left(\frac{1}{2}\frac{{n}^{{\prime\prime} }}{n}-\frac{1}{4}{\left(\frac{{n}^{{\prime} }}{n}\right)}^{2}+\frac{1}{r}\frac{{n}^{{\prime} }}{n}\right)-\frac{l(l+1)}{{r}^{2}}\right]g(r)=0.\end{array}\end{eqnarray}$
This equation shows that the functions g(r) is not affected by the conformal factor Ω(r).Substituting equation (7 ) into equation (22 ), we obtain
$\begin{eqnarray}\begin{array}{l}\frac{{{\rm{d}}}^{2}g(r)}{{\rm{d}}{r}^{2}}+\left\{{\omega }^{2}\frac{{(1+M/2r)}^{6}}{{(1-M/2r)}^{2}}-{\rm{i}}\omega p\left[\frac{2M{(M+2r)}^{2}}{{r}^{2}{(M-2r)}^{2}}-\frac{{(M+2r)}^{2}}{2{r}^{3}}\right]\right.\\ \left.\,-\,(4{p}^{2}-1)\frac{4{M}^{2}}{{({M}^{2}-4{r}^{2})}^{2}}-\frac{l(l+1)}{{r}^{2}}\right\}g(r)=0.\end{array}\end{eqnarray}$
An exact solution of this equation would be very difficult to obtain. However, as we shall see, the most important features of the solution can be found by the uniqueness theorem.4. Radial functions near and far the event horizon
It is clear that the coefficients in equation (23 ) are single-valued analytic functions in the interval (M/2, ∞). The uniqueness theorem tells us that there is only one analytic function satisfying the differential equation and the given boundary conditions. The clear inference is that the asymptotic behavior of the solutions of the equation is therefore dominated by the solutions of the asymptotic equation. For this reason, below we study the asymptotic equation of equation (23 ) and its solutions.
4.1. Solutions in the near zone
In the vicinity of the event horizon, equation (23 ) should have the form
$\begin{eqnarray}\frac{{{\rm{d}}}^{2}g}{{\rm{d}}{r}^{2}}+\left[{\omega }^{2}+\frac{{a}_{1}}{r-M/2}+\frac{{a}_{2}}{{(r-M/2)}^{2}}\right]g=0,\end{eqnarray}$
where $\begin{eqnarray}{a}_{1}=-\frac{2{a}_{2}}{M},\quad {a}_{2}=\frac{1}{4}-{(p+{\rm{i}}4M\omega )}^{2}.\end{eqnarray}$
It should be noted that if we restrict our consideration to the region near the event horizon, we can safely neglect the ω2 term inside the square brackets in equation (24 ). This is because whether or not the ω2 term is present, the form of the solution to equation (24 ) near the outer event horizon remains unchanged—specifically, all solutions reduce to equation (48 ) in section 5 . In other words, retaining the ω2 term does not affect the solution’s behavior near the event horizon.
Note that as r → ∞, the limiting form of equation (24 ) becomes $\frac{{{\rm{d}}}^{2}g}{{\rm{d}}{r}^{2}}+{\omega }^{2}g=0$. This is identical to the equation derived from the equation (23 ) under the same limit. The subtlety in retaining the ω2 term stems from the fact that equation (24 ) serves not only as an approximation of equation (23 ) near the event horizon, but also as its asymptotic form when r → ∞.
If we introduce a new variable 24 ) in the vicinity of the horizon takes the form
$\begin{eqnarray}x=r-\frac{M}{2},\end{eqnarray}$
in terms of x, the radial equation ( $\begin{eqnarray}\frac{{{\rm{d}}}^{2}g}{{\rm{d}}{x}^{2}}+\left[{\omega }^{2}+\frac{{a}_{1}}{x}+\frac{{a}_{2}}{{x}^{2}}\right]g=0.\end{eqnarray}$
Again, we transform to new variable and function. First we put
$\begin{eqnarray}z={\rm{i}}2\omega x,\end{eqnarray}$
and then we define $\begin{eqnarray}g(r)={{\rm{e}}}^{-{\rm{i}}\omega x}{x}^{\frac{1}{2}(1+\sqrt{1-4{a}_{2}})}u.\end{eqnarray}$
Thus equation (27 ) transforms to 30 ) becomes the standard degenerate hypergeometric equation:
$\begin{eqnarray}\begin{array}{l}z\frac{{{\rm{d}}}^{2}u}{{\rm{d}}{z}^{2}}+(1+\sqrt{1-4{a}_{2}}-z)\frac{{\rm{d}}u}{{\rm{d}}z}\\ \,-\,\left[\frac{1}{2}(1+\sqrt{1-4{a}_{2}})+\frac{{\rm{i}}{a}_{1}}{2\omega }\right]u=0.\end{array}\end{eqnarray}$
Clearly, by setting $\begin{eqnarray}\gamma =1+\sqrt{1-4{a}_{2}}=1+2({\rm{i}}4M\omega +p),\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}\alpha & = & \frac{1}{2}(1+\sqrt{1-4{a}_{2}})+\frac{{\rm{i}}{a}_{1}}{2\omega }\\ & = & \frac{1}{2}(1+\sqrt{1-4{a}_{2}})-\frac{{\rm{i}}{a}_{2}}{M\omega }\\ & = & \frac{1}{2}\left[1+2({\rm{i}}4M\omega +p)\right]+\frac{{\rm{i}}}{M\omega }\left[{({\rm{i}}4M\omega +p)}^{2}-\frac{1}{4}\right],\end{array}\end{eqnarray}$
equation ( $\begin{eqnarray}z\frac{{{\rm{d}}}^{2}u}{{\rm{d}}{z}^{2}}+(\gamma -z)\frac{{\rm{d}}u}{{\rm{d}}z}-\alpha u=0,\end{eqnarray}$
As is well known, if γ is not an integer, then the general solution of the degenerate hypergeometric equation has the form $\begin{eqnarray}u={C}_{1}F(\alpha ,\gamma ;z)+{C}_{2}{z}^{1-\gamma }F(\alpha +1-\gamma ,2-\gamma ;z),\end{eqnarray}$
where $\begin{eqnarray}F(\alpha ,\gamma ;z)=\displaystyle \sum _{k=1}^{\infty }\frac{{(\alpha )}_{k}}{k!{(\gamma )}_{k}}{z}^{k}.\end{eqnarray}$
Here F is called the degenerate hypergeometric function or the Kummer function.Therefore, if $1+\sqrt{1-4{a}_{2}}$ is not an integer, then the general solution of equation (27 ) has the form
$\begin{eqnarray}\begin{array}{rcl}g & = & {{\rm{e}}}^{-{\rm{i}}\omega x}{x}^{\frac{1}{2}(1+\sqrt{1-4{a}_{2}})}\left[{C}_{1}F\left(\frac{1}{2}(1+\sqrt{1-4{a}_{2}})\right.\right.\\ & & \left.-\,\frac{{\rm{i}}{a}_{2}}{M\omega },1+\sqrt{1-4{a}_{2}};{\rm{i}}2\omega x\right)\\ & & +\,{C}_{2}{x}^{-\sqrt{1-4{a}_{2}}}F\left(\frac{1}{2}(1-\sqrt{1-4{a}_{2}})\right.\\ & & \left.\left.-\,\frac{{\rm{i}}{a}_{2}}{M\omega },1-\sqrt{1-4{a}_{2}};{\rm{i}}2\omega x\right)\right],\end{array}\end{eqnarray}$
where C1 and C2 are arbitrary constants.4.2. Solutions in the far zone
Our this step is to find the behavior of g(r) at large r. In this case, the differential equation (23 ) becomes
$\begin{eqnarray}\frac{{{\rm{d}}}^{2}g}{{\rm{d}}{r}^{2}}+\left({\omega }^{2}+\frac{{b}_{1}}{r}+\frac{{b}_{2}}{{r}^{2}}\right)g=0,\end{eqnarray}$
where $\begin{eqnarray}{b}_{1}=2\omega (2\omega M+{\rm{i}}p),\quad {b}_{2}=\frac{15}{2}{\omega }^{2}{M}^{2}-l(l+1).\end{eqnarray}$
Using a method similar to that of equation (27 ), we transform equation (37 ) through the following changes37 ) into
$\begin{eqnarray}z={\rm{i}}2\omega r,\end{eqnarray}$
$\begin{eqnarray}g(r)={{\rm{e}}}^{-{\rm{i}}\omega r}{r}^{\frac{1}{2}(1+\sqrt{1-4{b}_{2}})}u.\end{eqnarray}$
This converts equation ( $\begin{eqnarray}\begin{array}{l}z\frac{{{\rm{d}}}^{2}u}{{\rm{d}}{z}^{2}}+(1+\sqrt{1-4{b}_{2}}-z)\frac{{\rm{d}}u}{{\rm{d}}z}\\ -\,\left[\frac{1}{2}(1+\sqrt{1-4{b}_{2}})+\frac{{\rm{i}}{b}_{1}}{2\omega }\right]u=0.\,\end{array}\end{eqnarray}$
If we set $\gamma =1+\sqrt{1-4{b}_{2}}$ and $\alpha =(1\,+\sqrt{1-4{b}_{2}})/2+{\rm{i}}{b}_{1}/(2\omega )$, equation (41 ) becomes the standard degenerate hypergeometric equation (33 ) as well. A natural conclusion is that, if $1+\sqrt{1-4{b}_{2}}$ is not an integer, the general solution of equation (37 ) can be expressed in the form
$\begin{eqnarray}\begin{array}{rcl}g & = & {{\rm{e}}}^{-{\rm{i}}\omega r}{r}^{(1+\sqrt{1-4{b}_{2}})/2}\left[\Space{0ex}{2.5ex}{0ex}{D}_{1}F\left(\Space{0ex}{2.5ex}{0ex}-p+{\rm{i}}2\omega M\right.\right.\\ & & \left.+\,\frac{1}{2}(1+\sqrt{1-4{b}_{2}}),1+\sqrt{1-4{b}_{2}};{\rm{i}}2\omega r\right)\\ & & +\,{D}_{2}{r}^{-\sqrt{1-4{b}_{2}}}F\left(\Space{0ex}{2.5ex}{0ex}-p+{\rm{i}}2\omega M\right.\\ & & \left.\left.+\,\frac{1}{2}(1-\sqrt{1-4{b}_{2}}),1-\sqrt{1-4{b}_{2}};{\rm{i}}2\omega r\right)\right],\end{array}\end{eqnarray}$
where D1 and D2 are arbitrary constants.The existence and uniqueness of the solution of differential equation ensure that the behavior of the solutions of equation (23 ) near the event horizon and far from the event horizon can be described by equations (36 ) and (42 ), respectively.
When solutions to a differential equation are not constrained by boundary conditions, the independent particular solutions of equations (24 ) and (37 ) exhibit a common characteristic as r → ∞: one particular solution contains the factor eiωr, while the other contains e−iωr. The particular solutions of equations (24 ) and (37 ) sharing the same exponential factor represent approximations of one independent analytic particular solutions of equation (23 ) in the near-region and far-region respectively. In other words, as r → ∞, the two particular solutions in equations (30 ) and (34 ) containing the factor eiωr approximate a single independent particular solution of equation (23 )—one localized near the event horizon and the other in the far region. Similarly, the two particular solutions in equations (36 ) and (42 ) with e−iωr approximate to another independent particular solution of equation (23 ), with analogous spatial distributions. Thus, retaining the ω2 term provides an effective method for matching near-region and far-region solutions.
5. Quasibound frequencies
The uniqueness theorem for second-order differential equations states that if all coefficients are single-valued and analytic within a region, then given two specified boundary conditions, there exists a unique single-valued analytic solution in that domain. The theoretical existence of this solution is guaranteed regardless of whether mathematical techniques are used to explicitly find it. Therefore, the spectrum is determined solely by the boundary conditions. Any method that artificially imposes additional constraints at non-boundary locations to determine the spectrum would contradict the differential equation’s uniqueness theorem.
The primary focus is on ensuring consistency between the solutions in the near-horizon region, far-horizon region, and the exact analytic solution across the entire domain, rather than on the consistency of these solutions in other regions. This is because, in those regions, the near-horizon and far-horizon solutions differ significantly from the exact analytic solution. Consequently, discussing the properties of these solutions in other regions holds no theoretical or practical value.
Specifically, requiring the near-region solution (36) and far-region solution (42) to coincide at an intermediate point in the overlapping region effectively imposes an extra constraint. This approach may initially appear reasonable for determining quasibound frequencies, but it violates the uniqueness theorem for second-order differential equations. Any method that artificially introduces additional conditions within the domain (excluding boundaries) to determine the quasibound frequencies lack rigorous justification; the resulting wave function may fail to satisfy both boundary conditions simultaneously. This occurs because a general solution to a second-order differential equation contains only two arbitrary constants, making it fundamentally impossible for a single solution to satisfy three independent conditions.
As mentioned above, equations (36 ) and (42 ) are consistent with the behavior of the solutions of equation (23 ) at small r − M/2 and larger r, respectively. The asymptotic behavior of the solutions of equation (23 ) is therefore dominated by equations (36 ) and (42 ). Thus we can study the quasibound states and frequencies by means of equations (36 ) and (42 ).
It is widely known that the quasinormal-mode boundary conditions are that the wave is purely ingoing on the event horizon and purely outgoing at spatial infinity, while quasibound-state boundary conditions are that the wave is purely ingoing on the event horizon and tend to zero at spatial infinity, that is to say that the complete radial function, (Ωr)p−sRp, must satisfy the following boundary conditions [16–20]
$\begin{eqnarray}{({\rm{\Omega }}r)}^{p-s}{R}_{p}\sim \,\left\{\begin{array}{cc}{{\rm{e}}}^{-{\rm{i}}\omega {r}_{* }}, & {r}_{* }\to -\infty ,\\ 0, & {r}_{* }\to \infty ,\\ \end{array}\right.\end{eqnarray}$
where r* is called the tortoise coordinate. It is determined by the equation $\begin{eqnarray}{\partial }^{\mu }v{\partial }_{\mu }v=0,\quad {\partial }^{\mu }u{\partial }_{\mu }u=0.\end{eqnarray}$
Here v and u are the Eddington-Finkelstein null coordinates, which take the form $\begin{eqnarray}v=t+{r}_{* },u=t-{r}_{* }.\end{eqnarray}$
By substituting equation (45 ) into equation (44 ) and using metric (1), we derive the exact form of the tortoise coordinate for Schwarzschild equivalent mediums:
$\begin{eqnarray}{r}_{* }=r+\frac{{r}_{H}^{2}}{r}-\frac{1}{2\kappa }\mathrm{ln}\frac{r}{{r}_{H}}+\frac{1}{\kappa }\mathrm{ln}\left|\frac{r-{r}_{H}}{{r}_{H}}\right|.\end{eqnarray}$
Note that κ = 1/(4M) is the surface gravity (not to be confused with the spin coefficient). When r → rH, we have $\begin{eqnarray}{r}_{* }\approx \frac{1}{\kappa }\mathrm{ln}\left|\frac{r-{r}_{H}}{{r}_{H}}\right|.\end{eqnarray}$
Near the event horizon r* → −∞, r → M/2, the complete radial function has the following asymptotic behavior at the exterior event horizon,
$\begin{eqnarray}\begin{array}{rcl}{({\rm{\Omega }}r)}^{p-s}{R}_{p} & = & {C}_{1}{\left(r-\frac{M}{2}\right)}^{2p+{\rm{i}}4M\omega }+{C}_{2}{\left(r-\frac{M}{2}\right)}^{-{\rm{i}}4M\omega }\\ & = & {C}_{1}{{\rm{e}}}^{(2p\kappa +{\rm{i}}\omega ){r}_{* }}+{C}_{2}{{\rm{e}}}^{-{\rm{i}}\omega {r}_{* }},\end{array}\end{eqnarray}$
where all the remaining constants are included in C1 and C2. The boundary condition at r = M/2 requires that C1 = 0, and hence we have $\begin{eqnarray}g={C}_{2}{x}^{1-\frac{\gamma }{2}}F(\alpha +1-\gamma ,2-\gamma ;{\rm{i}}2\omega x)\left({\rm{for}}\,{\rm{small}}\,r-\frac{M}{2}\right).\end{eqnarray}$
Once we have determined the solution near the event horizon, a natural question arises: which of the two particular solutions in equation (42 ) corresponds to the given particular solution in equation (36 )? To answer this question, we retain the first term in square brackets in equation (24 ). When r → ∞, the limiting function of equation (49 ) contains the factor eiωr, and similarly, the limiting function of the second particular solution in equation (42 ) also contains this factor as r → ∞. Therefore, we conclude that the second particular solution in equation (42 ) and equation (49 ) correspond to the particular solution of equation (23 ) that contains the factor eiωr in the asymptotic limit of r → ∞.
To match equation (49 ), we should have that the value of D1 in equation (42 ) is zero, then one finds that for large values of r,
$\begin{eqnarray}\begin{array}{rcl}g & = & {D}_{2}{r}^{(1-\sqrt{1-4{b}_{2}})/2}{{\rm{e}}}^{-{\rm{i}}\omega r}F\left(\Space{0ex}{3.5ex}{0ex}1-p+{\rm{i}}2\omega M\right.\\ & & \left.-\,\frac{1+\sqrt{1-4{b}_{2}}}{2},1-\sqrt{1-4{b}_{2}};{\rm{i}}2\omega r\right).\end{array}\end{eqnarray}$
To obtain expressions for the quasibound frequencies we need the r → ∞ (or r* → ∞) asymptotic form of the complete radial function (Ωr)p−sRp(r) which is given as
$\begin{eqnarray}{({\rm{\Omega }}r)}^{p-s}{R}_{p}(r)\sim {D}_{2}{r}^{-(s+p+1)}{{\rm{e}}}^{{\rm{i}}\omega r}.\end{eqnarray}$
It should be noted that equation (51 ) has the form eiωr/r when p = −s, which is an outgoing spherical wave we do not want. As the physical requirement is that the imaginary parts of the quasibound frequencies are negative (due to the time factor e−iωt being finite as t → ∞), this leads to an infinite value of equation (51 ) in the limit r → ∞, and thus does not satisfy the boundary condition (43) at spatial infinity. On the other hand, the Kummer function in equation (50 ) has the asymptotic form $F\sim {r}^{-p+(1-\sqrt{1-4{b}_{2}})/2}{{\rm{e}}}^{{\rm{i}}2\omega r}$ for very large r. The Kummer function’s contribution to equation (51 ) evidently includes the factor ei2ωr. Therefore, there is only one way to wiggle out of this: for quasibound states, the Kummer function must be forcibly truncated and become Laguerre polynomials. Thus, asymptotic form of the complete radial function is the product of a polynomial and e−iωr [this can also be seen directly from equation (50 )], which becomes zero as r → ∞, thereby satisfying the boundary condition.
The condition for the Kummer’s series in equation (50 ) to become the Laguerre polynomials is 52 ) are unphysical. Obviously all quasibound frequencies result in damped modes. The degrees of damping depend on the values of l and N. Underdamped cases occur only for $2l+1\gt \sqrt{15/7}N$, which exhibit oscillatory behavior and the damping absent when N = 0. However, the case N = 0 must be excluded because of the boundary condition at r → ∞. As shown in figure 1, a finite number of underdamped modes can arise for a given value of l.
$\begin{eqnarray}1-\frac{1}{2}(1+\sqrt{1-4{b}_{2}})-p+{\rm{i}}2M\omega =-k,\end{eqnarray}$
where k = 0, 1, 2, 3… . This gives remarkably simple formulae for the frequencies of the quasibound states (quasibound frequencies): $\begin{eqnarray}\frac{7}{2}M\omega =\left\{\begin{array}{cc}-{\rm{i}}N+\frac{1}{4}\sqrt{14{(2l+1)}^{2}-30{N}^{2}}, & \rm{if}\,2l+1\gt \sqrt{15/7}N;\\ -{\rm{i}}[N\pm \frac{1}{4}\sqrt{-14{(2l+1)}^{2}+30{N}^{2}}], & \rm{if}\,N\lt 2l+1\lt \sqrt{15/7}N;\\ -{\rm{i}}[N+\frac{1}{4}\sqrt{-14{(2l+1)}^{2}+30{N}^{2}}], & \rm{if}\,2l+1\leqslant N.\end{array}\right.\end{eqnarray}$
where N = 2k − 2p + 1 = 1, 2, 3,… . Note that the other solutions of equation (
Figure 1. Quasibound frequencies in units of 1/(3.5M) for $2l+1\gt \sqrt{15/7}N$. The absolute value of ${\rm{Im}}(3.5M\omega )$ is an odd integer for bosons (B) and an even integer for fermions (F). All frequencies on the same curve have the same quantum number l. The curves towards to the top left of the graph correspond to larger l values, which are 0, 1, 2, 3, 4, 5, 6, and 7. |
When $2l+1\lt \sqrt{15/7}N$, ω takes only discrete imaginary values, which leads to an exponentially decreasing function of distance (also time). The behavior is called the overdamped modes. There are two branches of overdamped modes for $N\lt 2l+1\lt \sqrt{15/7}N$ which is of special interest. As similar to the underdamped modes, for a given value of l, only a finite number of overdamped modes with the smaller ω of two branches can arise as shown in figure 2.
Figure 2. Quasibound frequency spectrum with the smaller absolute value of ω in units of 1/(3.5M) for $N\lt 2l+1\lt \sqrt{15/7}N$. N is an odd integer for bosons (B) and an even integer for fermions (F). All frequencies on the same curve have the same quantum number l. From left to right, l = 0, 1, 2, 3, 4, 5, 6, and 7. |
In this section, we have found solutions that match in the near-horizon and far-horizon regions. According to the uniqueness theorem for differential equations, the exact analytical solution must coincide with our derived solution in both the near-horizon and far-horizon regions. Therefore, the quasibound frequencies (53) obtained using the boundary conditions are correct, as it is guaranteed by the uniqueness of the solution to the differential equation.
6. Discussion and conclusion
It is well known that the wave equations for different massless spin particles are rather different. Contrary to conventional wisdom, here we have found a single statement, equation (14 ), governing them, which shows that they have such a similar structure in Schwarzschild equivalent mediums. Equation (14 ) separates, and the solutions can be written as products of the time-dependent function, the angular function, and the radial function. The time-dependent function is e−iωt, and the angular function is the spin-weighted spherical harmonics. The actual shape of the metric functions affects only the radial function. In particular, the function of Ω(r) in the solutions is a product factor of the radial function. Therefore, for a wave packet, the conformal factor affects the ‘envelope’, but does not affect the ‘ripples’.
We have investigated the radial equation in the near-horizon and the far zones, and found that the radial equation can be transformed to a degenerate hypergeometric equation in any zone. Hence the solutions are governed by Kummer functions. Using the boundary condition (43), we have obtained the quasibound frequencies, which are given by equation (53 ). It is important to remember that p = ±s. Thus equation (53 ) admits two kinds of particles: bosons, for which N is an odd integer, and fermions, for which N is an even integer. As one can see directly from figures 1 and 2, the quasibound frequencies of different bosons are identical at the same N and l, and the same is true of different fermions, but a quasibound frequency for bosons can never equal a quasibound frequency for fermions. This fact provides a theoretical basis for simulation between particle-waves of same statistical properties.
Due to the fast development of metamaterials and transformation optics, one can realize beam shaping within dielectric slab samples with predesigned refractive index varying so as to create curved space environment for electromagnetic waves. As mentioned in the introduction, artificial optical materials have been proposed to study the various aspects of curved spacetimes. However, simulation of gravitational waves in curved spacetimes remains a challenge. Here our results provide a theoretical basis for electromagnetic waves are used to mimic gravitational waves. In particular, artificial optical materials have been realized to mimic Schwarzschild spacetime [41, 55], therefore the transformation optics approach can be used to test our predictions.
Appendix A: Spin coefficients and Weyl scalars
The 12 spin coefficients are defined by [56–58]
$\begin{eqnarray}\begin{array}{rcl}\kappa & = & { \triangledown }_{\nu }{l}_{\mu }{m}^{\mu }{l}^{\nu },\quad \lambda =-{ \triangledown }_{\nu }{n}_{\mu }{\bar{m}}^{\mu }{\bar{m}}^{\nu },\\ \sigma & = & { \triangledown }_{\nu }{l}_{\mu }{m}^{\mu }{m}^{\nu },\quad \nu =-{ \triangledown }_{\nu }{n}_{\mu }{\bar{m}}^{\mu }{n}^{\nu },\\ \rho & = & { \triangledown }_{\nu }{l}_{\mu }{m}^{\mu }{\bar{m}}^{\nu },\quad \tau ={ \triangledown }_{\nu }{l}_{\mu }{m}^{\mu }{n}^{\nu },\\ \mu & = & -{ \triangledown }_{\nu }{n}_{\mu }{\bar{m}}^{\mu }{m}^{\nu },\quad \pi =-{ \triangledown }_{\nu }{n}_{\mu }{\bar{m}}^{\mu }{l}^{\nu },\\ \alpha & = & \frac{1}{2}({ \triangledown }_{\nu }{l}_{\mu }{n}^{\mu }{\bar{m}}^{\nu }-{ \triangledown }_{\nu }{m}_{\mu }{\bar{m}}^{\mu }{\bar{m}}^{\nu }),\\ \beta & = & \frac{1}{2}({ \triangledown }_{\nu }{l}_{\mu }{n}^{\mu }{m}^{\nu }-{ \triangledown }_{\nu }{m}_{\mu }{\bar{m}}^{\mu }{m}^{\nu }),\\ \gamma & = & \frac{1}{2}({ \triangledown }_{\nu }{l}_{\mu }{n}^{\mu }{n}^{\nu }-{ \triangledown }_{\nu }{m}_{\mu }{\bar{m}}^{\mu }{n}^{\nu }),\\ \varepsilon & = & \frac{1}{2}({ \triangledown }_{\nu }{l}_{\mu }{n}^{\mu }{l}^{\nu }-{ \triangledown }_{\nu }{m}_{\mu }{\bar{m}}^{\mu }{l}^{\nu }).\end{array}\end{eqnarray}$
Here, lμ, nμ, mμ and ${\bar{m}}_{\mu }$ are the Newman–Penrose null tetrad. The tetrad consists of two real null vectors, lμ and nμ, and a pair of complex null vectors, mμ and ${\bar{m}}_{\mu }$, which satisfies the orthonormal conditions, ${l}_{\mu }{n}^{\mu }=-{m}_{\mu }{\bar{m}}^{\mu }=1$, and ${l}_{\mu }{l}^{\mu }={n}_{\mu }{n}^{\mu }={m}_{\mu }{m}^{\mu }={\bar{m}}_{\mu }{\bar{m}}^{\mu }=0$. The indexes are raised and lowered using global metric gμν, which in terms of null vectors can be expressed as, ${g}_{\mu \nu }=2{l}_{(\mu }{n}_{\nu )}-2{m}_{(\mu }{\bar{m}}_{\nu )}$.The Newman–Penrose tetrad for Schwarzschild equivalent mediums can be chosen as
$\begin{eqnarray}\begin{array}{rcl}{l}_{\mu } & = & \frac{{{\rm{\Omega }}}^{4}}{n}\left(\frac{1}{n},-1,0,0\right),\\ {n}_{\mu } & = & \frac{1}{2{{\rm{\Omega }}}^{2}}(1,n,0,0),\\ {m}_{\mu } & = & -\frac{{\rm{\Omega }}r}{\sqrt{2}}(0,0,1,{\rm{i}}\sin \theta ).\end{array}\end{eqnarray}$
Using the metric (5) and null tetrad (A2), the spin coefficients can be written as $\begin{eqnarray}\begin{array}{rcl}\kappa & = & \sigma =\nu =\lambda =\pi =\tau =0,\\ \rho & = & -\frac{{{\rm{\Omega }}}^{2}}{n}\left(\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}}+\frac{1}{r}\right),\quad \mu =-\frac{n}{2{{\rm{\Omega }}}^{4}}\left(\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}}+\frac{1}{r}\right),\\ \varepsilon & = & \frac{{{\rm{\Omega }}}^{2}}{n}\left(2\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}}-\frac{{n}^{{\prime} }}{n}\right),\quad \gamma =-\frac{n}{2{{\rm{\Omega }}}^{4}}\frac{{{\rm{\Omega }}}^{{\prime} }}{{\rm{\Omega }}},\\ \alpha & = & -\beta =-\frac{1}{2\sqrt{2}{\rm{\Omega }}r}\cot \theta .\end{array}\end{eqnarray}$
The five Weyl scalars are defined by [56–58]
$\begin{eqnarray}\begin{array}{rcl}{\psi }_{0} & = & -{C}_{\mu \nu \rho \sigma }{l}^{\mu }{m}^{\nu }{l}^{\rho }{m}^{\sigma },\quad {\psi }_{1}=-{C}_{\mu \nu \rho \sigma }{l}^{\mu }{n}^{\nu }{l}^{\rho }{m}^{\sigma },\\ {\psi }_{2} & = & -\,\frac{1}{2}{C}_{\mu \nu \rho \sigma }({l}^{\mu }{n}^{\nu }{l}^{\rho }{n}^{\sigma }-{l}^{\mu }{n}^{\nu }{m}^{\rho }{\bar{m}}^{\sigma }),\\ {\psi }_{3} & = & -\,{C}_{\mu \nu \rho \sigma }{\bar{m}}^{\mu }{n}^{\nu }{l}^{\rho }{n}^{\sigma },\quad {\psi }_{4}=-{C}_{\mu \nu \rho \sigma }{\bar{m}}^{\mu }{n}^{\nu }{\bar{m}}^{\rho }{n}^{\sigma },\end{array}\end{eqnarray}$
where Cμνρσ is the Weyl tensor, which satisfies $\begin{eqnarray}\begin{array}{rcl}{R}_{\mu \nu \rho \sigma } & = & {C}_{\mu \nu \rho \sigma }+\frac{1}{2}({g}_{\mu \rho }{R}_{\nu \sigma }-{g}_{\mu \sigma }{R}_{\nu \rho }-{g}_{\nu \rho }{R}_{\mu \sigma }\\ & & +\,{g}_{\nu \sigma }{R}_{\mu \rho })+\frac{1}{6}({g}_{\mu \sigma }{g}_{\nu \rho }-{g}_{\mu \rho }{g}_{\nu \sigma })R.\end{array}\end{eqnarray}$
Appendix B: Spin field equations
The Weyl neutrino satisfies the wave equation [59]
$\begin{eqnarray}{{\rm{\nabla }}}_{A{A}^{{\prime} }}\,{P}^{A}=0,\end{eqnarray}$
where PA is the two-component spinor, ${{\rm{\nabla }}}_{A{A}^{{\prime} }}$ is the symbol for covariant spinor differentiation. Equation (B1 ) can be written in the Newman–Penrose formalism in the form
$\begin{eqnarray}\begin{array}{rcl} & & (D+\varepsilon -\rho ){P}^{0}+(\bar{\delta }+\pi -\alpha ){P}^{1}=0,\\ & & (\delta +\beta -\tau ){P}^{0}+({\rm{\Delta }}+\mu -\gamma ){P}^{1}=0.\end{array}\end{eqnarray}$
In type-D spacetime, the equation can be decoupled into [58]
$\begin{array}{l} {[(D+\bar{\varepsilon}-\rho-\bar{\rho})(\Delta-\gamma+\mu)} \\ -(\delta-\bar{\alpha}-\tau+\bar{\pi})(\bar{\delta}-\alpha+\pi)] \Phi_{+1 / 2}=0 \\ {[(\Delta-\bar{\gamma}+\mu+\bar{\mu})(D+\varepsilon-\rho)} \\ -(\bar{\delta}+\bar{\beta}-\pi-\bar{\tau})(\delta+\beta-\tau)] \Phi_{-1 / 2}=0 \end{array}$
with Φ+1/2 = P1, and Φ−1/2 = −P0.
$\begin{array}{l} {[(D-\varepsilon+\bar{\varepsilon}-2 \rho-\bar{\rho})(\Delta-2 \gamma+\mu)} \\ -(\delta+\bar{\pi}-\bar{\alpha}-\beta-2 \tau)(\bar{\delta}+\pi-2 \alpha)] \Phi_{+1}=0, \\ {[(\Delta+\gamma-\bar{\gamma}+2 \mu+\bar{\mu})(D+2 \varepsilon-\rho)} \\ -[\bar{\delta}-\bar{\tau}+\bar{\beta}+\alpha+2 \pi)(\delta-\tau+2 s \beta) \Phi_{-1}=0 . \end{array}$
$\begin{array}{l} {[D-2 \varepsilon+\bar{\varepsilon}-3 \rho-\bar{\rho})(\Delta-3 \gamma+\mu)} \\ \left.-(\delta+\bar{\pi}-\bar{\alpha}-2 \beta-3 \tau)(\bar{\delta}+\pi-3 \alpha)-\psi_{2}\right] \Phi_{+3 / 2}=0 \\ {[(\Delta+2 \gamma-\bar{\gamma}+3 \mu+\bar{\mu})(D+3 \varepsilon-\rho)} \\ \left.-(\bar{\delta}-\bar{\tau}+\bar{\beta}+2 \alpha+3 \pi)(\delta-\tau+3 \beta)-\psi_{2}\right] \Phi_{-3 / 2}=0 . \end{array}$
$\begin{array}{l} {[(D-3 \varepsilon+\bar{\varepsilon}-4 \rho-\bar{\rho})(\Delta-4 \gamma+\mu)} \\ \left.-(\delta+\bar{\pi}-\bar{\alpha}-3 \beta-4 \tau)(\bar{\delta}+\pi-4 \alpha)-3 \psi_{2}\right] \Phi_{+2}=0 \\ {[(\Delta+3 \gamma-\bar{\gamma}+4 \mu+\bar{\mu})(D+4 \varepsilon-\rho)} \\ \left.-[\bar{\delta}-\bar{\tau}+\bar{\beta}+3 \alpha+4 \pi](\delta-\tau+4 \beta)-3 \psi_{2}\right] \Phi_{-2}=0 \end{array}$