Oscillations of Highly Magnetized Non-rotating Neutron Stars

Highly magnetized neutron stars are promising candidates to explain some of the most peculiar astronomical phenomena, for instance, fast radio bursts, gamma-ray bursts, and superluminous supernovae. Pulsations of these highly magnetized neutron stars are also speculated to produce detectable gravitational waves. In addition, pulsations are important probes of the structure and equation of state of the neutron stars. The major challenge in studying the pulsations of highly magnetized neutron stars is the demanding numerical cost of consistently solving the nonlinear Einstein and Maxwell equations under minimum assumptions. With the recent breakthroughs in numerical solvers, we investigate pulsation modes of non-rotating neutron stars which harbour strong purely toroidal magnetic fields of $10^{15-17}$ G through two-dimensional axisymmetric general-relativistic magnetohydrodynamics simulations. We show that stellar oscillations are insensitive to magnetization effects until the magnetic to binding energy ratio goes beyond 10%, where the pulsation mode frequencies are strongly suppressed. We further show that this is the direct consequence of the decrease in stellar compactness when the extreme magnetic fields introduce strong deformations of the neutron stars.


Introduction
Neutron stars (NSs) are compact objects formed by core-collapse supernovae. Due to field amplification in the violent formation processes, most NSs are endowed with strong magnetic fields of 10 11−13 G 11 . In some extreme cases, the magnetars can harbour even stronger magnetic fields of 10 14−16 G, about 1000 times stronger than usual pulsars (for comparison, the magnetic field of a sunspot is 10 3 G 12 ). Younger magnetars may carry even higher magnetic fields since they have been subjected to dissipative processes for shorter times 13 .
These extreme magnetic fields affect the structure and evolution of NSs. For instance, strong magnetic fields can deform NSs 13,14 . A direct consequence of structural deformations of NSs could be significant gravitational wave emissions [15][16][17] . The geometry of magnetic fields of the NSs is a crucial factor governing the physics of NSs. However, the field configuration inside the NS is unknown. Studies of equilibrium models with simple field configurations suggest that a purely toroidal field makes NSs prolate [18][19][20] while a purely poloidal field forces the stars to become oblate [21][22][23] . Nevertheless, these simple geometries are expected to be unstable [24][25][26][27][28] . Numerical simulations suggest that the magnetic fields of the NSs are rearranged rapidly due to these instabilities, leading to a mixed configuration of toroidal and poloidal fields, which is roughly axisymmetric [29][30][31][32] . This mixed geometry is usually called twisted torus.
The major difficulties in studying magnetized NSs come from the non-linear nature of Einstein equations, and with Maxwell equations fully coupled, analytical calculations are generally impossible. Hence, numerical computations are inevitable to solve all the involved physics with a minimum number of assumptions. Until recently, due to breakthroughs in generalrelativistic magnetohydrodynamics (GRMHD) simulations, dynamical studies of magnetized NSs have become possible, e.g.
A novel approach to compute strongly magnetized equilibrium models is recently presented and demonstrated by an open-source code XNS 9, 13, 66-75 . Moreover, the GRMHD code Gmunu 10, 76-78 allows us to robustly evolve NSs in dynamical spacetime even with extreme magnetic fields of 10 15−17 G. With these powerful tools in hand, we are now in a much better position to systematically investigate the oscillation modes of magnetized NSs.
In this work, we numerically study the oscillations of highly magnetized non-rotating axisymmetric (two-dimensional) NSs. Specifically, we first construct 12 equilibrium models with different magnetic to binding energy ratios H /W using XNS, including one non-magnetized reference model named 'REF' and 11 magnetized models in ascending order of H /W named T1K1, T1K2, ..., T1K11 (Methods). Next, we utilize Gmunu to perturb and evolve the equilibrium models in dynamical spacetime, where we try three different initial fluid perturbations for excitation of stellar oscillations, namely = 0, = 2, and = 4 perturbations (Methods). After that, we perform a Fourier analysis of the simulation results to examine how the eigenfrequencies of oscillation modes vary with H /W of the NS (Methods), and we discuss possible reasons behind our results.

