Asteroseismic g-mode period spacings in strongly magnetic rotating stars

Strong magnetic fields are expected to significantly modify the pulsation frequencies of waves propagating in the cores of red giants or in the radiative envelopes of intermediate- and high-mass main-sequence stars. We calculate the g-mode frequencies of stars with magnetic dipole fields which are aligned with their rotational axes, treating both the Lorentz and Coriolis forces non-perturbatively. We provide a compact asymptotic formula for the g-mode period spacing, and universally find that strong magnetism decreases this period spacing substantially more than is predicted by perturbation theory. These results are validated with explicit numerical mode calculations for realistic stellar models. The approach we present is highly versatile: once the eigenvalues $\lambda$ of a certain differential operator are precomputed as a function of the magnetogravity and rotational frequencies (in units of the mode frequency), the non-perturbative impact of the Coriolis and Lorentz forces is understood under a broad domain of validity, and is readily incorporated into asteroseismic modeling.

In recent years, there has been a large amount of progress in developing asteroseismology as a probe of strong internal magnetic fields, particularly through their effects on the gravity (g) modes which propagate in radiative regions.Such fields likely have important consequences for the transport of angular momentum within evolved stars (Mathis & de Brye 2012;Fuller et al. 2019;Aerts et al. 2019).On the red giant branch, g modes propagate in the radiative core, which may possess strong magnetic fields left over from efficient core convective dynamos on the main sequence (Fuller et al. 2015;Stello et al. 2016b).In these cases, magnetism may have a significant effect on the frequency spectrum: by measuring these frequency patterns, Li et al. (2022b) strongly constrain both the rotational periods and field strengths (≳ 30 kG) as well as their geometries for a modest sample ★ E-mail: nrui@caltech.edu of red giants.Even stronger magnetic fields ≳ 100 kG are commonly invoked to explain the observed suppression of dipole (ℓ = 1) and quadrupole (ℓ = 2) oscillation modes in red giants (e.g., García et al. 2014;Stello et al. 2016a,b).Specifically, mode suppression is expected to occur in the non-perturbative "strong magnetogravity" regime (Fuller et al. 2015;Lecoanet et al. 2017;Rui & Fuller 2023), when a mode's frequency  is sufficiently close to the critical frequency In Equation 1,   =   / √︁ 4 is the radial component of the Alfvén velocity,  is the radial coordinate, and  is the Brunt-Väisälä (buoyancy) frequency, given by where  is the adiabatic index.Equivalently, mode suppression occurs at some frequency  crit when the magnetic field is at least comparable to some critical field (Fuller et al. 2015): Complementarily, main-sequence pulsators of intermediate mass (≳1.3 ⊙ ) have radiative, rather than convective, envelopes, and their g modes extend to their surfaces where they can be observed directly (Aerts 2021).Examining the slowly pulsating B-type (SPB) star HD 43317, Lecoanet et al. (2022) demonstrate that observations of mode suppression in main sequence (MS) pulsators may place meaningful constraints on their internal magnetism.The detection of g modes in the vicinity of magnetic suppression suggests MS pulsators as a separate platform for testing the effect of strong magnetic fields on propagating gravity waves.
In the absence of effects such as magnetic fields or rotation, successive g modes are evenly spaced in period by a constant g-mode period spacing   , which can be estimated as where the integral is over the part of the radial cavity within which g modes propagate ( < ).However, both rotation and magnetism leave distinctive signatures on the period spacing, both by lifting the degeneracy between modes of different  (by breaking the spherical symmetry of the system) and by introducing period dependence (Bouabid et al. 2013;Van Reeth et al. 2016;Dhouib et al. 2022).The period spacing as a function of period   =   () is therefore a valuable measure for rotational and magnetic effects (Van Beeck et al. 2020;Henneco et al. 2021).Characterizing this observable non-perturbatively is the primary focus of this work.
Our paper proceeds as follows.Section 2 presents the problem statement and motivates the asymptotic treatment of magnetism and rotation.Section 3 derives the differential operator which governs the horizontal structure of magnetic gravito-inertial modes.In Section 4, we numerically calculate this operator's eigenvalues, which enter directly into an asymptotic formula for the period spacing.In Section 5, we solve the radial oscillation problem directly, including both magnetism and rotation while partially relaxing the asymptotic assumption.Section 6 presents the results of such calculations for models of red giants,  Doradus ( Dor), and SPB pulsators.Finally, Section 7 concludes.The reader seeking our observational predictions is guided to Equations 33 and 40 (for the asymptotic period spacing) and the discussion in Section 6.
The central goal of this work is to calculate the period spacing pattern in the simultaneous presence of rotation and an axisymmetric radial magnetic field in a non-perturbative way.The work proceeds under the "traditional approximation of rotation and magnetism" defined in Section 2.1 (which restricts attention purely to the radial field).

The traditional approximation of rotation and magnetism (TARM)
Pure, low-frequency gravity waves follow the dispersion relation when  ≫ .Therefore, their wavenumbers are primarily radial, with their radial wavenumbers exceeding their horizontal wavenumbers by ratios In the presence of restoring forces other than buoyancy or pressure (e.g., Coriolis forces, magnetic tension), the dispersion relation will be modified from Equation 5.However, for modes which still have g mode character (i.e., buoyancy is still a significant restorative force),   / ℎ will still be comparable to / ≫ 1. Throughout, we restrict our attention to modes whose wavenumbers are primarily radial: this is a crucial assumption of our work.
This approximation underlies the standard analytic treatment of gravity waves in rotating stars.The qualitative behavior of lowfrequency gravito-inertial waves can be seen in the dispersion relation in the fully Jeffreys-Wentzel-Kramers-Brillouin (JWKB) limit: see, e.g., Bildsten et al. (1996) and Lee & Saio (1997).Because low-frequency g modes have primarily radial wavenumbers,  ≈   ∼ (/) ℎ ≫  ℎ , the radial part of the rotation vector ì Ω dominates in Equation 7. It is thus both convenient and accurate for many purposes to assume ì Ω ≈ Ω cos  r, i.e., to neglect the horizontal part of ì Ω.Given its usefulness, this assumption is appropriately known as the "traditional approximation of rotation" (TAR).
We emphasize that the TAR is only valid when   ≫  ℎ .It is therefore applicable when  ≪  and 2Ω≪ , with the interpretation that stratification is the dominant restorative force in the radial direction (such that the Coriolis force is only important in the horizontal directions).The utility of this approximation is that it allows the (buoyancy-driven) vertical dynamics to be decoupled from the (Coriolis-driven) horizontal dynamics.Because of this useful feature, the TAR has also found extensive use in geophysics (e.g., Eckart 1960;Longuet-Higgins 1968).However, if either of the hypotheses of the TAR above are not satisfied, the traditional approximation should be abandoned (Dintrans et al. 1999;Dintrans & Rieutord 2000;Gerkema & Shrira 2005;Ballot et al. 2010;Mathis et al. 2014).
When assumed, the TAR implies that the pressure perturbation varies in the horizontal directions according to the Laplace tidal equation: In the non-rotating limit ( → 0), the Laplace tidal equation approaches the usual generalized Legendre equation, for which  = ℓ(ℓ + 1) and the eigenfunctions are associated Legendre polynomials.Here,  = cos  is the colatitude and  = 2Ω/ is the spin parameter.When computing mode frequencies the TAR, the effect of rotation is thus simply to replace ℓ(ℓ + 1) with .
To handle the effect of a strong dipole magnetic field, Rui & Fuller (2023) borrow intuition from the TAR.The full JWKB dispersion relation for magnetogravity waves is where ì   = ì  0 / √︁ 4 0 is the Alfvén velocity, e.g., Unno et al. (1989) andFuller et al. (2015).Analogously with the rotational argument, we see that the radial part of ì  dominates, and it suffices for a dipole magnetic field to assume ì  ≈  0 cos  r.The pressure perturbation then follows where  =     /.The interpretation of this approximation is that the fluid is sufficiently stratified that buoyancy is the only important restorative force in the radial direction (i.e., the Lorentz force need only be included in the horizontal directions).As in the TAR, including magnetism in a calculation of mode frequencies under this approximation simply involves replacing ℓ(ℓ + 1) with a suitably computed  when solving the radial problem.We note the similar forms of Equation 8 (for rotation) and Equation 10 (for magnetism).However, unlike the singularities in Equation 8 (around which the eigenfunctions are smooth), the singularities in Equation 10 are of significantly different character, and imply sharp fluid features corresponding to resonances with Alfvén waves (Rui & Fuller 2023).For the frequency-shift analysis conducted in this work, this property of the singularities in Equation 10 motivates restriction to solutions for which  < 1 (so that the Alfvén resonances are not on the domain).
In this work, we generalize both the traditional approximation of rotation and its magnetic analogue to incorporate both effects: in other words, we consider only the effects of the radial components of both the rotation vector and magnetic field.Equivalently, we include only the horizontal components of the Coriolis and Lorentz forces.Hereafter, we refer to this joint approximation as the traditional approximation of rotation and magnetism (TARM).

