Oscillations of highly magnetized non-rotating neutron stars

Neutron stars are compact remnants with magnetic fields typically 10^11–10^13 G; magnetars can reach 10^14–10^16 G, potentially even higher in young objects. Such extreme fields can deform neutron stars and may power fast radio bursts, gamma-ray bursts, and superluminous supernovae. Stellar pulsations, excited by events such as core-collapse supernovae or giant flares, can emit gravitational waves and probe the neutron-star equation of state. While oscillation modes of non-magnetized stars are well studied, accounting for strong magnetic fields in full general relativity is challenging due to the nonlinear, coupled Einstein–Maxwell equations. Prior studies often used Newtonian frameworks or the Cowling approximation, which can overestimate mode frequencies by up to a factor of two. The research question is: how do strong, purely toroidal magnetic fields affect the oscillation spectra of non-rotating neutron stars in full dynamical spacetime? Leveraging recent advances in equilibrium modeling (XNS) and GRMHD evolution (Gmunu), the study aims to quantify magnetization effects on mode frequencies and link them to structural changes, especially compactness.

Magnetic field geometry crucially affects neutron-star structure: purely toroidal fields tend to make stars prolate while purely poloidal fields tend to make them oblate. However, such simple configurations are generally unstable, and simulations suggest rearrangement into mixed (twisted torus) configurations. Oscillation modes of non-magnetized stars have been extensively analyzed via perturbation theory and dynamical simulations, including with and without spacetime evolution. Magnetic effects have been studied in Newtonian and relativistic Cowling frameworks; these works often find frequency shifts scaling approximately with B^2 and generally small for moderate fields. However, the Cowling approximation can lead to significant errors (up to a factor of two) in mode frequencies. Recent tools such as XNS enable self-consistent general-relativistic equilibrium models with strong fields, and GRMHD codes like Gmunu allow evolution with dynamical spacetime, enabling a more accurate assessment of magnetization effects at high field strengths.

Equilibrium models: The XNS code is used to construct 12 axisymmetric, non-rotating neutron-star equilibria under the conformally flat condition, with a polytropic equation of state p = K ρ^γ (K = 1.6 × 10^5 cm^3 g^−1 s^−2, γ = 2) and a purely toroidal magnetic field specified by a polytropic-like prescription with toroidal magnetization index m = 1. All models share the same rest mass M0 = 1.68 M⊙. One non-magnetized reference model (REF) and 11 magnetized models (T1K1–T1K11) are arranged by increasing magnetic-to-binding energy ratio H/W. Stellar properties (central density, gravitational mass, circumferential and equatorial radii, polar-to-equatorial radius ratio, H/W, and maximum internal field Bmax up to ~6 × 10^17 G) are tabulated. Initial perturbations: To excite oscillations, three velocity perturbations are applied separately to each equilibrium: (i) ℓ = 0 perturbation in the radial component v_r: δv_r = a sin(π r/r_s(θ)) with amplitude a = 0.001 c; (ii) ℓ = 2 perturbation in v_θ: δv_θ = a sin(π r/r_s) sinθ cosθ with a = 0.01; (iii) ℓ = 4 perturbation in v_θ: δv_θ = a sin(π r/r_s) sinθ cosθ (3 − 7 cos^2θ) with a = 0.01. The radial sine factor has nodes at the center and surface to avoid boundary artifacts, while the angular dependence follows spherical harmonics. GRMHD simulations: Using Gmunu with full dynamical spacetime (no Cowling approximation for main results), each of the 12 models is evolved for 10 ms in 2D axisymmetry, yielding 36 simulations (three perturbations per model). Ideal GRMHD is assumed. Atmosphere floor density is set at 10^−10 of the initial maximum; the threshold ratio for atmosphere is 0.99. The computational domain is 0 ≤ r ≤ 60, 0 ≤ θ ≤ π with base resolution Nr × Nθ = 64 × 16, 4 AMR levels (effective 512 × 128). Refinement is guided by the relativistic potential Φ ≡ 1 − α with max Φ = 0.2; grids are coarsened near the center to satisfy CFL constraints. Frequency and eigenfunction extraction: At 361 spatial points inside the star, time series of the initially perturbed velocity component are FFT-transformed. Prominent spectral peaks, verified to be consistent across points, are fitted by parabolic interpolation to estimate eigenfrequencies f_eig and full-width-at-half-maximum (FWHM) as uncertainties. If needed, FFTs are integrated along radial lines to sharpen peaks. Eigenfunctions are reconstructed as spatial maps of FFT amplitude at the measured f_eig, serving to identify mode character across models and perturbations. Additional simulations under the Cowling approximation are performed for comparison (reported in Supplementary material), corroborating known biases of the approximation.

