© The Authors 2026

1. Introduction

The problem of the transport of angular momentum is one of the major obstacles in modern stellar physics. Helioseismic measurements of the internal rotation of the Sun (Schou et al. 1998; Chaplin et al. 1999) and subsequent spatial asteroseismology missions with CoRoT (Baglin et al. 2006) and Kepler (Borucki et al. 2010) led to a revolution in the field. Internal rotation rates were measured for main-sequence stars (Kurtz et al. 2014; Benomar et al. 2015; Li et al. 2020), red giant stars (Beck et al. 2012; Mosser et al. 2012b; Deheuvels et al. 2014; Triana et al. 2017; Gehan et al. 2018; Kuszlewicz et al. 2023), and white dwarfs (Hermes et al. 2017). These results, when compared to the predictions of current stellar models for main-sequence stars (Ouazzani et al. 2019) and red giants (Marques et al. 2013; Eggenberger et al. 2012; Ceillier et al. 2013) show that core rotation rates are largely overestimated by the models, and this demonstrates that additional physical mechanisms are needed to properly explain and understand the evolution of rotation in stars. It was suggested that angular momentum could be transported by internal gravity waves (Talon & Charbonnel 2005) or by magnetic fields (Eggenberger et al. 2005) in stellar interiors.

Magnetic fields have been shown to be promising candidates for the efficient redistribution of angular momentum (e.g. Cantiello et al. 2014; Jouve et al. 2015; Fuller et al. 2019; Meduri et al. 2024). They have been observed at the surface of stars for a long time (Landstreet 1992; Donati & Landstreet 2009), where they are generally linked to dynamo effects in convective envelopes (Donati & Landstreet 2009). Although magnetic fields are expected to be found in stellar interiors, detecting such fields was not possible until recently. Magnetic fields were measured for the very first time in the cores of red giant stars (Li et al. 2022, 2023; Deheuvels et al. 2023; Hatt et al. 2024) through their effects on the frequencies of so-called mixed oscillation modes. In red giants, non-radial modes are mixed modes, which means that they behave as pressure (p) modes in the envelope and as gravity (g) modes in the core, offering an opportunity to probe the deepest layers of the stars. Weak magnetic fields have very recently been discovered in the core of a slowly rotating γ Dor star (Takata et al. 2026) through study of asymmetric dipolar g-mode multiplets and linking their existence to predominantly toroidal fields. Another study associated the nine components of a quadrupolar (l = 2) low-frequency multiplet of a β Cep pulsator with the presence of a magnetic field in the core and measured its strength (Vandersnickt et al. 2025).

The effects of magnetic fields on oscillation modes were studied and explored prior to their discovery. Magnetic fields, if weak enough, can be interpreted as small perturbations of the oscillation equations (Gough & Thompson 1990). The induced perturbations on mode frequencies were computed in the specific case of a dipolar field aligned with the rotation axis for pure g modes (Hasan et al. 2005) and red giant mixed modes (Gomes & Lopes 2020; Mathis et al. 2021; Bugnet et al. 2021). The case of inclined dipolar fields was addressed by Loi (2021). The generalised case was addressed in Li et al. (2022), where the only assumption made was that the toroidal component of the field does not completely dominate over the radial component. These studies showed that weak magnetic fields cause a departure from the asymptotic expressions of g-mode frequencies in the form of a frequency shift that strongly depends on the eigenfrequency of the mode. This shift also creates asymmetries in the rotational multiplets of non-radial mixed modes, depending on the geometry of the field. Li et al. (2022, 2023) exploited this property of asymmetry to make the first detections of magnetic fields in red giants using dipolar (l = 1 mixed modes). The authors measured the intensities of the fields in the core and placed constraints on their topology. Deheuvels et al. (2023) detected strong near-critical magnetic fields through their impact on the regularity of the g-mode pattern. Hatt et al. (2024) studied low-luminosity, high signal-to-noise ratio (S/N) red giants to find magnetic field signatures in mixed mode triplets. Today, we count 48 direct detections of magnetic fields in red giant stars. Other studies have proposed magnetism as an explanation for the existence of a category of red giants with ‘suppressed’ dipolar modes, that is, dipolar modes with unexpectedly low amplitudes (Mosser et al. 2012a; García et al. 2014b). It has been suggested that these stars harbour core magnetic fields that exceed the critical field above which magneto-gravity waves cannot propagate (Fuller et al. 2015; Stello et al. 2016; Cantiello et al. 2016). However, the interpretation on the nature and origin of suppressed modes is under debate (Mosser et al. 2017). Super-critical fields have also been invoked to account for oscillation mode suppression in a B-type main-sequence pulsator (Lecoanet et al. 2022).

Li et al. (2022) suggest that magnetic fields, if unaccounted for, can bias the measurement of g-mode parameters, such as the asymptotic period spacing (ΔΠ₁) and the gravity offset (ϵ_g). The latter parameter is linked to the properties near the turning points of the g-mode cavity (Pinçon et al. 2019). The gravity offset has been shown to have a well-determined value for all red giant stars at ϵ_g = 0.28 ± 0.08 (Mosser et al. 2018) such that a deviation from this value is a priori not expected. The largest catalogue of seismic measurements of red giants to date was made by Li et al. (2024). They measured the rotation and other seismic parameters (including ϵ_g) for 2495 stars but without taking magnetic effects into account. This makes for an interesting opportunity to test the effect magnetic fields have on the measure of ϵ_g according to the suggestion of Li et al. (2022).

In this study, we investigate the use of the gravity offset as an indicator for magnetic fields using the catalogue of Li et al. (2024). In Sect. 2 we explain the effect the magnetic field has on the measurement of ϵ_g. In Sect. 3 we present the data used in this study and the steps taken to extract mixed mode frequencies from the oscillation spectra. Then, in Sect. 4 we present the asymptotic model we used to describe mode frequencies, accounting for both rotation and magnetic fields. This model is then adjusted to the data using Bayesian inference. For our magnetic detections, we measured the intensity of the field. Our results are compared to previous studies in Sect. 5. Lastly, in Sect. 6 we explore the distributions of different seismic parameters for the total sample of magnetic detections.

2. Gravity offset as a probe of core magnetic fields

Magnetic fields have proven hard to detect in preceding studies. Out of around 2000 targets, only 48 have confirmed detections of magnetic field. Increasing the number of detections is vital to the study of the origin and prevalence of magnetic fields. Past studies have established selection criteria to increase the chance to detect magnetic fields, such as looking for triplet asymmetries (Li et al. 2023), or searching for irregularities in the measurements of the asymptotic period spacing ΔΠ₁ in Deheuvels et al. (2023) for near-critical fields.

As presented earlier, Li et al. (2022) suggested that weak magnetic fields, if not accounted for, can bias the measurement of g-mode parameters such as ΔΠ₁ and ϵ_g. To illustrate this, we generated a set of g-mode periods using the asymptotic expression P_g = ΔΠ₁(n_g + 1/2 + ϵ_g). We then added the effects of rotation and a moderate magnetic field following a perturbative approach (this step is explained in more detail in Sect. 4.1). The resulting mode periods are represented in the left panel of Fig. 1 in the shape of an échelle diagram folded with the true period spacing ΔΠ₁. The two main effects of magnetic fields appear clearly. Since magnetic frequency shifts strongly depend on frequency (∼ω⁻³, see Mathis et al. 2021; Li et al. 2022), low-frequency modes are more affected, resulting in a bending of the ridges to the left at low-frequency in the échelle diagram. Secondly, the shift differs for m = 0 and m = ±1 components, which creates an asymmetry of the triplet. If magnetic fields are ignored, the m = 0 modes are fitted as P′_n = ΔΠ₁′(n_g + 1/2 + ϵ_g′), where the apparent ΔΠ₁′ is chosen so that the m = 0 ridge is as close as possible to a straight vertical line in the échelle diagram (red dashed line in Fig. 1). The result is shown in an échelle diagram folded with the apparent period spacing ΔΠ₁′ in the right panel of Fig. 1. For a weak field, the value of ΔΠ₁′ differs from ΔΠ₁ by only a few tenths of a second, which is not enough to make these measurements abnormal. On the contrary, the value of ϵ_g is much more affected, and ϵ_g′ can significantly differ from the expected value of 0.28 ± 0.08. As seen in the right panel of Fig. 1, the blue area represents the zone in which we would expect to find most modes considering the empirical value of ϵ_g. This is not the case here where the measured ϵ_g′ is subject to significant distortions (it corresponds to ϵ_g′=0.42 in the example shown in Fig. 1).

Fig. 1. Échelle diagrams of mock g-mode spectra perturbed by rotation and magnetic field. Left: Spectrum fitted with an asymptotic expression accounting for magnetic and rotational perturbations, folded by the true ΔΠ1. Right: Spectrum fitted with an asymptotic expression for rotation only, folded by the measured ΔΠ1′. The red lines represent the slope of the m = 0 component, which is purely vertical on the right. The blue area represents the area where m = 0 modes are expected to be found using ϵg = 0.28 ± 0.08.

