Zeeman Splitting of Torsional Oscillation Frequencies of Magnetars

Abstract

Magnetars form a special class of neutron stars possessing superstrong magnetic fields and demonstrating power flares triggered by these fields. Observations of such flares reveal the presence of quasi-periodic oscillations (QPOs) at certain frequencies; they are thought to be excited in the flares. QPOs carry potentially important information on magnetar structure, magnetic field, and mechanisms of magnetar activity. We calculate frequencies of torsional (magneto-elastic) oscillations of the magnetar crust treating the magnetic field effects in the first order of perturbation theory. The theory predicts the splitting of non-magnetic oscillation frequencies into Zeeman components. Zeeman splitting of the torsional oscillation spectrum of magnetars was suggested, clearly described and estimated by Shaisultanov and Eichler (2009), but their work has not been given considerable attention. To extend it, we suggest the technique of calculating oscillation frequencies, including Zeeman splitting at not too strong magnetic fields for arbitrary magnetic field configuration. Zeeman splitting enriches the oscillation spectrum and simplifies the theoretical interpretation of observations. We calculate several low-frequency oscillations of magnetars with a pure dipole magnetic field in the crust. The results qualitatively agree with the low-frequency QPOs detected in the hyperflare of SGR 1806–20 and in the giant flare of SGR 1900+14.

1. Introduction

It is well known that neutron stars are the most compact stars and contain superdense matter in their interiors (e.g., [1]). These stars have massive liquid cores and thin envelopes. The core contains matter of a supranuclear density; its equation of state (EOS) and other properties are still not certain and remain the fundamental problem of physics and astrophysics. The envelope above the core is also important; its thickness is ∼1 km, and its mass is ∼$0.01M_{\odot }$. The envelope consists mostly of electrons and atomic nuclei. In addition, at densities, $\rho$, larger than the neutron drip density, ${\rho }_{\mathrm{drip}}\approx 4.3\times {10}^{11}$ g cm${}^{-3}$, quasi-free neutrons appear. The atomic nuclei in the envelope form Coulomb crystals (e.g., [2]), which melt near the very surface layers of the star. The solidified envelope is called the crust; one distinguishes the outer ($\rho <{\rho }_{\mathrm{drip}}$) and inner ($\rho >{\rho }_{\mathrm{drip}}$) crust. The density of matter at the crust–core interface is ${\rho }_{\mathrm{cc}}\sim 1.4\times {10}^{14}$ g cm${}^{-3}$.

This paper studies torsional oscillations of the star due to the elasticity of the crystalline lattice. In a non-magnetic star, such pure torsion oscillations are confined in the crust. The foundation of the theory was laid by Hansen and Chioffi [3], Schumaker and Thorne [4], and McDermott et al. [5] in the 1980s. Later, the theory was developed in numerous publications (see, e.g., [6,7,8,9,10,11,12,13,14,15] and the references therein).

The theory attracted great attention after the discovery of quasi-periodic oscillations (QPOs) in spectra of soft-gamma repeaters (SGRs) (see [16,17,18,19,20,21]). SGRs belong to a class of magnetars. They are neutron stars possessing very strong magnetic fields $B\sim {10}^{15}$ G (e.g., [22,23,24]) and demonstrating flaring activity that was probably triggered by these fields. The activity is accompanied by the processes of enormous energy release in the form of hyperflares, superflares and ordinary flares. The QPOs have been discovered at the fading stages of such flares. Seismic activity of SGRs was predicted by Duncan [25] in 1998.

The detected QPO frequencies range from a few tens of Hz to several kHz. This is the range typical for theoretical torsion frequencies of non-magnetic neutron stars. Magnetar QPOs are separated into low-frequency (≲ a few hundred Hz) and high-frequency ones. The discovery of magnetar QPOs gave hope to developing reliable methods of exploring magnetar structure and evolution by comparing observations with elaborated seismological models.

The seismology of magnetars has been studied by many authors (e.g., Levin [26,27], Glampedakis et al. [28], Sotani et al. [29], Cerdá-Durán et al. [30], Colaiuda et al. [31], Colaiuda and Kokkotas [32,33], van Hoven and Levin [34,35], Gabler et al. [36,37,38,39,40,41], Passamonti and Lander [42], Link and van Eysden [43]). They considered different possibilities. The first idea was treating magnetar QPOs by ordinary torsion oscillations of neutron star crust nearly unaffected by magnetar fields (e.g., [8,9,10,11,12,14,15,29] ). Other models were based on magneto-elastic oscillations of the crust, where the elasticity of a crystalline lattice was combined with the elasticity of magnetic fields lines; such oscillations can spread outside the crust in the form of Alfvén waves (e.g., [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44]). In some cases, the presence of the lattice becomes fully unimportant; then, the oscillations are of purely Alfvénic type and can easily spread outside the crust.

Let us emphasize the remarkable paper by Shaisultanov and Eichler [45], who predicted the effect of Zeeman splitting of torsional oscillations in magnetar magnetic fields. Although the effect is clear and the paper is physically transparent, it has not been given much attention. It is our aim here to fill this gap, emphasize the importance of the effect, and extend the consideration of Shaisultanov and Eichler [45] (although our technique will be somewhat different, as we discuss later). Here, we focus on low-frequency magneto-elastic oscillations confined mainly to the crust, including Zeeman splitting. The paper is organized as follows.

In Section 2, we suggest calculating magneto-elastic oscillation frequencies using the first-order perturbation theory with respect to the magnetic terms. The basic equations are presented in Section 2, neglecting relativistic effects (for simplicity). Pure torsional oscillations are outlined in Section 3. The first-order perturbation approach is described in Section 4 for any $\mathit{B}$-field geometry. It demonstrates Zeeman splitting [45] of pure torsional oscillation frequencies into Zeeman components, greatly enriching the oscillation spectrum. In Section 5, we use the perturbation approach in the case of fundamental low-frequency torsional oscillations in a poloidal and axially symmetric $\mathit{B}$-field. In Section 6, we modify the equations by including relativistic effects. In Section 7, we consider the pure dipole $\mathit{B}$-field configuration in a neutron star crust. In Section 8, these results are used to sketch possible interpretations of detected low-frequency QPOs in afterglows of the hyperflare of SGR 1806–20 as well of the giant flare of SGR 1900+14. Finally, in Section 9, we summarize the results and compare our approach with those available in the literature.

