© The Authors 2026

1. Introduction

Magnetic fields can be found at all stages of stellar evolution, during the pre-main sequence (Villebrun et al. 2019), the main sequence (Donati & Landstreet 2009; Vidotto et al. 2014), the red giant phase (Aurière et al. 2015; Cantiello et al. 2016), and in compact stars such as white dwarfs (Ferrario et al. 2015) or neutron stars (e.g., Thompson & Duncan 1993; Raynaud et al. 2020). Magnetic fields are expected to impact stellar evolution by redistributing angular momentum, and they are considered as one of the most promising candidates to solve the long-standing problem of angular momentum transport in stars.

The seismology of subgiants and red giants have recently provided invaluable constraints to make progress on this question. These stars harbor oscillation modes with a mixed character, behaving both as pressure modes (p modes) in the envelope and as gravity modes (g modes) in the core. Mixed modes enabled measurement of internal rotation along all stages of the post-main-sequence evolution, i.e., for subgiants (Deheuvels et al. 2012, 2014, 2020), red giant branch stars (Beck et al. 2012; Mosser et al. 2012b; Gehan et al. 2018; Li et al. 2024; Hatt et al. 2024), and core-He burning giants (Deheuvels et al. 2015; Mosser et al. 2024). At all stages, these measurements have shown that the cores of post-main-sequence stars spin much more slowly than predicted by purely hydrodynamical processes of angular momentum transport (Ceillier et al. 2013; Marques et al. 2013), demonstrating the need for an additional mechanism redistributing angular momentum. Magnetic fields could efficiently transport angular momentum either directly by applying a magnetic torque on the stellar core (Maeder & Meynet 2014) or through turbulence associated with magnetohydrodynamic instabilities, such as the Tayler instability (Tayler 1973; Spruit 2002; Fuller et al. 2019; Petitdemange et al. 2023; Barrère et al. 2023) or the magneto-rotational instability (Rüdiger et al. 2015; Jouve et al. 2015; Meduri et al. 2024). The lack of observational constraints on internal magnetic fields has been a major obstacle to exploring the hypothesis of a magnetic transport of angular momentum.

Mixed modes were recently used to obtain the first observational constraints on the strength and topology of internal magnetic fields (Li et al. 2022). The signature of magnetic fields on mixed mode frequencies was first derived using a perturbative approach in the special case of a dipolar field with various radial profiles either aligned on the rotation axis (Gomes & Lopes 2020; Mathis et al. 2021; Bugnet et al. 2021) or inclined (Loi 2021). Li et al. (2022) obtained a more general expression of magnetic frequency shifts for magnetic fields with an arbitrary configuration, and they detected these signatures in three Kepler red giants. Overall, magnetic fields have now been detected in 71 red giants using Kepler data (Li et al. 2022, 2023; Deheuvels et al. 2023; Hatt et al. 2024; Villate et al. 2026). The detected fields have intensities ranging between about 10 and 600 kG in the deep core, and they show a wide diversity of field geometries. We also note the recent detection of a mainly toroidal magnetic field in the deep interior of a main-sequence F star (Takata et al. 2026). Complementarily, the existence of a population of red giants with suppressed dipole and quadrupole mixed modes (Mosser et al. 2012a; García et al. 2014; Stello et al. 2016) has been interpreted as the consequence of strong core magnetic fields that exceed the critical strength (B_c) above which magneto-gravity waves can no longer propagate (Fuller et al. 2015). When these waves reach regions with super critical fields, they are thought to be refracted upward into outgoing slow magnetic waves (Lecoanet et al. 2017; Loi & Papaloizou 2017; Rui & Fuller 2023; David et al. 2026) and eventually dissipated, which would cause mode suppression. However, the detection of partially suppressed dipole modes that retain a g-like character (Mosser et al. 2017) seems incompatible with this interpretation and remains to be explained.

Magnetic frequency shifts strongly depend on the frequency (ω) of the mode (they vary as ω⁻³) so that strong fields significantly modify the regular period spacing of high-order g modes at low frequency (Li et al. 2022; Bugnet 2022). Such departures from regular g-mode period spacing were detected by Deheuvels et al. (2023) in 11 Kepler red giants and attributed to strong near-critical core magnetic fields. All of these stars show a single component per dipole triplet, identified as the m = 0 component. In this work, we identify eight Kepler red giants whose dipole modes show a distortion of the g-mode pattern that is very similar to the stars studied in Deheuvels et al. (2023), but these stars show l = 1 doublets. We investigate the magnetic origin of these features.

In Sect. 2, we extract mode frequencies from the oscillation spectra of the eight red giants. We then show in Sect. 3 that the most obvious identification of the detected components as m = ±1 yields poor agreement with the observations, and we explore the identification as m = 0 and m = 1 components. To produce the observed distortion in the g-mode pattern, the magnetic fields need to be near-critical so that their effects on mode frequencies cannot be treated with a perturbative approach. Rui & Fuller (2023) and Rui et al. (2024) developed a formalism inspired from the traditional approximation of rotation in order to take magnetic effects into account in a non-perturbative manner. In Sect. 4, we adapt this formalism to our case. In Sect. 5, we show that measuring the seismic signature of near-critical fields can yield constraints on the radial profile of the radial magnetic field B_r. We apply this formalism to the targets of our sample in Sect. 6, and we discuss the results in Sect. 7.

2. Red giants with strongly distorted doublets

In Li et al. (2024), we measured g-mode properties and internal rotation rates for red giants in the catalog of Yu et al. (2018). In the process of this work, we identified a subsample of red giants showing strongly distorted g-mode patterns.

2.1. Stretched échelle diagrams

While high-radial-order g modes are nearly regularly spaced in period, mixed modes are not. Mosser et al. (2015) proposed to apply a “stretching” to the mixed mode periods in order to make them equally spaced by ΔΠ₁. For this purpose, stretched periods τ are computed by solving the differential equation dτ = dP/ζ, where ζ is defined as the fraction of the mode kinetic energy that is enclosed in the g-mode cavity (ζ tends to 1 for pure g modes, and 0 for pure p modes). In the nonmagnetic case, we expect stretched periods of same l to align on a vertical ridge when plotting them in an échelle diagram folded with the period spacing ΔΠ₁. For eight Kepler red giants, the resulting stretched échelle diagrams show two strongly distorted ridges, as shown in Fig. 1. At this stage, no evidence has yet been found of such a distortion in red giants that show more than two components per l = 1 mode.

Fig. 1. Stretched period échelle diagram for KIC 7458743 obtained from 4-yr Kepler time series. Only peaks with an S/N above eight are shown. The size of the symbols are proportional to the peak amplitude. For clarity, peaks in the vicinity of l = 0 and l = 2 modes are omitted. The échelle diagram is folded with an apparent period spacing of ΔΠ1meas = 76.5 s, which produces the best visual alignment of the ridges.

2.2. Extraction of mode frequencies

We performed a full seismic analysis of the eight Kepler red giants of our sample. The light curves were taken from Kepler Light Curves Optimized For Asteroseismology (KEPSEISMIC) and downloaded from Mikulski Archive for Space Telescopes (MAST). KEPSEISMIC light curves were corrected from outliers, jumps, and drifts through the implementation of the Kepler Asteroseismic Data Analysis and Calibration Software (KADACS, García et al. 2011). The gaps shorter than 20 days were filled following in-painting techniques described in García et al. (2014) and Pires et al. (2015). The power density spectra (PSD) were then obtained using the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982). The background of the PSD is modeled as a sum of the contributions from the granulation (two Harvey-like profiles, Harvey 1985, with distinct characteristic timescales) and the photon noise (white noise). To obtain an estimate of the background, we fit this model completed with a contribution from the oscillations (modeled as a Gaussian function) to the observed PSD with a maximum likelihood estimation (MLE) method.

To identify significant peaks in the PSD, we searched for peaks that exceed the threshold x_p corresponding to a false-alarm probability of p = 0.05 over the whole frequency range where oscillation modes are expected (ν_max ± 5 Δν). These peaks can be used to plot stretched échelle diagrams (see Fig. 1). The list of significant modes was then extended to include additional peaks from the PSD that align on the two clear ridges identified in the stretched échelle diagram, while having a power exceeding eight times the background level (see Mosser et al. 2015). Using the identified significant peaks, we obtained a list of guessed frequencies for dipole mixed modes. We then proceeded to fit a model of the oscillation spectrum to the observed PSD. Each mode was modeled as a Lorentzian function with free central frequency ν_n, l, m, line width Γ_n, l, m and height H_n, l, m. We used an MLE approach to fit this model to the PSD. The list of measured frequencies are given in Appendix A for all eight stars. The results are shown in the shape of stretched échelle diagrams in Fig. 2. As mentioned above, the ridges are not straight, so that we can only determine an apparent value of the period spacing ΔΠ₁^meas that provides a rough visual alignment of the ridges. The values of ΔΠ₁^meas for each star are given in Table 1.

Fig. 2. Stretched period échelle diagrams of the eight Kepler red giants of our sample. The échelle diagrams were folded using apparent period spacings ΔΠ1meas that offer good visual alignment of the ridges. To guide the eye, the ridges are highlighted with different colors (red for the most distorted ridge, blue for the other).

2.3. Distorted doublets

Figure 2 confirms that the eight stars under study all have l = 1 doublets that show strong distortion in the stretched échelle diagram. For clarity, two distinct colors have been used to highlight the ridges, but at this stage this is only to guide the eye. In the following sections, no assumption has been made about modes belonging to one ridge or the other. Several observations can be made about the stretched échelle diagrams:

• In all cases the two ridges either cross or get closer to each other in the high-frequency part of the échelle diagram.

• The intensity of the curvature seems different for the two ridges (in Fig. 2, we represent the ridge showing the strongest curvature in red).

• The three stars for which the ridges are the most curved (KIC 3326392, KIC 9227589, and KIC 9697268) show hints of mode suppression, the most distorted ridge (the red ridge) vanishing at low frequencies, while the less distorted one (the blue one) remains.

• There seems to be a difference in mode height between the two ridges. This can be seen by eye for KIC 7458743 in Fig. 1. Our MLE fits confirm that the oscillation modes of the red ridge have heights that are on average 1.9 times larger than those of the blue ridge for this star.

2.4. Location in the Δν − ΔΠ₁ plane