• Six dominant modes are identified: fundamental quasi-radial (ℓ = 0) mode F and first overtone H1; fundamental quadrupole (ℓ = 2) 2f and first overtone 2p1; fundamental hexadecapole (ℓ = 4) 4f and first overtone 4p1. Eigenfunctions match the expected angular structure and node counts.
• Threshold behavior with magnetization: For H/W ≤ 10^−2, eigenfrequencies are nearly unchanged across models despite Bmax up to ~10^17 G. For H/W ≥ 10^−1, all mode frequencies decrease markedly with increasing H/W, and higher-order ℓ = 4 modes become suppressed or masked by ℓ = 2 modes in the most magnetized cases.
• Representative frequencies (kHz, rounded) illustrate the trend (Table 1): REF has F ≈ 1.32, H1 ≈ 3.95, 2f ≈ 1.63, 2p1 ≈ 3.75, 4f ≈ 2.50, 4p1 ≈ 5.00. As magnetization increases, e.g., T1K10 shows F ≈ 0.49, H1 ≈ 1.33, 2f ≈ 0.78, 2p1 ≈ 1.30, 4p1 ≈ 1.84 (4f undetected); T1K11 shows H1 ≈ 0.77, 2f ≈ 0.60, 2p1 ≈ 0.99 with F and ℓ = 4 modes undetermined due to masking or data quality.
• Near-degeneracy between H1 and 2p1 emerges at high H/W as both are pushed to low frequencies.
• Compactness link: Stellar compactness M/R_circ remains almost constant for H/W ≤ 10^−2, but decreases sharply for H/W > 10^−1. Eigenfrequencies for all identified modes decrease quasi-linearly with M/R_circ, establishing a quasilinear correlation between f_eig and compactness in magnetized stars.
• Physical interpretation: Strong toroidal fields induce significant prolate deformation, reducing compactness; this structural change dominates over minor B^2-type shifts and suppresses oscillation frequencies when H/W ≳ 0.1.

The study demonstrates that magnetization effects on neutron-star pulsations are negligible below H/W ≈ 10^−2 but become substantial for H/W ≳ 10^−1, where mode frequencies drop and higher-ℓ modes are suppressed. This behavior directly correlates with a reduction in compactness due to magnetic deformation under strong toroidal fields, extending known compactness–frequency relations from non-magnetized to highly magnetized stars. Compared to Newtonian or Cowling-based studies that reported minor, approximately B^2-proportional frequency shifts, fully general-relativistic, self-consistent equilibria and dynamical spacetime evolutions reveal that large magnetic energies lead to non-perturbative structural changes that dominate the modal spectrum. The Cowling approximation is confirmed to overestimate frequencies by factors up to ~2. Astrophysically, strong internal toroidal fields at 10^15–10^17 G may arise in proto-neutron stars or merger remnants; the resulting oscillations could produce gravitational waves detectable with next-generation observatories (KAGRA, ET, NEMO). The findings underscore the need to include magnetic geometry and compactness changes when interpreting potential gravitational-wave asteroseismology signals from magnetized neutron stars.

This work provides a systematic, dynamical-spacetime GRMHD investigation of oscillations in non-rotating neutron stars with strong purely toroidal magnetic fields. It identifies a magnetization threshold in H/W beyond which oscillation frequencies are strongly suppressed and establishes a quasilinear correlation between eigenfrequencies and compactness M/R_circ in magnetized stars. These results refine expectations for gravitational-wave asteroseismology of magnetized neutron stars and highlight the structural role of extreme magnetic fields. Future research should address: (i) more realistic magnetic geometries (purely poloidal and twisted torus) with resistive GRMHD to treat exterior magnetospheres; (ii) rotation, to reflect observed neutron-star spins; and (iii) realistic equations of state to enhance constraints on dense matter from oscillation spectra.

• Magnetic geometry is restricted to purely toroidal fields, which are generally unstable in 3D; stability here is aided by 2D axisymmetry, potentially suppressing relevant instabilities.
• Simulations are axisymmetric and non-rotating; rotational effects on modes are not included.
• Ideal GRMHD without physical damping is used; resistive effects and crustal physics are not modeled, and exterior magnetosphere is not evolved.
• A simple polytropic equation of state (γ = 2) is employed; realistic microphysics is not included.
• Some eigenfrequencies are undetermined in highly magnetized cases due to mode masking or data quality, particularly for ℓ = 4 modes and F in select models.
• Evolution time is limited (10 ms), which may limit spectral resolution for very low-frequency features.