2. The Approach

Consider elastic oscillations of a neutron star crust containing a magnetic field B. At B = 0, these oscillations are pure torsional and occur due to elastic properties of Coulomb crystals of atomic nuclei in the crust. In a magnetic field, Alfvén waves, owing to the elasticity of magnetic field lines, modify the elastic properties of matter and lead to combined magneto-elastic oscillations where both elasticities contribute. We will restrict ourselves to the case of not too strong field B so that the magnetic effects can be treated as a perturbation. In this case, we retain the term ‘torsional oscillation’ for brevity.

We start with a simplified problem in which we describe the star, neglecting relativistic effects; we include them later. Otherwise, we adopt standard assumptions for studying magneto-elastic oscillations based on the linearization procedure. All the quantities involved are divided into unperturbed ones and small perturbations. The non-perturbed star is taken as stationary and spherically symmetric; ordinary spherical coordinates (r, $\theta$, and $\varphi$) are most convenient for specifying position vectors $\mathit{r}$. Unperturbed quantities, like mass density $\rho$, depend only on r. The magnetic field B is assumed to be not too high and treated as a perturbation. It is split further as $\mathbf{B}=\mathit{B}+{\mathbf{B}}_1$, where $\mathit{B}$ is the static (given) magnetic field of the star, while ${\mathbf{B}}_1$ is a smaller field induced by the oscillations. The field $\mathit{B}$ leads to small stationary perturbations, which do not interfere with oscillatory perturbations in the adopted linear approximation. Here, we focus on oscillatory perturbations.

Let $\mathbf{v}(\mathit{r},t)$ be a velocity of a matter element, and $\mathbf{u}(\mathit{r},t)$ be a shift of its position with respect to equilibrium position $\mathit{r}$ in a non-perturbed star; t is time. For an oscillation mode with frequency $\omega$ in the linear approximation one has

$$\mathbf{v}(\mathit{r},t)={\mathrm{e}}^{\mathrm{i}\omega t}\mathit{v}\left(\mathit{r}\right),\mathbf{u}(\mathit{r},t)={\mathrm{e}}^{\mathrm{i}\omega t}\mathit{u}\left(\mathit{r}\right),\mathbf{v}=\mathrm{i}\dot{\mathbf{u}},{\mathbf{B}}_1(\mathit{r},t)={\mathrm{e}}^{\mathrm{i}\omega t}{\mathit{B}}_1\left(\mathit{r}\right),$$

(1)

where italic vectors $\mathit{v}$, $\mathit{u}$ and ${\mathit{B}}_1$ are (generally complex and stationary) amplitudes of perturbations; dot denotes the time derivative.

It is well known (e.g., [28]) that the torsional (magneto-elastic) oscillations can be treated as incompressionable ($\nabla ·\mathit{v}=0$; $\nabla ·\mathit{u}=0$). Then, matter elements oscillate along spheres with fixed radius r and density $\rho \left(r\right)$, with $v_r=0$ and $u_r=0$. The density perturbations are negligible. The equation of motion for a matter element reduces to the stationary equation

$$\rho {\omega }^2\mathit{u}={\mathbf{T}}_{\mu }+{\mathbf{T}}_B,$$

(2)

where the force terms ${\mathbf{T}}_{\mu }$ and ${\mathbf{T}}_B$ (minus elastic forces per unit volume) come from the elasticity of the crystalline lattice and magnetic field lines, respectively. Elastic deformations of the lattice are determined by the shear modulus $\mu \left(r\right)$ taken in the approximation of an isotropic solid. Any i-th Cartesian component of ${\mathbf{T}}_{\mu }$ is given by

$${\mathbf{T}}_{\mu i}=-\frac{\partial {\sigma }_{ik}}{\partial x_k},$$

(3)

where ${\sigma }_{ik}=\mu (\partial u_i/\partial x_k+\partial u_k/\partial x_i)$ is the shear stress tensor.

The magnetic force term ${\mathbf{T}}_B$ is obtained in the standard manner. From the Maxwell induction equation (without dissipation), in our case, one has ${\mathit{B}}_1=\mathrm{curl}(\mathit{u}\times \mathit{B})$ and

$${\mathbf{T}}_B=\frac{1}{4\pi }\mathit{B}\times \mathrm{curl}\left[\mathrm{curl}(\mathit{u}\times \mathit{B})\right].$$

(4)

Equation (2), supplemented with (3) and (4), allows one to find oscillation eigenfrequencies $\omega$ and eigenvectors $\mathit{u}\left(\mathit{r}\right)$, which we discuss later in more detail. It is important that the ${\mathbf{T}}_{\mu }$ and ${\mathbf{T}}_{\mathrm{B}}$ terms are linear in the same shift vector $\mathit{u}$. This couples the lattice and $\mathit{B}$-field elasticities. The vector $\mathit{u}$ (with $u_r=0$) plays the role of the wave function of the problem.

Let us mention one feature of Equation (2). Namely, let us multiply (2) by the complex conjugated shift vector ${\mathit{u}}^{\ast }$ and integrate it over the star. In this way, we obtain the equality

$${\omega }^2={\omega }_{\mu }^2+{\omega }_B^2,$$

(5)

where

$${\omega }_{\mu }^2=\frac{\int \mathrm{d}V{\mathit{u}}^{\ast }{\mathbf{T}}_{\mu }}{\int \mathrm{d}V\rho |u^2|},{\omega }_B^2=\frac{\int \mathrm{d}V{\mathit{u}}^{\ast }{\mathbf{T}}_B}{\int \mathrm{d}V\rho |u^2|}.$$

(6)

Therefore, the true oscillation frequency $\omega$ is formally expressed via two auxiliary frequencies, ${\omega }_{\mu }$ and ${\omega }_B$. In the absence of the magnetic field, one has $\omega ={\omega }_{\mu }$, which refers to a pure torsional oscillation due to the elasticity of the lattice. In the absence of the lattice, one gets $\omega ={\omega }_B$, which describes oscillations due to Alfύen waves. If the lattice and the magnetic field are present at once, both auxiliary frequencies depend generally on crystal and magnetic elasticities. Both integrals in the expression for ${\omega }_{\mu }^2$ contain positive definite self-conjugated operators and bilinear combinations of $\mathit{u}$ and ${\mathit{u}}^{\ast }$, which guarantees that ${\omega }_{\mu }^2>0$. It may be not true for the integral in the nominator of ${\omega }_B^2$; in principle, ${\omega }_B^2$ may be a complex number.