Deheuvels et al. (2022) showed that red giants collapse on a clear sequence in the Δν − ΔΠ₁ plane (the “degeneracy sequence”) when the electron gas becomes degenerate in the core. The strongly magnetic stars of Deheuvels et al. (2023) had a measured period spacing ΔΠ₁^meas that was too low, so that these stars appeared below the degeneracy sequence in the Δν − ΔΠ₁ plane, where no stars are expected. Accounting for the magnetic field enabled to determine the true period spacing ΔΠ₁, which placed these stars back on the degeneracy sequence, as expected (Deheuvels et al. 2023). In Fig. 3, we place our stars in the Δν − ΔΠ₁ plane (blue star-shaped symbols) using the apparent period spacing ΔΠ₁^meas, and their position can be compared to the degeneracy sequence (well identified with the gray crosses). The degeneracy sequence is mass-dependent, higher-mass stars having a lower ΔΠ₁ for a given Δν. For each star of the sample, we determined the equation of the degeneracy sequence corresponding to its stellar mass (taken from seismic scaling relations, see Yu et al. 2018). For this purpose, we selected stars among the catalog of Li et al. (2024) that have the same mass within 1σ errors. We fitted a third-degree polynomial to the relation ΔΠ₁(Δν) for these stars, giving the degeneracy sequence for the target mass. The difference z between the measured period spacing ΔΠ₁^meas and the one predicted from the degeneracy sequence is given for each star in Table 1. We find that all our targets are located from 1.6 s to 9.8 s below the degeneracy sequence. The typical scatter around the degeneracy sequence does not exceed 4 s. We thus conclude that six out of eight of our targets can be considered significantly below the degeneracy sequence, confirming the similarity between our targets and those of Deheuvels et al. (2023).

Fig. 3. Location of the stars under study in the (Δν-ΔΠ1) plane. The apparent period spacings of g modes ΔΠ1meas are shown as blue star symbols, and the red star symbols indicate the asymptotic period spacing ΔΠ1 measured by accounting for magnetic fields in a non-perturbative manner (Sect. 6). Gray crosses correspond to all red giants studied in Li et al. (2024).

2.5. Origin of ridge distortion

In Deheuvels et al. (2023), we identified two potential sources for the ridge distortions: buoyancy glitches (corresponding to sharp variations in the profile of the Brunt-Väisälä frequency, with length scales shorter than the mode wavelength), or internal magnetic fields. We ruled out buoyancy glitches, which would produce very different modulation in the ridges (see Appendix A of Deheuvels et al. 2023) and fail to account for the abnormally low measured period spacings. These arguments also hold in the present case. Moreover, with buoyancy glitches, we would expect the two ridges of the doublet to show similar deviations. Indeed, to first order, the modes of each doublet have the same eigenfunctions, so that they should be affected by a buoyancy glitch in the same way. The fact that we observe different curvatures for the two ridges in all our stars makes is even less likely that buoyancy glitches are causing the observed ridge distortion. In the following sections we explore the magnetic hypothesis.

3. Mode identification

If the magnetic field is not axisymmetric (that is, if it does not have an axis of symmetry that coincides with the rotation axis), it splits every l = 1 mode into three components, each of which being also split in three components by the effect of rotation. In principle, a nonaxisymmetric field can thus produce nine components per multiplet (Gough 1990; Loi 2021). Although some of these components are likely to have much smaller amplitudes than the others, we still expect to have more than three components per multiplet if the effects of nonaxisymmetry of the field are nonnegligible. The fact that only two components per multiplet are observed suggests that the nonaxisymmetry effects are indeed weak. We thus discuss the mode identification under this assumption.

3.1. Identification as m = ±1 components

The most natural mode identification for an l = 1 doublet is that we are seeing the m = ±1 components. This occurs when the inclination angle between the rotation axis and the line of sight is close to 90°, that is, the star is seen equator-on. In this case, the m = 0 component geometrically cancels. For a predominantly radial magnetic field, the mode frequencies depends on |m|, so that the m = −1 and m = +1 are affected in the same way by the magnetic field. Only rotational effects can thus produce different curvatures for the two ridges. If the core rotation is slow, the rotational splitting is small and the m = ±1 ridges are nearly parallel, contrary to the observations. A faster core rotation can make the two ridges cross (Gehan et al. 2018), as we observe in most of our stars. However, for all eight stars, we found it impossible to satisfactorily fit the shape of the two ridges and the frequency at which they cross. In Fig. 4, we show the example of KIC 7458743, for which the two ridges cross around a frequency of ν_cross ≈ 150 μHz. The core rotation needed to produce a crossing of the two ridges at ν_cross is expressed as ⟨Ω⟩_g/(2π) = ν_cross²ΔΠ₁ (Gehan et al. 2018, the presence of a magnetic field is expected to modify this expression, but its effect is negligible for KIC 7458743). This gives a core rotation of ⟨Ω⟩_g/(2π)≈1.78 μHz (which would make this star one of the fastest rotators in the Kepler catalog, Li et al. 2024). Using the method that is described in Sect. 6, we fit an asymptotic expression of mixed mode frequencies including the effects of rotation and magnetic fields in a non-perturbative way to one of the observed ridges, taking ⟨Ω⟩_g/(2π) = 1.78 μHz. The results are shown in Fig. 4 in a stretched échelle diagram. As can be seen, one of the ridges is very well reproduced, but the second one is completely off. The difference of slopes between the two ridges at the frequency of the crossing is very different in the model and in the observations. If we fit the two ridges simultaneously, considering the core rotation rate and the magnetic field properties as free parameters, the best solution is obtained when the ridges cross around 130 μHz and the quality of the fit is poor. We reach the same conclusion for all eight stars of the sample, we find it impossible to reproduce the observations if we identify the ridges as the m = ±1 components.

Fig. 4. Stretched échelle diagram of KIC 7458743. Filled blue circles show the observed mode frequencies. Open red circles show asymptotic mixed mode frequencies including rotation and magnetic field effects, assuming that the observed ridges correspond to m = ±1 components. The core rotation was determined to ensure a crossing of the ridges at νcross = 150 μHz, and magnetic field properties were adjusted to reproduce the m = −1 ridge (see Sect. 3.1).

Moreover, the identification as m = ±1 components fails to account for the difference in mode amplitudes between the two ridges, as observed for KIC 7458743 (the m = ±1 components have the same average over the observed disk, so they are expected to have the same amplitude). It does not explain either the suppression of only one of the two ridges, as observed in three stars of the sample (since the magnetic field affects m = ±1 components in the same way, they should be identically suppressed). For all these reasons, we rule out the identification of the ridges as m = ±1 components.

3.2. Identification as m = 0 and m = 1 components

The alternative identification is that one of the two ridges is the m = 0 component. This offers much more flexibility to fit the observed ridges because in general the magnetic field has a different impact on m = 0 and m = ±1 components. The inclination angle between the stellar rotation axis and the line of sight would need to be intermediate, so that both m = 0 and m = 1 components are visible. This configuration could account for the fact that the heights of the two ridges differ in some cases, for example for KIC 7458743. However, it raises the question of why we do not observe three ridges. Two hypotheses can be made to potentially account for this.

First, the core rotation could be too slow to separate the m = ±1 components. This would occur if the core rotation rate is smaller than the typical uncertainties to which mode frequencies can be estimated in the Kepler data. Typically, these uncertainties reach values of about 10 nHz for the most g-like modes near ν_max, which is close to the frequency resolution with 4 years of data (ν_res = 7.9 nHz). In this case, these stars would be exceptionally slow rotators compared to the bulk of Kepler red giants. Indeed, among the 2006 red giants in the catalog of Li et al. (2024), 90% of stars have core rotation rates between 380 nHz and 1170 nHz, and only 8 stars have core rotation rates slower than 100 nHz. This catalog also comprises 489 red giants for which no significant rotational splitting was detected, which can be due to a low inclination angle, but also possibly to a very slow core rotation. One possible explanation for a slow core rotation in our targets would be that the strong internal field produces a very efficient transport of angular momentum that enforces solid-body rotation in the star. We note that near solid-body rotation was found by Li et al. (2024) for six red giants with slow-spinning cores. In Li et al. (2024), envelope rotation rates significantly above zero were measured for only 12% of the catalog (243 stars out of 2006), the other stars having envelope rotation rates compatible with zero. The slowest measured envelope rotation rate is 16 nHz, barely above the frequency resolution. This means that 88% of the stars in the catalog by Li et al. (2024) have envelope rotation below 16 nHz. If a mechanism enforces solid-body rotation in these stars, the core could rotate at rates that are undetectable even with 4 years of Kepler data.

Another possible explanation would be that through a yet unknown mechanism, when the magnetic field becomes near-critical (∼B_c), mode suppression could affect only the m = +1 component or only the m = −1 component. In the following sections, we work under the hypothesis that the two detected ridges correspond to m = 0 and m = 1 components, which is the only way we found of reproducing the observations.

4. Non-perturbative approach

The magnetic fields causing the strong deviation from classical g-mode pattern in our stars are near-critical, so that the effects of magnetic fields on the oscillations need to be addressed with a non-perturbative approach (Rui et al. 2024). Before presenting the details of our non-perturbative approach, adapted from Rui & Fuller (2023) (Sects. 4.2–4.6), we briefly recall a few aspects of the perturbative treatment of magnetic effects on oscillation mode frequencies in Sect. 4.1. This is useful when comparing with the results of the non-perturbative treatment in the next sections.

4.1. Perturbative approach

A perturbative approach was used early on to estimate the effects of a magnetic field on solar p modes (Gough 1990), on g modes in SPB stars (Hasan et al. 2005), and later on mixed modes in red giants, in the special case of a dipolar field with various radial profiles (Gomes & Lopes 2020; Mathis et al. 2021; Bugnet et al. 2021; Loi 2021). Li et al. (2022) generalized the expression of magnetic frequency shifts to magnetic fields with an arbitrary configuration. This study only assumes that the azimuthal field B_φ does not completely dominate over the radial field B_r (it is valid provided B_φ/B_r ≪ N/ω, where N is the Brunt-Väisälä frequency and ω is the mode angular frequency, with N/ω having typical values on the order of 10² in red giants).

Li et al. (2022) showed that the magnetic frequency shifts depend on only two parameters. The first parameter is the average frequency shift ω_B of all the components within the multiplet, which is given by the expression

$$\begin{array}{c}\hfill {\omega }_B=\frac{\mathcal{I}}{{\mu }_0{\omega }^3}⟨B_{\mathrm{r}}^2⟩,\end{array}$$(1)

where ℐ is a term depending on the internal structure of the reference model. The term ⟨B_r²⟩ is an average of B_r², expressed as

$$\begin{array}{c}\hfill ⟨B_{\mathrm{r}}^2⟩=\frac{1}{4\pi }{\displaystyle {\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}}K(r){\displaystyle {\int }_0^{2\pi }{\int }_0^{\pi }}{B}_{\mathrm{r}}^{2}sin\theta \mathrm{d}\theta \mathrm{d}\phi \mathrm{d}r,\end{array}$$(2)

