Spin-1/2 Particles under the Influence of a Uniform Magnetic Field in the Interior Schwarzschild Solution

Abstract

The relativistic wave equation for spin-1/2 particles in the interior Schwarzschild solution in the presence of a uniform magnetic field is obtained. The fully relativistic regime is considered, and the energy levels occupied by the particles are derived as functions of the magnetic field, the radius of the massive sphere and the total mass of the latter. As no assumption is made on the relative strengths of the particles’ interaction with the gravitational and magnetic fields, the relevance of our results to the physics of the interior of neutron stars, where both the gravitational and the magnetic fields are very intense, is discussed.

1. Introduction

The behavior of charged quantum particles inside a uniform magnetic field in Minkowski spacetime displays quantized energy levels—the so-called Landau levels—that were first derived by Rabi based on the relativistic Dirac equation [1] and by Landau based on the Schrödinger equation [2] (see e.g., Ref. [3] for a textbook introduction). This discovery led to many physical applications in condensed matter physics, such as the explanation of the Landau diamagnetism [2], of the Shubnikov–de Haas effect [4], of the de Haas-van Alphen effect [5] and of the integer quantum Hall effect [6]. Remarkably, the Landau levels have even found astrophysical applications. In fact, the physics of neutron stars and other highly magnetized stellar objects is expected to be governed to some extent by the discrete nature of the energy levels occupied by the charged fermions moving under the influence of the intense magnetic field of such astrophysical objects [7,8,9,10,11,12,13,14,15,16]. When we recall that such highly compact astrophysical objects are also sources of very intense gravitational fields, it becomes extremely important to investigate the fate of Landau levels under the influence of a combination of a magnetic field and the gravitational field of a spherical object. Such a study has recently been conducted in Refs. [17,18,19] based on the Klein-Gordon and Schrödinger equations by considering the gravitational interaction of the particle with a massive sphere as a small perturbation compared to the interaction of the particle with the magnetic field. Although relying on the Schrödinger equation and treating gravity perturbatively as done in these studies might be sufficient for laboratory applications [20,21,22,23], when seeking astrophysical applications, implementing a general-relativistic and a non-perturbative approach and appealing to the relativistic Dirac equation are highly recommended.

The Dirac equation in a curved spacetime caused purely by an intense magnetic field–known as the Melvin universe [24]—has been solved exactly in Ref. [25] where the quantized energy levels were extracted. Contemplating curved spacetimes caused purely by a magnetic field is motivated by the well-known fact that the strong magnetic fields of neutron stars could reach up to ${10}^{12-15}$ G on the surface of magnetars [26], and up to ${10}^{17}$ G in their interior (see, e.g., Refs. [27,28] and the references therein). On the other hand, neutron stars are made of neutrons, as well as of a plasma of electrons, protons and muons, the energies of which are affected by such strong magnetic fields. The aim of the present paper is therefore to include the general-relativistic effect of the intense gravitational field of such stellar objects by studying the dynamics of charged fermions moving inside a massive spherical object in the presence of a uniform magnetic field. To the best of our knowledge, such a study has not appeared in the literature before. Unlike the purely magnetic spacetime considered in Ref. [25], however, we consider in the present work the curved spacetime caused purely by the distribution of matter inside the core of the stellar objects rather than by the uniform magnetic field. The contribution of the latter to spacetime curvature will thus be considered negligible. The justification behind such an assumption will be discussed in Section 2.

The remainder of this paper is structured as follows. In Section 2, we derive the Dirac equation in a general static and spherically symmetric curved spacetime in the presence of a uniform magnetic field. We discuss the various symmetries of the equation and the complications one faces when attempting to solve the equation exactly. In Section 3, we apply the equation obtained in Section 2 to the case of the interior Schwarzschild solution. We then extract the quantized energy levels of charged fermions moving along the equator inside a massive spherical object by keeping only the leading terms of the differential equation. In Section 4, we briefly discuss the ways to overcome the limitations we encountered in solving our equation in Section 3. In Section 5, we apply our results to the study of the magnetization of neutron stars. We conclude this paper with a short discussion section in which we summarize our findings.

2. Curved-Spacetime Dirac Equation in a Uniform Magnetic Field

Since our aim in this paper is to consider the Dirac equation in the static and spherically symmetric interior Schwarzschild solution, let us start by deriving in this section the equation for a general static and spherically symmetric curved spacetime in the presence of a uniform magnetic field. For that purpose, let us write our spacetime metric in the following general form:

$$\mathrm{d}{s}^2=-{\mathcal{A}}^2\left(r\right)c^2\mathrm{d}{t}^2+{\mathcal{B}}^2\left(r\right)\mathrm{d}{r}^2+r^2(\mathrm{d}{\theta }^2+{sin}^2\theta \mathrm{d}{\varphi }^2),$$

(1)

where the radial functions $\mathcal{A}\left(r\right)$ and $\mathcal{B}\left(r\right)$ are everywhere regular except, maybe, at the origin $r=0$. In addition, let us assume the constant and uniform magnetic field $\mathbf{B}$ to be parallel to the z-direction. A convenient gauge for the vector potential would then be $A_{\mu }=(0,0,0,\frac{1}{2}B{r}^2{sin}^2\theta )$. The Dirac equation of a particle of mass m and of charge e (that we take to be the negative charge of the electron for definiteness) moving inside a curved spacetime and minimally coupled to a Maxwell field $A_{\mu }$ reads [29,30],

$$i\hslash e_a^{\mu }{\gamma }^a\left({\partial }_{\mu }+\frac{1}{8}{\omega }_{\mu }^{bc}[{\gamma }_b,{\gamma }_c]-\frac{ie{A}_{\mu }}{\hslash }\right)\psi =mc\psi .$$

(2)

The quantities $e_a^{\mu }$ are the inverse of the spacetime tetrads $e_{\mu }^a$, whereas ${\gamma }^a$ are the gamma matrices, ${\omega }_{\mu }^{bc}$ is the spin connection and $[{\gamma }_b,{\gamma }_c]$ is the commutator of the gamma matrices. The spin connection is antisymmetric, ${\omega }_{\mu }^{bc}=-{\omega }_{\mu }^{cb}$, and it is expressed in terms of the tetrad fields and the Christoffel symbols ${\Gamma }_{\mu \nu }^{\lambda }$ of the spacetime by ${\omega }_{\mu }^{bc}=e_{\nu }^b{\partial }_{\mu }{e}^{c\nu }+{\Gamma }_{\mu \nu }^{\lambda }{e}_{\lambda }^b{e}^{c\nu }$.

On the other hand, two choices for the four curved-spacetime gamma matrices ${\gamma }^{\mu }=e_{\mu }^a{\gamma }^a$ are possible. One may choose either to work with the tetrads’ axes parallel to the r, $\theta$ and $\varphi$ coordinates, or choose to orient the tetrads’ axes parallel to some rectangular coordinate system [31]. When combined with the standard representation of the constant flat-spacetime gamma matrices ${\gamma }^a$, given explicitly by

$${\gamma }^0=\left(\begin{array}{cc}\mathrm{𝟙}& 0\\ 0& -\mathrm{𝟙}\end{array}\right),{\gamma }^k=\left(\begin{array}{cc}0& {\sigma }_k\\ -{\sigma }_k& 0\end{array}\right),$$

(3)