Magnetization effects on oscillations of NSs
In total, six dominant oscillation modes are observed in our numerical study, namely the fundamental quasi-radial ( = 0) mode F and its first overtone H 1 , the fundamental quadrupole ( = 2) mode 2 f and its first overtone 2 p 1 , as well as the fundamental hexadecapole ( = 4) mode 4 f and its first overtone 4 p 1 (we follow the notations in Ref. 59 ). Each mode is predominantly excited under the initial perturbation with the corresponding index, and each eigenfunction qualitatively agrees with the spherical harmonic in the corresponding perturbation function, as shown in Fig. 1. The measured eigenfrequencies of the six modes in the 12 different NS models are summarized in Table 1, where the undetermined eigenfrequencies denoted by 'N/A' in different columns stem from different reasons below. For the column of F mode, the missing eigenfrequencies are due to unsatisfactory data quality in Gmunu simulations of T1K8 and T1K11 models under = 0 perturbation. On the other hand, for the columns of 4 f and 4 p 1 modes, some eigenfrequencies are missing because the hexadecapole ( = 4) modes are masked by the quadrupole ( = 2) modes and are no longer the dominant modes in Gmunu simulations of the most magnetized models under = 4 perturbation. To better illustrate the results in Table 1, we plot in Fig. 2 the eigenfrequencies f eig of the six modes as functions of the magnetic to binding energy ratio H /W of the NS model.
We have observed an H /W threshold for stellar magnetization to start affecting the oscillations of NSs. For NSs with H /W 10 −2 , stellar oscillations are insensitive to magnetization effects. This can be seen from Table 1 that f eig of every oscillation mode is nearly the same for the first six models (REF -T1K5) even though these models span a few orders of magnitude in H /W and can achieve a maximum field strength of 10 15−17 G; this can also be seen from Fig. 2 that the data points at H /W ∼ 0 show a nearly horizontal trend. On the other hand, for NSs with H /W 10 −1 , stellar oscillations are significantly suppressed by stronger magnetization. Refer to the data points at H /W > 10 −1 in Fig. 2, f eig decreases with H /W in general, and all the oscillation modes are pushed towards the low-frequency region, leading to the near-degeneracy of H 1 and 2 p 1 modes. Moreover, as afore-explained about the undetermined eigenfrequencies, = 4 perturbation excites the quadrupole ( = 2) modes preferentially over the expected hexadecapole ( = 4) modes in the most magnetized models, hinting at suppression or even disappearance of higher-order oscillation modes in a more magnetized NS for H /W 10 −1 . To summarize, magnetization effects start to hinder stellar oscillations if H /W of the NS passes the threshold somewhere between 10 −2 and 10 −1 .

Compactness as an underlying factor
The magnetization effects on NS oscillations discussed above may be understood by studying the compactness M/R circ of the NS, where M is the gravitational mass and R circ is the circumferential radius. As shown Ref. 79 , the eigenfrequencies of the fundamental quasi-radial and quadrupole modes are related to the stellar compactness for non-magnetized NSs, and we suspect this correlation also holds for highly magnetized NSs. Thus, based on our NS models, we plot in Fig. 3 the compactness M/R circ against the magnetic to binding energy ratio H /W . We find that M/R circ remains nearly unchanged for H /W 10 −2 but decreases dramatically for H /W > 10 −1 , which agrees with the trends of f eig (H /W ) shown in Fig. 2 and indeed reveals a correlation between eigenfrequencies of oscillation modes and stellar compactness. We also plot in Fig. 4 f eig against M/R circ . For all the modes, f eig decreases together with M/R circ in an almost linear way. Therefore, we found a quasilinear relation between f eig and M/R circ for magnetized NSs. The complete physical interpretation of our results is that a strong toroidal field

2/19
can cause deformation of the NS 13 and alter the stellar compactness, so the propagation of seismic activities inside the NS is affected. In consequence, the eigenfrequencies of oscillation modes are correspondingly modified.