The decomposition of frequencies $\omega$, similar to (5), was studied earlier by Schumaker and Thorne [4] (also see, e.g., [46,47]) for pure torsional oscillations at $\mathit{B}=0$ with the aim to separate contributions of elastic shear forces and centrifugal forces in ${\omega }_{\mu }^2$. It is not a surprise that similar decomposition appears in the case of magneto-elastic oscillations.

3. Pure Torsional Oscillations

This case has been studied in many publications starting from [3,4,5] and developed further (e.g., [6,7,8,9,10,11,12,13,14,15] and references therein).

Here, we present a brief summary of these results, which we will need later. Oscillation modes are characterized by three numbers. They are: (i) multipolarity $\mathsf{\ell }=$2, 3, …; (ii) azimuthal number m (which is an integer varying from $m=-\mathsf{\ell }$ to $m=\mathsf{\ell }$); and (iii) the number $n=$0,1, 2, …of radial nodes of wave functions $\mathit{u}(r,\theta ,\varphi )$.

The functions ${\mathit{u}}_{n\mathsf{\ell }m}$ and eigenfrequencies ${\omega }_{\mu n\mathsf{\ell }}$ can be found from Equation (2) with ${\mathbf{T}}_B$ = 0. In view of the spherical symmetry of our non-perturbed star, the quantization axis z can be arbitrary. Accordingly, the eigenfrequencies do not depend on m, meaning that any frequency is degenerate; it is one and the same for (2ℓ + 1) different state vectors $\mathit{u}$ with fixed n and ℓ but different m. The consequences of this fact will be important for our analysis.

In the given formulation, the solution for $\mathit{u}$ can be presented as

$$\begin{array}{ccc}\hfill u_{\varphi }(r,\theta ,\varphi )& =& r{Y}_{n\mathsf{\ell }}\left(r\right)exp\left(\mathrm{i}m\varphi \right)\frac{\mathrm{d}{P}_{\mathsf{\ell }}^m}{\mathrm{d}\theta },\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill u_{\theta }(r,\theta ,\varphi )& =& r{Y}_{n\mathsf{\ell }}\left(r\right)exp\left(\mathrm{i}m\varphi \right)\frac{\mathrm{i}{P}_{\mathsf{\ell }}^m}{sin\theta },\hfill \end{array}$$

(8)

where $P_{\mathsf{\ell }}^m(cos\theta )$ is an associated Legendre polynomial, and $Y_{n\mathsf{\ell }}\left(r\right)$ is a convenient dimensionless radial wave function, which satisfies the equation

$$Y_{n\mathsf{\ell }}^{″}+\left(\frac{4}{r}+\frac{{\mu }^{\prime }}{\mu }\right)Y_{n\mathsf{\ell }}^{\prime }+\left[\frac{\rho }{\mu }{\omega }_{\mu n\mathsf{\ell }}^2-\frac{(\mathsf{\ell }+2)(\mathsf{\ell }-1)}{r^2}\right]Y_{n\mathsf{\ell }}=0.$$

(9)

Here, the primes mean derivatives with respect to r. Since pure torsional oscillations are localized in the crystalline matter ($r_1\le r\le r_2$), Equation (9) is solved with the boundary conditions $Y^{\prime }\left(r_1\right)=Y^{\prime }\left(r_2\right)=0$ (of zero traction). The solution gives the eigenfrequencies ${\omega }_{\mu n\mathsf{\ell }}$. It is Equation (9) which contains information on the microphysics of neutron star matter. The angular dependence of the wave vectors $\mathit{u}$ is standard (universal for all stars).

4. The Simplest Iterative Solution

Now we suggest the simplest iterative procedure to include the effects of stellar magnetic field $\mathit{B}$. We treat the shift vector ${\mathit{u}}_{n\mathsf{\ell }m}$, given by (7) and (8), as the zero-order solution of (2), and we treat the term containing ${\mathbf{T}}_B$, as a perturbation. Let us restrict ourselves by the first-order perturbation. The problem is similar to finding a first-order correction to a non-perturbed degenerate energy level of a quantum-mechanical system in case the perturbation removes degeneracy and splits the energy into a series of sublevels (e.g., [48]).

In this case, it is sufficient to use only ($2\mathsf{\ell }+1$) zero-order wave functions corresponding to a fixed non-perturbed oscillation frequency ${\omega }_{n\mathsf{\ell }m}$. Such functions are often called the initial zero-order wave functions. The set of these functions is not unique: any linear superposition

$$\tilde{\mathit{u}}=\sum _{m=-\mathsf{\ell }}^{+\mathsf{\ell }}{\alpha }_m{\mathit{u}}_{n\mathsf{\ell }m}$$

(10)

with constant transformation coefficients ${\alpha }_m$ also describes an oscillatory state with the same ${\omega }_{n\mathsf{\ell }m}$.

The iteration procedure prescribes to calculate the perturbation matrix T${}_{m{m}^{\prime }}$ (of dimension $(2\mathsf{\ell }+1)\times (2\mathsf{\ell }+1)$) on the basis of the initial wavefunctions. In our case, such matrix elements are given by

$${\mathrm{T}}_{m^{\prime }m}={\int }_{\mathrm{crust}}\mathrm{d}V{\mathit{u}}_{n\mathsf{\ell }{m}^{\prime }}^{\ast }{\mathbf{T}}_B\left({\mathit{u}}_{n\mathsf{\ell }m}\right).$$

(11)

Since we deal with zero-order wave functions, which are localized in the crust, the integration domain has to be restricted by the crust alone.

At the next stage, one should use Equation (10) and introduce the basis of ‘true’ zero-order wave functions ${\tilde{\mathit{u}}}_{\nu }$ (labeled by an index $\nu$), which diagonalize ${\mathrm{T}}_{m^{\prime }m}$. Let us denote the diagonal values by ${\tilde{\mathrm{T}}}_{\nu }$. Their number is $(2\mathsf{\ell }+1)$; some of them can be equal if the perturbation does not completely remove degeneracy. Then, the oscillation frequency for the ‘true’ state $\nu$ is