where $\mathrm{𝟙}$ is the $2\times 2$ unit matrix and ${\sigma }_k$ (for $k=1,2,3$) are the three Pauli matrices, the curved-space gamma matrices ${\gamma }^{\mu }$ take on a simpler form in the first case than in the second [31]. The spinor wavefunction $\psi$ in the two cases will be different, but physical quantities as well as the radial wave equation will all be the same as the two representations are simply related by a unitary transformation [31]. However, to deal with the mixed cylindrical and spherical symmetries imposed, respectively, by the magnetic and gravitational fields, it is easier to adopt the spherical-coordinates representation of the gamma matrices, in which the three Pauli matrices take the form

$${\sigma }_r=\left(\begin{array}{cc}cos\theta & sin\theta e^{-i\varphi }\\ sin\theta e^{i\varphi }& -cos\theta \end{array}\right),{\sigma }_{\theta }=\left(\begin{array}{cc}-sin\theta & cos\theta e^{-i\varphi }\\ cos\theta e^{i\varphi }& sin\theta \end{array}\right),{\sigma }_{\varphi }=\left(\begin{array}{cc}0& -i{e}^{-i\varphi }\\ i{e}^{i\varphi }& 0\end{array}\right).$$

(4)

We now choose the spacetime tetrads for the metric (1) to be $e_{\mu }^a=\mathrm{diag}(\mathcal{A},\mathcal{B},r,rsin\theta )$. With these tetrads, we easily evaluate the nonvanishing components of the contracted spin connection with the commutator of the gamma matrices to be:

$$\begin{array}{cc}\hfill \frac{1}{8}{\omega }_0^{ab}[{\gamma }_a,{\gamma }_b]& =\frac{{\mathcal{A}}^{\prime }}{2\mathcal{B}}{\gamma }^0{\gamma }^1,\frac{1}{8}{\omega }_2^{ab}[{\gamma }_a,{\gamma }_b]=\frac{1}{2\mathcal{B}}{\gamma }^1{\gamma }^2,\hfill \\ \hfill \frac{1}{8}{\omega }_3^{ab}[{\gamma }_a,{\gamma }_b]& =\frac{sin\theta }{2\mathcal{B}}{\gamma }^1{\gamma }^3+\frac{cos\theta }{2}{\gamma }^2{\gamma }^3.\hfill \end{array}$$

(5)

Plugging these, together with the vector potential $A_{\mu }$, into Equation (2) and then multiplying the whole equation on the left by $\beta ={\gamma }^0$ and using the definition ${\alpha }^k=\beta {\gamma }^k$ of the Dirac alpha matrices, the resulting equation reads,

$$\left[\frac{i{\partial }_t}{\mathcal{A}c}\psi +\frac{i{\alpha }^1}{\mathcal{B}}\left({\partial }_r+\frac{1}{r}+\frac{{\mathcal{A}}^{\prime }}{2\mathcal{A}}\right)+\frac{i{\alpha }^2}{r}\left({\partial }_{\theta }+\frac{cot\theta }{2}\right)+\frac{i{\alpha }^3}{rsin\theta }\left({\partial }_{\varphi }-\frac{ie{A}_{\varphi }}{\hslash }\right)\right]\psi =\frac{mc}{\hslash }\beta \psi .$$

(6)

We are interested in this paper in stationary states of energy E, for which we have ${\partial }_t\psi =-\frac{i}{\hslash }E\psi$. On the other hand, we may simplify further Equation (6) by setting $\psi =\Psi /\left(r\sqrt{\mathcal{A}sin\theta }\right)$. With these two replacements, the resulting equation in $\Psi$ takes the form

$$\left[\frac{E}{\mathcal{A}\hslash c}+\frac{i{\alpha }^1}{\mathcal{B}}{\partial }_r+\frac{i{\alpha }^2}{r}{\partial }_{\theta }+\frac{i{\alpha }^3}{rsin\theta }\left({\partial }_{\varphi }-\frac{ieB}{2\hslash }{r}^2{sin}^2\theta \right)-\frac{mc}{\hslash }\beta \right]\Psi =0.$$

(7)

This is the general equation obeyed by fermions under the influence of the uniform magnetic field.

Unlike the well-known case of the Dirac equation in a spherically symmetric potential, i.e., in a central potential (see, e.g., Ref. [32]), our Equation (7) may not be separated in the radial variables r and the angular variables $\theta$ and $\varphi$. This is due to the mixture of the spherical symmetry imposed by the gravitational field and the cylindrical symmetry imposed by the vertical uniform magnetic field $\mathbf{B}$. One might then be tempted to switch to cylindrical coordinates $(\rho ,\varphi ,z)$ to benefit from the cylindrical symmetry pertaining to the magnetic field. Unfortunately, even doing so would not make Equation (7) easier to solve. In fact, by making the substitutions $\rho =rsin\theta$ and $z=rcos\theta$, Equation (7) would contain operators of the form ${({\rho }^2+z^2)}^{-\frac{1}{2}}z{\partial }_{\rho }$ and ${({\rho }^2+z^2)}^{-\frac{1}{2}}\rho {\partial }_z$. It is then evident that even if we neglect the z-variation of $\Psi$ by assuming ${\partial }_z\Psi \approx 0$, i.e., by neglecting the vertical momentum of the particle, we would still be left with an equation that is not separable in the variables $\rho$, z and $\varphi$.

It seems then that the only choice left is to set $z=0$ as well. Physically, this is due (in the absence of any other interaction of the particle) to the unstable nature of any position $z\ne 0$, caused by the transverse gravitational attraction that pulls the particle away from any plane other than the equator. This observation can also be reached by referring to the variable $\theta$ in Equation (7). In fact, the only way of removing the $\theta$-dependence in that equation is to set $\theta =\frac{\pi }{2}$, i.e., to restrict ourselves to particles moving along the equatorial plane. This is because setting $\theta ={\theta }_0\ne \frac{\pi }{2}$ would prevent us from letting $\rho =rsin{\theta }_0$ vary and have at the same time planar motion. If ${\theta }_0$ is set to a constant other than $\frac{\pi }{2}$, then for $\rho$ to change, r must change too, which makes the particle leave the plane $\theta ={\theta }_0$. Only along the equatorial plane ${\theta }_0=\frac{\pi }{2}$ does $\rho =rsin\frac{\pi }{2}=r$ become a variable describing planar motion.

We will come back to the possibility of considering a plane other than the equator in Section 4. Before setting $\theta =\frac{\pi }{2}$, however, we should first extract the general single differential equation from the coupled Dirac spinors.

3. The Quantized Energy Levels

To extract a quantization condition, if any, on the energy E from Equation (7), we need first to decompose the Dirac spinor into a pair of 2-spinors as follows

$$\Psi (r,\theta ,\varphi )=\left[\begin{array}{c}\Phi (r,\theta ,\varphi )\\ \Theta (r,\theta ,\varphi )\end{array}\right].$$

(8)

Substituting this ansatz into Equation (7), the latter splits into the following coupled first-order differential equations:

$$\begin{array}{c}\hfill \left(\frac{E}{\mathcal{A}\hslash c}-\frac{mc}{\hslash }\right)\Phi +\left[\frac{i{\sigma }_r}{\mathcal{B}}{\partial }_r+\frac{i{\sigma }_{\theta }}{r}{\partial }_{\theta }+\frac{i{\sigma }_{\varphi }}{rsin\theta }\left({\partial }_{\varphi }-\frac{ieB}{2\hslash }{r}^2{sin}^2\theta \right)\right]\Theta =0,\\ \hfill \left(\frac{E}{\mathcal{A}\hslash c}+\frac{mc}{\hslash }\right)\Theta +\left[\frac{i{\sigma }_r}{\mathcal{B}}{\partial }_r+\frac{i{\sigma }_{\theta }}{r}{\partial }_{\theta }+\frac{i{\sigma }_{\varphi }}{rsin\theta }\left({\partial }_{\varphi }-\frac{ieB}{2\hslash }{r}^2{sin}^2\theta \right)\right]\Phi =0.\end{array}$$