Li et al. (2024) measured mixed mode parameters for 2006 red giants with great precision. The authors also measured core and envelope rotation rates for red giants showing at least two components per dipole triplet. However, these measurements did not take magnetic fields into account. In Fig. 2, we show the measured values of ϵ_g as a function of the large separation Δν of p-modes. Clearly, most stars are found in the expected range for ϵ_g, but a smaller population of red giants have ϵ_g that clearly deviate from the empirical value (these stars are highlighted in red in Fig. 2). Following the behaviour explained above, these stars are interesting candidates for hosting a magnetic field.

Fig. 2. Relationship between ϵg and Δν in the Li et al. (2024) catalogue. Red dots represent stars selected for this study, in contrast to the grey dots. Solid blue lines represent the delimitation of the bulk distribution at ±σ, and the dashed blue line is made from the median of every 118 stars slices.

To select these candidates, we first calculated the distribution of ϵ_g in slices of Δν, each slice containing 118 stars. The standard deviation and median were then calculated for each slice and used to delimit the bulk of the full distribution. The stars that lie outside of ±1σ of the delimitation were taken as part of the sample. Also, in this study, we chose to focus on low-luminosity red giants, using as a constraint Δν >  10 μHz. Indeed, for higher-luminosity giants, the oscillation spectra become more complex, the rotational splitting being generally larger than the frequency separation between mixed modes. Higher-luminosity stars will be treated in a later study. In total, we selected 218 low-luminosity red giants with abnormal values of ϵ_g for analysis.

3. Analysis of oscillation spectra

3.1. Seismic data and pre-processing

In this study, we used the full Kepler (Borucki et al. 2010) dataset, corresponding to nearly 4-year long high-precision photometry time series. The oscillation spectra were downloaded from the Mikulski Archive for Space Telescopes (MAST). We used the calibrated KEPSEISMIC¹ data, high-pass filtered at 20 days (García et al. 2011, 2014a; Pires et al. 2015). We then fitted to the observed spectra a model of the background in order to determine the S/N of the power spectral density (PSD) of each star. The background was fitted using two Harvey laws for the granulation, a white noise component for photon noise, and a Gaussian function corresponding to the pressure mode bump. The fit was done with a Maximum Likelihood estimation (MLE), using apollinaire (Breton et al. 2022). Initial guesses for this estimation were chosen based on the results of previous fits for similar stars.

3.2. Extraction of mode frequencies

The signature of a magnetic field is found either by detecting asymmetry in the dipole triplet (only if three components are visible) or by measuring deviations from a constant g-mode period spacing. This translates into curved ridges in the stretched échelle diagram (similarly to Fig. 1), where stretching designates a transformation of mixed mode periods to account for the variation of period spacing (see Mosser et al. 2015). In the latter case, the detection of a field crucially depends on modes in the low-frequency part of the observed spectrum because magnetic effects are largest at low frequency. However, these modes also have lower S/N. For this reason, special care is needed to (i) maximise the probability of detection of lower-frequency modes, and (ii) limit as much as possible the false-alarm probability, in order to avoid false detections of magnetic fields. Many seismic studies assume an empirical detection threshold of eight times the background noise level (e.g. Mosser et al. 2015). In this work, we developed a more robust detection criterion.

We first elaborated a basic, conservative S/N threshold over which the signal is considered as a detection. To do so, we considered each bin in the oscillation spectrum as a realisation of a random variable X following a χ² distribution with two degrees of freedom (e.g. Appourchaux et al. 1998). Here two possible outcomes exist: either the realisation of X is due to noise (hypothesis H₀), or it is signal (hypothesis H₁). We want to define a threshold x_p in the PSD over which peaks are considered to be due to signal. The false-alarm probability, known as the p-value, corresponds to the probability of having at least one peak due to noise that exceeds x_p among the N considered frequency bins, that is, P( ∃i ∈ [​[1, N]​] ;  X_i >  x_p | H₀)=p. For a given p-value, assuming that all the bins are independent from one another (this is a good approximation because the duty cycle of Kepler is above 90%), the threshold x_p can be expressed as (e.g. Deheuvels 2010)

$$\begin{array}{c}\hfill x_p=-ln[1-{(1-p)}^{1/N}].\end{array}$$(1)

We consider modes in the frequency range ν_max ± 5Δν. For 4-yr-long Kepler datasets (with a frequency resolution of ∼7 nHz), this selects N ∈ [14 000; 25 000] bins for each star. For a p-value of 5%, we obtain x_p ≈ 12. This threshold guarantees a 95% chance of having no peak due to noise exceeding x_p in the considered frequency range.

For stochastically excited modes, the mode lifetime will generally be shorter than the duration of the observations, so that these modes are spread over several bins of the spectrum. We can concentrate the power of a mode by smoothing the spectrum, redistributing power into fewer bins, to maximise detection probability (Appourchaux 2004). The width of the smoothing window, hereafter ‘boxcar width’, is usually taken as a multiple of the mode width. Two courses of action are possible when we smooth the spectrum. Either re-binning the spectrum over the N bins of the boxcar, or keeping every bin, in which case we lose the independence between bins and the threshold cannot be computed analytically. Here, we chose the latter because the first option yields a significantly worse detection probability. We thus used a Monte Carlo approach to infer the threshold in S/N corresponding to a false alarm probability of 5% for each boxcar width. We generated series of pure noise (random 2-d.o.f. χ²) with the same number of bins as the real spectrum, and we applied a smoothing over the chosen boxcar. This procedure was repeated over 1000 iterations. Using this, the threshold is given by the height in S/N that only 5% of the realisations exceed.

We then searched for the boxcar width that maximises our chances of detecting the mode. For this purpose, we generated synthetic spectra containing a single mode, to which we applied noise following a 2-d.o.f. χ² distribution. This procedure was repeated over 10 000 iterations. For different mode energies (proportional to the product of mode height and width), we measure the detection probability for several boxcar widths. Using this, we were able to determine that the probability of detection is at its peak when the boxcar width is about three times the mode width (see Fig. 3). We used this boxcar size in our subsequent analysis.

Fig. 3. Probability of detection of a mode versus the boxcar width, in units of the mode width. The colour of the lines indicates the mode energy, defined as E ∝ Γh, where Γ is the mode width and h is mode height. The vertical dash-dotted line corresponds to a boxcar width of 3Γ.

In our case, the lifetimes of mixed modes vary depending on the ratio of mode inertia stored in the acoustic cavity of the buoyancy cavity. Modes that are more g-like have longer lifetimes and have narrow mode profiles in the PSD, with a width close to the frequency resolution, while p-dominated modes span several bins of the spectrum. The width of a mixed mode can be expressed as (Mosser et al. 2015)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{l=1}={\mathrm{\Gamma }}_{l=0}(1-\zeta ),\end{array}$$(2)

where ζ is the ratio of mode inertia stored in the buoyancy cavity to the total mode inertia, which varies with frequency and can be computed using the measured frequencies of mixed modes (Goupil et al. 2013). Here, we compute the expected trapping of the modes using parameter ζ, determined using the parameters in Li et al. (2024). Using this we smooth the PSD with a boxcar width of 3Γ_l = 1. We then use the first Monte Carlo scheme described earlier to compute the threshold at each boxcar width, with respect to its corresponding frequency. This results in a threshold that varies with the smoothing of the spectrum. After removing the radial and quadrupolar modes from the spectra, we show an example of this smoothing and threshold in Fig. 4, which presents a portion of the spectrum of KIC 9999911. Clearly, modes are more smoothed in p-mode regions (∼142.5 to 145 μHz) than in g-mode regions (the rest of the spectrum), where it is not smoothed. We can also notice how the variable threshold allows for the detection of weaker p-mode components (we detect one component at ∼143.5 μHz that was not detected with a flat threshold).

Fig. 4. Spectrum of KIC 9999911. Top: Unprocessed spectrum. Bottom: Pre-processed spectrum with variable width smoothing and radial and quadrupolar modes cut out. Red arrows indicate the position of detected modes using the red dashed-dot threshold.

Using this new threshold, we were able to obtain a list of reliable l = 1 and l = 3 mode frequencies guesses from the spectra. In practice, several peaks sometimes belong to the same mode and obtaining one frequency guess per mode requires an additional step, which is detailed in Appendix A.

The last step consists of extracting the mode frequencies from the PSD using the guesses obtained above. For this purpose, we modelled each mode as a Lorentzian profile with free central frequency, height, and width. The background was fixed to the estimate that was made in Sect. 3.1, using the contributions from granulation and the photon noise. The model of spectrum was fitted to the PSD for each overtone with detected modes using maximum likelihood estimation (MLE) similarly to what was done e.g. in Deheuvels et al. (2020). The errors were then estimated by computing the covariance matrix as the inverse of the hessian matrix.