$${\omega }_{n\mathsf{\ell }\nu }^2={\omega }_{\mu n\mathsf{\ell }}^2+\frac{{\tilde{T}}_{\nu }}{{\int }_{\mathrm{crust}}\mathrm{d}V\rho \left|u^2\right|}.$$

(12)

This means the splitting of ${\omega }_{\mu n\mathsf{\ell }}$ by the magnetic field that is traditionally treated as the Zeeman effect. This equation is in accord with the more general Equation (6) in which ${\omega }_B^2$ has to be identified with the last term in (12).

This procedure of studying the magnetic splitting of torsional oscillation frequencies is computationally simple for any static magnetic field $\mathit{B}\left(\mathit{r}\right)$ configuration, being mostly reduced to calculating and diagonalizing the matrix (11). The procedure gives a much richer spectrum of eigenfrequencies than the theory of non-magnetic torsional oscillations.

The disadvantage of the procedure is evident. It is certainly restricted by the not-too-high magnetic fields—the magnetic splitting of oscillation frequencies has to be smaller than the zero-order torsion oscillation frequencies ${\omega }_{\mu n\mathsf{\ell }}$ themselves. The suggested procedure needs only the magnetic field $\mathit{B}$ located in the crust; only zero-order wavefunctions are involved. Undoubtedly, sufficiently high $\mathit{B}$ breaks down this approximation; then, penetration of Alfvén perturbations to the neutron star core and possibly to the magnetosphere may become essential. At large $\mathit{B}$, the theory should be modified, but the imprints of the degeneracy problem need to be studied anyway. Another disadvantage of the iterative scheme is that, although its first step is feasible, higher-order iterations seem too complicated to be useful.

5. Fundamental Torsional Oscillations in Poloidal Axially Symmetric Magnetic Fields

By way of illustration, we apply the formalism of Section 4 to the case of fundamental torsional oscillations in the presence of a poloidal axially symmetric magnetic field configuration. In this case

$$B_r=B_r(r,\theta ),B_{\theta }=B_{\theta }(r,\theta ),B_{\varphi }=0,$$

(13)

where z is the magnetic axis. Since $\nabla ·\mathit{B}=0$, the field components $B_r$ and $B_{\theta }$ are related by

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{B}_r\right)+\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(B_{\theta }sin\theta )=0.$$

(14)

The main problem is to calculate the matrix elements T${}_{m{m}^{\prime }}$ from Equation (11) with the wavefunctions (8) and (7), and to diagonalize the matrix. Luckily, a careful inspection of the integral in (11) shows that for the magnetic field (13), the matrix T${}_{m{m}^{\prime }}$ is diagonal on the basis of the initial zero-order wavefunctions because the integration over $\varphi$ reduces to

$${\int }_0^{2\pi }\mathrm{d}\varphi exp(\mathrm{i}m-\mathrm{i}{m}^{\prime })=2\pi {\delta }_{m{m}^{\prime }},$$

(15)

where ${\delta }_{m{m}^{\prime }}$ is the Kronecker delta. Accordingly, the initial zero-order functions are identical to the true zero-order functions. Then, the magnetically split eigenfrequencies can be labeled by $\nu =m$ and readily given by Equation (12) with ${\tilde{\mathrm{T}}}_{\nu }={\mathrm{T}}_{mm}$. This solution can be rewritten in the form

$${\omega }_{n\mathsf{\ell }m}^2={\omega }_{\mu n\mathsf{\ell }}^2+{\omega }_{Bn\mathsf{\ell }m}^2,{\omega }_{Bn\mathsf{\ell }m}^2=\frac{{\mathrm{T}}_{mm}}{{\int }_{\mathrm{crust}}\mathrm{d}V\rho \left|u^2\right|}.$$

(16)

In addition, due to the symmetry with respect to the equatorial plane ($\theta =\pi /2$), one has T${}_{mm}=$T${}_{-m-m}$. Accordingly, the splitting of Zeeman states with azimuthal numbers, m and $-m$, is equal, so that the magnetic sublevels can be labeled by $m=0$, 1,…ℓ (with $\mathsf{\ell }+1$ different sublevels in total). The sublevel $m=0$ is nondegenerate, while those with $m>0$ are still degenerate twice. The latter degeneracy can be further removed by a more complicated $\mathit{B}$-field geometry.

Moreover, we restrict ourselves by considering fundamental torsional oscillations; they are those with no radial nodes. Then, we always have $n=0$, so we drop this index for simplicity. Therefore, our zero-order oscillation frequencies ${\omega }_{\mu }={\omega }_{\mu \mathsf{\ell }}$ depend only on ℓ. At low ℓ, these frequencies are known to be small enough to explain the lowest frequencies of QPOs observed in the afterglow of magnetar flares (see Section 8 for more details). For fundamental torsional oscillations with very good accuracy, the radial wave function Y in Equations (7) and (8) is independent of r: $Y_{\mathsf{\ell }}\left(r\right)\approx Y_0$, see [46,47]. Since the numerators and denominators in Equation (6) are proportional to $Y_0^2$, this quantity just drops out of consideration.

Then, ${\omega }_{B\mathsf{\ell }m}^2$ in Equation (16) can be rewritten as

$${\omega }_{B\mathsf{\ell }m}^2=\frac{\frac{1}{4\pi }{\int }_{\mathrm{crust}}\mathrm{d}V{I}_B}{{\int }_{\mathrm{crust}}\mathrm{d}V\rho r^2\left(P^{\prime 2}+\frac{m^2}{{sin}^2\theta }{P}^2\right)},$$

(17)

where $P=P_{\mathsf{\ell }}^m\left(\theta \right)$, and