(9)

It is clear from these coupled differential equations that given the $\rho$-dependence of all the coupling terms, extracting one of the 2-spinors from one equation and substituting it into the other would result in a very complicated single equation. This is unlike what happens with the standard Dirac equation for the case of the central potential of the hydrogen atom [32]. It would be natural then at this stage to attempt to solve this system of equations by applying the approach developed in Refs. [33,34,35], which consists of applying a gauge transformation followed by a unitary transformation to turn one of the coupling terms in one of the two equations into a constant. Indeed, such a method considerably simplifies the resulting final single equation one arrives at. Unfortunately, such an approach cannot be consistently applied in our case because that approach is specifically designed only for a certain type of interactions of the particle with the gauge field, which is not the case for our present physical system1. The other possibility would then be to attempt to rather adopt the approach developed in Ref. [36] for the Dirac equation in curved spacetime. Unfortunately, that approach is not helpful to us either as it becomes rapidly involved even when the particle is not coupled to a gauge potential. For this reason, such an approach has been applied in Ref. [36] only in $1+1$ dimensions. We are thus left only with the old strategy for dealing with coupled differential equations of the type (9).

Extracting the 2-spinor $\Theta$ from the second line of Equation (9) and plugging it into the first, we obtain the following single equation for the 2-spinor $\Phi$:

$$\begin{array}{cc}\hfill & \frac{{\Phi }_{,rr}}{{\mathcal{B}}^2}+\left(\frac{2}{\mathcal{B}r}-\frac{{\mathcal{B}}_{,r}}{{\mathcal{B}}^3}+\frac{{\mathcal{A}}_{,r}E}{\mathcal{A}{\mathcal{B}}^2\left[E+\mathcal{A}m{c}^2\right]}\right){\Phi }_{,r}\hfill \\ \hfill & +\frac{{\Phi }_{,\theta \theta }}{r^2}+\left(\frac{i\left[\mathcal{B}-1\right]{\sigma }_{\varphi }}{\mathcal{B}{r}^2}+\frac{cot\theta }{r^2}+\frac{i{\sigma }_{\varphi }}{r}\frac{{\mathcal{A}}_{,r}E}{\mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}\right){\Phi }_{,\theta }\hfill \\ \hfill & +\frac{{\Phi }_{,\varphi \varphi }}{r^2{sin}^2\theta }+\left(\frac{i\left[1-\mathcal{B}\right]{\sigma }_{\theta }}{\mathcal{B}{r}^{2}sin\theta }-\frac{ieB}{\hslash }-\frac{i{\sigma }_{\theta }}{rsin\theta }\frac{{\mathcal{A}}_{,r}E}{\mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}\right){\Phi }_{,\varphi }\hfill \\ \hfill & +(\frac{E^2}{{\mathcal{A}}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}-\frac{e^2{B}^2{r}^2{sin}^2\theta }{4{\hslash }^2}-\frac{\left[\mathcal{B}+1\right]eB{\sigma }_{\theta }sin\theta }{2\hslash \mathcal{B}}+\frac{{\sigma }_{r}eBcos\theta }{\hslash }\hfill \\ \hfill & -\frac{{\mathcal{A}}_{,r}EeBr{\sigma }_{\theta }sin\theta }{2\hslash \mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]})\Phi =0.\hfill \end{array}$$

(10)

With this result, we can now proceed to restrict the motion of the particle to the equatorial plane by setting $\theta =\frac{\pi }{2}$. Therefore, we also need to perform the replacement $rsin\theta =r\equiv \rho$ everywhere in this equation. On the other hand, for the $\theta$-derivatives, we use the following identities valid for motion restricted to a plane: ${\partial }_{\theta }=rcos\theta {\partial }_{\rho }$ and ${\partial }_{\theta }^2=-rsin\theta {\partial }_{\rho }+r^2{cos}^2\theta {\partial }_{\rho }^2$. Therefore, for $\theta =\frac{\pi }{2}$ we should discard the partial derivative ${\partial }_{\theta }$ and perform the replacement ${\partial }_{\theta }^2\to -\rho {\partial }_{\rho }$. The above equation then becomes

$$\begin{array}{cc}\hfill & \frac{{\Phi }_{,\rho \rho }}{{\mathcal{B}}^2}+\left(\frac{2-\mathcal{B}}{\mathcal{B}\rho }-\frac{{\mathcal{B}}_{,\rho }}{{\mathcal{B}}^3}+\frac{{\mathcal{A}}_{,\rho }E}{\mathcal{A}{\mathcal{B}}^2\left[E+\mathcal{A}m{c}^2\right]}\right){\Phi }_{,\rho }\hfill \\ \hfill & +\frac{{\Phi }_{,\varphi \varphi }}{{\rho }^2}+\left(\frac{i\left[\mathcal{B}-1\right]{\sigma }_z}{\mathcal{B}{\rho }^2}-\frac{ieB}{\hslash }+\frac{i{\sigma }_z}{\rho }\frac{{\mathcal{A}}_{,\rho }E}{\mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}\right){\Phi }_{,\varphi }\hfill \\ \hfill & +\left(\frac{E^2}{{\mathcal{A}}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}-\frac{e^2{B}^2{\rho }^2}{4{\hslash }^2}+\frac{\left[\mathcal{B}+1\right]eB{\sigma }_z}{2\hslash \mathcal{B}}+\frac{{\mathcal{A}}_{,\rho }EeB\rho {\sigma }_z}{2\hslash \mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}\right)\Phi =0.\hfill \end{array}$$

(11)

We have used here the fact that according to Equation (4), the spherical-coordinates Pauli matrix ${\sigma }_{\theta }$ reduces to $-{\sigma }_z$ for $\theta =\frac{\pi }{2}$ (as expected), where ${\sigma }_z$ is the usual third Pauli matrix in Cartesian coordinates. On the other hand, according to expressions (4) of the Pauli matrices, we see in Equation (11) that all the terms are diagonal and that the equation does not involve the angular variable $\varphi$ explicitly. Therefore, without loss of generality, we can write the two independent solutions ${\Phi }_{+}(\rho ,\varphi )$ and ${\Phi }_{-}(\rho ,\varphi )$ of the equation in a basis made of the spin-eigenstates along the z-direction with eigenvalues $s=\pm 1$, respectively. We thus introduce arbitrary radial functions $f_{+}\left(\rho \right)$ and $f_{-}\left(\rho \right)$ such that

$${\Phi }_{+}\left(\rho \right)=e^{i\ell \varphi }\left[\begin{array}{c}{f}_{+}\left(\rho \right)\\ 0\end{array}\right],{\Phi }_{-}\left(\rho \right)=e^{i\ell \varphi }\left[\begin{array}{c}0\\ f_{-}\left(\rho \right)\end{array}\right].$$

(12)

Imposing the periodic condition along the equator, $\Phi (\rho ,\varphi +2\pi )=\Phi (\rho ,\varphi )$, we learn that the quantum number ℓ must be an integer number, which we take here to be nonnegative in accordance with the negative sign we chose for the charge e. Plugging these ansatzes into Equation (11) yields the following second-order differential equation in $\rho$:

$$\begin{array}{cc}\hfill & \frac{f_s^{″}}{{\mathcal{B}}^2}+\left(\frac{2-\mathcal{B}}{\mathcal{B}\rho }-\frac{{\mathcal{B}}^{\prime }}{{\mathcal{B}}^3}+\frac{{\mathcal{A}}^{\prime }E}{\mathcal{A}{\mathcal{B}}^2\left[E+\mathcal{A}m{c}^2\right]}\right)f_s^{\prime }\hfill \\ \hfill & +(\frac{E^2}{{\mathcal{A}}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}+\frac{eB\ell }{\hslash }-\frac{e^2{B}^2{\rho }^2}{4{\hslash }^2}-\frac{\mathcal{B}{\ell }^2+\left[\mathcal{B}-1\right]s\ell }{\mathcal{B}{\rho }^2}+\frac{\left[\mathcal{B}+1\right]eBs}{2\hslash \mathcal{B}}\hfill \\ \hfill & -\frac{{\mathcal{A}}^{\prime }Es\left[2\hslash \ell -eB{\rho }^2\right]}{2\hslash \rho \mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]})f_s=0,\hfill \end{array}$$

(13)

Here, and henceforth, a prime denotes a derivative with respect to the variable $\rho$, and the function $f_s\left(\rho \right)$ stands for the two cases $f_{\pm }\left(\rho \right)$ of $s=\pm 1$, respectively. Next, we introduce the following ansatz:

$$f_s\left(\rho \right)={\left(\frac{\mathcal{B}}{\mathcal{A}}E+\mathcal{B}m{c}^2\right)}^{\frac{1}{2}}{F}_s\left(\rho \right),$$

(14)

for some radial function $F_s\left(\rho \right)$. Please note that only the convergence of the radial function $F_s\left(\rho \right)$ is required henceforth, for the argument of the square root in this ansatz converges for both large and small values of $\rho$ as we shall see below after introducing the explicit components of the spacetime metric. Plugging now this ansatz into Equation (13), the latter takes the following form:

$$\begin{array}{c}\frac{F_s^{″}}{{\mathcal{B}}^2}+\frac{2-\mathcal{B}}{\mathcal{B}\rho }{F}_s^{\prime }+[\frac{E^2}{{\mathcal{A}}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}+\frac{eB\ell }{\hslash }-\frac{e^2{B}^2{\rho }^2}{4{\hslash }^2}-\frac{\mathcal{B}{\ell }^2+(\mathcal{B}-1)s\ell }{\mathcal{B}{\rho }^2}\hfill \\ +\frac{(\mathcal{B}+1)eBs}{2\hslash \mathcal{B}}-\frac{{\mathcal{A}}^{\prime }Es(2\hslash \ell -eB{\rho }^2)}{2\hslash \rho \mathcal{A}\mathcal{B}(E+\mathcal{A}m{c}^2)}-\frac{1}{2{\mathcal{B}}^2}{\left(\frac{{\mathcal{A}}^{\prime }E}{\mathcal{A}[E+\mathcal{A}m{c}^2]}-\frac{{\mathcal{B}}^{\prime }}{\mathcal{B}}\right)}^{\prime }\\ \hspace{1em}\hspace{1em}\hspace{1em} -{\left(\frac{{\mathcal{A}}^{\prime }E}{2\mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}-\frac{{\mathcal{B}}^{\prime }}{2{\mathcal{B}}^2}\right)}^2-\frac{2-\mathcal{B}}{2\mathcal{B}\rho }\left(\frac{{\mathcal{A}}^{\prime }E}{\mathcal{A}\left[E+\mathcal{A}m{c}^2\right]}-\frac{{\mathcal{B}}^{\prime }}{\mathcal{B}}\right)]F_s=0.\hfill \end{array}$$

(15)

Finally, the following ansatz,

$$F_s\left(\rho \right)=G_s\left(\rho \right)exp\left[\int \frac{{(1-\mathcal{B})}^2}{2\rho }\mathrm{d}\rho \right],$$

(16)

for an arbitrary radial function $G_s\left(\rho \right)$, allows us to convert Equation (15) into the following more useful form:

$$\begin{array}{cc}\hfill & G_s^{″}+\frac{G_s^{\prime }}{\rho }+[\frac{E^2}{{\mathcal{A}}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}+\frac{eB\ell }{\hslash }-\frac{e^2{B}^2{\rho }^2}{4{\hslash }^2}-\frac{{\ell }^2}{{\rho }^2}+\frac{(1-\mathcal{B})s\ell }{\mathcal{B}{\rho }^2}+\frac{(1+\mathcal{B})eBs}{2\hslash \mathcal{B}}\hfill \\ \hfill & -\frac{{\mathcal{A}}^{\prime }Es\left[2\hslash \ell -eB{\rho }^2\right]}{2\hslash \rho \mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}-\frac{1}{2{\mathcal{B}}^2}{\left(\frac{{\mathcal{A}}^{\prime }E}{\mathcal{A}\left[E+\mathcal{A}m{c}^2\right]}-\frac{{\mathcal{B}}^{\prime }}{\mathcal{B}}\right)}^{\prime }-{\left(\frac{{\mathcal{A}}^{\prime }E}{2\mathcal{A}\mathcal{B}\left[E+\mathcal{A}m{c}^2\right]}-\frac{{\mathcal{B}}^{\prime }}{2{\mathcal{B}}^2}\right)}^2\hfill \\ \hfill & -\frac{2-\mathcal{B}}{2\mathcal{B}\rho }\left(\frac{{\mathcal{A}}^{\prime }E}{\mathcal{A}\left[E+\mathcal{A}m{c}^2\right]}-\frac{{\mathcal{B}}^{\prime }}{\mathcal{B}}\right)-\frac{{\mathcal{B}}^{\prime }(1-\mathcal{B})}{{\mathcal{B}}^2\rho }-\frac{{(1-\mathcal{B})}^4}{4{\mathcal{B}}^2{\rho }^2}]{\mathcal{B}}^2{G}_s=0.\hfill \end{array}$$

(17)

It is clear that Equation (17) would not be easy to solve exactly. However, we can clearly see now which terms of the equation could safely be dropped out even without making any prior assumption on the orders of magnitude of the various physical quantities involved, such as the energy of the particle and the strengths of the gravitational and magnetic interaction terms. In Section 4, we shall come back to this point to discuss the order of magnitude of the correction that would have been brought to our final result for the energy levels of the particle if we had kept those extra terms.

Now, while Equation (17) seems analytically challenging to solve even when keeping only the dominant terms, the fact that the metric components ${\mathcal{A}}^2$ and ${\mathcal{B}}^2$ of the interior Schwarzschild solution depend on ${\rho }^2$ instead of $\rho$ will greatly simplify our task as we shall see shortly. For a spherical body of mass M, of radius R and of uniform density, the interior Schwarzschild solution describing the gravitational field inside the body is given by the metric (1), with [37]

$$\mathcal{A}\left(\rho \right)=\frac{3}{2}{\left(1-\frac{r_s}{R}\right)}^{\frac{1}{2}}-\frac{1}{2}{\left(1-\frac{r_s{\rho }^2}{R^3}\right)}^{\frac{1}{2}},\mathcal{B}\left(\rho \right)={\left(1-\frac{r_s{\rho }^2}{R^3}\right)}^{-1/2},$$

(18)

where $r_s=2GM/c^2$ is the Schwarzschild radius of the massive body. For convenience, we introduce the dimensionless parameter $\eta =\frac{3}{2}{(1-\frac{r_s}{R})}^{\frac{1}{2}}$ and the inverse length squared $\lambda =\frac{r_s}{R^3}$. Please note that by choosing this metric, we have ignored the possible geometric effect of the magnetic field on the spacetime. We will come back to this point in Section 4 as well.