Discussion
In this work, we systematically investigate how a strong purely toroidal magnetic field with a field strength of 10 15−17 G affects the oscillations of non-rotating NSs via two-dimensional axisymmetric simulations. We carefully extract the eigenfrequencies of the excited oscillation modes and construct the corresponding eigenfunctions from the simulated data. We have found that stellar oscillations are insensitive to magnetization effects for NSs with magnetic to binding energy ratio H /W 10 −2 , even though the maximum magnetic field strength B max can reach O(10 17 ) G in the star. However, stellar oscillations are suppressed significantly by stronger magnetization if H /W 10 −1 . This behaviour can be understood by the decrease of stellar compactness due to strong magnetic fields. We show that the compactness has the same dependence on H /W as the eigenfrequencies and demonstrate that the correlation between eigenfrequencies and compactness exists not only in non-magnetized NSs 79 but also in highly magnetized NSs.
We compare our results with previous Newtonian studies, e.g. Ref. 45,48 . These studies considered either perturbative or self-consistent MHD to construct the equilibrium models in the Newtonian regime. Both approaches found that the magnetic distortion and frequency shift in oscillation modes due to toroidal fields are minor corrections approximately proportional to B 2 (or roughly H /W in this work). However, in our GRMHD simulations, the equilibrium models are constructed by solving self-consistent general-relativistic magnetohydrostatic equations in the code XNS 9, 13, 66-75 . When H /W 10 −1 , the magnetic deformations are far from small corrections, and thus the stellar compactness is significantly reduced. Therefore, the effect of decreasing compactness dominates and results in the suppression of oscillation modes. Besides, we compare our results with those under the Cowling approximation and we corroborate what has been shown in the literature 58,59 , namely that the Cowling approximation can lead to errors of factors of 2 (see Supplementary information).
The strongest magnetic field strength of 10 17 G in this work is not expected to be observed in the exterior of ordinary pulsars and magnetars. Nevertheless, since the toroidal fields are enclosed inside the NSs, this ultra-high field could exist in the interior regions. Moreover, such field strength could also be generated during the formation of a proto-NS 13 , and binary neutron star mergers 80 . The excited oscillation modes in these scenarios are potential sources for gravitational waves, and they could be detected with the next-generation detectors, such as the Kamioka Gravitational Wave Detector (KAGRA) 81 , the Einstein Telescope (ET) 82 , and the Neutron Star Extreme Matter Observatory (NEMO) 83 .
This work presents the first step to understanding how magnetic fields with different geometries affect the oscillations of NSs. Since stellar models with purely toroidal fields are generally unstable 25 , the instability is only suppressed due to the restriction to 2D axisymmetry in this work. Therefore, a natural extension considers strong purely poloidal fields and the more realistic twisted torus configuration. Since these field configurations extend to the regions outside NSs, an accurate and robust resistive GRMHD solver could be used to model these regions. This solver has already been implemented into Gmunu 78 for future studies. In addition to different configurations of magnetic fields, rotation should also be taken into account to work towards a more realistic problem, as the observed NSs are suggested to be rotating. Furthermore, introducing realistic EoSs is essential since one of the most important purposes of oscillation studies is to probe the structure and the EoSs of NSs.

Equilibrium models
Equilibrium models of NSs are constructed by the code XNS 9, 13, 66-75 . XNS is a branch of the X-ECHO code 66 developed to compute equilibrium models of highly magnetized axisymmetric NSs with rotations. Different magnetic field configurations 13 , uniformly and differentially rotating profiles 66 , and polytropic and non-polytropic tabulated equations of state 74 are admitted. XNS enforces the 3 + 1 formulism, the conformal flatness condition, and the assumption of axisymmetric and stationary space-time so that the line element can be written as where α(r, θ ) is the lapse function, ψ(r, θ ) is the conformal factor, and β φ (r, θ ) is the shift vector (β φ = 0 for non-rotating configurations). We assume a polytropic EoS p = Kρ γ for the stellar fluid, where p is the pressure, K is the polytropic constant, ρ is the density, and γ is the adiabatic index; as well as a polytropic expression B φ = α −1 K m (ρhϖ 2 ) m for the toroidal field, where K m is the toroidal magnetization constant, h is the specific enthalpy, ϖ 2 = α 2 ψ 4 r 2 sin 2 θ , and m ≥ 1 is the toroidal magnetization index. Although the field configuration of an isolated NS is expected to be a mixture of toroidal and poloidal fields, it is important first to assess how a simpler field geometry would affect the oscillations of NSs before we move on to the more complicated Twisted Torus case.
In total, 12 equilibrium models are computed with XNS, where one of them is a non-magnetized reference model named 'REF', and the remaining 11 models are magnetized. All the 12 models share the same rest mass M 0 = 1.68 M , and the same K = 1.6 × 10 5 cm 5 g −1 s −2 and γ = 2 in the fluid EoS. The 11 magnetized models have the same m = 1 but different values of K m in the B φ expression, and they are arranged in ascending order of magnetic to binding energy ratio H /W , where the one with the lowest ratio is named 'T1K1', and the one with the second-lowest ratio is named 'T1K2', so on and so forth. ('T1' specifies the toroidal magnetization index being 1 and 'K' stands for K m ). The detailed properties of all the 12 models are summarized in Table 2.