4. Measurement of internal magnetic fields

4.1. Asymptotic expression of mixed mode frequencies with rotation and magnetic fields

We used the asymptotic expression for dipolar mixed mode frequencies in the oscillation spectrum, assuming the validity of the WKB approximation. Shibahashi (1979) and Unno (1989) have shown that the mixed mode frequencies ν can be approximated by solving

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

where q is the coupling factor; θ_p, m(ν) and θ_g, m(ν) are the p and g-mode phases, respectively; and the subscript m refers to the dependence in the azimuthal order. Following Mosser et al. (2015), we express θ_p, m and θ_g, m as

$$\begin{array}{c}\hfill {\theta }_{p,m}(\nu )=\pi \frac{\nu -{\nu }_{p,m}}{\mathrm{\Delta }\nu (n_p)},{\theta }_{g,m}=\frac{\pi }{\mathrm{\Delta }{P}_g(n_g)}(\frac{1}{\nu }-P_{g,m})-\frac{\pi }{2},\end{array}$$(4)

with Δν(n_p)=ν_p, m(n_p + 1)−ν_p, m(n_p) as the local frequency large separation at radial order n_p and ΔP_g = P_g, m(n_g + 1)−P_g, m(n_g) as the local period spacing at radial order n_g. Frequencies ν_p, m are linked to the pure dipolar unperturbed p mode frequencies ν_p⁽⁰⁾, which can be expressed using an asymptotic expression of p-modes as Mosser et al. (2015)

$$\begin{array}{c}\hfill {\nu }_p^{(0)}=[n_p+\frac{1}{2}+{ϵ}_p+\frac{\alpha }{2}{(n_p-n_{\max })}^2]\mathrm{\Delta }\nu -\delta {\nu }_{01},\end{array}$$(5)

with ϵ_p as the p mode offset. Here, α and δν₀₁ are parameters introduced by the second-order asymptotic expression. The parameter δν₀₁ characterises the small separation between l = 0 and 1 modes, and n_max = ν_max/Δν. The periods P_g, m are linked to the pure dipolar unperturbed g-mode periods P_g⁽⁰⁾:

$$\begin{array}{c}\hfill P_g^{(0)}=\mathrm{\Delta }{\mathrm{\Pi }}_1(n_g+1/2+{ϵ}_g).\end{array}$$(6)

Equations (4) and (6) appear to be slightly different from the ones written by Mosser et al. (2015), but they are strictly equivalent and provide the same spectrum of mixed modes. However, with our expressions, g-dominated mixed mode frequencies tend to g-mode frequencies when q tends to zero.

The frequencies ν_p, m and periods P_g, m contain the contributions of rotation and magnetic field as perturbations to the model, which we formulate here. Since red giants are slow rotators, rotation can indeed be treated as a first-order perturbation. In this case, dipolar modes are split into three components, forming a symmetrical triplet. The frequency shifts induced by rotation between each component can be expressed as ν_R, p = ⟨Ω⟩_p/(2π) and ν_R, g = ⟨Ω⟩_g/(4π) for p and g modes, respectively, and are related to the average rotation rate ⟨Ω⟩ of their respective cavities (see Goupil et al. 2013).

Magnetic fields can also be treated as first-order perturbations if the field strength B stays small compared to the critical field B_c. Several studies investigated special cases of magnetic fields with specific geometries (Gomes & Lopes 2020; Mathis et al. 2021; Bugnet et al. 2021; Loi 2021). The generalised case was treated in Li et al. (2022), where the only assumption is that the toroidal component of the field does not completely dominate over the radial component. These expressions are valid as long as B_ϕ/B_r ≪ N/Ω ∼ 100 in red giants, where N is the Brunt-Väisälä frequency and Ω is the rotation of the star. In the most general case, when the effects of non-axisymmetry of the field are strong enough, each component of the dipolar mode may be divided in 3 components, with a total of 9 components. Li et al. (2022) showed that, when the rotational splitting is smaller than the magnetic shift such that b = ν_B/ν_R <  1, non-axisymmetric effects become negligible, and frequencies behave as if the field were axisymmetric. In this study, we make the hypothesis that this property is always verified and we assess the validity of that assumption in Sect. 6.5.

In our context, the frequencies of p modes ν_p, m and the periods of g modes P_g, m perturbed by rotation and magnetic field can be written as

$$\begin{array}{c}\hfill {\nu }_{p,m}={\nu }_p^{(0)}-m{\nu }_{\mathrm{R},p}\mathrm{and}{P}_{g,m}={(\frac{1}{P_g^{(0)}}-m{\nu }_{\mathrm{R},g}+f_m(a){\nu }_B)}^{-1},\end{array}$$(7)

where the term f_m(a)ν_B represents the magnetic perturbation. Magnetic fields were shown to affect mostly g-mode oscillations (Mathis et al. 2021), which is why we did not consider a magnetic perturbation to p-mode frequencies. Following Li et al. (2022), we can express the magnetic perturbation term of Eq. (7) as

$$\begin{array}{c}\hfill {\nu }_B={\nu }_{\mathrm{mag}}{\left(\frac{{\nu }_{\max }}{\nu }\right)}^3\mathrm{and}{f}_m(a)=1-a(1-\frac{3}{2}|m|),\end{array}$$(8)

where a is the asymmetry parameter, linked to the latitudinal structure of the field (Li et al. 2022). The value of a varies between −0.5 and 1, for a field concentrated at the equator and at the poles, respectively. The parameter ν_mag is the average magnetic frequency shift at ν_max, which is linked to the intensity of the magnetic field by the expression

$$\begin{array}{c}\hfill {\nu }_{\mathrm{mag}}=\frac{\mathcal{I}}{2\pi {\mu }_0{(2\pi {\nu }_{\max })}^3}⟨B_r^2⟩.\end{array}$$(9)

Here, ℐ is an integral depending on the core structure, and ⟨B_r²⟩ is an average of the radial component of the magnetic field squared weighted by the kernel function K_r(r). It is defined as

$$\begin{array}{c}\hfill ⟨B_r^2⟩={\displaystyle {\int }_{r_1}^{r_2}}{K}_r(r)\overline{B_r^2}\mathrm{d}r,\end{array}$$(10)

where $\overline{B_r^2}$ is the average of the field in the horizontal direction. Also, as K_r(r) peaks in the hydrogen burning shell (HBS) around the core, ⟨B_r²⟩ has a high sensibility to the field in this region (see extended Figure 1 from Li et al. 2022).

In practice, the visibility of mode components impacts the measurement of the magnetic perturbation. Indeed, in the best-case scenario of a triplet, both ν_mag and a are measurable. But in other cases, only the product of f_m(a)ν_mag is measurable. In the cases of a doublet or singlet, we measure

$$\begin{array}{cc}& {\stackrel{\sim }{\nu }}_{\mathrm{mag}}=f_{\pm 1}(a){\nu }_{\mathrm{mag}}=(1+a/2){\nu }_{\mathrm{mag}}\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\stackrel{\sim }{\nu }}_{\mathrm{mag}}=f_0(a){\nu }_{\mathrm{mag}}=(1-a){\nu }_{\mathrm{mag}}\hfill \end{array}$$(12)

where the tilde denotes the observable parameter. Since a is bounded between −1/2 and 1, ν_mag lies between $\frac{2}{3}{\stackrel{\sim }{\nu }}_{\mathrm{mag}}$ and $\frac{4}{3}{\stackrel{\sim }{\nu }}_{\mathrm{mag}}$ for doublets, and above $\frac{2}{3}{\stackrel{\sim }{\nu }}_{\mathrm{mag}}$ for singlets.

To obtain the mixed mode frequencies, we had to solve Eq. (3) for each azimuthal order m while considering the perturbations of the rotation and magnetic field introduced by Eq. (7). In the next subsection, we detail the process used to adjust these modelled frequencies to the observations.

4.2. Bayesian inference

We fitted the observed frequencies of radial and dipolar modes with the asymptotic model described in Sect. 4.1. The free parameters of the model are those describing the p-mode spectrum (Δν, ϵ_p, α and d₀₁), the dipolar g-mode spectrum (ΔΠ₁ and ϵ_g), the coupling parameter q, the mean rotation rates in the p and g cavities (ν_R, p and ν_R, g) and the core magnetic field (ν_mag and a). We fix ν_max at the value measured by Yu et al. (2018). These parameters are gathered in a parameter vector denoted θ. The modelled frequencies are generated using Eqs. (3) to (6) over a given range of n_p and n_g. This range is tailored to cover the whole observed spectrum. In order to estimate the parameters, we use a Bayesian inference approach. Following Bayes’ theorem, the posterior distribution of the parameters θ given a dataset y – which here consists of a set of observed frequencies – is

$$\begin{array}{c}\hfill p(\mathit{\theta }|\mathit{y})=\frac{p_I(\mathit{\theta })p(\mathit{y}|\mathit{\theta })}{p(\mathit{y})},\end{array}$$(13)