$$\begin{array}{ccc}\hfill I_B& =& -(B_{\theta }^2-B_{\theta }{B}_{\theta }^{\prime }cot\theta -B_r{B}_{\theta }cot\theta -r{B}_r{B}_{\theta ,r}cot\theta )P^{\prime 2}\hfill \\ & & -(B_{\theta }{B}_{\theta }^{\prime }+B_r{B}_{\theta }+r{B}_r{B}_{\theta ,r})P^{\prime }{P}^{″}-B_{\theta }^2{P}^{\prime }{P}^{‴}\hfill \\ & & +\left[B_{\theta }^2{P}^{\prime 2}-(B_{\theta }{B}_{\theta }^{\prime }+B_{\theta }^{2}cot\theta +2B_r{B}_r^{\prime }+B_r{B}_{\theta }+r{B}_r{B}_{\theta ,r})P{P}^{\prime }\right.\hfill \\ & & +\left.B_r(B_{\theta }^{\prime }+B_{\theta }cot\theta +r{B}_{\theta ,r}^{\prime }+r{B}_{\theta ,r}cot\theta -B_r^{″}+B_r^{\prime }cot\theta )P^2\right]\frac{m^2}{{sin}^2\theta }.\hfill \end{array}$$

(18)

Here the primes denote derivatives with respect to $\theta$. The subscript r after the comma denotes the derivative with respect to r. Equation (18) is obtained from (11) by substituting (7) and (8), taking $Y\left(r\right)=Y_0$, and integrating over $\varphi$ using (15). We have also used (14) to rearrange some terms. The denominator of (17) can be further simplified using the well known (e.g., [49]) property of associated Legendre polynomials $P_{\mathsf{\ell }}^m(cos\theta )$:

$${\int }_0^{\pi }sin\theta \mathrm{d}\theta \left(P^{\prime 2}+\frac{m^2}{{sin}^2\theta }{P}^2\right)\equiv \Xi (\mathsf{\ell },m)=\frac{2\mathsf{\ell }(\mathsf{\ell }+1)(\mathsf{\ell }+m)!}{(2\mathsf{\ell }+1)(\mathsf{\ell }-m)!}.$$

(19)

Then, we can rewrite (17) as

$${\omega }_{B\mathsf{\ell }m}^2=\frac{\frac{1}{4\pi }{\int }_{r_1}^{r_2}\mathrm{d}r{r}^2{\int }_0^{\pi }sin\theta \mathrm{d}\theta I_B}{\Xi (\mathsf{\ell },m){\int }_{r_1}^{r_2}\mathrm{d}r\rho r^4}=\frac{\frac{1}{4\pi }{\int }_{\mathrm{crust}}\mathrm{d}V{I}_B}{\Xi (\mathsf{\ell },m){\int }_{\mathrm{crust}}\mathrm{d}V\rho r^2}.$$

(20)

The integral in the denominator is exactly the same as for pure torsional oscillations at $\mathit{B}=0$. Thus, the problem of calculating Zeeman splitting of the zero-order oscillation frequency reduces to the integration of the known function $I_B$ over r and $\theta$ in the nominator. This integral is determined by the magnetic field configuration (13) and can be taken for any assumed poloidal axially symmetric $B_r$ and $B_{\theta }$. Otherwise, the integration is insensitive to microphysics in the neutron star crust. Although we have studied fundamental torsional modes ($n=0$), the results are easily generalized to modes with radial nodes ($n>0$).

According to the above expressions, ${\omega }_{B\mathsf{\ell }m}^2$ is quadratic in B. Then, the true oscillation frequency has the form $\omega \left(B\right)=\sqrt{{\omega }^2\left(0\right)+\alpha B^2}$, where B is a characteristic $\mathit{B}$-field value, and $\alpha$ is some constant. Since our consideration is strictly valid at ${\omega }^2\left(0\right)\gg \alpha B^2$, one gets $\omega \left(B\right)\approx \omega \left(0\right)+\frac{1}{2}\alpha B^2/\omega \left(0\right)$, meaning a small correction to $\omega \left(0\right)$ quadratic in B in the guaranteed applicability range. However, taking into account the generic rule (5), we will often retain the square root, which can formally exceed the applicability limit; at ${\omega }^2\left(0\right)\ll \alpha B^2$, one would get $\omega \left(B\right)\propto B$ although the validity of this expression is unclear.

Similar square-root expression for $\omega \left(B\right)$ has been used in various calculations of magneto-elastic oscillation frequencies (e.g., [29,50]), which neglected Zeeman splitting (note that in some cases, the dependence $\alpha B^2$ was declared to be inaccurate at large B, e.g., [38]). The factor $\alpha$ has often been used as a convenient parameter to fit the calculated $\omega \left(B\right)$. Our approach gives an exact prescription on how to determine $\alpha$, including ($m\ne 0$) or disregarding ($m=0$) Zeeman splitting at not too high B.

Let us mention that in the pioneering work on the Zeeman effect in magnetars, Shaisultanov and Eichler [45] used an elegant formalism of spherical vectors (e.g., [51]) to calculate (in our notations) T${}_{mm}$. They assumed a constant $\mathit{B}$ directed along the z-axis and a constant $\mu \left(r\right)/\rho \left(r\right)$ in the magnetar crust, and demonstrated the main features of Zeeman splitting in magnetars. In the case of more sophisticated magnetic field configurations, using the formalism of spherical vectors would be more complicated.

6. Relativistic Effects

So far, we have neglected relativistic effects (to make our theoretical sketch physically transparent), but for neutron stars, these effects are quantitatively important.

The first effect is due to the Special Theory of Relativity. Since the neutron star matter is essentially relativistic, one should replace (e.g., [4]) the mass density $\rho$ in the equations of motion by the inertial mass density $\rho +\mathrm{P}/c^2$, where P is the pressure.

Other effects are those due to general relativity (GR). The spacetime in and around a neutron star is curved. For our problem, it would be sufficient to use the relativistic Cowling approximation and the standard metric for a non-perturbed spherically symmetric star

$$\mathrm{d}{s}^2=-exp(2\Phi \left(r\right))\mathrm{d}{t}^2+exp(2\Lambda \left(r\right))\mathrm{d}{r}^2+r^2(\mathrm{d}{\theta }^2+{sin}^2\theta \mathrm{d}{\varphi }^2).$$

(21)

Here, t is a Schwarzschild time (for a distant observer); r is a radial coordinate that has the meaning of circumferential radius; $\theta$ and $\varphi$ are ordinary spherical angles which are not affected by GR here; $\Phi \left(r\right)$ and $\Lambda \left(r\right)$ are the two metric functions to be determined from the Tolman–Oppenheimer–Volkoff equations (e.g., [2,4]).