Initial perturbations to excite oscillations
Consulting a similar study on rotating non-magnetized NSs done by Ref. 59 , we try the following three types of initial fluid perturbations for exciting oscillations in the equilibrium models.
First, we have the = 0 perturbation on the r-component of the three-velocity field, where r s (θ ) locates the surface of the NS, and the perturbation amplitude a (in unit of c) is chosen to be 0.001. Second, we have the = 2 perturbation on the θ -component of the three-velocity field, where a is chosen to be 0.01. Lastly, we have the = 4 perturbation on the θ -component of the three-velocity field, where a is again set to be 0.01. All the three perturbation functions comprise a sine function of r and the θ -part of a spherical harmonic with the corresponding index. The sine function of r has its nodes at the centre and on the surface of the NS to avoid initial perturbations on sensitive boundaries of the problem and minimize any potential numerical errors. On the other hand, spherical harmonics are a natural choice for exciting oscillations on a sphere-like object. Moreover, for the higher-order = 2 and = 4 perturbations, the perturbation amplitude a has to be larger to induce any observable oscillations.

Simulations
Simulations are performed with our code Gmunu 10, 76-78 . For each of the 12 equilibrium models, we execute Gmunu three times, once for each initial perturbation function. Hence, 12 × 3 = 36 simulations are carried out in total. In all the 36 simulations, the models are evolved over a time span of 10 ms with the polytropic EoS p = Kρ γ , under the same setting as in the computation of equilibrium models, namely, γ = 2 and K = 110. The lowest allowed rest mass density ('atmosphere') is set to be ρ atmo = ρ max (t = 0) × 10 −10 , and the ratio of ρ atmo to threshold density ρ thr is ρ atmo /ρ thr = 0.99. For completeness, we also perform simulations under the Cowling approximation (see Supplementary information), with other settings unchanged.
The two-dimensional computational domain covers 0 ≤ r ≤ 60, 0 ≤ θ ≤ π with the resolution N r × N θ = 64 × 16 where each block has 8 2 cells, thus allowing 4 AMR level (an effective resolution of 512 × 128). The grid refinement used in this study is identical to the GR simulations in Ref. 77 . In particular, we define a relativistic gravitational potential Φ := 1 − α. As Φ is almost proportional to M/R, we can use Φ −1 as a measure of the characteristic length-scale 77 . For any Φ larger than the maximum potential Φ max (which is set as 0.2 in this work), the block is set to be the finest. While for the second-finest level, the same check is performed with a new maximum potential which is half of the previous one, so on and so forth. To avoid the rigorous Courant-Friedrichs-Lewy (CFL) condition at the centre of the star, the grids are enforced to be coarsened for keeping r∆θ ∼ ∆r when r is smaller than 0.5. (Unless otherwise specified, all quantities in this subsection are in dimensionless units c = G = M = 1.)

Extraction of eigenfrequencies and eigenfunctions
We analyze the data from a Gmunu simulation in the following three steps. For the first step, we extract the time evolutions of the initially perturbed component of the three-velocity field at 361 (r, θ )-points in the NS model and compute the Fast Fourier

4/19
Transform (FFT) of the temporal data at each (r, θ )-point. Hence, 361 FFT spectra, plots of magnitude of the complex FFT in the frequency domain, are obtained altogether. According to Ref. 59 , our initial perturbation amplitudes are small enough such that the overall evolution of the input model in a Gmunu simulation can be described as a superposition of a few global oscillation modes. We verify this by observing that the FFT spectra obtained at different spatial points show discrete peaks and agree well on the peak positions.
For the second step, we extract the eigenfrequencies of the excited oscillation modes. Usually, the FFT spectrum at a spatial point where the initial perturbation function has a large magnitude can reveal FFT peaks loud enough for further analysis (e.g. at (r, θ ) (r e /2, π/2), (r e /2, π/4), and (r e /2, 2π/15) for = 0, = 2, = 4 perturbations respectively). Nevertheless, occasionally, we may have to integrate the FFT spectra along a radial line for sharper FFT peaks (along θ = π/2, π/4, and 2π/15 for = 0, = 2, = 4 perturbations respectively). Since our study here is in the ideal GRMHD regime with no physical damping of the oscillations, we apply parabolic interpolation instead of Lorentzian fitting to the peaks in the single-point or integrated FFT spectrum for simplicity (see Fig. 5 as an example). We then take the interpolated peak positions as the measured eigenfrequencies f eig and the full-width-at-half-maximums (FWHMs) of the parabolic interpolations as the uncertainties in eigenfrequency extraction.
For the third step, we extract the eigenfunctions of the excited oscillation modes. According to Refs. 59,84 , the eigenfunction of a mode is correlated to the spatial map of FFT amplitude at the eigenfrequency of the mode, where FFT amplitude is the magnitude of the FFT multiplied by the sign of its real part. Using our FFT data computed at the 361 points, we spatially map the FFT amplitude at the frequency to which the measured eigenfrequency is the closest in the discretized frequency domain of our FFT analysis for simplicity. The eigenfunction visualized by such a spatial map can serve as a unique trademark to help us identify the same oscillation mode excited in different Gmunu simulations so that we can investigate the dependence of eigenfrequency f eig of a particular mode on the magnetic to binding energy ratio H /W of the input model.
In the end, we can obtain the curves of f eig (H /W ) for different oscillation modes to examine the magnetization effects on oscillations of NSs. Lastly, we determine the correspondence between the modes found in our study and the modes in the literature by comparing the eigenfrequencies at zero magnetic energy, f eig (H /W = 0), of the modes we found here with the mode frequencies previously reported for a non-magnetized non-rotating NS model with a similar gravitational mass 59 .    Table 2. Stellar properties of the 12 equilibrium models constructed by the XNS code. All numerical values are rounded off to two decimal places. ρ c is the central density, M is the gravitational mass, R circ is the circumferential radius, r e is the equatorial radius, r p /r e is the ratio of polar radius r p to equatorial radius r e (a purely toroidal field elongates the NS along the z-axis), H /W is the ratio of total magnetic energy H to total binding energy W , and B max is the maximum field strength achievable inside the star. These quantities are defined according to Ref. 13 . All the 12 models share the same rest mass M 0 = 1.68 M , polytropic constant K = 1.6 × 10 5 cm 5 g −1 s −2 , and adiabatic index γ = 2. The 11 magnetized models also share the same toroidal magnetization index m = 1. It can be seen that the eigenfunctions of the higher-order quadrupole ( = 2) and hexadecapole ( = 4) modes have more nodes in the θ -direction compared to the quasi-radial ( = 0) modes, while the eigenfunction of each first overtone has more nodes in the r-direction compared to its fundamental mode. Furthermore, each eigenfunction qualitatively agrees with the spherical harmonic in the corresponding perturbation function.