Assumptions, conventions, and scope
In addition to assuming that   ≫  ℎ , we adopt the JWKB approximation in the radial direction only, i.e., we assume that the equilibrium structure and field of the star vary on length scales much larger than the radial wavelength (the "asymptotic" regime).Because such length scales are typically ∼ , this assumption is usually justified, although it may be violated in the presence of sharp compositional gradients which are known to produce mode-trapping phenomena (e.g., Miglio et al. 2008;Pedersen et al. 2018;Michielsen et al. 2019).In Section 5, we solve for the full radial dependence of the wavefunction without directly assuming that the radial wavenumber is large.However, under the TARM, we perform this calculation using a precomputed grid of horizontal eigenvalues  (see Section 4) which does make this assumption.Therefore, the calculation described in Section 5 is expected to partially, but not fully, capture non-JWKB effects in the radial direction.
We index branches by the angular degree ℓ and order .In particular, a mode is said to have some value of ℓ and  when the horizontal dependence becomes the spherical harmonic  ℓ when both the field and rotation are smoothly taken to zero.We caution that, while we may refer to some mode as having some degree ℓ in a rotating and/or magnetized star,  ℓ is not the correct horizontal dependence, and the eigenvalues are no longer ℓ(ℓ + 1).For the angular order , we adopt the sign convention used by Lee & Saio (1997) and Rui & Fuller (2023) that  > 0 ( < 0) corresponds to retrograde (prograde) modes.Additionally, without loss of generality, we consider throughout the case where  > 0 and  > 0 (which appear in, e.g., Equations 8 and 10, respectively), i.e., positive (negative) azimuthal order  corresponds to retrograde (prograde) modes.In this problem, the sign of  is irrelevant, and the effect of a sign change in  can be fully compensated by changing the sign convention of .
In the presence of (solid-body) rotation, it is important to distinguish the mode frequency in the inertial frame (which is observable) from the mode frequency in the frame co-rotating with the star (in which the effect of rotation appears as a Coriolis force).Hereafter, we use  ( ω) to denote the mode frequency in the inertial (co-rotating) frame.Hence, we calculate the oscillation modes directly with respect to ω, but convert to  for observational purposes.
We restrict our attention to a magnetic field whose radial part has a dipolar horizontal dependence.However, our results are not sensitive to the radial dependence of the field (as long as it is not very steep), or the geometry of the horizontal field components (as long as they are not much larger than the radial component).While Section 3.2 makes no additional assumptions about the field than those listed above, Section 5 requires a radial magnetic field profile  0 =  0 ().For this work, we adopt the Prendergast magnetic field geometry (Prendergast 1956).For our purposes, it suffices to specify the radial component of the magnetic field: where  1 () = sin / 2 −cos / is the first spherical Bessel function and  is the radius of the star.Although Kaufman et al. (2022) have recently shown that the Prendergast geometry is likely unstable over timescales relevant to stellar evolution, we adopt it simply as a closed-form model for a large-scale, dipole-like field, and we expect our findings to be insensitive to the exact radial dependence of the field.Following, e.g., Kaufman et al. (2022), we take Λ ≈ 5.76346 and  ≈ 1.31765, corresponding to the normalized, lowest-energy field solution with a vanishing surface field.Hereafter,   should be understood to refer to the radial component of the core magnetic field amplitude, although it is typically expected that the radial and horizontal components of the field are comparable.We expect all of the chief results of this work to be robust to magnetic field geometry, as long as  0ℎ / 0 ≲ / and   ≫ d ln  0 /d.We specialize to the case where magnetism is not strong enough to suppress the modes (although we explore the mode frequencies right up to this limit).While the suppression mechanism of magnetogravity waves is not fully understood, suppression may occur when magnetogravity waves refract upwards at some critical  =   to infinite wavenumber (Lecoanet et al. 2017(Lecoanet et al. , 2022;;Rui & Fuller 2023) or are damped out by phase-mixing processes once resonant with Alfvén waves ( > 1) in a manner similar to that described by Loi & Papaloizou (2017).Therefore, we restrict the scope of our calculations to the case where  =     / < 1 and  <   .Under these circumstances, the effects of magnetism on g modes should be well-modeled by our method.
For demonstrative purposes, we restrict most of our attention in this work to the dipole (ℓ = 1) and quadrupole (ℓ = 2) modes, although our calculations do not assume this, and it is not more complicated to extend this analysis to higher ℓ.Low-degree g modes suffer the least from geometric cancellation and are thus the easiest g modes to observe (there are no radial g modes).For simplicity, we assume modes are adiabatic, and neglect perturbations to the gravitational potential (i.e., we adopt the Cowling approximation).The general result that the perturbative theory underestimates the impact of magnetism on the period spacings for the dipole modes (Section 6) is also expected to hold for the quadrupole modes, although the asymmetry in the frequency shifts of different multiplets is known to behave differently (cf.Section 6.1 of Bugnet et al. 2021).

ANALYTIC FORMULATION
In this Section, we derive an expression for the horizontal equation obeyed by low-frequency g modes under the simultaneous influence of uniform (or weak differential) rotation and a dipolar magnetic field (Section 3.1).Under the TARM, the eigenvalues associated with these normal modes can be easily translated to an asymptotic expression for the period spacing (Section 3.2).

Fluid equations for gravity modes
In the presence of gravity, magnetic tension and pressure, and Coriolis forces, the linearized momentum equation is where ì  is the fluid displacement, subscript 0 and primes denote equilibrium and perturbed quantities respectively,  is the gravitational acceleration, and is the magnetic field perturbation.Equations 12 ignore the centrifugal force, and apply a Cowling approximation to neglect perturbations in .Under the TARM,   → −  when acting on a perturbation, and the magnetic tension term in Equation 12 thus becomes 1 4 where |  cos | is the radial component of the Alfvén velocity, with the angular dependence explicitly factored out.
In spherical coordinates and applying the traditional approximation, the momentum equation becomes where we have assumed harmonic time dependence,   →  ω, and used axisymmetry to take   → .Magnetic tension dominates over magnetic pressure in the asymptotic regime, and so the latter is ignored in Equations 15.
For adiabatic oscillations, the pressure and density  and  are related by where D/D = / + ì  • ∇ denotes the advective derivative.Equation 16 can be linearized to where   = √︁   0 / 0 is the sound speed.For gravity waves, the first term dominates, so that Finally, the fluid perturbation must satisfy the equation of continuity, so that where we have applied the Boussinesq approximation (Proctor & Weiss 1982).Now, the horizontal momentum equations give a linear system of equations for   and   in terms of  ′ : where  =     / ω (Rui & Fuller 2023) and again  = 2Ω/ ω (Lee & Saio 1997) are the dimensionless parameters governing the effects of magnetism and rotation on the horizontal eigenfunctions.Equations 20 can be solved to obtain where  = cos .Likewise, the radial component of the momentum equation (Equation 15a) can be solved to yield Substituting Equations 17, 21, and 22 into the continuity equation (Equation 19), we obtain where the differential operator L ,, is given by The operator L ,, further reduces to the standard Laplace tidal operator (e.g., Lee & Saio 1997) when  = 0 (no magnetism), and to the magnetic operator discussed by Rui & Fuller (2023) when  = 0 (no rotation).Hereafter, we define the "eigenvalues"  of L ,, as constants admitting solutions  () to i.e., the eigenvalues of L ,, in the "standard" sign convention are −.
When  =  = 0 (i.e., no magnetism or rotation), L ,, [  ()] reduces further still to the standard Laplacian operator on a sphere, where solutions to the associated boundary value problem are the spherical harmonics, indexed by integers ℓ,  with eigenvalues ℓ(ℓ + 1).Equation 23 where kℎ ≡ √ / is an effective horizontal wavenumber, which incorporates the effects of rotation and magnetism.By analogy with the spherically symmetric case, we may define an effective degree such that  = ℓ  (ℓ  + 1).In the TARM, oscillation modes are calculated by replacing ℓ with ℓ  throughout the entire star, in the same manner as is done in the standard TAR.

Asymptotic period spacing
In the absence of rotation and magnetism, gravity modes obey the dispersion relation ω = ± ℎ /  ∝  −1  .In the asymptotic regime (where    → ∞), this implies that adjacent g modes (with relative radial orders   = 1, and   ∼   /) are spaced uniformly in the mode period .In this Section, we derive an expression for the asymptotic period spacing   for g modes.We note that further departures from the asymptotic formula are expected when the stellar structure varies over a comparable radial scale to the wavefunction, or when there is mode mixing.
Before proceeding, we review a fundamental difference between the inclusion of uniform rotation and magnetism.For rotation, the fluid equations are solved by eigenfunctions whose shapes are solely parameterized by the spin parameter  = 2Ω/ ω, which can be calculated using stellar model parameters and the mode frequency, i.e., without knowledge of   .Observed spin parameters for intermediatemass g-mode pulsations range from  ≃ 0.1 to  ≃ 30 (Aerts et al. 2017).However, for magnetism, the parameter which controls the shapes of the eigenfunctions,  =     / ω, does depend on   (which varies mode-to-mode and with  in a complicated way).Fortunately, Equation 26 can also be rewritten where the parameter  (described by Rui & Fuller 2023) is given by This parameter is the squared ratio of the magnetogravity frequency   (Equation 1) to the mode frequency ( ∼  2  / ω2 ) and, conveniently, can be computed in terms of the stellar model and ω alone.By computing the horizontal eigenfunctions  as a function of  and then inverting Equation 28,  can be found as a function of .
To compute the period spacing in the co-rotating frame, we first observe that the radial phase   across the gravity mode cavity is where we have used Equation 26, and the integral is over the region of the star where  <  and  <  ℎ   .In Equation 30,   is the radial order, and   is a (here unimportant) phase offset.Adjacent modes (with   = 1) will thus have where we have neglected the frequency dependence of the bounds of the buoyancy integral.
Because  ∝ ω−1 ∝ P and  ∝ ω−2 ∝ P2 , Combining Equations 31 and 32 and solving for  P gives This approaches the well-known, zero-field, zero-rotation asymptotic formula in the relevant limit (Equation 4), as well as Equation 4of Bouabid et al. (2013) which was derived for the purely rotational case.
Equation 33 requires the calculation of ( ln / ln )  and ( ln / ln )  , where subscripts denote fixed variables with respect to the partial derivative.In Section 4.1, we compute  and its derivatives on a discrete, rectangular grid of  and .While ( ln / ln )  is easy to calculate numerically via a finite difference formula (since fixing  is straightforward), computing ( ln / ln )  is slightly trickier because it is harder to fix .Via Equation 28, we see that Using the identity that we obtain We use Equation 37 in our numerical calculation of   .
In the inertial frame, the observed frequencies  are related to ω under our sign convention by so that the periods  and P in the inertial and co-rotating frames are related by where  rot is the rotation period.The (asymptotic) period spacing measured by an observer is thus given by Thus, the inclusion of either rotation or magnetism will also leave distinct imprints on   as a function of mode period: understanding these signatures is crucial for extracting these properties from   .

NUMERICAL SOLUTIONS OF THE HORIZONTAL PROBLEM
In preceding sections, we have introduced an analytic formulation of the magnetorotational pulsation problem.However, applying the TARM to concrete predictions of oscillation spectra requires robust numerical solutions for the horizontal eigenvalues .We describe our numerical procedure for this calculation in this Section.

Numerical collocation scheme
Rui & Fuller (2023) calculate numerical solutions to the horizontal problem (Equation 25) in the nonrotating case ( = 0) by introducing a small artificial dissipation and using a relaxation scheme.While this method satisfactorily treats numerical pathologies associated with a singularity at critical latitudes  ± = ±1/ for large fields, it is computationally inefficient.Relatedly, because the coefficients of Equation 24 vary quickly across  ± , unreasonably large dissipation coefficients must be assumed to avoid needing prohibitively high resolution near those latitudes.
The more general form of Equation 25 that we consider here is still of the Boyd type (e.g., Boyd 1981), but now has solutions, and singular points, indexed by two parameters,  and .In particular, Equation 25 produces four additional singular points, obeying two of which may lie within the solution domain even for fields too weak to resonate with a given oscillation mode (i.e.,  < 1).We therefore seek an alternative solution strategy that is robust to the presence of such regular singular points.For  < 1 −  2 , no singularities lie on the domain, and it suffices to perform standard Chebyshev collocation on the real line (e.g., Wang et al. 2016).However, the collocation procedure must be modified somewhat to work for  > 1 −  2 .We note that since the Sturm-Liouville linear operator in Equation 24 is analytic, it may be treated as defining an ordinary differential equation on the complex plane.Solutions to the standard Sturm-Liouville problem on the real line coincide with those of this analytically continued problem, restricted to the real line.Thus, we may construct numerical solutions to the analytically continued problem on a contour on the complex plane, chosen to match the boundary conditions of the real problem on the interval  ∈ [−1, +1].Eigenvalues of the analytically continued problem will not depend on this choice of contour.Thus, the contour may be chosen to avoid the singular points that we have described above, and therefore to improve the numerical conditioning (e.g.stiffness) of the problem.We refer the reader to, e.g., Boyd (1985) for a more detailed examination of this procedure, and nature of the resulting solutions.
We find the standard collocation procedure to be sufficient for  = 0 and  = ±2 for any values of  ∈ [0, 2] and  ∈ [0, 1).However, solutions for the  = ±1 modes under this procedure are numerically badly behaved for  > 1 −  2 .In these cases, we perform a complex coordinate transformation from  to  given by and then solve the resulting problem using Chebyshev collocation on the interval  ∈ [−1, 1].This contour is chosen to share endpoints with the original real interval, while being tangent with the real line from the || > 1 (rather than || < 1) direction.
The eigenvalue  depends on the relationship between the mode frequency ω, rotational frequency Ω (via  = 2Ω/ ω), and magnetogravity frequency   (via  =  2  / ω2 ).Therefore, once  is computed for a given pair of  and , we retroactively compute  = / √ ), and regard  as being a function of  and .Because this procedure only produces values of  below some critical  crit =  crit () (corresponding to the maximum field which permits propagating magnetogravity waves), we excise two families of solutions: the Alfvén resonant ones for which  > 1 (which are expected to experience phase mixing, e.g., Loi & Papaloizou 2017), and those which lie on the "slow" branch described by Rui & Fuller (2023) (which are expected to approach infinite wavenumber).Within this work, we consider both such solutions to be "suppressed": we do not otherwise make claims about the degree of suppression or the mode frequencies of suppressed modes.
Figure 1 shows values of  computed for all dipole (ℓ = 1) and quadrupole (ℓ = 2) modes.In particular, for the zonal ( = 0) and retrograde ( > 0) modes calculated here, the critical magnetic field needed to cause mode suppression decreases with increasing rotation rate.This is because, for these branches,  increases relatively strongly with  (for  = 0,  ∝  2 ; Bildsten et al. 1996;Townsend 2003Townsend , 2020)).Therefore,  crit =  crit / √  decreases with .However, since larger rotation rates cause the prograde Kelvin modes (which have  = −ℓ) to attain larger horizontal scales ( decreases to a smaller constant value with , when  = 0), the critical field increases with increasing rotation rate.For the (ℓ, ) = (2, −1) case, the dependence of  on  and  is slightly more complicated, hence the non-monotonic behavior of the corresponding critical field with rotation rate.In any case, a straightforward prediction of this formalism is thus that different branches of modes should undergo suppression at different mode frequencies.Observational measurements of these critical periods may therefore impose strong constraints on the magnetic and rotational properties of the star.