where r_i and r_o are the inner and outer turning points of the g-mode cavity. The weight function, K(r), probes the g-mode cavity and sharply peaks in the H-burning shell. We stress that the average frequency shift ω_B strongly increases with decreasing mode frequency (it varies as ω⁻³). The second parameter that is involved in magnetic frequency shifts is the “asymmetry” parameter a_asym, which is a dimensionless quantity that depends on a horizontal average of B_r² weighted by the second-order Legendre polynomial P₂(cos θ) and expressed as

$$\begin{array}{c}\hfill {\displaystyle a_{\mathrm{asym}}=\frac{{\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}K(r){\displaystyle \int \int B_{\mathrm{r}}^2}{P}_2(cos\theta )sin\theta \mathrm{d}\theta \mathrm{d}\phi \mathrm{d}r}{{\displaystyle {\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}}K(r){\displaystyle \int \int B_{\mathrm{r}}^2}sin\theta \mathrm{d}\theta \mathrm{d}\phi \mathrm{d}r}.}\end{array}$$(3)

Essentially, the parameter ω_B gives a measure of the field strength and a_asym provides information on the field topology. The average quantity a_asym alone can be produced by an infinity of field geometries (Mathis & Bugnet 2023) but its value gives broad constraints on the shape of the field. The parameter a_asym indeed ranges between −1/2 (for a field completely concentrated on the equator) and 1 (for a field concentrated at the poles). For a purely dipolar field, a_asym takes values between −0.2 and 0.4. Das et al. (2024) showed that having access to magnetic frequency shifts for l = 2 modes could help characterize the field geometry.

If the effects of nonaxisymmetry of the field are negligible (that is, if B_r² is axisymmetric or if the magnetic splitting ω_B is smaller than the rotational splitting, see Li et al. 2022), then the magnetic field produces an additional frequency shift to the m = 0 and m = ±1 components, expressed as

$$\begin{array}{cc}& \delta {\omega }_g(m=0)=(1-a_{\mathrm{asym}}){\omega }_B\hfill \end{array}$$(4)

$$\begin{array}{cc}& \delta {\omega }_g(m=\pm 1)=(1+\frac{a_{\mathrm{asym}}}{2}){\omega }_B.\hfill \end{array}$$(5)

In previous works (Li et al. 2022; Deheuvels et al. 2023), ω_B was expressed as ω_B = δω₀(ω_max/ω)³, where δω₀ is the average magnetic frequency shift at the frequency of maximum power of the oscillations ω_max.

We attempted to fit expressions given by Eqs. (4) and (5) to the observed mode frequencies of the stars in our sample following the method that was used in Deheuvels et al. (2023), assuming that we detect the m = 0 and m = 1 components with negligible core rotation. For the four stars that show the most strongly curved ridges (KIC 3326392, KIC 5700274, KIC 9227589, and KIC 9697268), no satisfactory solution could be found. We managed to fit either one ridge or the other but not both of them simultaneously. For the other four stars, with less strongly distorted ridges, reasonable fits could be obtained using the perturbative approach. However, the magnetic frequency shift that was measured indicates magnetic fields that are close to the critical field strength, so that a non-perturbative approach is necessary (Rui et al. 2024).

4.2. Non-perturbative approach: Assumptions and analytical problem

We adapted the approach introduced by Rui & Fuller (2023). The authors include the Lorentz force in the oscillation equations using the Boussinesq approximation, which means that the effects of compressibility are taken into account only in the momentum equation. They assume that the oscillations are adiabatic and they make the Cowling approximation, which consists in neglecting the perturbation to the gravitational potential. We make the same assumptions here, also neglecting the effects of rotation in the oscillation equations (even if rotation were to be considered, it could be treated a posteriori as a perturbation because red giants are slow rotators). The magnetic field is assumed to be axisymmetric, which simplifies a lot the equations (the solutions vary as e^imφ, with m the azimuthal order). We assume that B_φ/B_r ≪ N/ω, so that horizontal field components can be neglected. Finally, the horizontal geometry of B_r(r, θ) is assumed to be the same in the whole g-mode cavity, so that it is expressed as B_r(r, θ) = B₀f(θ)g(r).

For gravity modes, the radial wavenumber is much larger than the horizontal wavenumber (k_r/k_h ∼ N/ω ∼ 10² in the cores of red giants), so that their displacements are nearly horizontal. Since we consider high-order g modes, the radial wavenumber is also much larger than the length scale of variations of structural quantities. For this reason, Rui & Fuller (2023) adopt the JWKB approximation in the radial direction (but not in the horizontal direction). In the oscillation equations, the eigenfunctions can then be written as the product of a plane-like wave function E(r) and a slowly varying envelope function. For instance, for the pressure perturbation, this yields

$$\begin{array}{c}\hfill p\prime (r,\theta )=\widehat{p}(r,\theta )E(r).\end{array}$$(6)

The phase function can be expressed as

$$\begin{array}{c}\hfill E(r)=exp(-i{\displaystyle \int k_{\mathrm{r}}}\mathrm{d}r).\end{array}$$(7)

In the oscillation equations, the radial derivative ∂/∂r can then be replaced by −ik_r, which is considered to be much larger than the radial gradients of envelope functions (${\partial }_{\mathrm{r}}\widehat{p}$, ${\partial }_{\mathrm{r}}\widehat{{\xi }_{\mathrm{r}}}$...).

Upon making the change of variable μ = cos θ and inserting the equation of motion into the continuity equation (see details in Appendix B), one finds that the pressure perturbation is solution of the eigenvalue problem

$$\begin{array}{c}\hfill \mathcal{L}\widehat{p}+\lambda \widehat{p}=0,\end{array}$$(8)

where

$$\begin{array}{cc}\hfill \mathcal{L}\widehat{p}& =\frac{\partial }{\partial \mu }\left[\frac{1-{\mu }^2}{1-b^2{f}^2}\frac{\partial \widehat{p}}{\partial \mu }\right]-\frac{m^2}{(1-{\mu }^2)(1-b^2{f}^2)}\widehat{p},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill \lambda & ={\left(\frac{r\omega k_{\mathrm{r}}}{N}\right)}^2,\hfill \end{array}$$(10)

and

$$\begin{array}{c}\hfill b=\frac{k_{\mathrm{r}}{B}_{0}g(r)}{\omega \sqrt{4\pi {\rho }_0}}.\end{array}$$(11)

In the general case, the coefficients of operator ℒ depend on both r and on μ. Thus the problem cannot be separated into a radial part and a horizontal part. However, with the JWKB approximation, the eigenfunctions of Eq. (8) depend only parametrically on r, as explained in Ogilvie & Lin (2004) (see also Mathis 2009), so that Eq. (8) can be solved in μ for each value of r. It is important to stress that the eigenvalue λ also depends parametrically on r. In the nonmagnetic case (b = 0), Eq. (8) reduces to the general Legendre equation. We then recover λ = l(l + 1) and the horizontal eigenfunctions correspond to the spherical harmonics Y_l^m(θ, φ).

The quantity b in Eq. (8) has been introduced by Rui & Fuller (2023) as b = k_rv_A, r/ω, where $\overrightarrow{v_{\mathrm{A}}}=\overrightarrow{B_0}/\sqrt{4\pi {\rho }_0}$ is the Alfvèn velocity. Because k_r ≫ k_h, b approximately corresponds to the ratio ω_A/ω, where ${\omega }_{\mathrm{A}}=\overrightarrow{k}·\overrightarrow{v_{\mathrm{A}}}$ is the Alfvèn frequency. The eigenvalue λ can also be written as λ = (b/a)², where

$$\begin{array}{c}\hfill a=\frac{N{v}_{\mathrm{A}}}{r{\omega }^2}.\end{array}$$(12)

The quantity a is proportional to the ratio ω_supp²/ω², where ω_supp is the frequency below which the magnetic field is expected to prevent the propagation of gravity waves (Rui & Fuller 2023). Although the details of how magnetogravity waves are suppressed by the magnetic field remain to be fully understood, it is thought that when they reach the critical frequency ω_supp, these waves refract into outgoing slow magnetic waves that approach infinite wavenumbers and dissipate (Lecoanet et al. 2017; Rui & Fuller 2023; David et al. 2026). Alternately, magnetogravity waves can be damped out by the phase-mixing process when they become resonant with Alfvèn waves (b > 1).

4.3. Profile of the magnetic field

In Rui & Fuller (2023) and Rui et al. (2024), the authors assumed that the field has a dipolar geometry (f(θ) = cos θ). For the radial dependence g(r), they adopted the Prendergast magnetic field geometry (Prendergast 1956). However, it should be noted that Rui et al. (2024) built their Prendergast field by imposing a vanishing surface field. The resulting field slowly decreases as a function of the radius, and it is nearly uniform over the g-mode cavity. Thus, the results presented by the authors approximately correspond to the case g(r) = constant. In this study, we explore different latitudinal and radial dependences for the field. In order to produce different levels of confinement of B_r in the core, we assumed a Gaussian-shape profile for g(r) = exp[−r²/(2σ_B²)]. The parameter σ_B measures the radial extension of the core field and is considered as a free parameter in our following analysis.

The latitudinal dependence f(θ) of B_r is a priori not known. In the perturbative treatment, it intervenes only in the calculation of the asymmetry parameter a_asym (Eq. (3)). In this study, we chose to build latitudinal profiles f(θ) that correspond to given values of a_asym. Even though the asymmetry parameter does not enter the calculation of magnetic frequency shifts in the non-perturbative approach, it is still interesting to relate f(θ) to a_asym because the effects of the magnetic fields on the oscillations are less strong in the high-frequency part of the spectrum, so we expect to approximately recover the solution of the perturbative approach. More importantly, we found that latitudinal profiles corresponding to the same value of a_asym produce similar results for all the purposes of our study (see Sect. 5.2 below). Naturally, an infinity of latitudinal profiles produce a given value of a_asym. We chose an f(θ) of the type

$$\begin{array}{c}\hfill f(\theta )=P_1(cos\theta )e^{\delta P_2(cos\theta )},\end{array}$$(13)