Let us now show the convergence of the square root in our ansatz (14). Please note that for very small values of $\rho$, the functions (18) approximate to $\mathcal{A}\sim \eta -\frac{1}{2}+\frac{\lambda }{4}{\rho }^2$ and $\mathcal{B}\sim 1+\frac{\lambda }{2}{\rho }^2$, from which we see that the leading-order $\rho$-dependent term in the argument of the square root in Equation (14) is proportional to ${\rho }^2$. This establishes thus the convergence of that square root for very small $\rho$. On the other hand, for $\rho \ge R$, the functions (18) lead, by continuity, to the metric components of the exterior Schwarzschild solution, for we have then $\mathcal{A}={\mathcal{B}}^{-1}={(1-\frac{r_s}{\rho })}^{\frac{1}{2}}$. Therefore, for very large values of $\rho$, the functions $\mathcal{A}$ and $\mathcal{B}$ in the argument of the square root in our ansatz (14) approximate to $\mathcal{A}\sim 1-\frac{r_s}{2\rho }$ and $\mathcal{B}\sim 1+\frac{r_s}{2\rho }$. This implies that the leading-order $\rho$-dependent term in the argument of the square root is proportional to $1/\rho$. This establishes thus the convergence of that square root for large values of $\rho$ as well, whence we deduce the convergence of that square root for all values of $\rho$.

Next, we use the functions (18) to expand $1/\mathcal{A}$ and $1/{\mathcal{A}}^2$, up to the first order in the parameter $\lambda$, as follows:

$$\begin{array}{cc}\hfill \frac{1}{\mathcal{A}}& \sim \frac{2}{(2\eta -1)}-\frac{\lambda {\rho }^2}{{(2\eta -1)}^2},\hfill \\ \hfill \frac{1}{{\mathcal{A}}^2}& \sim \frac{4}{{(2\eta -1)}^2}-\frac{4\lambda {\rho }^2}{{(2\eta -1)}^3}.\hfill \end{array}$$

(19)

Expanding only up to the first order in $\lambda$ is justified by the order of magnitude of $\lambda$ (even for neutron stars as we will see in Section 4). On the other hand, we also have $1/\mathcal{B}\sim 1-\frac{\lambda }{2}{\rho }^2$ and ${\mathcal{B}}^2\sim 1+\lambda {\rho }^2$. With these expansions, Equation (15) takes, up to the first order in the parameter $\lambda$, the following form:

$$\begin{array}{cc}\hfill & \frac{{\mathrm{d}}^2{G}_s}{\mathrm{d}{\rho }^2}+\frac{1}{\rho }\frac{\mathrm{d}{G}_s}{\mathrm{d}\rho }+[\frac{4E^2}{{(2\eta -1)}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}+\frac{eB(\ell +s)}{\hslash }+\frac{\lambda eBs}{4\hslash }-\frac{\lambda (2{\ell }^2+s\ell -2)}{2}\hfill \\ \hfill & -\frac{2\lambda (s\ell +1)}{(2\eta -1)[2+(2\eta -1)\frac{m{c}^2}{E}]}-\frac{{\ell }^2}{{\rho }^2}-(\frac{e^2{B}^2}{4{\hslash }^2}+\frac{\lambda m^2{c}^2}{{\hslash }^2}+\frac{8\lambda E^2[\eta -1]}{{\hslash }^2{c}^2{[2\eta -1]}^3}+\frac{\lambda eB[2\ell +s]}{2\hslash }\hfill \\ \hfill & +\frac{\lambda eBs}{\hslash (2\eta -1)[2+(2\eta -1)\frac{m{c}^2}{E}]}){\rho }^2-\frac{\lambda e^2{B}^2}{4{\hslash }^2}{\rho }^4]G_s=0.\hfill \end{array}$$

(20)

This equation has the form

$$G_s^{″}+\frac{G_s^{\prime }}{\rho }+\left(\epsilon -\frac{{\ell }^2}{{\rho }^2}-{\kappa }^2{\rho }^2-{\tau }^2{\rho }^4\right)G_s=0,$$

(21)

where a prime denotes a derivative with respect to $\rho$ and the constants $\epsilon$, ${\kappa }^2$ and ${\tau }^2$ are straightforwardly read off from Equation (20). Setting $G_s={\chi }_s/\sqrt{\rho }$ in Equation (21), the latter takes the form

$${\chi }_s^{″}+\left(\epsilon -\frac{{\ell }^2-\frac{1}{4}}{{\rho }^2}-{\kappa }^2{\rho }^2-{\tau }^2{\rho }^4\right){\chi }_s=0.$$

(22)

This is a Schrödinger equation with a centrifugal barrier term and an isotropic quartic anharmonic oscillator potential term. Please note that (i) in the case of a weak gravitational interaction compared to the magnetic one, i.e., for $\lambda \ll \left|e\right|B/\hslash$, and/or (ii) in the case of charged fermions moving very close to the center of the spherical mass, i.e., for $\rho \ll R$, the purely quartic term ${\tau }^2{\rho }^4$ can be dropped out. In that case, Equation (22) reduces to the one that has already been solved exactly in Ref. [17], where the quantized eigenvalues ${\epsilon }_{n,\ell }$ have been extracted. We will come back to these special cases in Section 3.1. Here, we will extract the eigenvalues from Equation (22) by keeping the quartic term.

It is shown in Ref. [38] that even though Equation (22) is not solvable exactly, one may still extract an approximate analytical expression for its eigenvalues using the Jeffreys–Wentzel–Kramers–Brillouin (JWKB) approximation. To the fourth order in the approximation, the result reads [38]

$${\epsilon }_{n,\ell ,s}={\left(2\tau \right)}^{\frac{2}{3}}\sum _{k=0}^8{\left[\frac{3\sqrt{\pi }\Gamma \left(\frac{3}{4}\right)}{\sqrt{2}\Gamma \left(\frac{1}{4}\right)}\left(2n+1+\ell \right)\right]}^{\frac{4-2k}{3}}{N}_k\left(\kappa ,\tau ,\ell \right),$$

(23)

where the explicit expressions of the nine terms $N_k(\kappa ,\tau ,\ell )$ are given in Appendix A. It is clear from this expression that the usual flat-space Landau levels of a charged fermion moving under the influence of a uniform magnetic field are dramatically altered. In fact, instead of the splitting of the levels that occurs due to a weak gravitational field outside [17] or inside a spherical mass [19], what happens here is a complete redistribution of the levels. Both the separation of the latter and their strengths are redefined for each principle quantum number n.

For the sake of illustration and simplicity, we display here explicitly the expression we obtain when keeping only the first term of the sum in Equation (23). By substituting the constants $\epsilon$, $\kappa$ and $\tau$, as well as the term $N_0=1$ (as given in Appendix A) into Equation (23), we find

$$\begin{array}{cc}\hfill E_{n,\ell ,s}& =\left(\eta -\frac{1}{2}\right)\{m^2{c}^4-\hslash c^{2}eB\left(\ell +s+\frac{\lambda s}{4}\right)+\frac{\lambda {\hslash }^2{c}^2(2{\ell }^2+s\ell -2)}{2}\hfill \\ \hfill & +\frac{2\lambda {\hslash }^2{c}^2(s\ell +1)}{(2\eta -1)[2+(2\eta -1)\frac{m{c}^2}{E}]}+{\hslash }^2{c}^2{\left[\left(2n+1+\ell \right)\frac{3\Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{1}{4}\right)}\right]}^{\frac{4}{3}}{\left(\frac{\pi eB\sqrt{\lambda }}{2\hslash }\right)}^{\frac{2}{3}}{\}}^{\frac{1}{2}}.\hfill \end{array}$$

(24)