Non-asymptotic numerical scheme
In the asymptotic regime, the perturbations vary with radius as ∼     , where   is given by Equation 30 (using the appropriate bounds).In other words, in this regime, the wavefunctions in the gmode cavity are expected to be sinusoidal with respect to a modified buoyancy radius Π given by which we define over the entire main radiative cavity (with respective inner and outer boundaries  1 and  2 ).This quantity is normalized such that Π ranges from 0 to 1.
However, the asymptotic assumption is violated in the proximity of sharp features in  (i.e., buoyancy glitches) when their characteristic widths are ≲  −1  .In such cases, the period spacing is expected to be modified from the asymptotic estimate in Equation 33.Sharp peaks in  are known to develop at the lower boundaries of the radiative envelopes of evolved MS stars (in which they cause periodic "dips" in the   pattern with ; Miglio et al. 2008;Pedersen et al. 2018), and similar buoyancy glitches have recently been observed asteroseismically in red giants (Cunha et al. 2015;Vrard et al. 2022), although their structure is very sensitive to the details of convective boundary mixing (e.g., Michielsen et al. 2021;Lindsay et al. 2023).The asymptotic assumption is also strongly violated for g modes  24) for the dipole (top) and quadrupole (bottom) modes.We plot  against the dimensionless parameters  =  2  / ω2 =    / ω2 and  = 2Ω/ ω, which govern the effects of magnetism and rotation, respectively.The eigenvalue  enters the asymptotic period spacing as in Equation 33.The turquoise lines show contours to the right of which the integrand in Equation 33deviates from a perturbative treatment by 10% and 50%.The pink hatched zones indicate  =     / > 1, i.e., modes which occupy these values of  and  at some layer within the star are likely to be suppressed.
. with low radial order, which may be observable in subgiants or some pulsators on the MS.
To model some of the non-asymptotic effects, we use a shooting method to solve the stellar pulsation equations under the assumption of adiabaticity, cast in the dimensionless form of Dziembowski (1971).This form of the pulsation equations is also employed by commonly used mode-solving codes such as GYRE (Townsend & Teitler 2013).Rotation and magnetism are implemented only by replacing the angular degree ℓ in the equations with an effective degree ℓ  , defined in Equation 27.Thus, we account only for the dynamical effects of rotation and magnetism, and neglect their indirect effects on stellar structure itself.Additionally, we emphasize that this "1.5D"approach still includes both rotation and magnetism asymptotically (similarly to the treatment of rotation in GYRE), and thus relies on the rotation and magnetic field profiles varying slowly in  compared to the wavefunctions themselves.In other words, while this procedure captures phenomena like wave-trapping due to peaks in , it does not accurately model the effects of sharp radial gradients in the magnetic field or rotation profiles, or coupling to, e.g., inertial modes.
In what follows,  and  denote the total mass and radius of the star, and  denotes the mass interior to radius .We solve the radial problem for where  = /, and ℓ  = ℓ  ( 1 ) is evaluated at the inner boundary.
In buoyancy coordinates, the perturbed time-independent oscillation equations then become where and Equations 45 reflect the -mode localization scheme of Ong & Basu (2020) as well as the Cowling approximation (neglecting perturbations to the gravitational potential).These approximations are made to restrict our attention to the effect of magnetism and rotation on pure g modes, and to avoid boundary condition-related numerical artifacts (see Section 2.2 of Ong & Basu 2020).Because the Cowling approximation is well-justified at high radial orders (where the TARM is valid), this approach should capture all of the robust predictions of our formalism.For the red giant model (Section 5.2), the resulting modes should be compared to the output of the stretching procedure typically used to extract   from solar-like oscillators (Mosser et al. 2015).
For our numerical shooting, we first integrate Equation 45 outwards from the stellar centre as an initial value problem to produce inner basis solutions which are consistent with the boundary conditions imposed there.In this work, we impose the boundary condition  1 = 0 (  = 0) on both boundaries.The solution vector evaluated at any intermediate point (here taken to be Π = 1/2) should thus be equivalent (up to linear dependence) when obtained by integrating from either boundary (starting from ì  = (0, 1)).These two solution vectors (obtained using a Radau integration scheme; Wanner & Hairer 1996) can then be formed into a 2 × 2 matrix whose determinant D ( ω) must vanish at a normal mode ω = ω * .
The adiabatic prescription of Equation 45 produces strictly real eigenvalues.To search for modes, we evaluate D ( ω) over some frequency grid.Between frequency grid points where D changes sign, we use a bisection algorithm to locate the roots of D. These oscillation modes ω are then converted to their values  in the inertial frame via Equation 39 (when there is rotation).