where δ ∈ ℝ. This is a large-scale field that is antisymmetric with respect to the equator. It has the advantage that it can be made to concentrate near the poles or near the equator depending on the value of δ. Indeed, it can be shown that a_asym → 1 for large positive values of δ and a_asym → −1/2 for large negative values of δ (despite the fact that the profile vanishes at the equator itself). When plugging the chosen expression of f(θ) in Eq. (3), we find that the function a_asym(δ) is continuous and strictly increasing over the interval ] − 1/2, 1[. Thus, for any target a_asym, we can find a unique value of δ so that the ratio in Eq. (3) is equal to a_asym (see Appendix C.1). We also imposed that the maximum value of |f| is 1. This way, we ensure that the mode becomes suppressed if either a of b exceeds unity somewhere in the star. The multiplication by P₁(cos θ) in the expression of f(θ) ensures that ∫₀^πf(θ)sin(θ) dθ vanishes. In these conditions, the chosen field can be made divergence-free, as expected from Maxwell’s equations, by adding an adequate θ-component B_θ(r, θ) (see Appendix C.3).

4.4. Numerical solution of the horizontal problem

For each radius r, Eq. (8) must be solved numerically. For this purpose, we use a relaxation method. Following Rui & Fuller (2023), we introduce variables 𝒫 = p′/ρ₀ω²r² and ${\mathcal{Z}}_{\theta }=\sqrt{1-{\mu }^2}{\xi }_{\theta }/r$, and the system of equations is then written as

$$\begin{array}{cc}\hfill \frac{\partial \mathcal{P}}{\partial \mu }& =-\frac{1-b^2{f}^2(\mu )}{1-{\mu }^2}{\mathcal{Z}}_{\theta }\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{Z}}_{\theta }}{\partial \mu }& =(\frac{b^2}{a^2}-\frac{m^2}{(1-{\mu }^2)[1-b^2{f}^2(\mu )]})\mathcal{P}.\hfill \end{array}$$(15)

We started from b = 0 (for which the eigenfunctions correspond to spherical harmonics), and we gradually increased b using the solution from the previous iteration as an initial guess for the next. For each value of b, this procedure provides the eigenfunctions (𝒫 and 𝒵_θ) and the parameter a. For a dipole field geometry, Rui & Fuller (2023) showed that a increases with b until the internal gravity wave branch connects with a slow magnetic wave branch, at which point a reaches a maximum value a_max. For larger values of b, a then decreases, which shows that the wave propagates back to larger radii inside the star and does not reach the deep core. For an axisymmetric dipolar field (a_asym = 0.4), Rui & Fuller (2023) found that in the case (l = 1, m = 0), the maximum value of a is reached when b → ∞ so that the wave is not refracted upward, but reaches infinite wavenumber and dissipates in the layer where a is maximal. In this study, we integrate the equations until either b reaches 1 or a(b) reaches a maximum.

In Fig. 5, we plot the numerical solutions to the horizontal problem obtained with the latitudinal profile given in Eq. (13) and various values of a_asym. We plot λ = (b/a)² as a function of a. For the case a_asym = 0.4 corresponding to an axisymmetric dipolar field, we recover the same a(b) relation as Rui & Fuller (2023). For negative asymmetry parameters, we find that λ increases with increasing a_asym for a given a, while the situation is reversed for a_asym > 0. The cause of mode suppression depends on the value of a_asym. The modes become suppressed because of resonance with Alfvèn waves (b exceeds unity) for a_asym ≲ −0.204 or a_asym ≳ 0.086 with m = 0 modes, and for a_asym ≳ 0.464 with m = ±1 modes. For other values of a_asym, mode suppression occurs because a reaches a_max.

Fig. 5. Numerical solution of the horizontal problem for different values of the asymmetry parameter aasym. We show λ as a function of a for (l = 1, m = 0) modes (top panel) and for m = ±1 modes (bottom panel). The colored dots indicate that mode suppression is reached, and the dashed black line corresponds to b = 1.

4.5. Numerical solution of the radial problem

The JWKB approximation in the radial direction yields the relation

$$\begin{array}{c}\hfill \omega =\frac{1}{\pi (n_g+{\epsilon }_g)}{\displaystyle {\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}}\frac{N}{r}\sqrt{\lambda }\mathrm{d}r,\end{array}$$(16)

where the expression of k_r was taken from Eq. (10). As mentioned above, λ depends parametrically on r, and its value needs to be determined for each layer of the g-mode cavity. For a given field geometry, Eq. (8) can be solved as explained in Sect. 4.4 to obtain a relation between a and b, and thus between a and λ = (b/a)². This means that λ(r) is known if we can determine a(r) for all radii within the g-mode cavity. The parameter a depends on equilibrium quantities – N(r), ρ(r) – and on the mode frequency ω (see Eq. (12)). Thus to calculate a, we needed a stellar model. The search for a stellar model representative of the studied star is addressed in Sect. 6.2, and here, we assume that we have access to such a model. The dependence of a on ω is an additional complication because the mode frequency is not known at this stage. We therefore proceeded iteratively, initializing ω to its nonmagnetic asymptotic frequency (which is obtained by setting λ = 2 in Eq. (16)) to obtain initial values for a(r). We interpolate the λ(a) relation to determine the value of λ(r) at each radius in the cavity. Equation (16) is then solved to obtain the asymptotic mode frequencies including the effects of the strong field. This process is iterated, the obtained mode frequencies being used to compute updated values of a(r), until convergence. If for some r in the cavity a(r) exceeds a_max or b(r) reaches unity, then we consider that the mode is suppressed and it is ignored.

4.6. Example

We here show an example of the oscillation spectrum obtained in the presence of a strong magnetic field. We consider a 1.1-M_⊙ stellar model evolved with the evolution code MESA (Paxton et al. 2011) until it reaches a radius of R_★ = 5.6 R_⊙. We assume a magnetic field B_r(r, θ) = B₀g(r)f(θ) as explained in Sect. 4.3, setting σ_B = 0.003 R_★ for the radial profile g(r). We choose an asymmetry parameter a_asym = 0.34 and build the corresponding function f(θ) (see Appendix C.1). In order to easily compare the results from the perturbative and non-perturbative approaches, it is convenient to relate the field intensity B₀ to the magnetic perturbative magnetic frequency shift δω₀. In this example, we considered a value of δω₀/(2π) = 2.2 μHz, which corresponds to ⟨B_r²⟩^0.5 = 200 kG. Using Eq. (1), B₀ is then obtained as

$$\begin{array}{c}\hfill B_0={\left[⟨B_{\mathrm{r}}^2⟩\frac{2}{{\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}K(r)g^2(r)\mathrm{d}r{\displaystyle {\int }_0^{\pi }}{f}^2(\theta )sin\theta \mathrm{d}\theta }\right]}^{0.5}.\end{array}$$(17)

We then calculated the frequencies of (l = 1, m = 0) and (l = 1, m = 1) modes by solving Eq. (16) (m = ±1 modes have the same frequencies because we ignored rotation here). The results are shown in the shape of a period échelle diagram in the left panel of Fig. 6 (red symbols).

Fig. 6. Period échelle diagrams of l = 1 g modes for the example case described in Sect. 4.6. Left: frequencies of m = 0 (circles) and m = 1 (triangles) modes computed for a field with a radial extension σB = 0.003 R★. Magnetic effects are treated either with a non-perturbative approach (red symbols), or a perturbative approach (gray symbols). Right: frequencies of m = 0 modes for values of σ = 0.0015 (red), 0.002 (green), and 0.003 (blue) and for a constant field (black).

For comparison, we also computed mode frequencies using a perturbative approach as explained in Sect. 4.1 with the same field configuration. The knowledge of δω₀/(2π) = 2.2 μHz and a_asym = 0.34 uniquely determines the magnetic shifts, given by Eqs. (4) and (5). The results are overplotted in Fig. 6 (gray symbols). As expected, for higher-frequency modes, the magnetic effects are weaker and the mode frequencies from the non-perturbative approach nearly overlap those obtained with the perturbative approach. At lower frequencies, the two sets of frequencies clearly separate and the perturbative approach is no longer valid. Also, evidence of suppression of m = 1 modes is seen at low frequencies.

5. Impact of the profile of B_r on mode frequencies

For weak fields, the perturbative approach holds, and the magnetic frequency shifts depend only on weighted averages of the radial field B_r in the radial and latitudinal directions (Eq. (2) and (3)). Two different fields producing the same weighted averages are thus expected to have the same impact on mode frequencies. In the non-perturbative approach, the solutions of Eq. (8) depend on the full radial profile g(r) and latitudinal profile f(θ) of B_r. This makes the situation more complex, but offers the potential opportunity to constrain the shape of the magnetic field.

5.1. Radial profile of B_r

Using the same example case as described in Sect. 4.6 (stellar model and magnetic field configuration), we computed magnetically distorted mode frequencies with different radial extensions of the core field. We considered values of σ_B of 0.0015 R_★, 0.002 R_★, and 0.003 R_★. For our stellar model, the H-burning shell is located at approximately 6 × 10⁻³ R_★ so that the amplitudes of these three fields in the H-burning shell corresponds to 0.04%, 1.3%, and 14.5% of the central amplitude, respectively. Finally, we considered a field that is constant over the whole cavity, similar to the field assumed by Rui & Fuller (2023). For each field, B₀ was chosen to produce the same value of ⟨B_r²⟩ (Eq. (17)). The left panel of Fig. 7 shows the considered fields. Since the weight function K(r) peaks in the H-burning shell, fields that are confined below this shell have much larger values of B₀. The period échelle diagrams of m = 0 modes for these different fields are shown in the right panel of Fig. 6. At higher frequencies, the modes approximately correspond to the perturbative solutions, as expected, but at lower frequencies, it is striking to see that the considered fields produce very different spectra.

Fig. 7. Impact of the radial profile of Br on the mode frequencies. We considered fields with radial extensions of σB = 0.0015 R★ (red lines), 0.002 R★ (green lines), and 0.003 R★ (blue lines) and a constant field (black lines). Left: radial profile of Br. Middle: profile of the parameter a(r) in the g-mode cavity for the mode closest to νmax (Eq. (12)). Right: profile of λ(r) deduced from the interpolation of the a(λ) relation for the same mode (see text). We show the departure from the value λ = 2 corresponding to the perturbative approach.

The departure from the perturbative solution (gray symbols) is maximal for the constant field. It decreases as σ_B decreases to reach a minimum for σ_B ≈ 2 × 10⁻³ R_★ and then increases again as σ_B further decreases. This behavior can be understood from the profile of a(r) and λ(r), which are shown in the middle and right panels of Fig. 7 for the mode that is closest to ν_max. The constant field has a maximal amplitude in the H-burning shell where the sensitivity to the field is the highest. It thus has large values of a (and λ values well above two) in this shell, which produces large non-perturbative effects, and causes mode suppression at low-frequency. When the field becomes more confined in the core, the relative amplitude of the field in the H-burning shell decreases, so that non-perturbative effects weaken. For σ_B = 2 × 10⁻³ R_★, the λ value is very close to two in the H-burning shell. Consequently, the ridges are less distorted than in the constant-field case, and mode suppression occurs at much lower frequency (see right panel of Fig. 6). When the field becomes even more confined in the core (σ_B = 1.5 × 10⁻³ R_★), its strength B₀ increases to maintain the average ⟨B_r²⟩. As a result, the values of a become large at the center, which increases non-perturbative effects and can cause mode suppression. This example demonstrates that measuring the mode frequencies of near-critical fields could give constraints on the radial extent of the core magnetic field, contrary to weaker fields.