where p_I(θ) is the prior distribution of θ, p(y|θ) the likelihood and p(y)=∫p_I(θ)p(y|θ)dθ is the global likelihood. The likelihood p(y|θ) is expressed

$$\begin{array}{c}\hfill p(\mathit{y}|\mathit{\theta })=\frac{1}{\sqrt{{(2\pi )}^{n}det\mathbf{\Sigma }}}exp[-\frac{1}{2}{[\mathit{y}-\mathit{\mu }(\mathit{\theta })]}^T{\mathbf{\Sigma }}^{-1}[\mathit{y}-\mathit{\mu }(\mathit{\theta })]],\end{array}$$(14)

where y = (ν_i, obs)_{i = 1, n} is a vector containing the n observed frequencies. The observed modes are assigned to the closest modelled modes. μ(θ)=(ν_i(θ))_{i = 1, n} is a vector containing the frequencies modelled with a set of parameters θ, and Σ = Σ_obs + Σ_mod a covariance matrix composed of two terms. The first one, Σ_obs, is the observational error covariance matrix. As the frequency measurements are independent, it boils down to a diagonal matrix Σ_obs = diag(σ_1, obs²,…,σ_n, obs²), where σ_i, obs is the error of ν_i, obs. The second term, Σ_mod, takes into account the errors of the model. Indeed, as the asymptotic expression presented in Sect. 4.1 are derived from simplified versions of the equations of stellar oscillation, we expect deviations from these expressions. To take them into account, Li et al. (2024) added a fixed error. For this work, we modelled Σ_mod as a diagonal matrix with terms

$$\begin{array}{c}\hfill {\sigma }_{i,\mathrm{mod}}^2={\zeta }_i{\sigma }_g^2+(1-{\zeta }_i){\sigma }_p^2\end{array}$$(15)

on its diagonal. We introduced two free parameters, σ_g and σ_p, which characterise the errors on g- and p-modes, respectively. The quantity ζ_i is the ratio of inertia in the g cavity over the total inertia for the mode with frequency ν_i. It is evaluated from the parameters q, ΔΠ₁ and the p-mode frequencies ν_p (see e.g. Eq. 4 from Gehan et al. 2018). The parameters σ_g and σ_p are estimated during the fitting process. This is a simplified model since the model errors of close frequencies are correlated and their distribution is not necessarily normal. Thus, we performed tests on realistic synthetic spectra to verify that the errors of parameters θ correctly take into account deviations from the asymptotic regime. Our approach appears to be quite conservative and tends to overestimate the errors.

The prior distributions p_I(θ) capture all our a priori knowledge (or lack thereof) on those parameters. Compiled in Table 1, the priors were chosen in order to cover the known range of the parameter values. If this range is not or badly known, a large and uninformative prior was chosen for that parameter. Four types of prior distributions were used: uniform over the interval [x_min, x_max], ‘uniform-periodic’, which is a uniform prior for which the resulting posteriors are the results of steps by the walkers modulo the period, ‘uniform-normal’, which is a uniform distribution with Gaussian edges of width σ_min and σ_max, and a modified Jeffrey prior, which is an uninformative, scale-invariant prior over [x_min, x_max], and smoothly transitions to a quasi-uniform distribution over [0, x_min] (see detailed descriptions in e.g. Handberg & Campante 2011). These appellations are abridged in Table 1 as ‘U’, ‘UP’, ‘UG’, and ‘J’, respectively.

For the parameters ϵ_p and ϵ_g, which have a periodic behaviour (see Sect. 2), we used a uniform-periodic distribution. For Δν, we used a uniform prior around the value measured by Li et al. (2024, noted Δν₀), with a width of 6 μHz to have a relative freedom for this parameter.

For parameter ΔΠ₁, we chose a normal-uniform prior centred on a value ΔΠ_1, 0, with a spread of 14 s, setting σ_max = 2 s and σ_min = 10 s, thus allowing for lower and higher values of ΔΠ₁. The value of ΔΠ_1, 0 was obtained from the degeneracy sequence, using the knowledge of Δν. In the case where the electron gas in the core of the star is degenerate, stars of mass M coalesce on a sequence in the (Δν, ΔΠ₁) plane, referred to as the degeneracy sequence (Deheuvels et al. 2022). In the case where the core is not fully degenerate, stars are above the degeneracy sequence (ΔΠ₁ >  ΔΠ_1, 0). There also exist cases where ΔΠ₁ <  ΔΠ_1, 0, which are attributed to mass transfer (Deheuvels et al. 2022). These cases justify the asymmetric spread of the prior of ΔΠ₁.

For the rotation, we chose priors in agreement with the measurements for core and envelope rotation in Li et al. (2024). For the magnetic shift ν_mag, we chose as an uninformative prior a modified Jeffrey prior with bounds 0.05 μHz and 10 μHz, which is scale-invariant. Then, because parameter a only takes values between −1/2 and 1, its uniform prior covers all of possible values. Lastly, the model error σ_g on g-mode frequency has a uniform prior between 0 and 0.1 μHz and the model error on p modes is normal-uniform between 0 and 0.3 μHz. These values were chosen based on model errors computed for red giants (see Li et al. 2024). It is also useful to mention that in the case of a star showing a doublet (i.e. two ridges in the stretched échelle diagram) or a singlet (i.e. one ridge), we cannot measure both ν_mag and a, so we fitted the product ${\stackrel{\sim }{\nu }}_{\mathrm{mag}}=f_m(a){\nu }_{\mathrm{mag}}$, as mentioned in Sect. 4.1. In the case of singlets, rotational parameters are not fitted.

To sample the posterior distributions, we used the ‘Asteroseismic Bayesian Inference by Monte Carlo Markov Chain’ (ABIM) code that we developed. It is a Fortran code parallelised with OpenMP directives. The sampling was performed with a Markov chain Monte Carlo method implementing the ‘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. We typically sample the posterior distributions with 20 parallel chains and 300 walkers for each chain. The initial positions of the walkers are randomly drawn from the prior distributions. The walkers are iterated over 13 000 steps, with the 3000 first steps discarded as burn-in to ensure that chains are stationary.

After the sampling process, we took care in assessing the completeness of the sampling for every parameter. One way to do this is to observe the resulting corner plots of the distributions. We also checked the correlation times, which are defined as the time in steps necessary for the autocorrelation function of the chain to diminish by a factor of e. A high correlation time τ_c means that at least τ_c consecutive samples are highly correlated, which signals a bad sampling of the parameter space. In our case, a typical value for a thorough sampling is around τ_c ≈ 30, which ensures to get about 100 000 uncorrelated sampling points.

An example for a favourable detection of magnetic field can be seen in Fig. 5, showing corner plot reduced to magnetic and g-mode parameters. We can clearly see that the posterior distribution of ν_mag drops to zero for low values of ν_mag, which are strongly excluded because of the negative asymmetry detected in this star. We can also notice that the distribution of ϵ_g spans across the interval of its prior and has large error-bars. This is most likely due to the small number of detected g modes which are directly used in the measure of ϵ_g. A representation of the stretched échelle diagram (thereafter s.e.d.) of KIC 4660930 can be seen in Fig. 6 along with the two following examples. There, we can observe the clear asymmetry of the component of the triplet, as well as a significant shift to the left of all three mode components at lower frequency due to the magnetic field. Another example, this time of a doublet, is shown in Fig. 7. Here, the asymmetry parameter cannot be measured because of the lack of the m = 0 mode component. Nonetheless, we can clearly observe the distribution of ν_mag to have a non-zero value because of the strong curvature in the s.e.d. A last example is shown in Fig. 8, where the detection is not favourable. Here we see a sharply decreasing distribution for ν_mag, with its maximum value at ν_mag = 0. We also notice a very large distribution for the asymmetry, almost reproducing the prior but favouring a = 0. The distribution of ν_mag is a strong indicator that this observation is highly compatible with a null magnetic field. This can also be seen in the s.e.d. where the detected frequencies are well explained without magnetic field: the triplet is symmetrical and no discernible shift is observable on the diagram. Among the 218 stars analysed, 13 were dismissed because of a low S/N ratio and one was removed because of its suppressed dipole modes that could not be analysed. One other star had to be removed from the sample because its ν_max was too close to the Nyquist frequency, causing significant aliasing at high frequency.

Fig. 5. Corner plot of KIC 4660930 reduced to only show g-mode and magnetic parameters.

Fig. 6. Stretched échelle diagrams of KIC 4660930, KIC 6222530, and KIC 3634488. Blue dots represent detected frequency, and red circles represent the maximum a posterior model of the determination.

Fig. 7. Corner plot of KIC 6222530 reduced to show only g-mode parameters and the magnetic shift νmag.

Fig. 8. Corner plot of KIC 3634488 reduced to show only g-mode and magnetic parameters.

4.3. Probability of magnetic field detection