Stellar models
We find the oscillation modes of stellar models produced using version r22.11.1 of the Modules for Stellar Experiments (MESA) code (Paxton et al. 2010(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019)).We incorporate realistic convective overshoot using exponential overmixing with scale height  ov   = 0.015  (where   is the local pressure scale height), with the overshooting region starting a distance 0.005  inside the convective zone.
The stellar profiles as well as the rotation periods and magnetic fields we assume for them are summarized in Table 1.In particular, we choose three snapshots from a 1.5 ⊙ model to assess the behavior of the period spacing on the early-MS (MS-1.5-young),late-MS (MS-1.5-evolved),and lower RG (RG-1.5),and two snapshots from a 6.0 ⊙ star on the early-MS (MS-6.0-young)and late-MS (MS-6.0-evolved).These models are chosen to be representative of  Dor (MS-1.5-young,MS-1.5-evolved), slowly pulsating B-type (SPB; MS-6.0-young,MS-6.0-evolved), and red giant solar-like (RG-1.5)oscillators.We solve for the dipole (ℓ = 1) oscillation modes over a realistic range of frequencies.For RG-1.5, we compute these frequencies with both rotation and magnetism, as well as in the absence of either, in order to benchmark the prediction of perturbation theory (subsection 6.1).For the main-sequence models, the mode frequencies are computed three times, including the effects of magnetism and rotation both separately and simultaneously (subsection 6.2).The magnetic field is chosen to be strong enough to exhibit the effects of strong magnetic modification and suppression of some branch of oscillation modes.The mode period/frequency ranges shown in Table 1 are given in the inertial frame.When relevant, we solve only for co-rotating frequencies ω ≳ 0 to avoid the pile-up of g modes close to ω = 0.
Our models do not take into account distortions of the stellar structure due to centrifugal forces and magnetic pressure.While these effects are unlikely to matter in most observed  Dor and SPB stars (Henneco et al. 2021), they are likely to be important in rapidly rotating p-mode pulsators (such as  Sct stars, e.g., Lignières et al. 2006).