5.2. Latitudinal profile of B_r

In this study, we found that latitudinal profiles f(θ) producing the same asymmetry parameter a_asym lead to very similar magnetic frequency shifts. To reach this conclusion, we considered several different profiles f(θ). Beside the profile given in Eq. (13), we explored profiles of the type f(θ) = P₁(cos θ)+BP₃(cos θ)+CP₅(cos θ). In this expression, B and C can be determined to produce a given value of a_asym, although solutions can be found for only a restricted interval (the value C ≈ 2.37 maximizes the reachable range of a_asym ∈ [ − 0.252, 0.832]). The value of a_asym can also be ensured by profiles of the type $f(\theta )=\sqrt{1+5a_{\mathrm{asym}}{P}_2(cos\theta )}$, yielding solutions only over the interval a_asym ∈ [ − 0.2, 0.4]. The analyses presented in the following sections of this paper were repeated with the latitudinal profiles presented above, and we found that our results are statistically unchanged. This suggests that non-perturbative effects depend on the latitudinal shape of B_r only through the weighted average given by a_asym, as for the perturbative approach. In the rest of this study, we used the latitudinal profile given by Eq. (13).

6. Application to red giants with near-critical magnetic fields

We then used the non-perturbative approach described in Sect. 4 to search for strong field configurations that could reproduce the observations for the eight stars of the studied sample. For this purpose, we built an asymptotic expression of mixed modes including the effects of magnetic fields following a non-perturbative approach, and we explored the parameter space.

6.1. Asymptotic model

To derive asymptotic expressions of magnetic mixed modes, we followed the approach that was first described in Li et al. (2022) and then used in Li et al. (2023) and Deheuvels et al. (2023). From the work of Shibahashi (1979) and Unno et al. (1989), asymptotic expressions of mixed modes of azimuthal order m are obtained by solving the equation

$$\begin{array}{c}\hfill tan{\theta }_{p,m}=qtan{\theta }_{g,m},\end{array}$$(18)

where q corresponds to the coupling strength between the p- and g-mode cavities, and θ_p, m, θ_g, m are phase terms that can be expressed as a function of asymptotic frequencies of p and g modes. Following Mosser et al. (2012c), we have

$$\begin{array}{cc}\hfill {\theta }_{p,m}& =\frac{\pi }{\mathrm{\Delta }\nu (n_p)}(\nu -{\nu }_{p,m})\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill {\theta }_{g,m}& =\frac{\pi }{\mathrm{\Delta }P(n_g)}(\frac{1}{\nu }-P_{g,m})-\frac{\pi }{2},\hfill \end{array}$$(20)

where Δν(n_p) and ΔP(n_g) are the local large separation of p modes and period spacing of g modes, and ν_p, m and P_g, m are the frequencies of pure p and periods of pure g modes, respectively¹.

The effects of the magnetic field on pressure modes can be neglected (Mathis et al. 2021), and since we considered the case of a negligible rotation, the asymptotic frequencies, ν_p, m, can be written as

$$\begin{array}{c}\hfill {\nu }_{p,m}=[n_{\mathrm{p}}+\frac{1}{2}+{\epsilon }_{\mathrm{p}}+\frac{\alpha }{2}{(n-n_{\max })}^2]\mathrm{\Delta }\nu -d_{01},\end{array}$$(21)

where α measures the second-order effects in the asymptotic development around ν_max (as expressed by Mosser et al. 2012c), ε_p is a phase offset, n_max = ν_max/Δν, and d₀₁ corresponds to the small separation δν₀₁ = [ν_p(0, n_p)+ν_p(0, n_p + 1)−2ν_p(1, n_p)]/2 when mode frequencies are expressed at first order.

For g modes, we include the effects of the magnetic field in a non-perturbative way using Eq. (16). For this purpose, a stellar model is needed. This is an additional difficulty compared to the perturbative approach, for which the magnetic shifts only depend on the average quantity ω_B (Eqs. (4) and (5)). To explore the parameter space, we computed one reference stellar model for each star (in Sect. 6.2, we explain the choice of a reference model and we discuss the validity of this approach). For this model, we can compute the asymptotic mode frequencies taking into account magnetic fields in a non-perturbative way by solving Eq. (16) for each mode. We can also compute nonmagnetic asymptotic frequencies ${\omega }_g^{(0)}$ by setting λ = 2 in Eq. (16). We then obtain fractional magnetic frequency shifts as

$$\begin{array}{c}\hfill \frac{\delta {\omega }_{\mathrm{mag},m}}{{\omega }_g^{(0)}}=\frac{{\displaystyle {\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}}\frac{N}{r}\sqrt{{\lambda }_m(r)}\mathrm{d}r}{{\displaystyle {\int }_{r_{\mathrm{i}}}^{r_{\mathrm{o}}}}\frac{N}{r}\sqrt{2}\mathrm{d}r}-1.\end{array}$$(22)

When exploring the parameter space, we assume that the fractional frequency shift of a given mode is that of the reference model given by the expression above. To obtain asymptotic frequencies of pure g modes, we add the magnetic frequency shift given by Eq. (22), leading to the expression

$$\begin{array}{c}\hfill P_{g,m}=P_g^{(0)}{(1+\frac{\delta {\omega }_{\mathrm{mag},m}}{{\omega }_g^{(0)}})}^{-1},\end{array}$$(23)

where $P_g^{(0)}=(n_g+1/2+{\epsilon }_g)\mathrm{\Delta }{\mathrm{\Pi }}_1$.

Overall, the free parameters of the asymptotic model are the properties of pure p modes (Δν, ε_p, α, d₀₁), those of pure g modes (ΔΠ₁, ε_g), the coupling strength q, and the properties of the magnetic field given by the parameters (δω₀, σ_B, a_asym). To characterize the field, we proceed as in Sect. 4.6. The free parameters δω₀ and σ_B are used to build the Gaussian-shaped radial profile g(r) of B_r (B₀ is obtained from δω₀ using Eqs. (1) and (17)), and the parameter a_asym is used to generate the latitudinal profile f(θ) from Eq. (13). For each set of parameters, the asymptotic frequencies of dipole mixed modes for m = 0 and m = 1 are calculated by solving Eq. (18). Finally, we must remind that asymptotic models are only approximate solutions of the oscillation equations. Departures from asymptotics are taken into account by including a normal model error with widths σ_p and σ_g for p and g modes, respectively. The quantities σ_p and σ_g are considered as free parameters when searching for optimal solutions.

6.2. Search for reference stellar models

To search for reference models for the eight stars, we used a grid of stellar models computed with MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), with masses varying between 0.9 and 2 M_⊙ (step 0.1 M_⊙), metallicities between −0.4 and 0.4 dex (step 0.1 dex), and a step overshoot was added beyond the convective core with α_ov = 0, 0.1, or 0.2. We used the OPAL equation of state (Rogers & Nayfonov 2002) and opacity tables. We used nuclear reaction rates from NACRE (Angulo et al. 1999). Microscopic diffusion was included by solving the equations of Burgers (1969) at each time step. The models were evolved until they reach 15 R_⊙, keeping track of the internal structure for all models computed along the evolution. For each star, we searched for the model of the grid that best reproduces the asymptotic large separation of p modes Δν, the frequency of maximum power of the oscillations ν_max, and the asymptotic period spacing of g modes ΔΠ₁ by minimizing the function

$$\begin{array}{c}\hfill {\chi }^2=\frac{{(\mathrm{\Delta }\nu -\mathrm{\Delta }{\nu }^{\mathrm{obs}})}^2}{{\sigma }_{\mathrm{\Delta }\nu }^2}+\frac{{({\nu }_{\max }-{\nu }_{\max }^{\mathrm{obs}})}^2}{{\sigma }_{{\nu }_{\max }}^2}+\frac{{(\mathrm{\Delta }{\mathrm{\Pi }}_1-\mathrm{\Delta }{\mathrm{\Pi }}_1^{\mathrm{obs}})}^2}{{\sigma }_{\mathrm{\Delta }{\mathrm{\Pi }}_1}^2}.\end{array}$$(24)

We note that ΔΠ₁ is in fact not known at this stage because the measured period spacing ΔΠ₁^meas is an underestimated measurement of ΔΠ₁ (see Sect. 2.3). To obtain a first estimate of ΔΠ₁, we used the results of the perturbative approach obtained in Sect. 4.1. As mentioned above, for the stars with the strongest ridge distortion, this approach cannot fit both ridges simultaneously. In these cases, we used the results of the perturbative approach applied to the ridge showing the weakest distortion alone. We checked a posteriori that it produces estimates of ΔΠ₁ that are close to the more proper measurements obtained with the non-perturbative approach in Sect. 6.3.

To check that our conclusions do not critically depend on the choice of a reference model, the fits presented in Sect. 6.3 were performed again using the next best models in the grid and we found that the results were statistically indistinguishable. Similarly, for one of the stars (KIC 5700274), the reference model was updated using the asymptotic period spacing ΔΠ₁ found in Sect. 6.3 and a new fit was performed using this model. Again, the posterior distribution of the parameters was found to be nearly identical.

6.3. Asymptotic fits

We fit our model to the observed frequencies by following a Bayesian approach. We used the Asteroseismic Bayesian Inference by MCMC (ABIM) code (see Villate et al. 2026, for more details). ABIM is a Fortran code parallelized with OpenMP directives. Sampling is performed with a Markov Chain Monte Carlo method implementing the so-called stretched move algorithm proposed by Goodman & Weare (2010) and including parallel tempering (e.g., Benomar et al. 2009), which improved the exploration of parameter spaces containing numerous local maxima, as is the case here. We typically sample the posterior distributions with 15 parallel chains and 256 walkers for each chain. The initial positions of the walkers are randomly drawn from the prior distributions. The walkers are iterated over 12 000 steps, with the 2000 first steps discarded as burn-in to ensure that chains are stationary. The priors for each parameter are provided in Appendix D. To calculate the likelihood, we compare each observed mode to the closest mode in the asymptotic model considered.

For each step of the MCMC, we need to compute the relation λ(a) corresponding to the considered value of a_asym, which requires solving Eq. (8) using a relaxation scheme, as described in Sect. 4.4. In order to reduce the computing time, we precomputed tables of λ(a) relations for all possible asymmetry parameters with a step of 10⁻³. Then, during the MCMC procedure, the value of λ is obtained by performing a bilinear interpolation for the considered values of a and a_asym.

6.4. Results

