Zeeman splitting of torsional oscillation frequencies of magnetars

Magnetars form a special class of neutron stars possessing superstrong magnetic fields and demonstrating power flares triggered likely 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 splitting of non-magnetic oscillation frequencies into Zeeman components. Zeeman splitting of 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 theoretical interpretation of observations. We calculate several low-frequency oscillations of magnetars with pure dipole magnetic field in the crust. The results qualitatively agree with low-frequency QPOs detected in the hyperflare of SGR 1806--20, and in the giant flare of SGR 1900+14.


Introduction
It is well known that neutron stars are 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 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.01 M ⊙ .The envelope consists mostly of electrons and atomic nuclei.In addition, at densities ρ larger than the neutron drip density ρ drip ≈ 4.3 × 10 11 g cm −3 , there appear quasi-free neutrons.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 (ρ < ρ drip ) and the inner (ρ > ρ drip ) crust.The density of matter at the crust-core interface is ρ cc ∼ 1.4 × 10 14 g cm −3 .
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 ∼ 10 15 G (e.g.[22][23][24]) and demonstrating flaring activity that was most probably triggered by these fields.The activity is accompanied by the processes of enormous energy release in the form of hyperflares, supeflares and ordinary flares.The QPOs have been discovered at 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 just 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.

arXiv:2312.10022v2 [astro-ph.HE] 15 Mar 2024
The discovery of magnetar QPOs gave hope to develop reliable methods of exploring magnetar structure and evolution by comparing observations with elaborated seismological models.
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 magnetoelastic oscillations confined mainly to the crust, including Zeeman splitting.The paper is organized as follows.
In Section 2 we suggest to calculate 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 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 to the case of fundamental low-frequency torsional oscillations in a poloidal and axially symmetric B-field.In Section 6 we modify the equations by including relativistic effects.In Section 7 we consider the pure dipole 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.

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 elasticity of magnetic field lines, modify elastic properties of matter and lead to combined magnetoelastic 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 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 linearization procedure.All the quantities involved are divided into unperturbed ones and small perturbations.The non-perturbed star is taken stationary and spherically symmetric; ordinary spherical coordinates (r, θ, and ϕ) are most convenient for specifying position vectors r.Unperturbed quantities, like mass density ρ, depend only on r.The magnetic field B is assumed to be not too high and treated as perturbation.It is split further as B = B + B 1 , where B is the static (given) magnetic field of the star, while B 1 is a smaller field induced by the oscillations.The field 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 v(r, t) be a velocity of a matter element, and u(r, t) be a shift of its position with respect to equilibrium position r in a non-perturbed star; t is time.For an oscillation mode with frequency ω in the linear approximation one has v(r, t) = e iωt v(r), u(r, t) = e iωt u(r), v = iω u, B 1 (r, t) = e iωt B 1 (r), where italic vectors v, u and B 1 are (generally complex and stationary) amplitudes of perturbations; dot denotes time derivative.It is well known (e.g., [28]) that the torsional (magneto-elastic) oscillations can be treated as incompressible (∇ • v = 0; ∇ • u = 0).Then matter elements oscillate along spheres with fixed radius r and density ρ(r), with v r = 0 and u r = 0.The density perturbations are negligible.The equation of motion for a matter element reduces than to the stationary equation where the force terms T µ and T B (minus elastic forces per unit volume) come from the elasticity of crystalline lattice and magnetic field lines, respectively.Elastic deformations of the lattice are determined by the shear modulus µ(r) taken in the approximation of isotropic solid.Any i-th Cartesian component of T µ is given by where σ ik = µ (∂u i /∂x k + ∂u k /∂x i ) is the shear stress tensor.The magnetic force term T B is obtained in the standard manner.From the Maxwell induction equation (without dissipation) in our case one has B 1 = curl(u×B) and Eq. ( 2), supplemented with (3) and (4), allows one to find oscillation eigenfrequencies ω and eigenvectors u(r) which we discuss later in more detail.It is important that the T µ and T B terms are linear in the same shift vector u.This couples the lattice and B-field elasticities.The vector u (with u r = 0) plays role of the wave function of the problem.
Let us mention one feature of Eq. (2).Namely, let us multiply (2) by the complex conjugated shift vector u * and integrate over the star.In this way we obtain the equality where Therefore, the true oscillation frequency ω is formally expressed via two auxiliary frequencies, ω µ and ω B .In the absence of the magnetic field, one has ω = ω µ , that refers to a pure torsional oscillation due to elasticity of the lattice.In the absence of the lattice, one gets ω = ω B , that describes oscillations due to Alf ven 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 ω 2 µ contain positive definite self-conjugated operators and bilinear combinations of u and u * which guarantees that ω 2 µ > 0. It may be not true for the integral in the nominator of ω 2 B ; in principle, ω 2 B may be a complex number.
The decomposition of frequencies ω, similar to (5), was studied earlier by Schumaker and Thorne [4] (also see, e.g., [46,47]) for pure torsional oscillations at B = 0 with the aim to separate contributions of elastic shear forces and centrifugal forces in ω 2 µ .It is not a surprise that similar decomposition appears in the case of magneto-elastic oscillations.
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 ℓ =2, 3, . . .; ii) azimuthal number m (which is an integer varying from m = −ℓ to m = ℓ); and iii) the number n =0,1, 2, . . . of radial nodes of wave functions u(r, θ, ϕ).
The functions u nℓm and eigenfrequencies ω µnℓ can be found from Eq. ( 2) with T B =0.In view of spherical symmetry of our non-perturbed star, the quantisation 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 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 u can be presented as where P m ℓ (cos θ) is an associated Legendre polynomial, and Y nℓ (r) is a convenient dimensionless radial wave function which satisfies the equation Here the primes mean derivatives with respect to r.Since pure torsional oscillations are localized in the crystalline matter (r 1 ≤ r ≤ r 2 ), Eq. ( 9) is solved with the boundary conditions Y ′ (r 1 ) = Y ′ (r 2 ) = 0 (of zero traction).The solution gives the eigenfrequencies ω µnℓ .It is Eq. ( 9) which contains information on microphysics of neutron star matter.The angular dependence of the wave vectors u is standard (universal for all stars).