Strong fields in red giant cores
Strong magnetic fields in red giant cores have two main asteroseismic manifestations.First, they may produce frequency shifts on the nonradial modes which tend to shift modes of all  in the same direction (as opposed to rotation, which creates a frequency multiplet).Measurements of such frequency shifts have recently been used to make inferences about the field strength and, in one case, even geometry (Li et al. 2022b(Li et al. , 2023)).Second, if the magnetic field is extraordinarily strong, the magnetic field is expected to suppress the amplitudes of dipole modes whose frequencies lie below some  crit ∼   (Fuller et al. 2015;Lecoanet et al. 2017;Rui & Fuller 2023).
Our red giant model (RG-1.5;described in Section 5.2) is chosen to mimic a star on the lower red giant branch (for which mixed modes are easiest to observe) with a typical rotation rate ( rot = 30 d).For a frequency of maximum power  max ≈ 300 Hz, we calculate all dipole modes within the frequency range  max /2 and 3 max /2 in the simultaneous presence of magnetism and rotation, using the scheme described in Section 5.1.The width of the adopted frequency range is comparable to the full width at half maximum value  env ≈ 100 Hz ≈  max /3 calculated using the scaling relation of Mosser et al. (2012a).The large central magnetic field   ≈ 820 kG is chosen such that the  = ±1 sectoral modes are suppressed at the lower frequency range, to show the effect of a strong field.Note that   refers to the maximum value of the radial component of the field at the center of the star, rather than some horizontally averaged version of this quantity.Therefore, this value of   corresponds to a horizontally averaged field  2 The middle panels of Figure 2 show mock period echelle diagrams corresponding to these calculations.
We additionally calculate the mode frequencies for the same stellar model in the absence of either rotation and magnetism, in order to test the perturbative formalism.At high frequencies (where both rotation and magnetism are perturbative), the mode frequencies are closely consistent with the perturbative frequency shifts derived by Li et al. (2022b) (the red crosses in Figure 2).However, at low frequencies close to suppression ( ≲ 220 Hz), the TARM and perturbative results deviate substantially, with the TARM results tending to predict much larger frequency shifts than the perturbative formulae.This effect becomes increasingly dramatic until, at  ≈ 170 Hz, the sectoral modes are totally suppressed (although the zonal  = 0 mode remains propagating, and is suppressed at a frequency below the chosen observed frequency range).Disagreement between the perturbative and TARM frequency shifts is fully expected: at or near suppression, the effects of magnetism are, by definition, highly nonperturbative.
To formally demonstrate consistency with the perturbative formulae at high mode frequencies, we can expand the operator in Equation 24 in  and  and perform a perturbation analysis.Corrections to the subsequent analysis enter at O ( 3 ,  2 ,  4 ) ∼ O (Ω 3 , Ω 4  ,  8  ).We obtain the following eigenvalue equation: where Model MS-1.5-youngMS-1.5-evolvedRG-1.5 MS-6.0-youngMS-6.0-evolvedTable 1.Summary of stellar models for which we calculate oscillation modes using the non-asymptotic scheme described in Section 5.1. (50) To find the effect of L ,, pert on the eigenvalues, we perform firstorder perturbation theory.If the dipole eigenvalues are given by where and  ′0  () are the unperturbed eigenfunctions (of L  0 ).We emphasize that this is a perturbative expansion on the space of latitudinal functions all of the same  (for the generalized Legendre operators), and not on the full space of spherical harmonics (as done by Li et al. 2022b).Degenerate perturbation theory is thus not necessary here, since the eigenvalues of the generalized Legendre operator for a given  do not repeat.Furthermore, while in principle corrections may enter in an expression at second-order in perturbation theory, the only relevant term ∝ − in L ,, pert shifts all of the eigenvalues of a given  equally, and thus does not induce a second-order perturbation in .
The unperturbed pressure perturbations are the associated Legendre polynomials, where we have normalized the functions to square-integrate to unity and ignored the overall (Condon-Shortley) phase.The integral in Equation 52 can therefore be evaluated to give To transform the independent variable  to  (which can directly be specified given a field and stellar profile), we note that up to the relevant order.Then The mode frequencies in the co-rotating frame are given by ω = 1 in the asymptotic regime, where   = (  +   ) (Tassoul 1980) is the total radial phase (note that  depends implicitly on ω in a complicated way).We again proceed in ignoring the frequency dependence of the bounds of the integral in Equation 58(which should formally only enclose the part of the main radiative cavity where  < ).We define the "buoyant average": Assuming that  ω ≪ ω0 (sufficient for the desired order of the expansion), we may expand Equation 58as for  = 0, and for  = ±1.Equations 60 and 61 can be solved to yield the following frequency shifts: We keep one higher order of the rotation rate than do Li et al. (2022b).We distinguish between ⟨Ω 2 ⟩  and ⟨Ω⟩ 2  in the above to allow for the possibility of weak differential rotation (e.g., Beck et al. 2012), which may distinguish between the two.However, in the case of uniform rotation (assumed throughout this work), ⟨Ω 2 ⟩  = ⟨Ω⟩ 2  = Ω 2 .In the inertial frame, these frequency shifts become where ω0 =  0 for the unperturbed modes.We have full consistency with the perturbation formulae of Li et al. (2022b) (their Equations 61 and 62, with  = 1).Note that the star-averaged quantity which they define to be   (≡  ,L22 ) is equal to  ,L22 = ⟨ 4  ⟩  /3 ω3 0 .We caution that both the direct role of the centrifugal force as a restorative force and its indirect impact on the stellar structure (e.g., Ballot et al. 2010) also enter at ∝ Ω 2 .Inclusion of these effects is likely necessary to accurately capture the second-order effects of rotation.
Our non-perturbative mode calculations imply a few straightforward predictions.First, as mentioned previously, the magnetic frequency shifts become substantially stronger than implied by a perturbative estimate.While the relative change in the period spacing is still small (  decreases by ≈ 10% before suppression), the frequency shifts still substantially modify the period echelle diagram.Conversely, if the period spacing pattern of a strongly magnetic red giant is fit using the perturbative formulae, the inferred magnetic field is likely to be a significant overestimate.For example, Deheuvels et al. (2023) claim the detection of a red giant (KIC 6975038) whose magnetic field (≈ 286 kG) significantly exceeds the critical field  crit by a factor ∼ 1.7.Under our formalism, a field near or exceeding  crit should efficiently damp magnetogravity waves, either through phase mixing or refraction to infinite wavenumbers.Indeed, Deheuvels et al. (2023) observe nearly total dipole suppression in the same star for only low-frequency modes ≲  max , consistent with  crit lying on their observed frequency range.Their results could potentially be brought into accord with ours if non-perturbative effects have caused an observational overestimate of the field by a factor of a few.
To characterize the severity of such systematic overestimates, we compute the dipole frequency shifts in the red giant model for a range of internal magnetic fields (by numerically solving Equation 58).For each order , we then calculate the internal magnetic field which would be needed to produce the same frequency shift in perturbation theory.Figure 3 shows that the magnetic field ⟨ 2  ⟩ 1/2 implied by perturbation theory can exceed the "true" value for fields which are almost strong enough to cause suppression.Specifically, we use ⟨ 2 ⟩ 1/2 to denote the field averaged over all angles and over the radial kernel (following Li et al. 2022b): where  () is given by Equation 67 in the asymptotic limit.While the errors accrued by the perturbative formulae in Figure 3 are relatively small and do not rise to a factor ∼ 1.7, the degree to which perturbation theory overestimates the field likely depends on the field geometry adopted and the exact structure of the star (via, e.g., how far up the red giant branch the star is).Moreover, it likely depends on the exact procedure used to extract the field.For example, Figure 3 shows magnetic field values inferred using only one azimuthal order at one frequency, but an inference using the whole oscillation spectrum may yield a different answer.In the future, the manner in which perturbation theory misestimates the field should be characterized in more detail as a function of these factors.Large relative errors in the inferred magnetic field may also appear at low fields end if second-order rotational effects are mistaken for magnetic shifts (top panel of Figure 3).
Second, Li et al. (2022b) and Li et al. (2023) measure the dipole asymmetry parameter, defined by This should not be confused with the parameter  =  2  / 2 defined in this work and by Rui & Fuller (2023).In the perturbative regime, they show that where  2 () = (3 2 −1)/2 is the second-order Legendre polynomial and  () is a radial kernel function given by In particular, when the horizontal dependence of  0 is given by (, ) (i.e., the horizontal geometry is radius-independent), the radial integral in Equation 66 can be eliminated, yielding In the special case of a dipole magnetic field whose axis is aligned with the rotational axis ((, ) = cos ), it can be seen that  asym = 2/5 = 0.4 in this expression.
In the bottom panel of Figure 2, we see that this expectation holds at high frequencies, but increases slightly to ≈ 0.5 at lower frequencies (near  crit ).While likely difficult to measure, a value of  asym that varies towards lower frequencies (coinciding with the inference of a large magnetic field from the frequency shifts) may be an independent signature of a near-critical field.This non-perturbative asymmetry effect is related to the different magnetic fields implied by perturbation theory's predictions for the frequency shifts of different azimuthal orders (Figure 3).
In stars with especially weak magnetic fields, it is in principle possible for the dipole asymmetry to be dominated by rotation, even if it is slow enough for perturbation theory to be applicable.From Equation 63and Equations 65, we have such that, for a uniform rotation rate Ω, ⟨Ω 2 ⟩  = ⟨Ω⟩ 2  = Ω 2 ), Equation 69possesses the limiting behavior When the magnetic asymmetry dominates (⟨ 4  ⟩  / ω2 0 ≫ ⟨Ω⟩ 2  ≃ Ω 2 ),  asym ≈ 2/5 = 0.40.However, when the Coriolis-induced rotational asymmetry dominates (Ω 2 ≫ ⟨ 4  ⟩  / ω2 0 ), we instead have  asym ≈ 1/25 = 0.04.We stress that this is a fully perturbative effect: it only deviates from the result of Li et al. (2022b) because it includes a single higher-order effect of rotation.The upshot is that, even when both rotation and magnetism are individually small,  asym ≠ 2/5 for aligned rotational/magnetic axes if the effect of rotation is relatively at least comparable to that of magnetism.We again caution that the centrifugal force (which is also relevant at this order in Ω) has been neglected-this likely implies that the rotation-dominated asymmetry does not exactly approach 1/25 but some other value.Inclusion of such effects (as done by, e.g., Mathis & Prat 2019;Dhouib et al. 2021a,b) is needed to properly predict the true rotation-dominated asymmetry value.Nevertheless, we expect the qualitative ability for rotation to dominate over magnetism in determining the dipole asymmetry to be robust.Solid symbols denote mode frequencies calculated using the TARM, whereas hollow symbols denote the lowest-order prediction of perturbation theory.Bottom: The dipole asymmetry parameter (Equation 65) plotted against unperturbed mode frequency.Li et al. (2022b) and Li et al. (2023) neglect the rotational asymmetry effect on the basis that the core rotation rates in the stars in their sample are typical (i.e., low): we hereafter check this explicitly.As a crude estimate, the magnetic asymmetry dominates the rotational asymmetry in a red giant core when ⟨ 4  ⟩  / 2 max Ω 2 ≫ 1.In the three stars investigated by Li et al. (2022b), ⟨ 4  ⟩  / 2 max Ω 2 ≳ 10 2 and their asymmetries are thus indeed very magnetically dominated.Most of the stars reported by Li et al. (2023) have values of ⟨ 4  ⟩  / 2 max in the tens or hundreds.However, this parameter reaches a minimum for KIC 8540034, for which ⟨ 4  ⟩  / 2 max Ω 2 ≈ 9.In this star, rotation may affect the asymmetry parameter for lowfrequency modes (note the frequency dependence of ⟨ 4  ⟩  / 2 Ω 2 ) In general, magnetic domination of the dipole asymmetry may not be the case for giants with either fast core rotation rates or weak fields, and we caution against using  asym alone to make an inference of the field geometry without checking this criterion explicitly. ⟩ 1/2 associated with perturbation theory.⟨ 2  ⟩ 1/2 refers to an angle-and radial kernel-averaged field, following the notation of Li et al. (2022b) (see Equation 64).Bottom: The internal magnetic field ⟨ 2  ⟩ 1/2 implied by perturbation theory using the frequency shift for some angular degree , plotted against the "real" value (given by our non-perturbative TARM formalism).The frequency shift is evaluated using   ≃   for a physically realistic radial order   = 70 ( max ≈ 150 Hz), roughly the bottom of the frequency range shown in the period echelle diagram in Figure 2.