For all eight stars, we were able to reproduce very well the observed mode frequencies with our asymptotic model. This can be seen in Fig. 8, which shows stretched échelle diagrams all eight stars We recall that no satisfactory fits had been obtained using a perturbative approach, or identifying the modes as m = ±1 components. We note that the non-perturbative approach as we have applied it in this study involves one more degree of freedom compared to the perturbative approach (the radial extent of the field σ_B). This increase in the degrees of freedom is fully justified because the perturbative approach (i) fails to fit the most distorted stars, and (ii) leads to field measurements that are close to the critical field strength, which makes the perturbative approach not valid (see Rui et al. 2024).

Fig. 8. Stretched échelle diagrams of three stars of our sample. Blue circles and horizontal blue lines correspond to observed modes and their 1σ uncertainties. Red empty circles correspond to the best asymptotic models including magnetic effects in a non-perturbative manner.

The optimal parameters of the asymptotic models are given in Table 2 and we show the posterior probability of a selected sample of parameters in Fig. 9 (see also Figs. D.1 and D.2). In particular, we obtained measurements of the asymptotic period spacing ΔΠ₁ that are unbiased by the magnetic field, contrary to the apparent period spacing ΔΠ₁^meas. We thus revised the location of the eight red giants in the (Δν, ΔΠ₁) plane (red star-shaped symbols in Fig. 3). As can be seen in this figure, the stars now clearly belong to the degeneracy sequence, as expected. This result is quite encouraging since we did not include any prior stating that ΔΠ₁ should be such that the stars lie on the degeneracy sequence.

Fig. 9. Corner plots of KIC 9227589 (left) and KIC 11408970 (right), restricted to the asymptotic period ΔΠ1 of dipolar g modes, and the parameters characterizing the core magnetic field: the magnetic shift δν0 = δω0/(2π), σB, and aasym. The corner plots for other stars are given in Appendix D.

6.4.1. Field strength

The measured values of δω₀ could be used to obtain estimates of the measured field strength. We expressed them as the weighted average ⟨B_r²⟩^0.5, for consistency with previous core field measurements. To relate δω₀ to ⟨B_r²⟩^0.5, we used Eq. (1) and the reference stellar models obtained in Sect. 6.2. We accounted for the uncertainties in the measurements of δω₀ (Table 2) and uncertainties related to the choice of a reference model were estimated by considering all stellar models from the grid yielding χ² < χ_min² + 9, where χ_min² is the value of χ² (see Eq. (24)) for the reference model. We obtained field strengths ranging from about 100–700 kG. These measurements are in line with those found by Deheuvels et al. (2023) for similar stars, except that in this previous study, only lower limits to the field strength were obtained (because only one component per triplet was detected), and a perturbative approach was used even though the fields are near-critical. In the last column of Table 2, we give the value of the field strength at the center B₀ for the median model of each star. We find that B₀ is always much larger than ⟨B_r²⟩^0.5. This is because the optimal magnetic fields are all confined well below the H-burning shell (see Sect. 6.4.3). For the four stars that yielded reasonable fits with the perturbative approach (KIC 7458743, KIC 10613561, KIC 10976412, and KIC 11408970), we found that the perturbative approach overestimates the field strength by 5–10%.

In the results presented here, we assumed that the core rotation is too weak to separate m = ±1 components (see Sect. 3.2). Our alternative explanation for the detection of two components only was that for some as-yet-unknown reason, only one of the m = ±1 components is visible. To check the impact of this alternative scenario on our results, we repeated our analysis of one of our stars (KIC 3326392) assuming a core rotation of ⟨Ω⟩_g = 0.5 μHz, which is typical for red giants (Hatt et al. 2024; Li et al. 2024). We found that the results are nearly unchanged, the measured field strength being modified by only 10 kG, which is within our quoted uncertainties.

6.4.2. Mode identification

In this study, we assumed that the two detected ridges correspond to m = 0 and m = ±1 ridges, but we made no assumption about which ridge corresponds to which m-component. With the perturbative approach, the nonrotating case holds an ambiguity. The oscillation spectrum corresponding to δω₀ and a_asym is identical to the spectrum obtained with δω′₀ = δω₀(2 − a_asym)/2 and a′_asym = 2a_asym/(a_asym − 2), the other parameters being unchanged. This degeneracy exists only for a_asym < 0.4, to ensure that a′_asym ∈ [ − 1/2, 1]. For strong fields, using the non-perturbative approach, this degeneracy is lifted. For all eight stars, the preferred solutions with positive asymmetries are clearly preferred, which corresponds to identifying the most distorted ridge of the m = ±1 component. For some stars, the negative asymmetry solutions appear as a small bump in the posterior distribution of a_asym (e.g., for KIC 11408970), but it remains negligible compared to the positive asymmetry solution. This is in line with the observation made for weaker magnetic fields that stars with positive a_asym are more frequent that stars with a_asym < 0.

6.4.3. Radial extent of the core magnetic field

As demonstrated in Sect. 5.1, we have the possibility to constrain the radial extent of the core magnetic field when the field is near-critical. For all eight stars, we obtained clear contraints on the width σ_B of the Gaussian function that was used to model the radial profile of B_r (see Table 2). The posterior distribution of σ_B can be seen in the corner plots shown in Figs. 9, D.1, and D.2. For four stars (KIC 5700274, KIC 9227589, KIC 9697268, and KIC 10976412), the posterior distribution of σ_B peaks on a single value. For KIC 3326392, KIC 10613561, KIC 11408970, and to a lesser extent KIC 7458743, we find a bimodality of the solutions. For instance, for KIC 3326392, the observations can be reproduced either with σ_B = 0.0018 R_★ or with σ_B = 0.0029 R_★. This phenomenon is linked to the behavior described in Sect. 5.1. For a given value of ⟨B_r²⟩, non-perturbative effects can be large either when the field is strong near the HBS, or when the field is confined in the central layers. These two possibilities give quite similar oscillation spectra, so that both solutions are possible. In practice, the two solutions lead to similar values of ⟨B_r²⟩ and a_asym.

For all eight targets, we found that the radial extent of the core field σ_B is significantly smaller than the radius of the HBS r_HBS (the ratio σ_B/r_HBS ranges from 30% to 50% for all eight stars). Consequently, the field strength in the HBS is less than 10% of its maximal intensity for all stars. If we consider magnetic fields that extend closer to the HBS, non-perturbative effects strongly increase, the ridges become too distorted, and mode suppression occurs at too-high frequency compared to the observations, so that no satisfactory solution can be found.

6.4.4. Mode suppression

For KIC 3326392, KIC 9227589, and KIC 9697268, seismic observations have shown that one ridge undergoes mode suppression at low frequency. Even though we did not include as a constraint in our fit that one ridge is suppressed, our best models indeed show mode suppression on the same ridge (the m = ±1 ridge according to Sect. 6.4.2) for all three targets. Moreover, we find that the suppression of the ridge in our optimal models occurs at a frequency that is very close to the lowest-frequency mode detected in the ridge. This can be seen in the stretched échelle diagrams shown in Fig. 8. For KIC 9227589 and KIC 9697268, the lowest-frequency detected mode coincides exactly with the last non-suppressed m = 1 mode in our optimal model. For KIC 3326392, mode suppression occurs at a frequency of 114.83 μHz, which is only slightly higher than the lowest-frequency observed modes (114.41 μHz). We find this agreement quite encouraging, although more thorough work will be needed to determine to what extent the frequency below which modes are suppressed can be used to infer magnetic field properties. We also checked the reason for mode suppression in the three stars. For KIC 9227589 and KIC 9697268, the solution for σ_B is unimodal and mode suppression occurs because a(r) exceeds the maximal value of a_max in the core. For KIC 3326392, the solutions are bimodal, and as could be expected, we find that a(r) exceeds a_max in the HBS for the solution with higher σ_B, and near the core for the solution with lower σ_B.

7. Discussion

7.1. Origin of the detected fields

One of the main scenarios is that the magnetic fields detected in red giants originate from dynamo action in main-sequence convective cores. Numerical simulations of convective cores in A-type stars show the presence of such fields, with intensities of a few tens of kiloGauss (Brun et al. 2005). After the end of the main sequence, these fields could relax into stable configurations (Becerra et al. 2022), although they might be reshaped by the core contraction (Gouhier et al. 2022). The ohmic diffusion timescale being longer than the star’s lifetime (Cantiello et al. 2016), we expect these fields to survive in the cores of red giants and to remain mostly confined in previously convective layers. As stars evolve in the post-main sequence phase, the HBS progresses outward in mass and eventually outreaches the maximal extent of the convective core during the MS. In this case, the core magnetic field does not extend to the HBS. This reduces the seismically measured field strength ⟨B_r²⟩, which has maximal sensitivity near the HBS. This also reduces non-perturbative effects for near-critical fields, as shown in Sect. 5.1, and it corresponds to the configuration that gives best agreement with the seismic observations for our eight targets (Sect. 6.4). The moment when the HBS leaves the previously convective core strongly depends on the stellar mass: it occurs sooner for lower-mass stars (M ∼ 1.2 M_⊙), which have small MS convective cores.

Interestingly, we observe that red giants with magnetically distorted doublets tend to have lower masses than other magnetic red giants. Indeed, seven stars (out of eight) have seismic masses that are between 1.12 and 1.24 M_⊙, in agreement with the masses of our best-fit models (lying between 1.1 and 1.2 M_⊙), while the average seismic mass of magnetic giants detected so far is 1.31 ± 0.15 M_⊙ (Li et al. 2022, 2023; Deheuvels et al. 2023; Hatt et al. 2024). The last star, KIC 9697268, has an estimated mass of 1.43 ± 0.09 M_⊙ from seismic scaling relations (Yu et al. 2018), but the period spacing ΔΠ₁ we measure for this star (see Table 2) makes it more compatible with lower-mass models (our best-fit model has a mass of 1.2 M_⊙). For all eight stars, we used our best-fit models to estimate the radius r_cc of the layer corresponding to the maximal extent of the convective core during the MS. We found that r_cc is below the HBS for all stars (r_cc/r_HBS is between 40% and 60%). This observation holds regardless of the amount of core overshooting that is added beyond the MS core (r_cc varies by only ∼10% when α_ov is increased from 0 to 0.2). In fact, we find that the radius r_cc is remarkably similar to the radial extent of the magnetic field σ_B that we measured in Sect. 6.4 (the ratio r_cc/σ_B is between 1 and 1.3 for all eight stars). The magnetic fields that we measure for the stars of our sample are thus consistent with the hypothesis that they originate from convective core dynamo.

We note that if our stars were more massive M ≳ 1.4 M_⊙, we would be in the inverse situation, the HBS would still be within the layers that were previously convective during the MS. In this case, we have shown that non-perturbative effects are much more severe, and mode suppression occurs at higher frequency. It is therefore tempting to suggest that we are detecting strongly distorted doublets in the stars of our sample precisely because they are low-mass stars, and that more massive stars at the same evolutionary state and with similar field intensity might have lost one or both components to mode suppression. A more thorough search for red giants with near-critical fields could help test this hypothesis.