In this section we explain the different steps leading to the measurement of core magnetic fields for the cases where a field was significantly detected.

Assessing the probability of detection of a magnetic field using the posterior distribution can be ambiguous. Stars for which non-magnetic models (i.e. with ν_mag ∼ 0 μHz) have been seldom visited during the Monte Carlo Markov chain sampling are good candidates for magnetic field detection, but we need to quantify this. For this, we compare a non-magnetic model M₀, where ν_mag is fixed to 0, and a magnetic model M₁ by computing the odd ratio of M₀ over M₁ expressed as

$$\begin{array}{c}\hfill O_{0,1}=\frac{p(M_0|\mathit{y})}{p(M_1|\mathit{y})}=\frac{p_I(M_0)}{p_I(M_1)}\frac{p(\mathit{y}|M_0)}{p(\mathit{y}|M_1)},\end{array}$$(16)

where y represents the dataset, p(M|y) represents the probability of the model given the data, p_I(M) is the model prior and p(y|M) represents the global likelihood, or evidence. Here we consider that both models are equiprobable. This implies that the first factor of Eq. (16) cancels, leaving what is known as the Bayes factor B_0, 1. To estimate B_0, 1, we computed the global likelihood for models M₁ and M₀ by taking advantage of properties on the parallel tempering scheme we used (as done e.g. by Benomar et al. 2009). However, in this particular case, we can simplify the computation by taking advantage of M₀ being included in M₁ (M₀ is indeed a special case of M₁ with ν_mag = 0).

We considered that in M₀ and M₁, all parameters except ν_mag constitute nuisance parameters and are represented by a vector a. These parameters share the same support 𝒜 and the same priors p_I, a for models M₀ and M₁. In addition, the support of ν_mag for model M₀ (which corresponds to {0}) is included in the support of model M₁. In this case, the expression of B_0, 1 can be simplified: using Bayes law and strategically marginalising the posterior distribution over vector a, accounting for the independence of the priors on a and ν_mag, B_0, 1 can be expressed as

$$\begin{array}{c}\hfill B_{0,1}={\displaystyle {\int }_{\mathcal{B}}}\frac{p_{I,{\nu }_{\mathrm{mag}},0}({\nu }_{\mathrm{mag}})}{p_{I,{\nu }_{\mathrm{mag}},1}({\nu }_{\mathrm{mag}})}{p}_{{\nu }_{\mathrm{mag}}}({\nu }_{\mathrm{mag}})\mathrm{d}{\nu }_{\mathrm{mag}},\end{array}$$(17)

where ℬ is the support of ν_mag and p_νmag(ν_mag) is the marginalised posterior distribution. In our case, ν_mag is fixed to 0 in model M₀. This means that p_{I, νmag, 0}(ν_mag) is a Dirac distribution at ν_mag = 0, reducing (17) to

$$\begin{array}{c}\hfill B_{0,1}=\frac{p_{{\nu }_{\mathrm{mag}}}(0)}{p_{I,{\nu }_{\mathrm{mag}},1}(0)},\end{array}$$(18)

because p_{I, νmag, 1}(0) is non-zero. Using the first equality in Eq. (16), we can express the probability of a favourable detection as

$$\begin{array}{c}\hfill P_1=\frac{1}{1+B_{0,1}}·\end{array}$$(19)

In practice, determining p_νmag(0) can be problematic. Ideally, it could be best determined for a continuous distribution, where p_νmag would be perfectly described by a function p(ν_mag). However, because of the nature of our sampling, p_νmag is a discrete distribution. In order to determine p_νmag(0), we can simply extrapolate the function p(ν_mag) using only the first bin of p_νmag. By varying the bin width h, we can use a Taylor series expansion to construct a system of equations to solve for p(0), which is effectively p_νmag(0). This method is further detailed in Appendix B. While this method is quite reliable when p_νmag has a non-zero value for ν_mag = 0, it is not the case when very few models were explored around zero (in the case of a very clear detection of a magnetic field). In that case, the values can at best be 0 or very low, as expected, and at worst can yield negative values. In the latter case, p_νmag(0) is simply set to zero.

We chose to accept a model as favourable if P₀ >  85%, with 85% being an arbitrary percentage. Among the 207 red giants that were investigated, we find 23 stars with a favourable probability of possessing a magnetic core, the list of which is given in Table 2. This number drops to 18 for a 90% probability and to 15 for a 98% probability.

4.4. Magnetic field intensity

The measurement of the field intensity requires the knowledge of the core structure, according to Eq. (9). We can obtain the value of integral ℐ by computing an evolution model that matches the seismic parameters that we measured using mixed modes. We computed a grid of models using the MESA 10108 release (Paxton et al. 2011, 2013, 2015) in which we searched for the best fitting model. The EOS used is the OPAL (Rogers & Nayfonov 2002) EOS. Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005). The used nuclear reaction rates are from NACRE (Angulo et al. 1999). This grid contains 324 models, varying in mass (0.9 to 2 M_s with steps of 0.1 M_s), metallicity (−0.4 to 0.4 dex with steps of 0.1 dex) and overshoot parameter (0, 0.1 and 0.2). The best fitting model is determined by calculating χ², taking into account ν_max, Δν, and ΔΠ₁ such that

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

where ‘mod’ refers to the value in the model grid, and ‘obs’ refers to the median and standard deviation of the posterior distribution for the corresponding parameter. Fitting for ΔΠ₁ ensures that the core structure is well reproduced and adding Δν and ν_max as constraints ensures that the mass of the model corresponds to the one computed using asteroseismic scaling relations. We managed to find good statistical fits to all our magnetic giants, except for one star, KIC 4350501. This star lies below the degeneracy sequence in the (Δν, ΔΠ₁) plane, which is impossible for single-star evolution. This star was already identified by Deheuvels et al. (2022) among stars that are below the degeneracy sequence, and the authors attributed this phenomenon to mass transfer occurring during the post-main-sequence. This explains why no model of our grid can satisfactorily fit KIC 4350501. Mass transfer modified the amount of mass stored in the envelope, but the structure of the core is unchanged. For this star, we thus computed a χ² using only the measured ΔΠ₁ to reproduce the core structure. We recall that when only m = 0 or m = ±1 are visible, the asymmetry parameter a is not measurable, so that we can only determine intervals for the value of ⟨B_r²⟩. In the case of a singlet, the lower bound of the field intensity is ${⟨B_r^2⟩}_{\min }=\frac{2}{3}\stackrel{\sim }{B_r^2}$, where

$$\begin{array}{c}\hfill \stackrel{\sim }{B_r^2}=\frac{{\stackrel{\sim }{\nu }}_{\mathrm{mag}}2\pi {\mu }_0{\omega }_{\max }^3}{\mathcal{I}},\end{array}$$(21)

whereas for doublets, the field strength ranges from ${⟨B_r^2⟩}_{\min }=\frac{2}{3}\stackrel{\sim }{B_r^2}$ to ${⟨B_r^2⟩}_{\max }=\frac{4}{3}\stackrel{\sim }{B_r^2}$. Using Eq. (9), we could measure the field intensity for the favourable detections. These values, given in Table 2, range from 34 kG to 227 kG, which is compatible with the values found in Li et al. (2022, 2023), Hatt et al. (2024). The error on the field intensity was estimated by scaling the distribution of ν_mag according to Eq. (9) and measuring the standard deviation. We can notice that the majority of the detections are below the critical field B_c, above which magneto-gravity waves no longer propagate in the core (Fuller et al. 2015; Stello et al. 2016; Rui & Fuller 2023). The ratio between field intensity and critical field intensity B/B_c varies between 0.18 and 0.68, with an exception for KIC 12168062 where B/B_c = 1.7. On that star, one can observe that low-frequency modes are suppressed, while higher-frequency modes retain a mixed mode behaviour. This star is further discussed in Sect. 6.6.

5. Comparison with previous studies

5.1. Gravity offset

To build our sample, we chose targets showing abnormal values of ϵ_g, on the ground that the presence of magnetic fields could bias their measurements. To verify this, we plotted the newly measured values of ϵ_g (taking magnetic fields into account) as a function of the previous measurements (ignoring magnetic fields) in Fig. 9. We notice that all the magnetic detections in this study have a measured ϵ_g that is compatible with the value expected for red giants (0.28 ± 0.08). This confirms that for these stars, the abnormal value of ϵ_g that was obtained in Li et al. (2024) was indeed caused by the neglect of magnetic fields. We also found that only 4.8% of our non-detections are more than 2σ away from the expected value. This means that almost all non-detections are compatible with the expected value, contrary to what was found in Li et al. (2024). This disagreement is caused by the uncertainties in the determination of ϵ_g, which are smaller than the ones we obtained in this study. The difference between the two can be attributed to the way the model error was taken into account. For the stars we studied, we obtained a mean model error of 24 nHz for g modes and 97 nHz for p-modes. Comparatively, Li et al. (2024) used 7 and 14 nHz respectively, which thus seems underestimated.