Supplementary Note 1: Comparison to the Cowling approximation
It is known that the Cowling approximation (i.e. fixed spacetime) overestimates the oscillation mode frequencies of NSs up to a factor of ∼ 2. Simulations with dynamical spacetime are essential for more accurate results. In this supplementary information, we show that the Cowling approximation is indeed inadequate under the circumstances described in this work. The corresponding visualizations of eigenfunctions of the excited oscillation modes are shown in Supplementary Figure 6. Supplementary Table 3 summarizes the measured eigenfrequencies of the six modes in the 12 different NSs with the undetermined eigenfrequencies due to unsatisfactory data quality denoted by 'N/A'. We plot in Supplementary Figure 7 the eigenfrequencies f eig of the six modes with and without Cowling approximation as functions of the magnetic to binding energy ratio H /W of the NS model. We further summarize the relative differences in f eig with the Cowling approximation with respect to those without the approximation in Supplementary Table 4 to better illustrate the discrepancies. We found a similar qualitative relation between f eig of the six modes and H /W for both simulations with and without the Cowling approximation. However, a clear quantitative difference in the values of f eig is observed. The largest discrepancy is found in the F-mode, where the eigenfrequencies in all NS models are overestimated by a factor of ∼ 2 under the Cowling approximation. The discrepancies are also significant for H 1 -mode with relative differences of 20 per cent in most cases and even up to ∼ 65 per cent in the NS model with the strongest magnetic field strength. The discrepancies are less severe for higher-order modes (i.e. = 2 and = 4 modes). The relative differences are less than 20 per cent for 2 f and 2 p 1 modes, while they are less than 10 per cent for 4 f and 4 p 1 modes.
Hence, imposing the Cowling approximation in our simulations also results in overestimating the correct eigenfrequency up to a factor of ∼ 2. Dynamical spacetime is still necessary for magnetized NS simulations to obtain more realistic oscillation mode frequencies, especially for the lower-order modes.  Table 3. Measured eigenfrequencies of the six dominant oscillation modes with the Cowling approximation in the 12 NS models, including the fundamental quasi-radial ( = 0) mode F and its first overtone H 1 , the fundamental quadrupole ( = 2) mode 2 f and its first overtone 2 p 1 , as well as the fundamental hexadecapole ( = 4) mode 4 f and its first overtone 4 p 1 , all predominantly excited under the perturbation with the corresponding l index. All eigenfrequencies are in kHz and rounded off to two decimal places. The undetermined eigenfrequencies in specific models are denoted by 'N/A'. The missing eigenfrequencies are due to unsatisfactory data quality in Gmunu simulations under the perturbations.  Table 4. Relative difference in per cent in the frequencies of the six dominant oscillation modes with the Cowling approximation with respect to those without the approximation. The models containing undetermined eigenfrequencies are denoted by 'N/A'. All relative differences are rounded off to two decimal places.