The simplest iterative solution
Now we suggest the simplest iterative procedure to include the effects of stellar magnetic field B. We treat the shift vector u nℓm , given by ( 7) and ( 8), as the zero-order solution of (2), and we treat the term containing T B , as perturbation.Let us restrict ourselves by the first-order perturbation.The problem is similar to finding 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ℓ + 1) zero-order wave functions corresponding to a fixed non-perturbed oscillation frequency ω nℓm .Such functions are often called the initial zero-order wave functions.The set of these functions is not unique: any linear superposition with constant transformation coefficients α m also describes an oscillatory state with the same ω nℓm .The iteration procedure prescribes to calculate the perturbation matrix T mm ′ [of dimension (2ℓ + 1) × (2ℓ + 1) ] on the basis of the initial wavefunctions.In our case, such matrix elements are given by 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 Eq. ( 10) and introduce the basis of 'true' zeroorder wave functions u ν (labelled by an index ν) which diagonalize T m ′ m .Let us denote the diagonal values by T ν .Their number is (2ℓ + 1); some of them can be equal if the perturbation removes degeneracy not completely.Then the oscillation frequency for the 'true' state ν is This means the splitting of ω µnℓ by the magnetic field that is traditionally treated as the Zeeman effect.This equation is in accord with the more general Eq. ( 6) in which ω 2 B 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 B(r) configuration, being mostly reduced to calculating and diagonalizing the matrix (11).The procedure gives 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 not too high magnetic fields: the magnetic splitting of oscillation frequencies has to be smaller than the zero-order torsion oscillation frequencies ω µnℓ themselves.The suggested procedure needs only the magnetic field B located in the crust; only zero-order wavefunctions are involved.Undoubtedly, sufficiently high 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 B the theory should be modified but the imprints of 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.

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 Here, z is the magnetic axis.Since ∇ • B = 0, the field components B r and B θ are related by The main problem is to calculate the matrix elements T mm ′ from Eq. ( 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 mm ′ is diagonal on the basis of the initial zero-order wavefunctions because the integration over ϕ reduces to where δ mm ′ 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 labelled by ν = m and readily given by Eq. ( 12) with T ν = T mm .This solution can be rewritten in the form In addition, due to the symmetry with respect to equatorial plane (θ = π/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 labelled by m = 0, 1,. . .ℓ (with ℓ + 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 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 that we drop this index, for simplicity.Therefore, our zero-order oscillation frequencies ω µ = ω µℓ depend only on ℓ.At low ℓ these frequencies are known to be small enough to explain lowest frequencies of QPOs observed in afterglow of magnetar flares (see Sect. 8 for more details).For fundamental torsional oscillations, with a very good accuracy, the radial wave function Y in Eqs. ( 8) and ( 7) is independent of r: Y ℓ (r) ≈ Y 0 , see [46,47].Since the numerators and denominators in Eq. ( 6) are proportional to Y 2 0 , this quantity just drops out of consideration.Then ω 2 Bℓm in Eq. ( 16) can be rewritten as where P = P m ℓ (θ), and Here, the primes denote derivatives with respect to θ.The subscript r after comma denotes derivative with respect to r. Eq. ( 18) is obtained from (11) by substituting ( 8) and ( 7), taking Y(r) = Y 0 , and integrating over ϕ 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 m ℓ (cos θ): Then we can rewrite (17) as The integral in the denominator is exactly the same as for pure torsional oscillations at B = 0. Thus, the problem of calculating Zeeman splitting of zero-order oscillation frequency reduces to the integration of known function I B over r and θ 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 θ .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, ω 2 Bℓm is quadratic in B. Then the true oscillation frequency has the form ω(B) = ω 2 (0) + αB 2 , where B is a characteristic B-field value, and α is some constant.Since our consideration is strictly valid at ω 2 (0) ≫ αB 2 , one gets ω(B) ≈ ω(0) + 1 2 αB 2 /ω(0), meaning a small correction to ω(0) 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 ω 2 (0) ≪ αB 2 one would get then ω(B) ∝ B although the validity of this expression is unclear.
Similar square-root expression for ω(B) 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 αB 2 was declared to be inaccurate at large B, e.g., [38]).The factor α has often been used as a convenient parameter to fit the calculated ω(B).Our approach gives exact prescription how to determine α including (m ̸ = 0) or disregarding (m = 0) Zeeman splitting at not too high B.
Let us mention that in the pioneering work on 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 constant B directed along the z-axis and constant µ(r)/ρ(r) in the magnetar crust, and demonstrated the main features of Zeeman splitting in magnetars.In case of more sophisticated magnetic field configurations, using the formalism of spherical vectors would be more complicated.

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 Special Theory of Relativity.Since the neutron star matter is essentially relativistic, one should replace (e.g.[4]) the mass density ρ in the equations of motion by the inertial mass density ρ + 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 ds 2 = − exp(2Φ(r)) dt 2 + exp(2Λ(r)) dr 2 + r 2 (dθ 2 + sin 2 θ dϕ 2 ).
Here, t is a Schwarzschild time (for a distant observer); r is a radial coordinate that has meaning of circumferential radius; θ and ϕ are ordinary spherical angles which are not affected by GR here; Φ(r) and Λ(r) 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 Φ(r) = −Λ(r) = 1 2 ln(1 − r g /r).Here, r 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 comoving the crust.
As demonstrated in [47], while calculating pure torsional vibration frequencies detected by a distant observer, it is quite accurate to use the metric functions constant though the crust, where x g * = 2GM * /(c 2 r * ), r * is any fiducial radial coordinate within the crust, and M * is the gravitational mass within a sphere of radius r * .One can easily show that this approximation, being applied to the redshifted frequency ω µℓ of a pure torsional fundamental mode, reduces to where the volume element dV = r 2 dr sin θ dθ dϕ.The factor (1 − x g * ) 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 ω µℓ calculated in full GR does not exceed a few per cent for any choice of r * (although choosing r * 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 ω 2 Bℓm .Then instead of (20), we will use In numerical calculations we take r * at the crust-core interface.In any case ( 24) is an approximation, which seems to be sufficient for a semi-quantitative analysis, although exact derivation of ω 2 Bℓm in full GR would be desirable.

The dipole magnetic field
Let us apply the above formalism to a pure dipole magnetic field in the neutron star crust, where B 0 is the field at the magnetic pole on the surface in the local reference frame.
In this case Eq. ( 18) reduces to Let us substitute this expression into (24).The integral over r in the nominator is taken analytically, and we obtain Note that r 1 = R cc is the radius at the crust-core interface, and r 2 is taken slightly lower than R for convenience of calculations; shifting r 2 → 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 θ in ζ ℓm .We have calculated them for ℓ=2, 3, 4, and 5.In all these cases, ζ ℓm is virtually exactly described by the expression where the coefficients c 0 (ℓ) and c 2 (ℓ) are listed in Table 1.
In this approximation ω 2 Bℓm is proportional to B 2 0 .In the range of strict validity of our approach ω Bℓ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 ω ℓ0 with ℓ = 2 is not affected by the magnetic field at all (c 0 (2) = 0 in our approximation, which seems true only for the dipole field), while other frequencies with ℓ > 2 increase as B 2 0 .would be instructive to check that c 0 (2) = 0 for the dipole magnetic field from previous calculations of ω 20 (B) but they are controversial: for instance, compare fig. 4 in [29] with figs. 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 ω ℓm acquires the term behaving as m 2 B 2 0 .For ℓ < 4 this term is positive, but for ℓ ≥ 4 it becomes negative.
The presented expressions allow one to calculate (evaluate) ω Bℓm and, hence, total eigenfrequencies ω ℓ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 oscillations 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 ν ℓm = ω ℓm /(2π) (expressed in Hz) at ℓ =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 contains also quasi-free neutrons, and there are quasi-free protons near the crust-core interface (ρ cc = 1.34 × 10 14 g cm −3 ).The cores of such stars are nucleonic, containing also electrons and muons.The crust and core are described by the same energydensity functional of nuclear interaction.The maximum mass of the star of this type is M = 2.27 M ⊙ .For instance, the stellar model with M = 1.4 M ⊙ has radius R = 12.60 km, and the radius of its crust-core interface is R 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 ν ℓ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 ⊙ , 1.4 M ⊙ , 1.8 M ⊙ , and 2.2 M ⊙ , respectively.The field B 0 ≤ 4 × 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 (ℓ + 1) Zeeman components.The components with 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 ℓ = 2 and 3, the frequencies ν ℓ0 are the lowest among respective Zeeman components (ν ℓm increases with m), while at higher ℓ they become the largest (ν ℓm decreases with m).This is reflected in inversion of sign of the coefficient c 2 (ℓ) in Table 1.
According to Figure 1, the dependence of ν ℓm on B 0 for stars of different masses is similar.The frequencies ν ℓ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 Sections 4 and 5).For ℓ ≤ 4 this seems to happen at B 0 ≳ 4 × 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.

SGR 1806-20 and SGR 1900+14
For illustrating the above results, let us sketch their use for 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.
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 (Fig. 1), in Figure 2 we analyze the possibility to simultaneously detect 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 interpretation.
Figure 2a shows the theoretical frequencies for the star of a typical mass 1.4 M ⊙ .By varying B 0 , we can choose such a B 0 -interval that is consistent with more detections.This interval, B 0 ≈ (3.2 − 3.4) × 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.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 Fig. 2b, where we take M = 2.2 M ⊙ .With B 0 ≈ (3.5 − 3.7) × 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 Fig. 2b seems more satisfactory than in Fig. 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.4 M ⊙ star and show one possible explanation, where the field strength B 0 ≈ (2.2 − 2.4) × 10 15 G is consistent with current estimates of magnetar fields.
The results of this section can be treated as very preliminary.They just 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 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 magnetic field.Then the effects of bulk viscosity µ 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 nearly closed form.
The difference of these solutions from those which were extensively studied previously is that they deal with not axially symmetric perturbations of matter velocities and magnetic field perturbations (although the non-perturbed B-field is axially symmetric).It brings into play with 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 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.
The exclusion was made by Shaisultanov and Eichler [45] who considered not 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 basing 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 B-field geometry.The approach reproduces Zeeman splitting [45] of pure torsional oscillation frequencies ω ℓn into a bunch of components, enriching thus 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 B-field.This case is especially simple.Finding the oscillation frequencies reduces 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 (≤ 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 in this way theoretical interpretation of detected QPOs.Needless to say that our restricted first-order iterative approach is far from being prefect and can be elaborated.First of all, one can consider other B-field geometries (e.g.[53]), not pure dipole, but dipole-like fields, magnetic quadrupole, toriodal configurations, or mixtures of these, etc. (see, e.g., [37,38]).Preliminary estimates show that for different 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 approximate description of magneto-elastic oscillations in GR (Section 6) by the exact description, in the spirit of Ref. [29].Moreover, one can extend the iterative approach to ordinary torsional oscillations, with finite amount of nodes of radial wave functions.Another direction of study could be to include delicate details of 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]).
Funding: This research was supported by the Russian Science Foundation (grant 19-12-00133 P)

5 Figure 1 .
Figure1.Fundamental torsional oscillation frequencies ν ℓm = ω ℓm /(2π) at four lowest ℓ (from 2 to 5) versus the dipole magnetic field B 0 at the magnetic pole.Calculations are performed for neutron star models with the BSk21 EOS at four masses; see panels (a)-(d).One can observe Zeeman splitting of any frequency in a non-magnetic star into (ℓ + 1) Zeeman components.The components with m = 0 are plotted by dashed lines.With increasing m (at fixed ℓ) the components either monotonically increase with respect to ν ℓ0 (for ℓ=2 and 3) or decrease (for ℓ=4 and 5).

B= 2 ℓ = 3 ℓ = 4 ℓ = 5 Figure 2 .
Figure 2. Calculated oscillation frequencies with ℓ ≤ 5 versus B 0 for two neutron star models, (a) M = 1.4 M ⊙ (same as in Fig.1b) and (b) M = 2.2 M ⊙ (same as in Fig.1d) compared with the low-frequency QPOs detected in the afterglow of the hyperflare of SGR 1806-20.The detected frequencies are plotted by horizontal dotted lines.More reliable detections (18, 26 and 30 Hz) are plotted by denser dotted lines.Less reliable ones (17, 21, 36 and 59 Hz) are displayed by more rarefied dotted lines.Vertical shaded strips show possible ranges of B 0 simultaneously consistent with some detections (see the text for details).

5 Figure 3 .
Figure 3. Theoretical oscillation frequencies with ℓ ≤ 5 versus B 0 for a neutron star model with M = 1.4 M ⊙ (same as in Fig.1b) compared with two low-frequency QPOs (28 and 53 Hz) detected in the afterglow of the giant flare of SGR 1900+14 (the horizontal dotted lines).Vertical shaded strip shows possible ranges of B 0 consistent with both detected frequencies.