Fig. 9. Comparison of ϵg between Li et al. (2024) and this study. Blue dots represent stars with no detected magnetic field and red dots represent stars with magnetic detections. The grey-shaded area represents the empirical value of ϵg at 0.28 ± 0.08 (Mosser et al. 2015). The dashed black line represents the 1:1 line.

Among the non-detections, we can identify two clumps at ϵ_g, L24 ≈ 0.15 and ϵ_g, L24 ≈ 0.38. These stars make up the tails of the distribution of ‘normal’ stars in Li et al. (2024), which is why we recover values compatible with the expected value in this study. We also see that almost no magnetic detection lies among these clumps, further reinforcing the idea that these stars should be considered as normal. The rest of the non-detections are more spread out and have larger error bars. This can be linked to difficulties in fitting the observed frequencies, where we recovered either multiple solutions or flat posteriors, which drives up the width of error bars.

5.2. Rotation rates

As a result of our fitting procedure, we also obtained estimates of the core and envelope rotation rates for stars with l = 1 triplets or doublets. To assess the robustness of our method, we show the measured core rotation rates as a function of previous measurements in Fig. 10. Most measurements fall on the 1:1 line, indicating that we generally recover the same rotation rates as in Li et al. (2024). We also noticed that a small fraction of measurements fall on the 1:2 line. In these stars, we detected a previously undetected component, effectively doubling the rotation rate. Not represented on Fig. 10, are cases where rotation was discovered in this study (Li et al. 2024 had detected a singlet, but we detect a doublet or a triplet), concerning 16 stars. We can also see that magnetic stars do not follow a particular tendency on this figure. We were also able to compare core rotation rates for seven stars in common with Hatt et al. (2024). In both studies, no magnetic field was detected in these stars. Our measurements are statistically compatible with those of Hatt et al. (2024) for these stars. Concerning envelope rotation, our measurements are compatible with those of Li et al. (2024), although our error bars are larger, for reasons explained above.

Fig. 10. Comparison of core rotation measured in Li et al. (2024) and in this study. The blue dots represent measurements for non-magnetic stars, and the orange dots represent magnetic stars with 1σ error bars. The dashed black line represents the 1:1 line, and the dashed-dotted line represents the 1:2 line.

6. Discussion

6.1. Mass distribution of magnetic red giants

The new detections in this study have an average mass of 1.375 ± 0.196 M_s. With the addition of these 23 new stars, the total number of detections is 71 magnetic red giants. This number allowed us to begin exploring the distributions of different stellar parameters (e.g. mass and rotation) for these stars and see how they compare to those of red giants in general. The distributions for all well-characterised red giants are taken from the catalogue of Li et al. (2024). Since all the studies that searched for magnetic fields in red giants so far have focused on low-luminosity stars, we decided to select only the stars from Li et al. (2024) for which Δν >  10 μHz, for consistency reasons.

By plotting the histograms of the mass distribution in Fig. 11, we find that the two distributions seem quite similar. We can see that the magnetic distribution follows the (much more populous) complete distribution. Computing a Kolmogorov-Smirnov test shows that the two distributions have a 50% chance of following the same probability law. If this tendency holds up with the growing number of magnetic detections, it shows that magnetic red giants are not peculiar from a stellar mass standpoint. The distribution of magnetic stars and non-magnetic stars being similar could also imply that more possible detections still lie in the non-detection distribution.

Fig. 11. Top: Normalised histogram of stellar mass comprised of data from Hatt et al. (2024), Li et al. (2024), Deheuvels et al. (2023), and this study. The distribution for magnetic detections is shown in orange, and the entire distribution of red giants is represented in blue. The distribution for dipole depressed modes in Stello et al. (2016) is shown on green. Bottom: Normalised histograms of the core rotation rate νR, g comprised of data from the same previous studies. The colours are the same as above.

We also represented the distribution for stars showing suppressed dipole modes in Stello et al. (2016), which are suspected of harbouring high-intensity magnetic fields. The stars represented here all have Δν >  10 μHz to be coherent with our sample in terms of evolutionary state. We can clearly see the difference between the two distributions, the suppressed-dipole-mode stars having much higher masses. One possible explanation for this discrepancy could be that the two populations are indeed different, and that the fields causing suppressed dipole modes are different from the fields we detect. This could also be caused by observational biases. Magnetic fields are directly detectable using seismology if they are strong enough to produce a detectable signature but not stronger than the critical field B_c. This gives a range in magnetic field strength for which fields are detectable. If this range becomes narrower with increasing stellar mass, this could explain the relative lack of more massive red giants with detected magnetic fields. This question will be addressed in a later study. The detection threshold decreases as stars ascend the RGB, so the prospects of detecting fields in more massive stars might improve when studying more luminous red giants.

6.2. Rotation rates

For the magnetic stars in this study, the core rotation can be measured in cases where either two or three mode components per triplet are visible, as reminded in Sect. 5.2. Here, rotation was measured in 12 magnetic stars, which can be added to the ones from Li et al. (2022, 2023, 13 stars) and the ones from Hatt et al. (2024, 24 stars). Here, we compare the distribution of core rotation rates for magnetic giants with the one made up of the remaining non-detections from this study, Li et al. (2024) and Hatt et al. (2024). However the size of the first distribution is fairly small, so it should not be over-interpreted.

In the lower panel of Fig. 11, we can see that the distribution of core rotation for magnetic stars is very similar to the distribution of core rotation for red giants in general, with a peak around 0.8 μHz. This peak contains mostly stars from Hatt et al. (2024) which can also be seen in the red giant global distribution. A Kolmogorov-Smirnov test using these two distributions indicates that they have high odds of following the same statistical law. It is quite striking that the core rotation rate does not seem to be different for stars that are found to host core magnetic fields. This could of course mean that the detected fields are not responsible for the efficient transport of angular momentum occurring inside red giants. However, this could also mean that magnetic fields are ubiquitous in red giant cores and that we have not detected them in most cases, either because they produce a seismic signature that is below the detection threshold or because the detection methods applied so far are ill-suited (for instance for strong non-axisymmetric fields).

6.3. Impact of evolution

We can also explore the distribution of magnetic stars with respect to evolution. In Fig. 12, we plot the magnetic field strength as a function of the mixed mode density 𝒩 = Δν/(ΔΠ₁ ν_max²), which serves as an indicator for evolution (Gehan et al. 2018). We confirm previous findings that the intensity of detected fields decreases with evolution (Deheuvels et al. 2023; Li et al. 2023).

Fig. 12. Variation of magnetic field intensity with mode density, used as a proxy for stellar evolution. Four different studies are represented in the plot. The arrows pointing upwards represent lower limits to the field intensity while red, while vertical error bars represent the lower and upper limits of the magnetic fields for stars that only show two l = 1 mode components. The grey dotted line at 𝒩 = 9.87 represents the upper limit in evolution reachable in this study (calculated using Δν = 10 μHz and the relevant values for ΔΠ1 and νmax obtained from the degeneracy sequence and seismic scaling laws).

This decreasing trend can seem counter-intuitive. Indeed, if the magnetic flux is conserved, the core contraction as the star ages should cause the field strength to increase. This may be explained through two different phenomena. First, as suggested by Deheuvels et al. (2023), the field intensity could be increasing with age, but would eventually reach the critical field strength B_c (which decreases with evolution), making the observation of stronger fields impossible. This is supported by the fact that our detections decrease similarly as B_c. Another potential explanation can be proposed, assuming that the detected fields are generated in main-sequence convective cores and survive on the RGB after relaxing into a stable configuration. During this phase, hydrogen is consumed in a shell (the HBS), which extends outwards as the star ages, leaving an inert helium core behind. Eventually, the HBS passes the limit of the main-sequence convective core. Since the seismic signature of magnetic field is mainly sensitive to the HBS, we expect the measured field strength to decrease with evolution when the HBS has left the previously convective core.

We can also see that the different studies shown cover different parts of the graph, such as Hatt et al. (2024) covering low-intensity, low-luminosity red giants, or Deheuvels et al. (2023) covering intense fields. This patchwork of detections exhibits the observational biases of these different studies. We can also notice a gap between the near-critical fields and the rest of the measurements. For now, no explanation was given for this gap. This representation also exhibits detection thresholds for field intensity. Indeed, Eq. (9) shows that ν_mag varies as (ν_max)⁻³. This means that for low-luminosity red giants with high ν_max, the magnetic intensity must be very high to observe significant shifts. In a similar way, evolved red giants with low ν_max must have field intensities below the critical field to be observable.

6.4. Asymmetry parameter a