This does, in fact, not look at all like the familiar Landau levels. It must be noted here that this formula does not reduce to the flat-space result when setting $\lambda =0$ because the expansion in terms of the JWKB integrals is about the pure quartic oscillator levels, not the other way around [38]. Therefore, to achieve a high accuracy in the Formula (23), ${\kappa }^2$ should not be larger than ${\tau }^2$. To properly recover the flat-space case, one needs to set $\lambda =0$ in Equation (20). Similarly, when the gravitational interaction is weaker than the magnetic interaction, one needs to extract the corresponding Landau levels by starting from Equation (20) as we shall do now.

3.1. Weak Quartic Term

When the gravitational interaction is weak compared to the magnetic interaction and/or when the charged particles move very closely to the center of the spherical mass, the quartic term ${\tau }^2{\rho }^4$ in Equation (21) can safely be dropped out so that the equation takes the form,

$$G_s^{″}+\frac{G_s^{\prime }}{\rho }+\left(\epsilon -\frac{{\ell }^2}{{\rho }^2}-{\kappa }^2{\rho }^2\right)G_s=0.$$

(25)

The solution to this equation is easily found to be given in terms of the confluent hypergeometric function as follows (see Ref. [17] for the details of the derivation):

$$G_s\left(\rho \right)=C{\rho }^{\ell }{e}^{-\frac{\left|\kappa \right|}{2}{\rho }^2}{}_1{F}_1\left(\frac{1+\ell }{2}-\frac{\epsilon }{4\left|\kappa \right|};\ell +1;\left|\kappa \right|{\rho }^2\right).$$

(26)

This solution is one of the two independent solutions of Equation (20) that is finite at the origin $\rho =0$. C is one of the constants of integration, and ${}_1{F}_1(a;b;\xi )$ is the confluent hypergeometric function of a variable $\xi$ for any values of the parameters a and b, except when b is a negative integer for which case the function has a simple pole [39]. Therefore, to guarantee a square-integrable character for the wavefunction, the following condition must be imposed for an arbitrary positive integer n:

$$\frac{1+\ell }{2}-\frac{\epsilon }{4\left|\kappa \right|}=-n,$$

(27)

for which case, the confluent hypergeometric function becomes indeed a polynomial of finite degree n. By substituting into this condition, the values of $\epsilon$ and $\kappa$ as read off from Equation (20), we arrive at the following condition involving the energy E of the charged fermion:

$$\begin{array}{cc}\hfill & \frac{4E^2}{{(2\eta -1)}^2{\hslash }^2{c}^2}-\frac{m^2{c}^2}{{\hslash }^2}+\frac{eB(\ell +s)}{\hslash }+\frac{\lambda eBs}{4\hslash }-\frac{\lambda (2{\ell }^2+s\ell -2)}{2}\hfill \\ \hfill & -\frac{2\lambda (s\ell +1)}{(2\eta -1)[2+(2\eta -1)\frac{m{c}^2}{E}]}=(2n+1+\ell )[\frac{e^2{B}^2}{{\hslash }^2}+\frac{4\lambda m^2{c}^2}{{\hslash }^2}+\frac{32\lambda E^2(\eta -1)}{{\hslash }^2{c}^2{(2\eta -1)}^3}\hfill \\ \hfill & +\frac{2\lambda eB(2\ell +s)}{\hslash }+\frac{4\lambda eBs}{\hslash (2\eta -1)(2+[2\eta -1]\frac{m{c}^2}{E})}{]}^{\frac{1}{2}}.\hfill \end{array}$$

(28)

We may solve this equation for E by keeping again only the first-order terms in the parameter $\lambda$, to arrive at the following quantization condition on the energy E of the particle:

$$\begin{array}{cc}\hfill E_{n,\ell ,s}& =\hslash c\left(\eta -\frac{1}{2}\right)\{\frac{m^2{c}^2}{{\hslash }^2}-\frac{eB(\ell +s)}{\hslash }-\frac{\lambda eBs}{4\hslash }+\frac{\lambda (2{\ell }^2+s\ell -2)}{2}\hfill \\ \hfill & +\frac{\lambda (s\ell +1)}{(2\eta -1)}\frac{\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}{1+\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}+(2n+1+\ell )[\frac{e^2{B}^2}{{\hslash }^2}+\frac{4\lambda m^2{c}^2}{{\hslash }^2}\hfill \\ \hfill & +\frac{2\lambda eB(2\ell +s)}{\hslash }+(\eta -1)\frac{8\lambda (m^2{c}^2+[2n+1-s]\hslash eB)}{{\hslash }^2(2\eta -1)}\hfill \\ \hfill & +\frac{2\lambda eBs}{\hslash (2\eta -1)}\frac{\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}{1+\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}{]}^{\frac{1}{2}}{\}}^{\frac{1}{2}}.\hfill \end{array}$$

(29)

We readily note that for $M=0$ (i.e., by setting $\eta =\frac{3}{2}$ and $\lambda =0$), the result (29) reduces to the usual quantized energy levels of a charged fermion inside a uniform magnetic field in Minkowski spacetime: $E_{n,s}=[m^2{c}^4+\hslash c^2(2n+1-s\frac{e}{\left|e\right|}){\left|e\right|B]}^{\frac{1}{2}}$. In addition, we also note, by setting $B=0$, that even in the absence of the magnetic field the gravitational field induces quantized energy levels. Furthermore, we neatly recognize a general-relativistic correction to the pure relativistic Landau levels arising from the magnetic field. Such a correction arises from the factor $(\eta -\frac{1}{2})$ multiplying the right-hand side.

We are now going to discuss the ways to overcome the limitations we encountered in this section and what modifications would be brought to the results we obtained so far.

4. Refinements and Prospective Extensions

Two limitations, we could not avoid when deriving our results in the previous section, were (i) the fact that we considered only the possible motion of the particle along the equator inside the massive sphere and (ii) the fact that we ignored the geometric effect that might be caused to spacetime by the magnetic field. The two other limitations we imposed on our analysis were (iii) the uniform mass density we assumed for the interior of the spherical mass and (iv) the first-order approximation in $\lambda$ we limited ourselves to when writing Equation (20) to be able to solve the latter exactly. Our aim in this section is to address these limitations and argue that they do not actually alter much the actual physical conclusions we would have reached had we not adopted those various approximations.

The first limitation consists of the restricted number of degrees of freedom we imposed on the charged fermion. By restricting the motion of the particle to the equatorial plane inside the massive sphere, we implicitly assumed the particle’s motion to be confined to such a plane due to the collective interaction of the particles among themselves. However, it is actually possible, for the same reason, to also consider the possibility of motion along a plane other than the equatorial plane. In fact, the vertical pull of gravity is then balanced by these other interactions so that only the horizontal variation of the gravitational potential matters. For such a case, however, we need to switch entirely to cylindrical coordinates. The interior Schwarzschild solution needs to be expressed in such coordinates and we need to set $r^2={\rho }^2+z_0^2$ for some vertical distance $z_0$ from the equator in the resulting differential equation that would be similar in form to Equation (10). As the differential operator ${\partial }_{\theta }$ would be replaced by $z_0{\partial }_{\rho }$ and the operator ${\partial }_{\theta }^2$ would be replaced by $z_0^2{\partial }_{\rho }^2-(\rho +z_0){\partial }_{\rho }$, the equation would still be a second-order differential equation in the single variable $\rho$. A first-order approximation in the parameter $\lambda$ might then be adopted. Although such an equation would certainly not be as easy to solve as Equation (20), we expect the energy of the particle to be only shifted by a constant while obeying a quantization condition similar to our results (23) and (29). A further refinement consists then of considering the case ${\partial }_z\Psi \ne 0$ and a non-constant z. Unfortunately, the fact that two variables would then be involved, no simple procedure could be expected. A separate future work should specifically be devoted to solving that case using cylindrical coordinates. Indeed, the problem of solving our very general Equation (10) without restricting the latter to the equatorial plane might be very important for other astrophysical applications besides neutron stars. We believe that the problem might very well be tackled by adapting to our case the existing numerical methods that have already proved fruitful elsewhere for solving such equations (see, e.g., the recent work [40] and the references therein).