7.2. Assumptions

We now discuss some of the assumptions that were made in this study.

7.2.1. Radial JWKB approximation in non-perturbative approach

Our non-perturbative approach relies on a number of assumptions, which were already made by Rui & Fuller (2023) and Rui et al. (2024) in their development. The use of a radial JWKB approximation in the unseparated problem would deserve to be validated, ideally using a complete non-perturbative treatment that includes the effects of the magnetic field. This is beyond the scope of this work, but we note that developments pursuing this goal are currently under way (Fort et al. 2024). One of the underlying assumptions is that the radial wavenumber k_r introduced in Eq. (7) depends only on r. We investigated this assumption using the ray tracing approach proposed by Loi (2020). The latter study integrates Hamilton’s equations in spherical coordinates using a stellar model and assuming a particular magnetic configuration (a Prendergast field). This procedure gives access to the profile of k_r for magnetic fields with arbitrary intensities. We reproduced the approach of Loi (2020) and found that k_r remains approximately a function of r only, even for field intensities that approach B_c. This would need to be shown for arbitrary magnetic configurations.

7.2.2. Profile of the magnetic field

In this study, we have assumed that the magnetic field is axisymmetric, mainly because only two components per dipole multiplet are visible. We note that if the core rotation of our stars is indeed very slow, as proposed in Sect. 3.2, then the ratio between the magnetic splitting and the rotational splitting is large. In this case, the alignment of the pulsations is expected to be dictated by the magnetic field rather than the rotation. If the magnetic field possesses an axis of symmetry, then we would indeed be in an axisymmetric configuration along the magnetic axis. So far, non-perturbative treatments of magnetic effects on oscillation mode frequencies were developed under the assumption of an axisymmetric field (Rui & Fuller 2023). When complete non-perturbative calculations are available, it will be interesting to investigate the non-axisymmetric case as well. In this study, we also separated the radial field component into B_r(r, θ) = B₀f(θ)g(r). This assumption seems reasonable because the oscillation modes are sensitive to the magnetic field only in the g-mode cavity, which has a very small extension in the star (typically ∼1% of the stellar radius) and most of their sensitivity is in the very thin H-burning shell.

8. Conclusion

In this paper, we detected strong deviations from the regular period spacing pattern of g modes in eight Kepler red giants exhibiting l = 1 doublets. The measured period spacings vary with frequency, and they appear abnormally low compared to what is expected for red giant branch stars with strong electron degeneracy in the core (Deheuvels et al. 2022). We ruled out buoyancy glitches as the cause of these features and investigated the magnetic origin. To produce such distortions in the g-mode pattern, the core magnetic field needs to be close to the critical strength, B_c, above which the propagation of magneto-gravity waves is inhibited (Fuller et al. 2015).

We took into account the effects of magnetic fields on oscillation modes using the non-perturbative approach developed by Rui & Fuller (2023) and Rui et al. (2024). We showed that the radial dependence of the radial component of the field B_r has a large impact on the intensity of the magnetic frequency shift. This is in contrast with results from the perturbative approach, whereby magnetic frequency shifts only depend on a weighted average of B_r² (Li et al. 2022). This demonstrates that information on the radial profile of B_r can be obtained for red giants that host near-critical core magnetic fields. On the contrary, the latitudinal dependence of B_r seems to impact magnetic frequency shifts only through the weighted average a_asym (Eq. (3)), as in the case of the perturbative approach.

One critical question in our work regards the identification of the two components detected per l = 1 multiplet. If these components are identified as m = ±1, as would be most natural, no satisfactory fit to the data can be found for any of the eight stars. Also, this identification fails to explain the difference of amplitude between the two components in KIC 7458743 or the suppression of one of the two ridges in KIC 3326392, KIC 9227589, and KIC 9697268. As an alternative, we explored the identification of the two components as m = 0 and m = 1. This assumption enabled us to satisfactorily account for all observational features. First, using a non-perturbative approach, we obtained very good fits to the two distorted ridges in all eight stars. Secondly, the period spacings, ΔΠ₁, obtained under this assumption brought the eight stars back to the degeneracy sequence where they belong in the (Δν, ΔΠ₁) plane (see Fig. 3). Finally, our best-fit solutions for KIC 3326392, KIC 9227589, and KIC 9697268 have a partially suppressed m = 1 ridge, in agreement with the observations, and the frequencies of mode suppression are very close to the observed ones. Despite these successes, it should be noted that we currently cannot satisfactorily explain why only two ridges (and not three) are observed. In this paper, we worked under the assumption that the core rotation is too slow to separate the m = ±1 ridges. So far, magnetic red giants have been found to have a similar core rotation compared to regular red giants (Li et al. 2022; Hatt et al. 2024; Villate et al. 2026). However, the stars studied here have stronger fields than the stars from these previous works, so it is conceivable that in the strong-field regime, the transport of angular momentum could be more efficient. Finally, we note that a toroidal field could act differently on the m = ±1 components (Dhouib et al. 2022), so it would be interesting to investigate the impact of a strong near-critical toroidal field on the oscillations modes of our eight stars in the future.

To reproduce the observations, we had to invoke core magnetic fields, B_r, that are confined well below the H-burning shell (Gaussian-shaped profiles for B_r with widths σ_B corresponding to 30–50% of the radius of the HBS). Fields extending all the way to the HBS induce stronger non-perturbative effects so that g-mode period distortion is too large, and mode suppression occurs at too a high frequency to reproduce the observations. If our interpretation is correct, this is the first time we can place constraints on the radial extent of the core magnetic field.

The eight stars under study have a low mass compared to other red giants. Our best-fit stellar models have masses between 1.1 and 1.2 M_⊙. They thus had a small convective core during the main sequence. Interestingly, the maximal size of the main-sequence convective core in our models approximately corresponds to our measurements of the radial extent, σ_B, of the core field. This means that the detected fields could be the remnant of fields generated by a dynamo effect in the main-sequence convective core. More massive stars have larger convective cores so that the HBS is still within the previously convective layers at the bottom of the red giant branch. As we have shown in this paper, non-perturbative effects are stronger in this case and lead to more mode suppression. This might partially explain why the mass distribution of red giants with direct detection of core magnetic fields peaks at a lower mass than the distribution of red giants with suppressed dipolar modes (see Villate et al. 2026). This calls for further investigation.

Acknowledgments

The authors thank the anonymous referee for comments and suggestions that improved the clarity of the manuscript. S.D. and J.B. acknowledge support from from the project BEAMING ANR-18-CE31-0001 of the French National Research Agency (ANR). This work has been supported by CNES, focused on the preparation of the PLATO mission. This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555.

References

Appendix A: Mode frequencies

We give the measured mode frequencies for all eight stars of the sample as obtained from our maximum likelihood estimation technique (see Sect. 2.2) in Tables A.1 and A.2.

Appendix B: Oscillation equations in non-perturbative framework

We here detail the main steps that lead to the expressions for g mode frequencies with magnetic field taken into account in a non-perturbative manner. The development essentially follows that of Rui & Fuller (2023). With the assumptions made in Sect. 4.2, the oscillation equations can be written as follows. The continuity equation reads:

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\partial }{\partial r}(r^2{\xi }_{\mathrm{r}})+\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta {\xi }_{\theta })+\frac{1}{rsin\theta }\frac{\partial {\xi }_{\phi }}{\partial \phi }=0\end{array}$$(B.1)

The momentum equations are written as:

$$\begin{array}{cc}& -{\rho }_0{\omega }^2{\xi }_{\mathrm{r}}=-\frac{\partial p^{\prime }}{\partial r}-{\rho }^{\prime }{g}_0\hfill \end{array}$$(B.2)

$$\begin{array}{cc}& -{\rho }_0{\omega }^2{\xi }_{\theta }=-\frac{1}{r}\frac{\partial p^{\prime }}{\partial \theta }+\frac{1}{4\pi }{B}_{\mathrm{r}}^2\frac{{\partial }^2{\xi }_{\theta }}{\partial r^2}\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& -{\rho }_0{\omega }^2{\xi }_{\phi }=-\frac{1}{rsin\theta }\frac{\partial p^{\prime }}{\partial \phi }+\frac{1}{4\pi }{B}_{\mathrm{r}}^2\frac{{\partial }^2{\xi }_{\phi }}{\partial r^2}\hfill \end{array}$$(B.4)

We assume the oscillations are adiabatic, so the energy equation can be replaced by the adiabaticity relation

$$\begin{array}{c}\hfill \frac{{\rho }^{\prime }}{{\rho }_0}=\frac{1}{{\mathrm{\Gamma }}_{1,0}}\frac{p^{\prime }}{p_0}+\frac{N^2{\xi }_{\mathrm{r}}}{g_0},\end{array}$$(B.5)

where Γ_1, 0 is the adiabatic index. For gravity waves, the first term in the right-hand-side of Eq. (B.5) can be neglected, so that ρ′≃ρ₀N²ξ_r/g₀.

All the eigenfunctions are written as the product of a plane-like wave function E(r) = exp(−i∫k_r dr) and an envelope function (for instance $\widehat{p}(r,\theta )$ for the pressure perturbation) that is a slowly varying function of the radius. The continuity equation then becomes

$$\begin{array}{c}\hfill -i{k}_{\mathrm{r}}\widehat{{\xi }_{\mathrm{r}}}-\frac{1}{r}\frac{\partial }{\partial \mu }\left(\sqrt{1-{\mu }^2}\widehat{{\xi }_{\theta }}\right)+\frac{1}{r\sqrt{1-{\mu }^2}}im\widehat{{\xi }_{\phi }}=0\end{array}$$(B.6)

Combining the radial equation of motion with the adiabaticity relation, we obtain

$$\begin{array}{c}\hfill -{\rho }_0{\omega }^2\widehat{{\xi }_{\mathrm{r}}}=i{k}_{\mathrm{r}}\widehat{p}-{\rho }_0{N}^2\widehat{{\xi }_{\mathrm{r}}}\end{array}$$(B.7)

Since ω ≪ N, the left-hand term of this equation can be neglected. The horizontal equations of motion yield

$$\begin{array}{cc}\hfill \widehat{{\xi }_{\theta }}& =-\frac{1}{{\rho }_{0}r{\omega }^2}\frac{\sqrt{1-{\mu }^2}}{1-b^2{f}^2}\frac{\partial \widehat{p}}{\partial \mu }\hfill \end{array}$$(B.8)

$$\begin{array}{cc}\hfill \widehat{{\xi }_{\phi }}& =\frac{1}{{\rho }_{0}r{\omega }^2}\frac{1}{\sqrt{1-{\mu }^2}(1-b^2{f}^2)}im\widehat{p}\hfill \end{array}$$(B.9)