Strong fields threading the envelopes of main-sequence pulsators
Stars with masses ≳1.3 ⊙ have radiative envelopes and convective cores on the main sequence.Therefore, such stars may pulsate in g modes which are directly detectable, without needing to be disentangled from p modes as in solar-like oscillators.Such oscillators are ubiquitous: as discussed previously, they include  Dor (AF-type) and SPB (B-type) variables.The pulsations are driven by coherent mechanisms such as convective flux blocking (in  Dors; Guzik et al. 2002;Dupret et al. 2004) and the  mechanism (in B-type pulsators; Gautschy & Saio 1993;Dziembowski et al. 1993).This is in contrast to the broadband, stochastic driving present in solar-like oscillators (Samadi et al. 2015).Crucially, in these pulsators, there is no guarantee that the measurable modes are complete over some observed frequency range.The selection mechanism for mode excitation is poorly understood, and the asteroseismic power spectra are often sparse.Observational studies of such pulsators thus typically apply a forward-modeling approach based on the identified modes (e.g., Aerts et al. 2018), which rely on good models for predicting observed oscillation spectra.
In this Section, we primarily focus on the period spacing pattern   =   () for modes of a given .This is a standard observable in the study of main-sequence pulsators.The period spacing pattern is known to encode the rotation rate of the star (through an overall slope; Bouabid et al. 2013;Ouazzani et al. 2017), as well as the presence of buoyancy glitches (e.g., Miglio et al. 2008).The Brunt-Väisälä ( ), rotational (Ω), and magnetogravity (  ) frequencies, plotted in relation to the range over which we solve for mode frequencies.Bottom: The period spacing    versus period  in the inertial frame for the dipole modes, in the magnetic, rotating, and magnetic and rotating cases.Predictions for the asymptotic period spacing for the  = 1 branch (using Equation 33) are shown in solid red.
We also show predictions for the asymptotic period spacing handling rotation non-perturbatively but magnetism only perturbatively (using Equation 71; dashed blue curves).
We first calculate the dipole oscillation modes for two  Dor-like models, one near the zero-age main sequence (MS-1.5-young)and one near the terminal-age main sequence (MS-1.5-evolved),shown respectively in Figures 4 and 5.The chief difference between these models is that the convective core in the latter model has had time to develop a large compositional gradient at the base of its radiative envelope: this produces a jump in  (see the top panel of Figure 5).Qualitatively, this sharp feature in  results in a trapping phenomenon which results in a period spacing   which oscillates as a function of mode period  (Miglio et al. 2008;Pedersen et al. 2018;Vanlaer et al. 2023).We adopt a fairly typical core rotation rate of 1.5 d to accentuate the effects of rotation (Van Reeth et al. 2016;Li et al. 2020).Unlike in the red giant model described in Section 6.1 (where realistic rotation rates are small, such that  ≪ 1), rotation in the MS models is fast enough to cause frequency splittings/shifts which are nonlinear with respect to Ω.
The lower panels of Figures 4 and 5 show   versus   for the young and evolved 1.5 ⊙ models, under the effects of magnetism and rotation individually as well as simultaneously.First, since rotation distinguishes between prograde and retrograde modes, the slope it imparts onto the period spacing pattern is different for the  = +1 and  = −1 modes.In contrast, the oscillation modes are not sensitive to the overall sign of the magnetic field, and thus magnetism affects the  = ±1 modes identically (but still differently than the  = 0 mode).Moreover, while rotation produces values of   which vary fairly linearly with , magnetism produces a curvature in the pattern, especially near suppression.This effect is similar to what was demonstrated by Dhouib et al. (2022) in the case of a purely toroidal field.In particular, when the maximum allowed value of  =  crit is determined by connection to the evanescent region (rather than the presence of a critical Alfvén latitude), the asymptotic expression (in Equation 33) predicts that   sharply approaches zero at  ≈   .This is because the term ∝ d ln /d ln  in Equation 33 diverges at radii where the main magnetogravity wave branch connects to the slow branch described by Lecoanet et al. (2017) and Rui & Fuller (2023).In reality, there is not likely to be an infinitely dense forest of modes, since the asymptotic formula is based on a linear approximation which is likely to break down close to suppression.Nevertheless, the curvature is conspicuous, especially for the young model, where the period spacing drops from its high-frequency value by ≃ 50% near the critical frequency.Moreover, this curvature is apparent even when rotation is included alongside magnetism, with the added feature that fast rotation can cause the  = +1 and  = −1 modes to become magnetically suppressed at very different frequencies.This curvature effect on the period spacing pattern is very different than those caused by inertial-mode coupling in main-sequence convective cores (which manifest as isolated "dips"; Tokuno & Takata 2022) and mode-trapping near strong compositional gradients outside of those cores (which manifests as "oscillations"; Miglio et al. 2008).This sharp curvature feature is not adequately captured by any low-order perturbative treatment of magnetism.To make comparison to the perturbative prediction generous, we expand Equation 33around  = 0, while treating rotation non-perturbatively (through the traditional approximation of rotation, cf.Van Beeck et al. 2020).The effect of magnetism then enters the period spacing earliest through  2 ∝  4  / 4 ∝  2 / 2 crit (as predicted by Cantiello et al. 2016).Specifically, defining   to be the eigenvalue calculated including rotation only, we have where we have used In addition to lacking the suppression phenomenon entirely, the perturbative prediction (shown for the  = +1 mode as the bluedashed lines in Figures 4 and 5) dramatically underestimates the magnetic curvature predicted by the full TARM-based formalism.To further demonstrate this point, in Figure 1, we show contours where the perturbative estimate misestimates the integrand of the integral in the asymptotic formula (Equation 33) by 10% and 50%, respectively.As expected, departure from the full TARM calculation becomes increasingly severe close to suppression.Non-perturbative effects must therefore be taken into account predicting the frequency spectrum close to  crit .For example, the magnetic "sawtooth" pattern in the period spacing pattern predicted by some authors (Prat et al. 2019(Prat et al. , 2020a;;Van Beeck et al. 2020) was derived using perturbation theory at low frequencies, and preliminary results suggest that this feature does not appear once magnetism is incorporated non-perturbatively (Dhouib et al., in prep.).An important observation is that the magnetically induced curvature in the period spacing pattern is more conspicuous in the young model than in the evolved one.This is because the relative magnetic frequency shifts are primarily determined by the quantity / ω0 (as shown in Section 6.1), which is maximized when as many layers of the star have   ∼ ω0 as possible.However, within our physical picture, the entire oscillation mode becomes suppressed when even a small layer of the star has ω0 ≲  crit ∼   .Because  accounts for most of the variation of   ∝ √  (the Prendergast field we adopt varies comparatively more slowly with radius),   is a much broader function of  in the young model versus in the evolved one, where it is peaked at the composition gradient at the lower boundary of the radiative envelope.Therefore, the young model reaches a larger maximum value of ⟨ 4  ⟩ 1/4  / ω0 than the evolved one, and furthermore in general attains large values of ⟨ 4  ⟩ 1/4  / ω0 over a wider frequency range.This heuristic explanation is even stronger for higher-order terms in the perturbative expansion, which involve buoyant integrals of higher powers of  4  / ω4 0 .The magnetic curvature is in principle detectable even in evolved main-sequence pulsators, as long as it can be deconvolved from other effects.It should be noted that the typical uncertainties in  Dor period spacings in Kepler are small, comparable to the marker sizes of Figures 4 and 5 (Van Reeth et al. 2015;Li et al. 2020).Moreover, because of the sensitivity of the magnetic curvature in the period spacing pattern to the compositional profile, strongly magnetized main-sequence pulsators may be a promising avenue for constraining mixing processes.However, in nonasymptotic cases where sharp features in the buoyancy profile are expected, the limitations of the TARM must carefully be considered.
For completeness, we also examine young (MS-6.0-young)and evolved (MS-6.0-evolved)SPB analogues, with masses 6 ⊙ (Figures 6 and 7, respectively).The qualitative features of the period spacing pattern are similar, except that the peak in   at the base of the radiative region in the young model (due to the peak in ) exceeds the value of   throughout the rest of the cavity.Therefore, for similar reasons as in the evolved  Dor model, the curvature in the period spacing pattern due to magnetism is not as prominent as in the young  Dor model.
For illustrative purposes, we calculate the critical mode frequencies for a variety of internal magnetic fields and rotation rates, using the MS-6.0-youngmodel.Figure 8 shows the critical mode period  crit = 2/ crit for the dipole and quadrupole prograde sectoral modes.Interestingly, although rotation is expected to make prograde modes suppress at higher frequencies in the co-rotating frame (see Figure 1), higher rotation rates actually cause modes to suppress at lower frequencies in the inertial frame.Simultaneous knowledge of the suppression frequency for one identified mode branch together with the rotation rate should be sufficient to make a model-dependent estimate of the magnetic field at the interior of the star.Alternatively, while potentially challenging, simultaneous measurement of the suppression frequencies for two identified mode branches may be able to put a constraint on both the internal magnetic field as well as rotation rate.Because the shapes of the contours in Figure 8 are largely determined by change-of-frame effects (vis-à-vis Equation 39), the latter method is most viable when the two mode branches have different azimuthal order .
Roughly ∼ 10% of massive dwarf stars possess significant (inclined dipolar) fossil fields up to tens of kilogauss at their surfaces (Grunhut et al. 2016;Shultz et al. 2019).Such fields may be strong enough in the interiors of such stars to suppress low-frequency gmode oscillations.Recently, Lecoanet et al. (2022) attributed missing low-frequency modes in the magnetic SPB star HD 43317 (observed with CoRoT; Buysschaert et al. 2017Buysschaert et al. , 2018) ) to magnetic suppression caused by a near-core radial field   ≃ 500 kG.As in our MS-6.0earlymodel, suppression in their model occurs when  crit >  in the compositional peak in  at the base of the radiative cavity (see their Figure 2).Moreover, Aerts et al. (2017) predict that core dynamos in B-type (AF-type) pulsators may produce strong magnetic fields 20-400 kG (0.1-3kG) where non-perturbative magnetic effects may be realized.Magnetic g-mode main-sequence stars thus appear to be natural environments to observe g modes which are non-perturbatively modified by magnetism.
Pulsators in the  Dor mass range may also possess influential magnetic fields (Aerts et al. 2021).Surface fields of hundreds to thousands of gauss are typical of the enigmatic family of rapidly oscillating Ap-type (roAp) stars (Hubrig et al. 2004), and the magnetic field is believed to play an important role in the (still not fully understood) driving mechanism of their high-overtone p-mode oscillations (Gautschy et al. 1998;Balmforth et al. 2001).It has been speculated (e.g., by Handler 2011) and claimed (Balona et al. 2011) that some roAp stars may also pulsate in g modes (on the basis of overlap between roAp and  Dor stars on the Hertzsprung-Russell diagram).However, this is far from certain.On the basis of non-adiabatic mode calculations, Murphy et al. (2020) argue that high-order g modes are likely to be very efficiently damped, possibly explaining the current lack of observed hybrid  Dor/roAp pulsators.However, if roAp stars containing high-order g modes do turn out to exist, they would serve as ideal laboratories for strong magnetogravity waves.Moreover, the understanding of high-order magnetogravity waves presented in this work may extend some insight into the behavior of low-order magnetic g modes (for which the asymptotic limit is not appropriate).