The second interesting refinement that should be discussed here is the possibility of adding to the metric components, given by Equation (18), the spacetime curvature caused by the magnetic field. The geometric effect of combining the external gravitational field sourced by a spherical mass with a uniform magnetic field B gives rise to the so-called Schwarzschild–Melvin spacetime [41]. Therefore, an analogue modification to the interior Schwarzschild metric components given by Equation (18) is also expected when including the spacetime curvature caused by B. However, since the corrections brought to the familiar exterior Schwarzschild solution by the magnetic field are of the order of ${ϵ}_{0}G{B}^2/c^2$, where ${ϵ}_0$ is the vacuum permittivity constant [42], we expect the correction that would be brought to the parameter $\lambda$ in our result (29) to be of the order of ∼${ϵ}_{0}G{B}^2/c^2$ as well. For a neutron star of mass ∼$1.4M_{\odot }$ and of radius 10 km, the Schwarzschild radius is about ∼$4.1$ km for which we find that $\lambda \sim 0.004$ km${}^{-2}$. On the other hand, for such a neutron star with an interior magnetic field as strong as ∼${10}^{17}$ G, we have ${ϵ}_{0}G{B}^2/c^2\sim 7\times {10}^{-9}$ km${}^{-2}$. The expected correction is thus quite insignificant.

The third limitation arises from having discarded in Equation (20) terms of the second order and higher in the parameter $\lambda$. We can justify now such an approximation by the fact that, as we just saw for a typical neutron star, $\lambda \sim 4\times {10}^{-3}$ km ${}^{-2}$. This entails that if we had kept those higher order terms the expected correction to be brought to the quantized energy levels (29) can be obtained using time-independent perturbation theory in the manner worked out in Refs. [17,18] for the case of Landau levels in the exterior Schwazschild solution.

Another refinement that should be taken into account is the fact that matter density inside neutron stars is not really uniform as we have assumed it in our calculations in the previous section. The interior matter density increases with increasing depth. However, the core of a neutron star makes up the largest part of the star and may be subdivided into outer and inner parts, each characterized by slowly varying mass densities with depth [43]. For this reason, although considering a uniform mass density is only a simplifying model for the interior of a neutron star, we expect that the conclusions that would be derived by assuming a slowly varying $\rho$-dependent mass density would not differ much from those obtained here. Unfortunately, an analysis based on the interior Schwarzschild solution with a $\rho$-dependent mass density is beyond the scope of the present paper and should be considered in future works as well. It must also be noted here that numerical methods would certainly be beneficial for such a task.

The last two extensions of this work we would like to mention here is the possibility of studying the effect of the interior Schwarzschild solution on (i) charged fermions when the magnetic field is nonuniform and (ii) on the neutron star matter when the latter is taken to be a superfluid [44]. For the former case, one needs to use again numerical methods (as in Ref. [45]) to solve the resulting equation, whereas for the latter case, one needs to consider, instead of the Dirac equation, the Gross–Pitaevskii equation coupled to the magnetic field inside the interior Schwarzschild solution.

5. Applications

Despite the limitations we discussed above, it is still interesting to apply our results to some systems that rely on the Landau energy levels of charged fermions. We will first examine the charged fermions inside a neutron star, and then we examine charged fermions inside a laboratory-made spherical mass.

5.1. Inside a Neutron Star

We propose to use here our results (23) and (29) to examine the effect of gravity on the magnetization $\mathcal{M}$ inside a neutron star. The magnetization $\mathcal{M}$ is given by the following general formula [9]:

$$\mathcal{M}=\sum _{i=e,\mu ,p}\left(\frac{\partial {ϵ}_i}{\partial B}-{\mu }_f^i\frac{\partial n_i}{\partial B}\right),$$

(30)

where the summation is over the electrons, muons and protons, respectively; ${ϵ}_i$ are the energy densities of those charged fermions and $n_i$ are their number densities. The magnetic field B in this formula is the original seed field (magnetizing field) that is thus gradually modified by the magnetization into $H=B+4\pi \mathcal{M}$, in CGS units [46]. A strong magnetization has been suggested by some authors to play a role even in giving rise to the original (primal) magnetic field itself [47]. However, despite considerable research on the topic, there is in the literature no consensus yet on the origin of such strong seed magnetic fields at the surface nor in the core of such compact stars. For a quick guide through the various proposals put forward in the literature, see, e.g., the reviews [48,49].

Considering the magnetizing field B to be constant in time and uniform over the radius of the star, as we do in this paper, is also only an approximation. However, given that the relaxation times of the star’s magnetic field decay is of the order of ∼${10}^3$ years [50], our present assumption of a constant magnetic field fits amply within the early stages of the magnetic evolution of the star. Similarly, assuming a slow variation of the magnetic field with distance inside highly magnetized stars as presented, for example, in Ref. [51] is still consistent with our present purpose of providing a preliminary study on the effect of gravity on the magnetic properties of such compact stellar objects.

First, recall that when ignoring the contribution of the gravitational field, the quantities ${ϵ}_i$ and $n_i$ are given by [9]

$$\begin{array}{cc}\hfill {ϵ}_i& =\frac{\left|e\right|B}{4{\pi }^2}\sum _s\sum _{n=0}^{n_{max}}\left({\mu }_f^i{k}_{f,n,s}^i+{\tilde{m}}_{n,s}^{i2}ln\left|\frac{{\mu }_f^i+k_{f,n,s}^i}{{\tilde{m}}_{n,s}^i}\right|\right),\hfill \\ \hfill n_i& =\frac{\left|e\right|B}{2{\pi }^2}\sum _s\sum _{n=0}^{n_{max}}{k}_{f,n,s}^i,\hfill \end{array}$$

(31)

where

$$\begin{array}{cc}\hfill {\tilde{m}}_{n,s}^{i2}& =m^{i2}+\left(2n+1-\frac{e}{\left|e\right|}s\right)\left|e\right|B,\hfill \\ \hfill k_{f,n,s}^{i2}& ={\mu }_f^{i2}-{\tilde{m}}_{n,s}^{i2}.\hfill \end{array}$$

(32)

The Fermi energies ${\mu }_f^i$ are fixed by the chemical potentials and the maximum integer $n_{max}$ in the summations represents the integer preceding the value of n for which $k_{f,n,s}^{i2}$ becomes negative. As the equations are more involved, we will work in this subsection in the natural units $\hslash =c=1$. Formula (30) then gives [9]

$$\mathcal{M}=\sum _{i=e,\mu ,p}\left[\frac{{ϵ}_i-{\mu }_f^i{n}_i}{B}+\frac{B}{2{\pi }^2}\sum _s\sum _{n=0}^{n_{max}}\left(n+\frac{1}{2}-\frac{s}{2}\right)ln\left|\frac{{\mu }_f^i+k_{f,n,s}^i}{{\tilde{m}}_{n,s}^i}\right|\right].$$

(33)

We will adapt these general formulas to our case. Since the formulas to be developed here are more relevant for very strong magnetic fields, we also extract from Equation (29) the following approximation for the quantized energy levels:

$$\begin{array}{cc}\hfill E_{n,\ell ,s}& =\left(\eta -\frac{1}{2}\right)\{m^2+\left(2n+2\ell +1-\frac{e}{\left|e\right|}s\right)\left|e\right|B-\frac{\lambda eBs}{4\hslash }+\frac{\lambda (2{\ell }^2+s\ell -2)}{2}+\frac{\lambda (s\ell +1)}{2\eta -1}\hfill \\ \hfill & +(2n+1+\ell )\left[\lambda (2\ell +s)+(\eta -1)\frac{4\lambda (2n+1-s)}{2\eta -1}+\frac{\lambda s}{2\eta -1}\right]{\}}^{\frac{1}{2}}.\hfill \end{array}$$

(34)

Please note that in the absence of gravity, i.e., when $\lambda =0$, this formula reduces to the first identity in Equation (32). Let us then start with the consequences of this formula which is valid either (i) for weak gravitational fields and/or (ii) for fermions very close to the center of the star. Since the terms proportional to $\lambda$ in this expression de not contain B, it is straightforward to compute the B-derivatives in Equation (30). The result for the magnetization remains thus exactly the same as in Equation (33), but with ${\tilde{m}}_{n,s}^{i2}$ replaced by $E_{n,\ell ,s}$ as given by Equation (34):

$$\mathcal{M}=\sum _{i=e,\mu ,p}\left[\frac{{ϵ}_i-{\mu }_f^i{n}_i}{B}+\frac{B(2\eta -1)}{4{\pi }^2}\sum _s\sum _{n=0}^{n_{max}}\left(n+\frac{1}{2}-\frac{s}{2}\right)ln\left|\frac{{\mu }_f^i+k_{f,n,s}^i}{E_{n,\ell ,s}^i}\right|\right].$$

(35)

We neatly see the general-relativistic correction factored out in the second term of the sum. More important, however, is the fact that the general expression of the magnetization remains unaltered in this case.

On the other hand, for the case of a non-negligible gravitational field the quartic term in Equation (20) should be taken into account, leading to the result (24) for $E_{n,\ell ,s}$ at the first-order term of the JWKB approximation. As the Landau levels are completely destroyed, we see that even at this order in the approximation the magnetization is completely different from the expression (35). Actually, even Formulas (31) and (32) do not hold anymore, for such formulas were derived by assuming the linear dependence (32) of ${\tilde{m}}_{n,s}^{i2}$ on the integer n. As expression (24) of $E_{n,\ell ,s}$ in terms of n displays no such linearity, and since that expression (24) is itself only an approximation for the highly nonlinear full expression (23), we conclude that the Landau-levels-based formalism for computing the magnetization $\mathcal{M}$ of neutron stars does not generally apply to the charged plasma in the core of highly dense neutron stars. The condition for the applicability of the formalism is to have the quartic anharmonic oscillator potential term in Equation (22) negligible, i.e., to have $\lambda \ll \left|e\right|B/\hslash$ or to apply the formalism only to fermions close to the vicinity of the star’s center. In other words, we conclude that only for a magnetic neutron star of mass M, of radius R and with a magnetic field B such that $\left|e\right|B{R}^3{c}^2\gg GM\hslash$, that the usual formalism generally applies to the charged matter inside the core of the star.

Magnetization in neutron stars affects both the equation of state of the stars’ matter and the stars’ structure and properties, such as their masses and transport properties. Indeed, both the pressure anisotropy and the softening of the equation of state of the stars matter is believed to be related to magnetization and the distribution of the charged fermions as they occupy the Landau levels. The equation of state determines, in turn, the mass-radius relation of neutron stars (see, e.g., the recent work [52] and the references therein).

It is believed that the net magnetization in neutron stars remains negligible until around ${10}^{16}$ G, beyond which it starts increasing and becomes oscillatory under the de Haas-van Alphen effect. Since we just saw that the Landau levels become dramatically altered only for non-negligible gravitational fields compared to the magnetic interaction, we conclude that such oscillations will not be much affected in typical stars. For weaker magnetic fields, however, the general-relativistic correction we found in Equation (35) for the magnetization might leave an observable signature on the equation of state of the star’s matter and, hence, on the mass-radius relation of the star as well. More work dedicated specifically to such an interesting investigation is required though.

5.2. At a Laboratory Level

Another application of our results is to consider physical systems at the laboratory level. In this case, an interesting system would be the free electrons moving inside a massive metallic sphere put inside a uniform magnetic field. For this case, the gravitational contribution is much weaker than the magnetic interaction, so Equation (29) is what we need to apply here. Expanding that expression up to the first order in $\lambda$, we obtain

$$\begin{array}{cc}\hfill E_{n,\ell ,s}& \approx \left(\eta -\frac{1}{2}\right)\{m^2{c}^4+\hslash c^2\left(2n+2\ell +1-\frac{e}{\left|e\right|}s\right)\left|e\right|B-\frac{\hslash c^2\lambda eBs}{4}+\frac{{\hslash }^2{c}^2\lambda (2{\ell }^2+s\ell -2)}{2}\hfill \\ \hfill & +\frac{\lambda {\hslash }^2{c}^2(s\ell +1)}{(2\eta -1)}\frac{\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}{1+\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}+(2n+1+\ell ){\hslash }^2{c}^2\lambda [\frac{2m^2{c}^2}{\hslash eB}+2\ell +s\hfill \\ \hfill & +(\eta -1)\frac{4(m^2{c}^2+[2n+1-s]\hslash eB)}{\hslash eB(2\eta -1)}+\frac{s}{2\eta -1}\frac{\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}{1+\sqrt{1+(2n+1-s)\frac{\hslash eB}{m^2{c}^2}}}]{\}}^{\frac{1}{2}}.\hfill \end{array}$$

(36)

These are the energy levels of the electrons inside the massive sphere. The $\lambda$-terms are the corrections brought to the familiar levels observed in the absence of gravity. This expression is valid for both relativistic and non-relativistic fermions, and it reduces to the familiar relativistic Landau levels in the absence of gravity. This result offers novel ways for exploiting the Landau levels of fermions in the presence of gravity at the laboratory level.

6. Summary

We have obtained the Dirac equation for a charged fermion in a general static and spherically symmetric curved spacetime in the presence of a static and uniform magnetic field. We applied the general equation we obtained to the case of fermions in the interior Schwarzschild solution describing the gravitational field inside a massive sphere of uniform density, and we derived the quantized energy levels of the charged particles. We found that our result reduces to the familiar flat-spacetime relativistic Landau levels of charged fermions inside a uniform magnetic field only when the gravitational interaction is weaker than the magnetic interaction and/or when the charged particles move very closely to the center of the spherical mass.

We then discussed, based on these results, the consequences on the physics of neutron stars. We found that the magnetization of the core of the latter would dramatically be altered for an extremely high gravitational field compared to the magnetic interaction and/or when not focusing solely on the fermions in the vicinity of the center of the star. We arrived at the conclusion that in the general case, the usual formalism for extracting the magnetization of neutron stars applies only when the latter satisfy the condition $\left|e\right|B{R}^3{c}^2\gg GM\hslash$. Although we arrived at such a conclusion by relying on a first-order approximation in the gravitational parameter $\lambda$, the generality of our results is not affected. In fact, such a first-order approximation remains valid for any relative strengths of the gravitational and magnetic interactions experienced by the charged particles inside typical neutron stars. Any extra correction to those levels would be obtained simply using time-independent perturbation theory.

In Section 4, we outlined a few refinements that might be brought to our model and the prospective outlook on future improvements that will allow us to go beyond the limitations we imposed on our present study. However, even with the simple model dealt with here, we have already glimpsed with new insights into the contribution of gravity in shaping the magnetization of highly compact neutron stars and into novel ways of exploiting the dynamics of charged fermions in the presence of gravity at the laboratory level.