Let M be the gravitational mass of the star and R be its circumferential radius. Outside the star ($r>R$), the metric (21) reduces to the Schwarzschild metric with $\Phi \left(r\right)=-\Lambda \left(r\right)=\frac{1}{2}ln(1-r_{\mathrm{g}}/r)$. Here, $r_{\mathrm{g}}=2GM/c^2$ is the Schwarzschild radius of the star, G is the gravitational constant, and c is the speed of light.

The space–time in a relatively thin and low-massive neutron star crust can be substantially curved, but the curvature is nearly constant there. Therefore, the expressions for the oscillation frequencies, obtained in previous sections, are expected to be rather accurate in a reference frame co-moving the crust.

As demonstrated in [47], while calculating the pure torsional vibration frequencies detected by a distant observer, it is quite accurate to use the metric functions constant through the crust,

$$\Phi =-\Lambda =\frac{1}{2}ln(1-x_{\mathrm{g}\ast }),$$

(22)

where $x_{\mathrm{g}\ast }=2G{M}_{\ast }/\left(c^2{r}_{\ast }\right)$, $r_{\ast }$ is any fiducial radial coordinate within the crust and $M_{\ast }$ is the gravitational mass within a sphere of radius $r_{\ast }$.

One can easily show that this approximation, being applied to the redshifted frequency ${\omega }_{\mu \mathsf{\ell }}$ of a pure torsional fundamental mode, reduces to

$${\omega }_{\mu \mathsf{\ell }}^2=(1-x_{\mathrm{g}\ast })\frac{{\int }_{\mathrm{crust}}\mathrm{d}V\mu }{{\int }_{\mathrm{crust}}\mathrm{d}V(\rho +\mathrm{P}/c^2)r^2},$$

(23)

where the volume element $\mathrm{d}V=r^2\mathrm{d}rsin\theta \mathrm{d}\theta \mathrm{d}\varphi$. The factor $(1-x_{\mathrm{g}\ast })$ can be naturally interpreted as due to the gravitational redshift. According to [47], the maximum deviation of the approximate frequencies (23) from the values of ${\omega }_{\mu \mathsf{\ell }}$ calculated in full GR does not exceed a few percent for any choice of $r_{\ast }$ (although choosing $r_{\ast }$ near the crust–core interface is slightly more accurate). We do not pretend to be very accurate here and accept this accuracy.

Moreover, we introduce similar relativistic corrections to the redshfited quantity ${\omega }_{B\mathsf{\ell }m}^2$. Then, instead of (20), we will use

$${\omega }_{B\mathsf{\ell }m}^2=(1-x_{\mathrm{g}\ast })\frac{\frac{1}{4\pi }{\int }_{\mathrm{crust}}\mathrm{d}V{I}_B}{\Xi (\mathsf{\ell },m){\int }_{\mathrm{crust}}\mathrm{d}V(\rho +\mathrm{P}/c^2)r^2}.$$

(24)

In numerical calculations, we take $r_{\ast }$ at the crust–core interface. In any case, (24) is an approximation, which seems to be sufficient for a semi-quantitative analysis, although the exact derivation of ${\omega }_{B\mathsf{\ell }m}^2$ in the full GR would be desirable.

7. The Dipole Magnetic Field

Let us apply the above formalism to a pure dipole magnetic field in the neutron star crust,

$$B_r=B_0{\left(\frac{R}{r}\right)}^{3}cos\theta ,B_{\theta }=\frac{B_0}{2}{\left(\frac{R}{r}\right)}^{3}sin\theta ,$$

(25)

where $B_0$ is the field at the magnetic pole on the surface in the local reference frame.

In this case, Equation (18) reduces to

$$\begin{array}{ccc}\hfill I_B& =& -B_{\theta }^2\left[P^{\prime 2}(1+3{cot}^2\theta )-3P^{\prime }{P}^{″}cot\theta +P^{\prime }{P}^{‴}\right.\hfill \\ & & \left.+\frac{m^2}{{sin}^2\theta }\left(-P^{\prime 2}-10P{P}^{\prime }cot\theta +8P^2{cot}^2\theta \right)\right].\hfill \end{array}$$

(26)

Let us substitute this expression into (24). The integral over r in the nominator is taken analytically, and we obtain

$${\omega }_{B\mathsf{\ell }m}^2=\frac{B_0^2{r}_2^3[{\left(r_2/r_1\right)}^3-1]}{12\pi {\int }_{r_1}^{r_2}\mathrm{d}r(\rho +\mathrm{P}/c^2)r^4}{\zeta }_{\mathsf{\ell }m},{\zeta }_{\mathsf{\ell }m}=\frac{1}{\Xi (\mathsf{\ell },m)}{\int }_0^{\pi }sin\theta \mathrm{d}\theta I_B.$$

(27)

Note that $r_1=R_{\mathrm{cc}}$ is the radius at the crust–core interface, and $r_2$ is taken slightly lower than R for the convenience of calculations; shifting $r_2\to R$ in the final expression would not affect the result; see, e.g., [46].

The Zeeman splitting is seen to be governed by the integrals over $\theta$ in ${\zeta }_{\mathsf{\ell }m}$. We have calculated them for ℓ = 2, 3, 4, and 5. In all these cases, ${\zeta }_{\mathsf{\ell }m}$ is virtually exactly described by the expression

$${\zeta }_{\mathsf{\ell }m}=c_0(\mathsf{\ell })+c_2(\mathsf{\ell })m^2,$$

(28)

where the coefficients $c_0(\mathsf{\ell })$ and $c_2(\mathsf{\ell })$ are listed in

Table 1. Coefficients $c_0(\mathsf{\ell })$ and $c_2(\mathsf{\ell })$, which describe ${\zeta }_{\mathsf{\ell }m}$ in Equation (28).

ℓ	2	3	4	5
$c_0$	0	2.6667	6.5455	11.487
$c_2$	0.66667	0.13333	–0.08182	–0.252

.

In this approximation, ${\omega }_{B\mathsf{\ell }m}^2$ is proportional to $B_0^2$. In the range of strict validity of our approach, ${\omega }_{B\mathsf{\ell }m}$ is also quadratic in $B_0$, as already stated in Section 5.