Future prospects
This work presents a non-perturbative formalism for calculating the g-mode oscillation frequencies of a magnetized and rotating star, including both effects asymptotically (i.e., applying the TARM).We have considered only with the case where the magnetic field is dipolar and aligned with the rotational axis.As test examples we have only applied it to red giant cores and g-mode pulsators on the main sequence.Here, we describe future possible directions of study in relation to the TARM formalism, and potential extensions.
This work represents a joint generalization of the traditional approximation of rotation (Lee & Saio 1997) and an analogous approximation for a purely dipolar magnetic field (Rui & Fuller 2023), in order to non-perturbatively incorporate the effects of both.Generalizations of the traditional approximation have, in the past, also incorporated centrifugal distortion (Mathis & Prat 2019;Dhouib et al. 2021a,b), differential rotation (Ogilvie & Lin 2004;Mathis 2009;Van Reeth et al. 2018;Dhouib et al. 2021b), and axisymmetric toroidal fields, both with constant Alfvén and rotation frequencies (Mathis & De Brye 2011) as well as with more general field geometries together with differential rotation (Dhouib et al. 2022).Based on observational demands (or theoretical intrigue), it is likely possible to add any combination of these effects to the operator L ,, defined in Equation 24.Although  would then be a function of more than two dimensionless parameters, such an approach would retain much of the advantage of non-perturbatively capturing complex rotational/magnetic effects while only interpolating over a precomputed eigenvalue grid.
Unlike Rui & Fuller (2023), this work has focused on the regime where suppression is not likely to occur, i.e., when there are no Alfvén resonances on the domain and where the slow magnetic branch has been ignored.We have ignored modes with these effects because their observational implications are unclear, but the behavior of the operator L ,, in this regime is an extremely rich mathematical problem with so far unexplored structure.Rui & Fuller (2023) find that solutions with  > 1 develop sharp fluid features at the Alfvénresonant critical latitudes, where processes such as phase-mixing are likely to efficiently damp the waves.In this regime, the magnetic operator in Equation 10 is of Boyd-type (Boyd 1981), and the interior singularities give dissipation an important role in determining the physically appropriate branch cut.The eigenvalues  for the  > 1 are thus not guaranteed to be real even in the formal limit where dissipation is taken to zero (and the numerical results of Rui & Fuller (2023) suggest that they are not).For reasons of scope, we have also ignored magneto-Rossby waves and magnetically stabilized gravity waves (Rui & Fuller 2023), which do not connect to any spherical harmonic in the limits ,  → 0. These, too, may conceal detectable predictions which are implied by the breakdown of positive-(semi)definiteness of L ,, .As such, our calculations also do not capture the coupling between magnetic g modes and magneto-inertial modes which propagate in the convective core of intermediate-mass main-sequence stars (within which dynamo-generated magnetic fields are expected; Brun et al. 2005;Featherstone et al. 2009).Coupling with inertial modes is known to result in isolated dips in the   - diagram at frequencies corresponding to those of inertial modes.This effect provides a seismic probe of the core rotation rates of such stars (Ouazzani et al. 2020;Saio et al. 2021;Tokuno & Takata 2022).In the future, it may be interesting to explore how this picture is modified by magnetism, and whether similar inference of the magnetic field in these convective cores is possible.We emphasize that coupling to (magneto-)inertial waves produces localized dip features in the period spacing pattern, and is very different than the global curvature in the pattern predicted by this work.
While we have only explicitly modeled analogues of  Dor and SPB stars, our analysis applies to any magnetized pulsator with pulsations of high radial order.This includes compact pulsators such as white dwarfs and hot subdwarfs.Since both of these species result from red giants whose envelopes have been lost (either in isolation or through binary evolution), it is natural to expect that they will retain the strong fields believed to cause dipole suppression in red giants.While a small handful of magnetized hot subdwarfs (100s of kG) are known (Pelisoli et al. 2022), white dwarfs with kilogauss surface fields are believed to make up a fourth of all white dwarfs (Cuadrado et al. 2004;Valyavin et al. 2006), and a number of magnetized white dwarfs with fields up to hundreds of megagauss have been discovered (Kepler et al. 2013;Bagnulo & Landstreet 2021).The latter fields are likely to be so strong that they outright suppress g mode oscillations altogether (Lecoanet et al. 2017).However, it may be possible for a white dwarf to have a field strong enough to significantly shift the frequencies of the g modes, without being not strong enough to suppress them outright.
While a dipolar field is expected at the surfaces of stars with fossil fields (Braithwaite & Nordlund 2006;Duez & Mathis 2010), that field need not be aligned with the rotation axis (Duez 2011;Keszthelyi 2023), and is unlikely to be dipolar at all if the field is generated by a dynamo.In the perturbative regime, Mathis & Bugnet (2023) recently characterized the frequency shifts associated to an inclined dipole field.Extending the TARM formalism to describe a nonaxisymmetric horizontal field dependence requires solving for the eigenvalues of families of two-dimensional differential operators over the sphere, rather than a one-dimensional one (as in L ,, ), and this analysis would need to be repeated for every different horizontal field dependence desired.Nevertheless, near suppression, departures in the frequency shifts from the perturbative theory are likely, and may be required for accurate magnetic field inference in this regime.
Finally, low-frequency propagating gravity waves are one of the best candidates for the strong angular momentum transport needed in stellar radiative zones to reproduce the observed internal rotation revealed in all types of stars by helio-and asteroseismology (e.g., Schatzman 1993;Zahn 1997;Charbonnel & Talon 2005;Aerts 2015;Rogers 2015;Pinçon et al. 2017;Neiner et al. 2020).The manner in which this wave-mediated angular momentum transport occurs can be significantly modified by the presence of a magnetic field.In general, the net angular momentum flux implied by this mechanism is given by the sum of the wave-induced Reynolds and Maxwell contributions to the stress tensor.The relevant gravity waves are precisely those which are most strongly affected by the combined action of rotation and magnetism (see, e.g., Mathis & de Brye 2012, in the case of weak, shellular differential rotation and a purely toroidal field with constant Alfvén frequency).Because our TARM-based formalism is relevant to exactly this kind of wave, its application to this problem is likely to yield insights into the rotational state and internal chemical mixing of rotating, magnetic stars.