In this study, we identified six magnetic red giants for which three components per l = 1 are visible, so that we were able to measure the asymmetry parameter a. Adding these measurements to those of previous studies, we now have 44 stars with measured asymmetries (Li et al. 2023; Hatt et al. 2024). Figure 13 shows the relation between the parameter a and the evolution proxy 𝒩. We count 7 stars whose asymmetry is compatible with zero, 27 stars with positive asymmetry and 10 with negative asymmetry. The interpretation of tendencies in Fig. 13 is complicated owing to the observational biases of the corresponding studies. Li et al. (2023) specifically searched for non-zero asymmetries. The study of Hatt et al. (2024) can detect fields with vanishing asymmetry parameter, but they are limited to stars with 𝒩 ≲ 8.

Fig. 13. Relation between parameter a and 𝒩.

From Fig. 13, one can notice that there is a larger number of stars with positive asymmetry than ones with negative asymmetry. This implies a higher probability for stars to form fields concentrated at the poles rather than ones concentrated at the equators. Also, we clearly see here that a sizeable portion of the detected fields seem incompatible with pure dipolar fields. Indeed, for dipolar fields a ranges from −0.2 and 0.4, and we see that a large proportion of stars are still outside these bounds. Thus, considering only dipolar fields as representative configurations for internal magnetic fields is not enough, and different configurations should be considered for accurate representations.

6.5. Potential non-axisymmetric effects

We assessed our assumption on the axisymmetry of the field by computing the parameter b = 2ν_B/ν_R, g. Li et al. (2022) showed that non-axisymmetric effects can be ignored if b <  1. We computed the value of b for all our magnetic stars at three different points in the spectrum: for the lowest-frequency detected mode, at ν_max, and for the highest-frequency detected mode. Since ν_B ∝ ν⁻³, b also follows this trend. This means that we are more likely to observe non-axisymmetric effects at lower frequencies. We also recall that in cases where only two components are visible, the value of ν_mag is bounded between $(2/3){\stackrel{\sim }{\nu }}_{\mathrm{mag}}$ and $(4/3){\stackrel{\sim }{\nu }}_{\mathrm{mag}}$, and hence the same intervals hold for b. In Table 3 these intervals are represented with brackets, along with the rest of the magnetic detections.

We notice that at least eight stars of this sample exhibit conditions in which non-axisymmetric effects could be observed. These effects would manifest by additional components, up to nine in total for dipole modes (see Loi 2021; Li et al. 2022, supplementary information). In these stars, no additional component appears at low frequency. For now, we thus do not detect signature of non-axisymmetric effects for the magnetic field, even though these stars could show some.

6.6. Mode suppression

In this study, we found evidence of a magnetic field in two stars showing signs of suppressed dipole modes The first star, KIC 12168062 (see Sect. 4.4), exhibits a strong visibility gradient of dipolar mixed modes. Indeed, the visibility of dipole-modes (measured relatively to the visibility of radial modes) goes from an average of 0.14 below ν_max to 1.49 above ν_max. This behaviour is shown in the top panel of Fig. 14, where we can clearly see that dipolar modes are almost absent from the first two large separations, but suddenly appear as they normally would in a star without suppressed dipoles. The second star, KIC 8687248, was already shown to have suppressed dipole modes in Mosser et al. (2017). In the bottom panel of Fig. 14 we can see how, again, the visibility of dipole modes is lower than that of radial modes, with a mean ratio of the dipolar to radial mode energy of about 0.51. We can also see that only the most p-dominated modes are visible in the spectrum. Despite this fact, we clearly found that these modes follow a mixed mode pattern disturbed by a core magnetic field. These two stars also share the two highest B/B_c ratios of our detections, which further links the effect of the critical field to the suppression of modes. This also means that the perturbative approach we used on these two stars likely overestimates the field intensity, and a non-perturbative would be required to better treat these stars (Rui et al. 2024).

Fig. 14. Top: Spectrum of KIC 12168062. Bottom: Spectrum of KIC 8687248. Black areas denote radial modes, while red areas denote quadruple modes, and orange signals indicate the presence of a visible l = 3 mode.

7. Conclusion

In this study we have investigated the use of the gravity offset ϵ_g as a means to identify magnetic red giants. We worked on 218 stars considered as abnormal in terms of ϵ_g in the catalogue of Li et al. (2024). We developed a method potentially applicable to all red giants to detect and measure magnetic fields. Since magnetic fields have stronger effects on low-frequency modes, it is crucial to reliably detect these modes. We applied significance tests that we designed in order to optimise the detection probability of the modes while guaranteeing a low false-alarm probability. We then adjusted an asymptotic expression of mixed modes to the frequencies of the detected modes, taking magnetic fields and rotation into account. For this purpose, we used a Bayesian inference, which allowed us to explore the posterior distributions for every parameter considered in the model and to assess the odds of detecting a magnetic field. We then used a grid of the stellar model to link the magnetic frequency perturbation to the magnetic field intensity.

Among the 218 stars selected, 23 exhibit magnetic fields with a high degree of confidence in the detection. We found intensities for the radial component of the magnetic field ranging from 34 kG to 262 kG, in line with measurements from previous studies. For these stars, we found values of ϵ_g that are now in agreement with the expected value for red giants, validating the suggestion of Li et al. (2022). For seven stars, we placed constraints on the topology of the field through the asymmetry parameter a. We were also able to measure the internal rotation for 12 of those newly detected magnetic giants.

Adding our new measurements to the previous studies measuring core magnetic fields (Li et al. 2022, 2023; Deheuvels et al. 2023; Hatt et al. 2024), there are now 71 detections of magnetic red giants, which allowed us to start investigating the distributions of the stellar parameters for magnetic giants. We confirm the decreasing trend of the radial field component with the evolution of the star, which was already noticed in previous studies. This trend is likely due to observational biases, but it could also result from the progression of the HBS with evolution, as we propose in this study. We find that the mass distribution of the stars in which we detect magnetic fields is so far indistinguishable from the mass distribution of the complete catalogue of seismic red giants. In contrast, the masses of the red giants with direct field detection are significantly lower than those of stars showing suppressed dipole modes at the same evolutionary stage (Stello et al. 2016). This discrepancy could mean that the two populations are truly distinct and that the fields producing mixed-mode suppression differ from the ones we detect directly. It could also be that the range of detectable field strengths (between the detection threshold and the critical field) becomes narrower as stellar mass increases, leading to fewer detections of more massive stars. We also find that the magnetic red giants have a distribution of core rotation rates that is strikingly similar to that of other red giants. This shows that the transport of angular momentum does not act in a different way for stars with detected magnetic fields compared to the rest of the red giants. There are two possible interpretations regarding the angular momentum redistribution in these stars: Either magnetic fields do not have a predominant role in the transport of angular momentum in red giants or the rest of the red giants also host core magnetic fields that have not yet been detected.

This study also brings new insight on field geometry. By studying the distribution of the asymmetry parameter a, we showed that positive asymmetries are more frequent than negative ones. Also, many detected fields have measured values of a that are incompatible with a dipolar geometry for B_r so that the assumption of dipolar fields, often made for simplicity, is not relevant. The effects of non-axisymmetry of the fields were also investigated. We showed that at least eight stars exhibit conditions favourable for the detection of non-axisymmetric effects. However, no such effects were discovered in these stars.

To conclude, this study represents an effort to bring new detections of core magnetic fields in which we have a high degree of confidence. These detections are another step towards building an exhaustive catalogue of magnetic red giants. In this work, we only studied stars with Δν <  10 μHz because more luminous stars have more complex spectra. In a forthcoming study, these stars will be explored using the same method described here. This new method, from the identification of mode frequencies to the measure of magnetic field intensity, is in principle applicable to any oscillation spectrum of red giants and can already be used on TESS spectra. In the near future, we aim to automate and use this method on the totality of the catalogue of Li et al. (2024) to look for core magnetic fields. Eventually, this method will be used on PLATO (Rauer et al. 2025) spectra also, further expanding the number of detected magnetic red giants.

Acknowledgments

This work has been supported by CANES, 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. Software: Python (Van Rossum & Drake 2009), numpy (Oliphant 2006; Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), corner (Foreman-Mackey 2016), astropy (Astropy Collaboration 2013, 2018, 2022).

References

Appendix A: Mode identification method

In this section we detail the process used to assign a unique frequency and width guess to each significant mode component in the spectrum. As a first step, we apply the variable threshold presented in Sect. 3 to the smoothed spectrum. This creates groups of consecutive points containing the peaks of the spectrum that exceed the threshold. Figure A.1 shows an example of such groups on a rotational triplet of KIC 4660930. There, we see that the m = 0 and m = −1 components are contained within one single group.

A second step consists of reducing the number of points in each group. This is done by selecting the local maxima in the group two consecutive times. Local maxima are defined as points where the two neighbouring points are lower than the central one. Selecting these maxima two consecutive times helps reduce the effects of noise in the spectrum, while conserving split modes that are found in one group. The result of this step is shown in Fig. A.2. We clearly see that this reduces the number of identifications to a manageable number, while preserving both m = 0 and m = 1 components.