At $m=0$, the presented expressions should be valid for describing previous results, where Zeeman splitting was neglected. Then, the lowest frequency ${\omega }_{\mathsf{\ell }0}$ with $\mathsf{\ell }=2$ is not affected by the magnetic field at all ($c_0\left(2\right)=0$ in our approximation, which seems true only for the dipole field), while other frequencies with $\mathsf{\ell }>2$ increase as $B_0^2$. It would be instructive to check that $c_0\left(2\right)=0$ for the dipole magnetic field from previous calculations of ${\omega }_{20}\left(B\right)$ but they are controversial. For instance, compare Figure 4 in [29] with Figures 1 and 2 in [44].

The number of oscillation frequencies with $m=0$ is much lower than the total number of frequencies with different m. If we fix ℓ and vary m, the true oscillation frequency ${\omega }_{\mathsf{\ell }m}$ acquires the term behaving as $m^2{B}_0^2$. For $\mathsf{\ell }<4$, this term is positive, but for $\mathsf{\ell }\ge 4$, it becomes negative.

The presented expressions allow one to calculate (evaluate) ${\omega }_{B\mathsf{\ell }m}$ and, hence, total eigenfrequencies ${\omega }_{\mathsf{\ell }m}$ of fundamental torsional (magneto-elastic) oscillations in the crust of any neutron star, whose model is given and the magnetic field is not too strong. If several oscillation frequencies are measured from one and the same star, one can try to choose a value of $B_0$ and a stellar model to explain a set of observed oscillations at once. This gives a method to constrain $B_0$ and the neutron star model.

For illustration, let us outline the main properties of cyclic fundamental oscillation frequencies ${\nu }_{\mathsf{\ell }m}={\omega }_{\mathsf{\ell }m}/\left(2\pi \right)$ (expressed in Hz) at $\mathsf{\ell }=$2,…5 and $m=$0,…ℓ as functions of $B_0$. For this aim, we take neutron star models composed of matter with the BSk21 equation of state (EOS). This is a typical EOS; its basic properties, as well as the properties of corresponding stellar models, are nicely described by analytic expressions by Potekhin et al. [52]. The crust of these stars consists of spherical atomic nuclei and electrons; the inner crust also contains quasi-free neutrons, and there are quasi-free protons near the crust–core interface (${\rho }_{\mathrm{cc}}=1.34\times {10}^{14}$ g cm${}^{-3}$). The cores of such stars are nucleonic, also containing electrons and muons. The crust and core are described by the same energy–density functional of the nuclear interaction. The maximum mass of the star of this type is $M=2.27M_{\odot }$. For instance, the stellar model with $M=1.4M_{\odot }$ has radius $R=12.60$ km, and the radius of its crust–core interface is $R_{\mathrm{cc}}=11.55$ km. More information on neutron star models with the BSk21 EOS and on pure torsional ($B=0$) oscillations of these stars can be found in [47].

The behavior of the oscillation frequencies ${\nu }_{\mathsf{\ell }m}$ with increasing $B_0$ is shown in Figure 1. The four panels (a), (b), (c) and (d) correspond to the four neutron-star masses 1.0 $M_{\odot }$, 1.4 $M_{\odot }$, 1.8 $M_{\odot }$ and 2.2 $M_{\odot }$, respectively. The field $B_0\le 4\times {10}^{14}$ G is insufficiently strong to produce a noticeable Zeeman splitting with respect to m, and the oscillation frequencies remain very close to those in non-magnetic stars. At higher $B_0$, the Zeeman splitting becomes progressively more pronounced. Each non-magnetic frequency with fixed ℓ splits into $(\mathsf{\ell }+1)$ Zeeman components. The components with the lowest $m=0$ are plotted by the dashed lines. These components have been studied in the literature; the real spectrum is seen to be much more complicated. At $\mathsf{\ell }=2$ and 3, the frequencies ${\nu }_{\mathsf{\ell }0}$ are the lowest among respective Zeeman components (${\nu }_{\mathsf{\ell }m}$ increases with m), while at higher ℓ, they become the largest (${\nu }_{\mathsf{\ell }m}$ decreases with m). This is reflected in the inversion of the sign of coefficient $c_2(\mathsf{\ell })$ in Table 1.

According to Figure 1, the dependence of ${\nu }_{\mathsf{\ell }m}$ on $B_0$ for stars of different masses is similar. The frequencies ${\nu }_{\mathsf{\ell }m}$ for more massive stars are smaller. At sufficiently high $B_0$, our iterative solutions should become inaccurate (because we extend the plots to very high B, see Section 4 and Section 5). For $\mathsf{\ell }\le 4$, this seems to happen at $B_0\gtrsim 4\times {10}^{15}$ G. With increasing ℓ, the Zeeman splitting gets wider so that splittings for different ℓ start to overlap, manifesting possible cases of (avoided) crossings and mixtures of states with different ℓ. At still higher ℓ, these effects are expected to occur at still lower $B_0$, resulting in a rich, densely spaced and complicated spectrum of frequencies, a good project for further studies.

8. SGR 1806–20 and SGR 1900+14

To illustrate the above results, let us sketch their use for a possible interpretation of QPO oscillation frequencies observed in X-ray afterglows of two flaring SGRs. Taken the simplicity of our model, we will not try to be accurate and seek the complete sets of solutions, but present a few possible cases chosen by eye.

We analyze the low-frequency QPOs detected in the hyperflare of SGR 1806–20 (in 2004) and the giant flare of SGR 1900+14 (in 1998); e.g., [16,17,18,19,20] and references therein.

In particular, the low-frequency QPOs in the spectrum of the hyperflare of SGR 1806–20 were detected at 18, 26, 30, 92 and 150 Hz (and with less confidence at 17, 21, 36, 59 and 116 Hz). Taking into account a restricted range of frequencies calculated here (Figure 1), in Figure 2 we analyze the possibility of simultaneously detecting seven QPOs at 18, 26, 30, 17, 21, 36 and 59 Hz from one and the same star with a dipole magnetic field in the crust (leaving more complete analysis for future studies). Naturally, the Zeeman splitting of oscillation frequencies greatly increases the probability of such an interpretation.

Figure 2a shows the theoretical frequencies for the star of a typical mass 1.4 $M_{\odot }$. By varying $B_0$, we can choose such a $B_0$-interval that is consistent with more detections. This interval, $B_0\approx (3.2-3.4)\times {10}^{15}$ G, is shown by a vertical bar. It allows one to be qualitatively consistent with the three detected frequencies: 26, 30 and 59 Hz. The other four frequencies are not explained in this way. In particular, the frequencies ≲ 20 Hz remain unexplained.