CONCLUSION
Rapidly evolving progress in observational magnetoasteroseismology demands refinements in our theoretical understanding of magnetic effects on stellar pulsations.In this work, we develop a formalism for incorporating the effects of an aligned dipole magnetic field into g mode calculations, valid for rapidly rotating stars.This method relies on an asymptotic treatment of magnetism and rotation (under a "traditional approximation of rotation and magnetism"), and can be partitioned into two main steps: (i) Calculate the eigenvalues  of the horizontal differential operator L ,, (Equation 24) as a function of the dimensionless magnetic and rotational parameters  =  2  / ω2 and  = 2Ω/ ω. (ii) In either an asymptotic mode formula (Equation 58) or a nonasymptotic numerical scheme (e.g., shooting; Section 5.1), include the effects of magnetism and rotation by replacing ℓ(ℓ+1) throughout the star with a suitably interpolated , calculated using the magnetic and rotational profiles.
These steps are done relatively independently of each other: once the eigenvalues  are computed once over a sufficiently large grid of  and  (for the desired ℓ and ), they do not need to be calculated again for any individual stellar model.Moreover, modifications to existing mode solving procedures are "minimal" in the sense of being localized to the interpolation of  and its substitution into the relevant equations.
As proofs of concept, we have computed the g modes in the cores of red giants as well as in the radiative envelopes of high-mass mainsequence stars.In both cases, strong magnetic fields tend to decrease the period spacing significantly more than is suggested by the perturbative theory, especially for low frequencies close to the critical frequency  crit ∼ √︁   /.This results in a curvature in the period spacing pattern which can in some cases be very conspicuous (e.g., Figure 4).Non-perturbative effects may also introduce asymmetry in the dipole frequency shifts which is not predicted by perturbation theory.
This regime is expected to be directly realized in the SPB star described by Lecoanet et al. (2022) and some of the red giants described by Deheuvels et al. (2023).Refined understanding of these effects is therefore prerequisite to perform accurate magnetic field inference using asteroseismology.

Figure 1 .
Figure 1.The eigenvalues  of the differential operator L ,, (Equation24) for the dipole (top) and quadrupole (bottom) modes.We plot  against the dimensionless parameters  =  2  / ω2 =    / ω2 and  = 2Ω/ ω, which govern the effects of magnetism and rotation, respectively.The eigenvalue  enters the asymptotic period spacing as in Equation33.The turquoise lines show contours to the right of which the integrand in Equation33deviates from a perturbative treatment by 10% and 50%.The pink hatched zones indicate  =     / > 1, i.e., modes which occupy these values of  and  at some layer within the star are likely to be suppressed.

Figure 2 .
Figure 2. Top: The Brunt-Väisälä ( ) and magnetogravity (  ) frequencies for the red giant model (RG-1.5),plotted in relation to the range over which we solve for mode frequencies.The rotational frequency Ω ≃ 2.4 Hz ( rot = 30 d) is below the bottom bound of this plot.Center: Period echelle diagram for the red giant's core g modes.The right panel zooms into the low frequency modes of the left panel, and folds on a different period for clarity.Solid symbols denote mode frequencies calculated using the TARM, whereas hollow symbols denote the lowest-order prediction of perturbation theory.Bottom: The dipole asymmetry parameter (Equation65) plotted against unperturbed mode frequency.

Figure 3 .
Figure 3.For the red giant model (RG-1.5):Top: The relative error on the inferred magnetic field ⟨ 2  ⟩ 1/2 associated with perturbation theory.⟨ 2  ⟩ 1/2 refers to an angle-and radial kernel-averaged field, following the notation ofLi et al. (2022b) (see Equation64).Bottom: The internal magnetic field ⟨ 2  ⟩ 1/2 implied by perturbation theory using the frequency shift for some angular degree , plotted against the "real" value (given by our non-perturbative TARM formalism).The frequency shift is evaluated using   ≃   for a physically realistic radial order   = 70 ( max ≈ 150 Hz), roughly the bottom of the frequency range shown in the period echelle diagram in Figure2.

Figure 4 .
Figure 4. Characteristic frequency profiles and mode frequencies for a young  Dor analogue (MS-1.5-young).Top: The Brunt-Väisälä ( ), rotational (Ω), and magnetogravity (  ) frequencies, plotted in relation to the range over which we solve for mode frequencies.Bottom: The period spacing    versus period  in the inertial frame for the dipole modes, in the magnetic, rotating, and magnetic and rotating cases.Predictions for the asymptotic period spacing for the  = 1 branch (using Equation33) are shown in solid red.We also show predictions for the asymptotic period spacing handling rotation non-perturbatively but magnetism only perturbatively (using Equation71; dashed blue curves).

Figure 8 .
Figure8.The critical period  crit against the rotation period  rot for a young SPB-like model (MS-6.0-young),for fixed values of the field near the compositional gradient at the base of the radiative envelope (which most easily experiences magnetic suppression). crit is given in the inertial frame.