We insert these expressions into the continuity equation, yielding

$$\begin{array}{c}\hfill \lambda \widehat{p}+\frac{\partial }{\partial \mu }\left(\frac{1-{\mu }^2}{1-b^2{f}^2}\frac{\partial \widehat{p}}{\partial \mu }\right)-\frac{m^2}{(1-{\mu }^2)(1-b^2{f}^2)}\widehat{p}=0,\end{array}$$(B.10)

where

$$\begin{array}{c}\hfill \lambda ={\left(\frac{r\omega k_{\mathrm{r}}}{N}\right)}^2={\left(\frac{b}{a}\right)}^2.\end{array}$$(B.11)

Appendix C: Magnetic field profile

C.1. Latitudinal profile of B_r

To build the latitudinal profile f(θ) introduced in Eq. (13), we need to determine the value of δ corresponding to a given a_asym = ∫₋₁¹f²P₂dμ/∫₋₁¹f²dμ. The denominator can be expressed as I₁e^−δ, where

$$\begin{array}{c}\hfill I_n={\displaystyle {\int }_{-1}^1}{\mu }^{2n}{e}^{3\delta {\mu }^2}\mathrm{d}\mu \end{array}$$(C.1)

To calculate these integrals numerically, we express them as:

$$\begin{array}{c}\hfill I_n=2{\displaystyle \sum _{k=0}^{+\infty }}\frac{{(3\delta )}^k}{k!(2k+2n+1)}.\end{array}$$(C.2)

and terms in the summation are added until they become smaller than 10⁻¹⁵. The numerator can be written as e^−δ(3I₂ − I₁)/2, so that

$$\begin{array}{c}\hfill a_{\mathrm{asym}}=\frac{3}{2}\frac{I_2}{I_1}-\frac{1}{2}.\end{array}$$(C.3)

The value of δ corresponding to a given a_asym is obtained by using a Newton-Raphson algorithm and making use of the following expression for the derivative da_asym/dδ:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{a}_{\mathrm{asym}}}{\mathrm{d}\delta }=\frac{9}{2}\frac{I_3}{I_1}-2{(a+\frac{1}{2})}^2.\end{array}$$(C.4)

C.2. Alternative expression for f(θ)

We can also build a latitudinal profile f(θ) for the radial field B_r in the shape f(θ) = P₁(cos θ)+BP₃(cos θ)+CP₅(cos θ). For a given value of the asymmetry parameter a_asym, the field must satisfy Eq. (3). This is achieved provided

$$\begin{array}{c}\hfill a_{\mathrm{asym}}=\frac{\frac{4}{15}+\frac{8}{105}{B}^2+\frac{20}{429}{C}^2+\frac{12}{35}B+\frac{40}{231}BC}{\frac{2}{3}+\frac{2}{7}{B}^2+\frac{2}{11}{C}^2}\end{array}$$(C.5)

We find that the value C ≈ 2.37 maximizes the interval of asymmetry parameters that are reachable. Solutions to the previous equation can be found for −0.252 ≲ a_asym ≲ 0.832. This does not cover the full range of possible values for a_asym, but gives satisfactory coverage for this study at least for the solutions with positive asymmetries. We then divide the obtained profile f(θ) by is maximum absolute value to satisfy the condition max(|f|) = 1.

C.3. Making the field divergence-free

We need to ensure that the magnetic field that we build satisfies the condition ∇ · B = 0. The field is axisymmetric, so that the azimuthal component B_φ(r, θ) can be arbitrarily chosen without affecting the divergence of B. We add a θ-component B_θ(r, θ) in such a way that

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\partial }{\partial r}(r^2{B}_{\mathrm{r}})+\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta B_{\theta })=0\end{array}$$(C.6)

This is achieved if B_θ takes the form

$$\begin{array}{c}\hfill B_{\theta }(r,\theta )=-\frac{B_0}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left[r^{2}g(r)\right]\frac{{\int }_0^{\theta }f(\theta \prime )sin\theta \prime \mathrm{d}\theta \prime }{sin\theta }\end{array}$$(C.7)

This solution is regular at the poles (μ = cos θ → ±1). Indeed, by expending μ = ±(1 − ε), one finds that the numerator is equivalent to e^δε while the denominator varies as $\sqrt{2\epsilon }$.

Appendix D: Priors and results for ABIM

D.1. Priors

When using the ABIM code to fit asymptotic expressions of mixed modes (including effects from rotation and magnetic fields) to the observations (Sect. 6.3), we impose a set of priors on the free parameters. The priors are chosen either uniform (U) over a range [x_min, x_max], “uniform-periodic” (UP), which is a uniform distribution for a periodic parameter (the steps of the walkers are calculated modulo the period), or “uniform-normal” (UG), which is a uniform distribution with Gaussian tails at the edges (σ_min and σ_max are the widths of the tails at both ends).

Our choices for the priors are summarized in Table D.1. For the large separation Δν, we used a uniform prior centered on the value measured by Yu et al. (2018) with an interval of ±1 μHz. The phase terms ε_p and ε_g have a periodicity of one. The priors on α and d₀₁ were designed based on measurements of these parameters for other similar red giants. The uniform prior on the period spacing ΔΠ₁ was chosen to include the apparent period spacing ΔΠ₁^meas (which is lower than the expected value) and extend 6 s above the degeneracy sequence (Deheuvels et al. 2022). For the coupling strength q, we chose a uniform distribution ranging up to 0.2, which covers the expected values on the red giant branch. The uniform prior on the magnetic frequency shift δν₀ (which is related to the quantity ⟨B_r²⟩) was chosen to extend well beyond the values measured by Deheuvels et al. (2023). The asymmetry parameter a_asym takes values between −0.5 and 1, which are in principle all reachable by the latitudinal profile that we chose (Eq. (13)). However, we must have δ → ±∞ to reach the edges of this range, which causes numerical problems. For this reason, we restricted the prior distribution to [ − 0.45, 0.9]. The radial extension of the field is parametrized by σ_B, for which we chose a uniform prior distribution in the range [0.001, 0.009] R_★. The upper limit of this interval corresponds to a field that is nearly constant within the g-mode cavity. Finally, the priors on σ_p and σ_g, which characterize the model errors, were chosen based on model errors computed for the red giants studied in Li et al. (2024).

D.2. Corner plots

In Fig. D.1 and D.1, we show corner plots similar to those of Fig. 9 for KIC 3326392, KIC 5700274, KIC 7458743, KIC 9697268, KIC 10976412, and KIC 11408970. For clarity, the corner plots are restricted to the asymptotic period ΔΠ₁ of dipolar g modes, and the parameters characterizing the core magnetic field: the magnetic shift δν₀ = δω₀/(2π), σ_B, and a_asym.

Fig. D.1. Same as Fig. 9 for KIC 3326392 (left) and KIC 5700274 (right).

Fig. D.2. Same as Fig. 9 for KIC 7458743 (top left), and KIC 9697268 (top right), KIC 10613561 (bottom left) and KIC 10976412 (bottom right).

All Tables

All Figures

Fig. 1. Stretched period échelle diagram for KIC 7458743 obtained from 4-yr Kepler time series. Only peaks with an S/N above eight are shown. The size of the symbols are proportional to the peak amplitude. For clarity, peaks in the vicinity of l = 0 and l = 2 modes are omitted. The échelle diagram is folded with an apparent period spacing of ΔΠ1meas = 76.5 s, which produces the best visual alignment of the ridges.

In the text

	Fig. 2. Stretched period échelle diagrams of the eight Kepler red giants of our sample. The échelle diagrams were folded using apparent period spacings ΔΠ1meas that offer good visual alignment of the ridges. To guide the eye, the ridges are highlighted with different colors (red for the most distorted ridge, blue for the other).
In the text

	Fig. 3. Location of the stars under study in the (Δν-ΔΠ1) plane. The apparent period spacings of g modes ΔΠ1meas are shown as blue star symbols, and the red star symbols indicate the asymptotic period spacing ΔΠ1 measured by accounting for magnetic fields in a non-perturbative manner (Sect. 6). Gray crosses correspond to all red giants studied in Li et al. (2024).
In the text

Fig. 4. Stretched échelle diagram of KIC 7458743. Filled blue circles show the observed mode frequencies. Open red circles show asymptotic mixed mode frequencies including rotation and magnetic field effects, assuming that the observed ridges correspond to m = ±1 components. The core rotation was determined to ensure a crossing of the ridges at νcross = 150 μHz, and magnetic field properties were adjusted to reproduce the m = −1 ridge (see Sect. 3.1).

In the text

	Fig. 5. Numerical solution of the horizontal problem for different values of the asymmetry parameter aasym. We show λ as a function of a for (l = 1, m = 0) modes (top panel) and for m = ±1 modes (bottom panel). The colored dots indicate that mode suppression is reached, and the dashed black line corresponds to b = 1.
In the text

Fig. 6. Period échelle diagrams of l = 1 g modes for the example case described in Sect. 4.6. Left: frequencies of m = 0 (circles) and m = 1 (triangles) modes computed for a field with a radial extension σB = 0.003 R★. Magnetic effects are treated either with a non-perturbative approach (red symbols), or a perturbative approach (gray symbols). Right: frequencies of m = 0 modes for values of σ = 0.0015 (red), 0.002 (green), and 0.003 (blue) and for a constant field (black).

In the text

Fig. 7. Impact of the radial profile of Br on the mode frequencies. We considered fields with radial extensions of σB = 0.0015 R★ (red lines), 0.002 R★ (green lines), and 0.003 R★ (blue lines) and a constant field (black lines). Left: radial profile of Br. Middle: profile of the parameter a(r) in the g-mode cavity for the mode closest to νmax (Eq. (12)). Right: profile of λ(r) deduced from the interpolation of the a(λ) relation for the same mode (see text). We show the departure from the value λ = 2 corresponding to the perturbative approach.

In the text

	Fig. 8. Stretched échelle diagrams of three stars of our sample. Blue circles and horizontal blue lines correspond to observed modes and their 1σ uncertainties. Red empty circles correspond to the best asymptotic models including magnetic effects in a non-perturbative manner.
In the text

	Fig. 9. Corner plots of KIC 9227589 (left) and KIC 11408970 (right), restricted to the asymptotic period ΔΠ1 of dipolar g modes, and the parameters characterizing the core magnetic field: the magnetic shift δν0 = δω0/(2π), σB, and aasym. The corner plots for other stars are given in Appendix D.
In the text

	Fig. D.1. Same as Fig. 9 for KIC 3326392 (left) and KIC 5700274 (right).
In the text

	Fig. D.2. Same as Fig. 9 for KIC 7458743 (top left), and KIC 9697268 (top right), KIC 10613561 (bottom left) and KIC 10976412 (bottom right).
In the text