A more attractive solution can be obtained by taking a more massive star, where torsional oscillation frequencies are smaller (e.g., [29,47]). This is shown in Figure 2b, where we take $M=2.2M_{\odot }$. With $B_0\approx (3.5-3.7)\times {10}^{15}$ G, we can now qualitatively explain the four detections at 18, 21, 26 and 59 Hz, out of six. The detected QPO at 17 Hz is quite close to QPOs at 18 and 21 Hz; it might be produced by nonlinear interactions of QPOs at 18 and 21 Hz, especially if the magnetic field configuration is more complicated than the pure dipole. However, the QPO at 30 Hz remains not explained. In any case, the explanation in Figure 2b seems more satisfactory than in Figure 2a, indicating that SGR 1806–20 may be a massive star. The obtained field strength is in line with current estimates of magnetar fields (e.g., [24]).

Figure 3 compares our theory with the detection of low-frequency QPOs in the giant flare of SGR 1900+14. The detected frequencies were 28, 53, 84 and 155 Hz; our current results can be used for explaining QPOs at 28 and 53 Hz. This problem of interpreting two detected points is much simpler than for SGR 1806–20. For instance, we take the $1.4M_{\odot }$ star and show one possible explanation, where the field strength $B_0\approx (2.2-2.4)\times {10}^{15}$ G is consistent with the current estimates of magnetar fields.

The results of this section can be treated as very preliminary. They illustrate the importance of Zeeman splitting of torsional oscillations of magnetars suggested by Shaisultanov and Eichler [45].

Let us stress once more that our approach is asymptotically accurate at not too strong $\mathit{B}$-fields. However, our attempts (in this section) to interpret magnetar QPOs indicate that real magnetar fields are just those, at which our approach can be already broken. Finding accurate values of breaking fields requires much effort.

The approach is simple because it is valid in the regime in which the elasticity of crystal matter is meant to be more important than that due to the magnetic field. Then, the effects of bulk viscosity $\mu$ dominate, and the magneto-elastic oscillations have much in common with pure torsional oscillations at $B=0$. In particular, the oscillations are mainly confined in the crystalline matter; in this way, they depend on the microphysics of the matter and the magnetic field configuration in the crust. This allows one to present solutions in a nearly closed form.

The difference between these solutions from those that were extensively studied previously is that they deal with non-axially symmetric perturbations of matter velocities and magnetic field perturbations (although the non-perturbed $\mathit{B}$-field is axially symmetric). It brings into play new degrees of freedom and enriches the oscillation spectrum. It would be an interesting project to extend previous numerical simulations (e.g., [37,38]) to this case.

With increasing B, the contribution of Alfvénic perturbations to the oscillations will increase (and finally dominate); these perturbations will spread out of the crystalline crust, and the presented approach will break. We fully appreciate many previous investigations (e.g., [30]) of magnetar oscillations based on Alfvén wave propagation over the entire star (not confined in the crust). Such waves can produce sophisticated oscillations of the star, which can be important for magnetar QPOs. Currently, this approach is mainly restricted by axially symmetric deformations. It would be interesting to extend such studies to non-axially symmetric ones.

9. Discussion and Conclusions

We have studied fundamental torsional (magneto-elastic) oscillations of neutron star crust in magnetars. The solution at $\mathit{B}=0$ is well known and corresponds to pure torsional oscillations of the stellar crust (Section 3). The solution for non-zero $\mathit{B}$ has been extensively studied since the beginning of the 2000s (e.g., [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44]) mostly numerically and for axially symmetric $\mathit{B}$-field configurations assuming axially symmetric velocity field of matter elements and magnetic field perturbations.

The exclusion was made by Shaisultanov and Eichler [45], who considered non-axially symmetric perturbations and predicted the effect of Zeeman splitting of torsional oscillation frequencies in magnetars. Unfortunately, their publication has not been given considerable attention.

We have tried to extend their consideration based on the same first-order perturbation theory with respect to the magnetic terms in the linearized oscillation equations. The equations are outlined in Section 2, neglecting relativistic effects. The properties of pure torsional oscillations are summarized in Section 3. The standard first-order perturbation approach is described in Section 4 for any $\mathit{B}$-field geometry. The approach reproduces Zeeman splitting [45] of pure torsional oscillation frequencies ${\omega }_{\mathsf{\ell }n}$ into a bunch of components, enriching the theoretical spectrum of magneto-elastic oscillations. We stress the well-defined nature of the solutions and their breakdown at very strong magnetic fields.

Section 5 applies the formalism of Section 4 to the case of fundamental oscillations (without nodes of radial wave functions) in a poloidal axially symmetric $\mathit{B}$-field. This case is especially simple. Finding the oscillation frequencies reduces it to taking 2D integrals (17) with well-defined integrands (18). In Section 6, we modified the equations to include relativistic effects. For illustration, Section 7 considers the simplest dipole field configuration in a neutron star crust. In Section 8, our results are used for sketching possible interpretations of detected low-frequency QPOs ($\le 70$ Hz) in afterglows of the hyperflare of SGR 1806–20 and of the giant flare of SGR 1900+14.

Let us emphasize that the Zeeman effect complicates the theory but greatly enriches the oscillation spectrum and simplifies the theoretical interpretation of detected QPOs. Needless to say, our restricted first-order iterative approach is far from being perfect and can be elaborated. First of all, one can consider other $\mathit{B}$-field geometries (e.g., [53])—not pure dipole but dipole-like fields, magnetic quadrupole, toroidal configurations or mixtures of these (see, e.g., [37,38]). Preliminary estimates show that for different $\mathit{B}$-field geometries, the Zeeman splitting can be qualitatively similar to that for the magnetic dipole but different in details. With the loss of magnetic field symmetry, the splitting will be more complete, removing residual degeneracy of oscillation frequencies. Additionally, one can replace the approximate description of magneto-elastic oscillations in GR (Section 6) with the exact description, in the spirit of Ref. [29]. Moreover, one can extend the iterative approach to ordinary torsional oscillations, with a finite amount of nodes of radial wave functions. Another direction of study could be to include delicate details of the microphysics of crustal matter, like advanced shear modulus, the effects of superfluidity, possible nuclear pasta phases and so on (as detailed, e.g., in [47]).