The third step tries to eliminate cases where consecutive selected peaks belong to the same mode component. This is done in two sub-steps. First, the distance between each peak within a group is checked. Then the distance between each group is checked. If this distance is smaller than three times the expected width of the dipolar mode at its frequency, then the peak with the smallest amplitude is removed. We found that using a minimal distance of three times the mode width works best to separate split modes. This step is effective in about 90% of cases, such as in Fig. A.3. There, the widths of the peaks do not overlap, thus preserving each component of the triplet, which is the best case scenario. In Fig. A.4 and A.5, we show typical caveats that arise in the remaining 10% of cases. In the first case, shown in Fig. A.4, we see another triplet of KIC 4660930 (a p-dominated mode). Again, the m = 0 and m = −1 components are within the same group, with an additional detection at around 156.65 μHz. However, the width of the modes overlaps, which means that only the guess with the highest intensity is retained. The second case, shown in Fig. A.5, corresponds to a high-S/N p-dominated mode. We see how noise in the left wing of the mode has caused an unwanted secondary detection. This detection is not removed because the red areas do not overlap. In most problematic cases, the modes involved are p-dominated modes because they have larger mode widths. These types of issues have to be addressed manually.

These few select cases show the limits of such a simple approach. In the global case, such cases are few and far between (around one or two modes per star) but still require additional vigilance and a thorough manual examination of the identifications after the fact. In a future study, we aim to reduce this kind of error by refining the detection algorithm.

Fig. A.1. Portion of the PSD of KIC 4660930 showing a mixed l = 1 mode. The raw spectrum is plotted in blue, the variable smoothed spectrum is plotted in orange, and the dash-dot grey line represents the significance threshold of the smoothed spectrum. The dot markers represent the points of the spectrum that exceed the threshold, and the colours represent the group to which they belong.

Fig. A.2. Same as Fig. A.1, but after reducing the number of peaks in each group (see text).

Fig. A.3. Same as Fig. A.1. The transparent red zones delimit an area of three times the l = 1 mode width, estimated using ζ. The red stem plots mark the points that have been selected as guesses for mode frequency.

Fig. A.4. Same as Fig. A.3, but for another p-dominated triplet in the PSD of KIC 4660930.

Fig. A.5. Same as Fig. A.3, but for a p-dominated singlet in the PSD of KIC 11087371.

Appendix B: Computing p_νmag(0)

As mentioned in Sect. 4.3, we need to compute the value of p_νmag(0) from a Monte Carlo Markov chain sampling. We built a density histogram P_h of the parameter ν_mag with a bin size h based on Scott’s rule (Scott 1979):

$$\begin{array}{c}\hfill h=3.49\frac{\sigma }{N^{1/3}},\end{array}$$(B.1)

where σ is the standard deviation of ν_mag and N is the number of samples. This rule is generally a good choice for high count histograms, which is the case here. The value of p(0) – we dropped the index ν_mag for visual clarity – can be extrapolated from the first bin of the histogram, P_h[1]. Since the first bin is the interval [0, h], P_h[1] is related to p through the relation

$$\begin{array}{c}\hfill P_h[1]\approx \frac{1}{h}{\displaystyle {\int }_0^h}p({\nu }_{\mathrm{mag}})\mathrm{d}{\nu }_{\mathrm{mag}}.\end{array}$$(B.2)

We expanded p around 0 using a Taylor series:

$$\begin{array}{c}\hfill p(x)\approx p(0)+xp\prime (0)+\frac{x^2}{2}p″(0)+O(x^3),\end{array}$$(B.3)

where the script ′ indicates the derivative. By combining Eqs. B.2 and B.3, and computing histograms with bin sizes h/2 and 2h, we obtain the following system of equations:

$$\begin{array}{c}\hfill \left(\begin{array}{c}{P}_h[1]\\ P_{h/2}[1]\\ P_{2h}[1]\end{array}\right)=\left[\begin{array}{ccc}1& 1& 1\\ 1& 1/2& 1/4\\ 1& 2& 4\end{array}\right]·\left(\begin{array}{c}p(0)\\ hp\prime (0)/2\\ h^{2}p″(0)/6\end{array}\right).\end{array}$$(B.4)

Inverting this system yields the approximated value of p(0).

All Tables

All Figures

Fig. 1. Échelle diagrams of mock g-mode spectra perturbed by rotation and magnetic field. Left: Spectrum fitted with an asymptotic expression accounting for magnetic and rotational perturbations, folded by the true ΔΠ1. Right: Spectrum fitted with an asymptotic expression for rotation only, folded by the measured ΔΠ1′. The red lines represent the slope of the m = 0 component, which is purely vertical on the right. The blue area represents the area where m = 0 modes are expected to be found using ϵg = 0.28 ± 0.08.

In the text

	Fig. 2. Relationship between ϵg and Δν in the Li et al. (2024) catalogue. Red dots represent stars selected for this study, in contrast to the grey dots. Solid blue lines represent the delimitation of the bulk distribution at ±σ, and the dashed blue line is made from the median of every 118 stars slices.
In the text

	Fig. 3. Probability of detection of a mode versus the boxcar width, in units of the mode width. The colour of the lines indicates the mode energy, defined as E ∝ Γh, where Γ is the mode width and h is mode height. The vertical dash-dotted line corresponds to a boxcar width of 3Γ.
In the text

	Fig. 4. Spectrum of KIC 9999911. Top: Unprocessed spectrum. Bottom: Pre-processed spectrum with variable width smoothing and radial and quadrupolar modes cut out. Red arrows indicate the position of detected modes using the red dashed-dot threshold.
In the text

	Fig. 5. Corner plot of KIC 4660930 reduced to only show g-mode and magnetic parameters.
In the text

	Fig. 6. Stretched échelle diagrams of KIC 4660930, KIC 6222530, and KIC 3634488. Blue dots represent detected frequency, and red circles represent the maximum a posterior model of the determination.
In the text

	Fig. 7. Corner plot of KIC 6222530 reduced to show only g-mode parameters and the magnetic shift νmag.
In the text

	Fig. 8. Corner plot of KIC 3634488 reduced to show only g-mode and magnetic parameters.
In the text

	Fig. 9. Comparison of ϵg between Li et al. (2024) and this study. Blue dots represent stars with no detected magnetic field and red dots represent stars with magnetic detections. The grey-shaded area represents the empirical value of ϵg at 0.28 ± 0.08 (Mosser et al. 2015). The dashed black line represents the 1:1 line.
In the text

	Fig. 10. Comparison of core rotation measured in Li et al. (2024) and in this study. The blue dots represent measurements for non-magnetic stars, and the orange dots represent magnetic stars with 1σ error bars. The dashed black line represents the 1:1 line, and the dashed-dotted line represents the 1:2 line.
In the text

Fig. 11. Top: Normalised histogram of stellar mass comprised of data from Hatt et al. (2024), Li et al. (2024), Deheuvels et al. (2023), and this study. The distribution for magnetic detections is shown in orange, and the entire distribution of red giants is represented in blue. The distribution for dipole depressed modes in Stello et al. (2016) is shown on green. Bottom: Normalised histograms of the core rotation rate νR, g comprised of data from the same previous studies. The colours are the same as above.

In the text

Fig. 12. Variation of magnetic field intensity with mode density, used as a proxy for stellar evolution. Four different studies are represented in the plot. The arrows pointing upwards represent lower limits to the field intensity while red, while vertical error bars represent the lower and upper limits of the magnetic fields for stars that only show two l = 1 mode components. The grey dotted line at 𝒩 = 9.87 represents the upper limit in evolution reachable in this study (calculated using Δν = 10 μHz and the relevant values for ΔΠ1 and νmax obtained from the degeneracy sequence and seismic scaling laws).

In the text

	Fig. 13. Relation between parameter a and 𝒩.
In the text

	Fig. 14. Top: Spectrum of KIC 12168062. Bottom: Spectrum of KIC 8687248. Black areas denote radial modes, while red areas denote quadruple modes, and orange signals indicate the presence of a visible l = 3 mode.
In the text

	Fig. A.1. Portion of the PSD of KIC 4660930 showing a mixed l = 1 mode. The raw spectrum is plotted in blue, the variable smoothed spectrum is plotted in orange, and the dash-dot grey line represents the significance threshold of the smoothed spectrum. The dot markers represent the points of the spectrum that exceed the threshold, and the colours represent the group to which they belong.
In the text

	Fig. A.2. Same as Fig. A.1, but after reducing the number of peaks in each group (see text).
In the text

	Fig. A.3. Same as Fig. A.1. The transparent red zones delimit an area of three times the l = 1 mode width, estimated using ζ. The red stem plots mark the points that have been selected as guesses for mode frequency.
In the text

	Fig. A.4. Same as Fig. A.3, but for another p-dominated triplet in the PSD of KIC 4660930.
In the text

	Fig. A.5. Same as Fig. A.3, but for a p-dominated singlet in the PSD of KIC 11087371.
In the text