© The Authors 2025

1. Introduction

In recent years, high-quality data from space-based missions such as Kepler (Borucki et al. 2010) have led to major advances in the research field of asteroseismology. Of particular interest are the oscillations of red giants (i.e., evolved low-mass stars), as their multipolar modes are sensitive to both the stellar envelope and the core. This is due to the mixed character of their oscillation modes, which means that they behave like a gravity mode in the stellar core and like a pressure mode in the outer layers (for a review, see Hekker & Christensen-Dalsgaard 2017). Analyzing the mixed modes can provide information about the interior of red giants, for example about the nuclear burning reaction and the evolutionary stage (Bedding et al. 2011; Mosser et al. 2011), the core rotation rate (e.g., Beck et al. 2012, 2014; Mosser et al. 2012a; Deheuvels et al. 2012, 2014; Gehan et al. 2018; Li et al. 2024), as well as the strength and topology of the internal magnetic field (Li et al. 2022, 2023; Deheuvels et al. 2023; Hatt et al. 2024).

Among the more than 16 000 red giants identified by Kepler (Yu et al. 2018), Mosser et al. (2012b) have identified a population of stars that have dipole modes with unusually low amplitudes, while their radial modes appear to be normal. Since the mode amplitude is a measure for the mode energy, this indicates that the stochastic mode excitation by turbulent convection works as usual in this particular group of red giants, and suggests an additional source of damping confined to the inner regions where the radial modes are least sensitive (e.g., García et al. 2014; Fuller et al. 2015; Coppée et al. 2024). The dipole modes are therefore suppressed. The red giants with suppressed dipole modes are often referred to as the “suppressed dipole mode stars”.

Fuller et al. (2015) proposed that the additional damping of the suppressed dipole mode stars is caused by a strong core magnetic field. Such a field could trap part of the total energy of the dipole modes within the core, thus removing it from the observable oscillations (the “magnetic greenhouse effect”). In the theoretical framework of Fuller et al. (2015), a dipole mode is suppressed if the frequency of the mode is lower than a characteristic frequency that depends on the strength of the internal magnetic field. For a dipole mode with a given frequency, this results in a step-like onset of the suppression as a function of the magnetic field strength. The strength of the internal magnetic field is thus decisive for whether a mode is suppressed or not. Red giants with normal dipole mode amplitudes are expected to have no or a weak internal magnetic field, while the suppressed dipole mode stars are expected to contain strong magnetic fields. Modes with a higher spherical degree can also be suppressed, although the onset of the suppression is happening at a different field strength and the effect of the suppression is less pronounced (Fuller et al. 2015; Stello et al. 2016a,b; Cantiello et al. 2016).

A key assumption of the mode suppression as described by Fuller et al. (2015) is that all multipole mode energy that leaks into the gravity mode cavity is lost. Crucially, this means that suppressed modes cannot be mixed, as this requires a gravity mode-like behavior in the inner regions of the star. In their prescription, the suppressed modes are pure pressure modes and their amplitudes are determined by the amount of energy that does not leak into the core (Fuller et al. 2015).

This scenario is challenged by Mosser et al. (2017), who show that a significant fraction of suppressed dipole mode stars have dipole modes with a mixed character. This implies that the damping mechanism does not fully suppress the energy leaking into the gravity mode cavity. Instead, a portion of the energy must be able to leave the gravity mode cavity and contribute to the observable oscillations. Therefore, the transition between no suppression and full suppression (i.e., the suppression as described by Fuller et al. 2015) cannot be a step function. A more gradual transition is needed to connect these two regimes. In fact, the suppressed dipole mode stars exhibit a variety of power spectral density (PSD) morphologies, of which the one showing dipole modes of mixed character is one (Coppée et al. in prep.).

The true nature of the suppression mechanism is still under investigation. Magnetism is often treated as a small perturbation to the restoring forces of the oscillations (e.g., Prat et al. 2019, 2020; Mathis et al. 2021; Bugnet et al. 2021; Loi 2021; Li et al. 2022; Bugnet 2022; Bhattacharya et al. 2024; Das et al. 2024). However, this method can only be used for weak magnetic fields and is not applicable if the magnetic field is strong enough to cause mode suppression (Rui et al. 2024). Another common assumption is that of a purely horizontal magnetic field (e.g., Mathis & de Brye 2011; Dhouib et al. 2022). In reality, on the other hand, the radial component of the field is expected to dominate the interaction with the gravity waves (e.g., Fuller et al. 2015; Mathis et al. 2021; Li et al. 2022).

Despite these difficulties, some authors have tackled the phenomenon of the magnetic mode suppression under different assumptions. Lecoanet et al. (2017) investigated the interaction of internal gravity waves with a strong magnetic field and find perfect refraction into upward propagating slow Alfvén waves if the magnetic field strength exceeds a certain threshold. These Alfvén waves are likely to dissipate as they travel outward. Loi & Papaloizou (2017) propose that the magneto-gravity modes might be dissipated by resonant interactions with standing torsional Alfvén modes. However, they assume that the internal magnetic field does not affect the structure of the spheroidal component of the modes, which is not the case for high field strengths (Loi & Papaloizou 2018). In a follow-up study, Loi & Papaloizou (2018) suggest that both the toroidal and spheroidal component of the modes are affected by the magnetic field, implying that the period spacing of the suppressed dipole modes stars should be larger than that of RGB stars with normal multipole mode amplitudes. Loi (2020a) numerically studied the onset of the suppression by applying the Hamiltonian ray tracing method to magneto-gravity waves in the core and find that there are ray trajectories that retain a gravity-dominated character even at high magnetic field strengths. Finally, Rui & Fuller (2023) analytically characterized the eigenfunctions of magneto-gravity modes in the presence of a dipole field and find that they do become evanescent in regions with a sufficiently strong magnetic field. However, they cannot exclude that a portion of the wave energy might escape the gravity mode cavity in a narrow equatorial band. Apart from red-giant branch (RGB) stars, the study of Lecoanet et al. (2022) suggests that the suppression of gravity modes by interaction with the magnetic field could also occur in the radiative layers of more massive main sequence (MS) stars. A universal solution for the global oscillation modes in the presence of an internal magnetic field has yet to be found.

In this work, we investigate the transition from no mode suppression to full mode suppression using the method of Loi (2020a). We updated their formalism to reduce the uncertainties and calculated more sets of rays for each oscillation mode to achieve a higher resolution of the transition as a function of the magnetic field strength. In addition, we considered different stellar masses, metallicities, and topologies of the internal magnetic field. We subsequently used these insights into the onset of the suppression mechanism to study the evolution of the mode suppression as the star climbs the RGB.

This paper is the first in a series of publications in which we will theoretically investigate the properties of the multipole modes of stars exhibiting magnetic mode suppression where part of the mode energy entering the gravity mode cavity can escape it. Here, we present a simple expression that describes the mode suppression as a function of the magnetic field strength as predicted by our Hamiltonian ray tracing analysis. In the future, we will relate the different morphologies observed in the PSD of suppressed dipole modes stars (Coppée et al., in prep.) to physical quantities and investigate the evolution of the multipole mode amplitudes on the RGB, as well as during the helium core-burning phase.

This article is structured as follows. In Sect. 2, we describe the simulation setup of the ray tracing analysis. In Sect. 3, we present the results of our parameter study and provide a simple expression for the magnetic mode suppression. We also study the evolution of the mode suppression along the RGB. Finally, we discuss our findings in Sect. 4 and summarize in Sect. 5.

2. Methods

2.1. Stellar models

We used version r23.05.1 of the publicly available stellar evolution code MESA (Modules for Experiments in Stellar Evolution; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) to compute four evolutionary tracks with varying stellar mass (M_* = 1 M_⊙ or 2 M_⊙) and metallicity (Z_init = 0.005 or 0.02). We chose ten numerical models on these tracks to serve as a background for the Hamiltonian ray tracing analysis (see Sect. 2.4). Table 1 lists the stellar mass and metallicity of all chosen models, as well as the corresponding stellar radius R_*, the central density ρ_cen, and the dynamical angular frequency ${\omega }_{\mathrm{dyn}}=\sqrt{G{M}_{\ast }/R_{\ast }^3}$ with G being the gravitational constant. We also show the critical magnetic field strength (Eq. (6)) of a dipole mode with the angular frequency ω = 10 ω_dyn in the stellar center (B_crit, cen) and in the hydrogen-burning shell (B_crit, H). The models 2A, 2B, and 2C have the same stellar mass, metallicity, and radius as the models A, B, and C of Loi (2020a) respectively and can be regarded as reproductions of these models. They are not identical to the models of Loi (2020a), because we used a different version of MESA and made small adjustments to the MESA inlist (see Appendix C). Nevertheless, they yield similar results for the ray tracing analysis. Figure 1 shows a Hertzsprung–Russell diagram of the four evolutionary tracks, in which the position of the models is indicated.

Fig. 1. Hertzsprung-Russell diagram of the four evolutionary tracks with different stellar masses and metallicities described in Sect. 2. The evolution of the MS is shown as a dotted line and the evolution after the MS as a solid line. We do not show the full evolution after the MS. The locations of the stellar models that served as the background for the ray tracing analysis are indicated by the various symbols. Star symbols indicate that the central field strength Bcen, the frequency ω, and the spherical degree l were varied during the ray tracing analysis of a particular model. Full dots indicate that Bcen and ω were varied and circles indicate that only Bcen was varied.

Loi (2020a) report some coarseness in their MESA models with a resolution similar to the grid scale. This was particularly evident in the profile of the buoyancy frequency N. It was likely caused by the receding convective core, which leaves small discontinuities in the profile of the chemical composition due to the finite timestep. To remove this coarseness, Loi (2020a) smoothed the profiles of the relevant output quantities and interpolated them onto a linear grid with a lower resolution than the grid of the MESA models. Although this post-processing removed the grid-scale coarseness, the resulting profiles had a relatively small number of grid points in the core region.

Instead of performing a similar smoothing procedure, we enabled the set_min_D_mix option in our MESA inlist and set the min_D_mix parameter to 0.5. This led to low artificial diffusive mixing everywhere in the star where the mixing was not already higher than this minimum value. The result is a smooth profile without coarsening, while retaining the original resolution of the MESA grid. The additional mixing was not efficient enough to significantly alter the evolution of the star. In Appendix A.1, we show the profile of the buoyancy frequency of the models 2A, 2B, and 2C.

2.2. Internal magnetic field

Stars with a mass M_* ≳ 1.1 M_⊙ are assumed to have a convective core while they are on the MS. The convective motion can give rise to a dynamo, which can generate a magnetic field. When the star evolves away from the MS and becomes a subgiant and later a red giant, the core is fully radiative and the potential dynamo can no longer operate. However, the magnetic field generated during the dynamo phase might stabilize into a large-scale magnetic field (e.g., Braithwaite 2008; Emeriau-Viard & Brun 2017; Villebrun et al. 2019). Since the magnetic diffusion timescale has been estimated to be small for magnetic fields contained in the radiative core (e.g., Cantiello et al. 2016), they could survive the post-MS evolution. The strength of these buried magnetic fields is expected to be relatively large on the RGB. The reason for this is that the radiative core harboring the magnetic fields contracts after the MS, which increases the field strength due to the conservation of magnetic flux (see Fuller et al. 2015; Cantiello et al. 2016; Bugnet et al. 2021).

Both analytical and numerical studies have shown that purely toroidal or poloidal magnetic fields are unstable (e.g., Tayler 1973; Markey & Tayler 1973; Braithwaite 2006, 2007). Instead, stable magnetic field configurations are expected to have both toriodal and poloidal components of significant strength (e.g., Tayler 1980; Braithwaite 2008; Akgün et al. 2013; Becerra et al. 2022) and to take the form of a twisted torus configuration (Braithwaite & Nordlund 2006; Braithwaite & Spruit 2017). A commonly used analytical prescription for an equilibrium axisymmetric twisted torus field was derived by Prendergast (1956) and later generalized by Duez & Mathis (2010) to the compressible case. This formalism is sometimes referred to as the “Prendergast field” and has several useful properties. First, the magnetic field can be easily constructed from analytic expressions that can be directly applied to stellar models. Second, the field resembles a dipole field, but also has a significant toroidal component. Third, the field strength does not diverge at the center of the star and has a finite maximum value. Finally, all components of the field vanish as the distance from the center approaches a pre-selected maximum radius. Note that although Prendergast fields satisfy the stability criterion established by various studies (e.g., Braithwaite 2009; Akgün et al. 2013; Becerra et al. 2022) in stably stratified stellar models, recent findings by Kaufman et al. (2022) indicate that these fields may not be stable throughout stellar evolution. For our purposes, this is not crucial because we demonstrate in Sect. 3.3 that our results are expected to be valid for any field topology with a dipole-like dependence on latitude and a radial component comparable or larger than the latitudinal and longitudinal component.

The Prendergast field is described in terms of the radial flux function Ψ(r) given by

$$\begin{array}{cc}& \mathrm{\Psi }(r)=\frac{\beta \lambda r}{j_1(\lambda r_{\mathrm{f}})}[f_{\lambda }(r,r_{\mathrm{f}}){\displaystyle {\int }_0^r}\rho {\xi }^3{j}_1(\lambda \xi )\mathrm{d}\xi \hfill \\ \hfill & +j_1(\lambda r){\displaystyle {\int }_r^{r_{\mathrm{f}}}}\rho {\xi }^3{f}_{\lambda }(\xi ,r_{\mathrm{f}})\mathrm{d}\xi ],\hfill \end{array}$$(1)

with f_λ(r₁, r₂)≡j₁(λr₂)y₁(λr₁)−j₁(λr₁)y₁(λr₂). Here, r is the distance to the center of the star, r_f is the maximum radial extend of the magnetic field, ρ is the density of the star, β is a scaling factor that controls the field strength, j₁ and y₁ are spherical Bessel functions, and λ can be determined by

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{r_{\mathrm{f}}}}\rho {\xi }^3{j}_1(\lambda \xi )\mathrm{d}\xi =0.\end{array}$$(2)

Finding values of λ that satisfy Eq. (2) is usually possible for subgiant and red giant models for a given r_f. In this work, we always chose the smallest possible value of λ for each stellar model. The value for r_f used for a model can be found in Table 2. The components of the magnetic field B = (B_r,B_θ,B_ϕ) are given by

$$\begin{array}{cc}\hfill B_{\mathrm{r}}& =\frac{2}{r^2}\mathrm{\Psi }(r)cos\theta ,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill B_{\theta }& =-\frac{1}{r}{\mathrm{\Psi }}^{\prime }(r)sin\theta ,\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill B_{\varphi }& =-\frac{\lambda }{r}\mathrm{\Psi }(r)sin\theta ,\hfill \end{array}$$(5)

which we used to insert a magnetic field into the radiative core of each of our MESA models before performing the ray tracing analysis.

The strength of the components of the magnetic field is scaled by adjusting β such that the absolute value of the component with the highest central field strength is equal to a preselected parameter B_cen at the center of the star. We refer to B_cen as the central field strength. It is specified in terms of the critical field strength

$$\begin{array}{c}\hfill B_{\mathrm{crit}}=\sqrt{{\mu }_0\rho }\frac{r}{N}\frac{{\omega }^2}{\sqrt{l(l+1)}},\end{array}$$(6)

where l is the spherical degree of the oscillation mode and μ₀ is the magnetic permeability. It marks the approximate field strength at which the magnetic field is strong enough to be dynamically relevant for the oscillations¹. In this work, we always express the field strength in the center B_cen as a factor that is multiplied by the central critical field strength B_crit, cen. The internal magnetic field strengths considered in this work are not high enough to significantly influence the evolution of the star, which means that the magnetic effects can be neglected in the MESA computations.

2.3. Hamilton’s equations

Hamiltonian ray tracing describes the propagation of a group velocity wave packet by assigning a position x = (r,θ,ϕ) and a wavevector k = (k_r,k_θ,k_ϕ) to the wave packet at any time t. The trajectory of the wave packet is described by Hamilton’s equations:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}={\mathrm{\nabla }}_{\mathbf{k}}H,\frac{\mathrm{d}\mathbf{k}}{\mathrm{d}t}=-\mathrm{\nabla }H,\end{array}$$(7)

where the Hamiltonian H is given by the dispersion relation ω (see Eq. (8)). In this work, we focus on the stellar core region, where the sound speed is so high that the dynamics of the oscillations are dominated by buoyancy and the interaction with a potential internal magnetic field (i.e., magneto-gravity waves). The dispersion relation of magneto-gravity waves is:

$$\begin{array}{c}\hfill {\omega }^2={\omega }_{\mathrm{A}}^2+{\omega }_{\mathrm{g}}^2={(\mathbf{k}·{\mathbf{V}}_{\mathrm{A}})}^2+{\left(\frac{k_{\perp }}{k}N\right)}^2.\end{array}$$(8)

Here, ω_A is the Alfvén wave frequency, ω_g is the gravity wave frequency, V_A = (V_Ar,V_Aθ,V_Aϕ) is the Alfvén velocity, $k=\sqrt{k_r^2+k_{\theta }^2+k_{\varphi }^2}$ is the wavenumber, and $k_{\perp }=\sqrt{k_{\theta }^2+k_{\varphi }^2}$ is the horizontal wavenumber.

Loi (2020a) formulated the Hamilton’s equations in the case of magneto-gravity waves in three-dimensional spherical polar coordinates, resulting in a system of six linear first-order ordinary differential equations. The derivation can be found in Loi (2020a). For completeness, we repeat the equations in Appendix A.2 (Eqs. (A.1)–(A.6)).

It is important to note that this approach relies on the Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) approximation in both the radial and horizontal direction, which requires that the wavelength is smaller than the length scale of the variation of the stellar structure. For non-magnetized red giant stars, this is typically a good assumption. However, the picture becomes more complicated when an internal magnetic field is present, since the horizontal wavenumber of the observable modes (i.e., modes with low l) has a length scale comparable to that of the horizontal variation of the magnetic field. The resulting limitations and our motivation for nonetheless using this approach are described in Sect. 4.4. Overall, we do not expect a strong impact on our main results.

2.4. Ray initialization and integration

We solved the magneto-gravity ray tracing equations (Eqs. (A.1)–(A.6)) by integrating them forward in time using a fourth-order Runge-Kutta scheme. The integration timestep was always set to $0.001{\omega }_{\mathrm{dyn}}^{-1}$. The total integration time t_int of each model is listed in Table 2. To initialize a ray, we selected a numerical stellar model which serves as the background. In addition, we chose the frequency of the wave packet ω, the spherical degree l, and the initial latitude θ₀. We also specified a polarization angle α with 0 ≤ α ≤ 2π, which defines the initial horizontal orientation of the ray.

Using these quantities, we determined the starting values for the six parameters to be integrated (i.e., r₀, θ₀, ϕ₀, k_r0, k_θ0, k_ϕ0) in the same way as Loi (2020a). The exact value of the initial radius r₀ does not matter as long as it is inside the radiative core and exterior to the maximum extend of the internal magnetic field r_f (Loi 2020a). Therefore, we did not test different values of r₀ for a given model. The value for r₀ used for each model is listed in Table 2. The initial latitude θ₀ can take values from 0 to π and the choice of the initial longitude ϕ₀ is irrelevant in the context of this work, since we only considered axisymmetric magnetic fields. We followed Loi (2020a) and set ϕ₀ = α. The starting values of the components of the wavevector are given by

$$\begin{array}{cc}\hfill k_{r0}& =\sqrt{k_0^2-k_{\perp 0}^2},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill k_{\theta 0}& =k_{\perp 0}sin\alpha ,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill k_{\varphi 0}& =k_{\perp 0}cos\alpha ,\hfill \end{array}$$(11)

with $k_{\perp 0}=\sqrt{l(l+1)/r_0}$ and k₀ = k_⊥0N(r₀)/ω.

We approximate the behavior of a global oscillation mode by calculating a set of 1200 individual rays and subsequently averaging over the initial latitude θ₀ and the polarization angle α. This is intended to mimic the arbitrariness of the propagation direction of waves traveling on a spherically symmetric background. Each set of rays is characterized by the parameters ω, l, and B_cen, and is assigned to a specific stellar model that serves as the background. The values used when varying frequency, the spherical degree, and the magnetic field strength are ω = (8,  10,  12) ω_dyn, l = (1,  2,  3), and B_cen = (0.1,  0.6,  0.7,  0.8,  1,  2,  3,  10) B_crit, cen, where the values of ω and l, which are used if they are not varied, are shown in bold. We chose θ₀ and α in such a way that the initial positions of the rays are evenly distributed on the surface of a sphere when θ₀ is considered as the latitude and α as the longitude.

The Runge-Kutta integration scheme does not inherently conserve ω (Loi 2020a)². However, since the Hamiltonian has no explicit time dependence, it is crucial that ω remains constant throughout the integration. While the integration routine conserves ω to machine precision when using infinitely differentiable background profiles (Loi 2020a), the finite grid of the MESA models introduces errors that accumulate and render the resulting trajectories of the rays unusable for further analysis. Loi (2020a) attempted to correct these errors by adjusting the value of k_r after each integration timestep such that ω remains constant. This was done by performing a Newton–Raphson iteration scheme for a maximum of five iterations to solve Eq. (A.7) (see Appendix A.3). Although this correction ensured the conservation of ω within 1% of the original value for most of their rays, the Newton-Raphson scheme did sometimes not converge, leading to uncertainties on the behavior of the global mode.

Here, we built on the method of Loi (2020a) and also adjusted k_r after each integration timestep to ensure that ω remains constant. However, instead of using a Newton–Raphson scheme, we directly calculated the roots of Eq. (A.7) after each timestep. We therefore faced no convergence problems. The roots can be complex. We then identified the root with the real part closest to our current estimate of k_r. The imaginary part of the wavenumber is typically small and corresponds to the damping or excitation of the oscillations, which are not included in the ray tracing formalism. We thus changed the current value of k_r to the real part of the selected root after each timestep. If the imaginary part was unequal to zero, Eq. (A.7) was not satisfied anymore and ω was not conserved. The possible discrepancy between the initial value of ω and the value of ω after each timestep depends on how much the imaginary part of the selected root deviated from 0. Therefore, the value of ω after each timestep is a sensitive measure of accuracy of the rays. If the value of ω was different to the original value by more than 0.0001% at any point during the integration, we did not consider that ray for further analysis³. We compare our method for correcting k_r and the Newton–Raphson scheme used by Loi (2020a) in Appendix A.5.

2.5. Trapped fraction

There are two possible outcomes for rays launched into the magnetized stellar core: They can be either reflected or trapped (Loi 2020a). In Fig. 2, we show examples for both a reflected and a trapped ray. Reflected rays propagate back and forth (in distance from the center) between two turning points and regularly enter and exit the area of the internal magnetic field. As shown in Fig. 2, the components of the wavevector, and thus also the wavenumber of a reflected ray, undergo periodic changes, but remain bounded in time. This is in contrast to the trapped rays. Once a trapped ray enters the area of the magnetic field, it never emerges and the wavenumber diverges. Since the rays visualize the paths of the energy flow of the oscillations, a trapped ray can be interpreted as a loss of the observable mode energy. To decide whether a ray is reflected or trapped, we used a classification scheme, which is described in Appendix A.4. It is an updated version of the ray classification scheme used by Loi (2020a). We compare our ray classification scheme with that of Loi (2020a) in Appendix A.5.

Fig. 2. Trajectory of a reflected and trapped ray (top row) and the corresponding wavevector components as a function of time (bottom row) for the model 2A with Bcen = Bcrit, cen, ω = 10 ωdyn, and l = 1. The reflected ray with θ0 = 128.9° and α = 3.7° is shown on the left, the trapped ray with θ0 = 98.7° and α = 90.9° is shown on the right. The dotted black circle in the top row represents the maximum extend of the magnetic field rf. Note that the limits of the axes do not coincide between left and right for the sake of clarity.

When all 1200 rays of a set have been classified as either reflected or trapped, or identified as unclassified due to numerical issues during the integration (see Appendix A.4), the trapped fraction f_T of that set of rays could be determined. It is an estimate of the fraction of energy in the gravity mode cavity that is lost from the observable global mode due to the interaction with the internal magnetic field. The trapped fraction can be calculated as follows:

$$\begin{array}{c}\hfill f_{\mathrm{T}}=\frac{1}{\#\mathrm{total}}(\#\mathrm{trapped}+\frac{\#\mathrm{unclassified}}{2}\pm \frac{\#\mathrm{unclassified}}{2}),\end{array}$$(12)

where #trapped is the number of trapped rays, #unclassified is the number of unclassified rays, and #total = 1200 is the total number of rays in the set. The estimated value of f_T corresponds to the case where half of the unclassified rays are actually reflected and the other half is trapped. The uncertainty of f_T introduced by the last term in Eq. (12) is numerical and corresponds to the fact that not all rays of a set could be identified as reflected or trapped. If these numerical difficulties did not exist, we would expect the unclassified rays to be either reflected or trapped. The minimum value of f_T permitted by the uncertainties therefore represents the case in which all unclassified rays are reflected, while the maximum value of f_T represents the case in which all unclassified rays are trapped.

3. Results

Before we discuss the results of our magneto-gravity ray tracing analysis, we summarize the changes that we have made to the procedure of Loi (2020a) described in Sect. 2. First, we used input stellar structure profiles with an improved radial resolution. This was achieved by smoothing the profiles directly in MESA by introducing a small amount of additional mixing (see Sect. 2.1). Second, we chose the initial positions of the rays so that they are evenly distributed on the surface of a sphere if θ₀ is considered as latitude and α as longitude (see Fig. A.2). Third, we used a different method to correct the wavenumber after each timestep during the integration (see Sect. 2.4 and Appendix A.3). Finally, we used an updated ray classification scheme (see Appendix A.4).

3.1. Numerical estimates of the trapped fraction

In Fig. 3, we show the trapped fraction of the models 2A, 2B, and 2C as a function of the magnetic field strength B_cen for different combinations of ω and l. Although we used a similar method as Loi (2020a) to compute f_T, the trapped fractions calculated in this work are smaller than those obtained by Loi (2020a) when comparing our model 2A (2B, 2C) to their model A (B, C). The reason for this is most likely a systematic difference in the propagation of the rays. We discuss these differences in Appendix A.6. In particular, we show that the trajectories of our rays match those of Loi (2020b), who also uses a slightly adapted version of the ray integration procedure presented in Loi (2020a).

Fig. 3. Trapped fraction as a function of the central field strength of the internal magnetic field for the models 2A (top), 2B (middle), and 2C (bottom). Different colors indicate different combinations of the frequency ω and the spherical degree l (see legend). The dotted lines mark the corresponding central critical field strength. The error bars are determined by the uncertainty of fT (see Eq. (12)). The line connecting the symbols is a fit to the trapped fraction presented in Sect. 3.2.

Overall, we recover three qualitative results reported by Loi (2020a):

• The trapped fraction f_T is negligible when B_cen is smaller than a certain threshold value. If B_cen is greater than this threshold, f_T increases with increasing B_cen.

• The threshold value of B_cen is lower for lower ω. If B_cen is greater than the threshold, f_T is higher at a certain B_cen when ω is lower.

• The threshold value of B_cen is lower for higher l. If B_cen is greater than the threshold, f_T is higher at a certain B_cen when l is higher.

Furthermore, we find that f_T always approaches 100% for values of B_cen much higher than the corresponding central critical field strength B_crit, cen. This applies to all models shown in Fig. 3, including the more evolved models 2B and 2C, and is in contrast to the trend with stellar evolution observed by Loi (2020a). They conclude that the maximum value of f_T appears to saturate at values less than 100% for their more evolved models. This discrepancy cannot be explained by the uncertainties of their f_T-estimate. However, as we show in Appendix A.5, it is possible that these uncertainties are underestimated due to numerical artifacts. This means that the value at which f_T saturates could actually be consistent with 100% also for their more evolved models.

As shown in Fig. 3, the exact value of f_T as a function of B_cen depends on the input stellar model, ω, and l. The same can be said about B_crit, cen, which also depends on these parameters (see Eq. (6)). In Fig. 4, we show f_T for the same models as in Fig. 3 with the central magnetic field strength expressed in terms of B_crit, cen. In short, this means that we aligned the dotted lines indicating B_crit, cen shown in Fig. 3. Figure 4 shows that the entire dependence of f_T on the stellar model, ω, and l is explained by the dependence of B_crit, cen on these parameters. Crucially, this means that f_T can be expressed as a function of B_cen/B_crit, cen only. Therefore, f_T can be directly estimated for a given B_cen if B_crit, cen is known.

Fig. 4. Same as Fig. 3, now as a function of the central field strength of the internal magnetic field divided by the corresponding central critical field strength. The upper limit of the horizontal axis is now 10. Note that Bcrit, cen is different for each model, as well as for each combination of ω and l (see Eq. (6)).

This also holds true if the mass and metallicity of the star are changed. We show plots similar to Fig. 4 for the models 1A, 1B, 1C, 1a, 1b, 2a, and 2b in Appendix B.

3.2. Fitted expression for the trapped fraction

The trapped fraction can be expressed as a function of the ratio between central field strength and the central critical field strength (see Fig. 4). It is therefore possible to fit an analytical expression to f_T that only depends on B_cen/B_crit, cen. Since B_crit, cen can be determined using Eq. (6), f_T is reduced to a function of B_cen for a given stellar model.

For simplicity, we only considered the model 2b with ω = 10 ω_dyn and l = 1 for the fitting procedure. We chose this combination of parameters because it yields small uncertainties for f_T and shows good agreement with the trapped fractions of the other models. The resulting fitted function describing the trapped fraction is as follows:

$$\begin{array}{cc}& f_{\mathrm{T}}^{(\mathrm{fit})}(B)=max[0,1-0.635{\left(\frac{B}{B_{\mathrm{crit}}}\right)}^{-1}\hfill \\ \hfill & +0.438{\left(\frac{B}{B_{\mathrm{crit}}}\right)}^{-2}-0.272{\left(\frac{B}{B_{\mathrm{crit}}}\right)}^{-3}].\hfill \end{array}$$(13)

According to this expression, the onset of the magnetic suppression takes place at B ≈ 0.628 B_crit. Here, we have omitted the subscript _cen because we show in the following section that Eq. (13) is also valid when considering the ratio between the field strength and the critical field strength at a location other than the center of the star. Note that the choice of the form of the function describing $f_{\mathrm{T}}^{(\mathrm{fit})}$ is purely empirical.

The fit to the trapped fraction is shown in Figs. 3 and 4, as well as Figs. B.1 and B.2, together with the numerical results from the ray tracing analysis. We discuss the applicability of Eq. (13) in the following section.

3.3. Applicability of the fitted trapped fraction

Since we have only considered Prendergast field configurations (see Sect. 2.2) for the calculation of the trajectories of the rays, Eq. (13) is strictly speaking only valid for exactly this type of internal magnetic field. Adopting other magnetic field topologies, such as a pure dipole field, led to numerical difficulties and likely requires changes to our integration method. This will be subject of a future work. However, we can deduce some characteristics of how the behavior of the rays is expected to change with a modification of the field topology from the results presented in this work.

In Fig. 5, we compare the trapped fraction obtained when using a Prendergast field with the trapped fraction when only the radial component of the Prendergast field is considered (i.e., B_θ = B_ϕ = 0), as well as with the radial component of an improvised magnetic field B_r^pla. The improvised field has a plateau and scales with the critical field strength at the hydrogen-burning shell B_crit, H. In addition, it has a dipole-like dependence on the latitude, similar to the Prendergast field (i.e., B_r^pla ∝ cos θ). For the Prendergast field, Fig. 5 shows that f_T does not change significantly when neglecting B_θ and B_ϕ, even though B_θ has the same maximum absolute strength as B_r. This indicates that B_r dominates the interaction of the magnetic field with the rays when the absolute values of B_θ and B_ϕ are roughly comparable or smaller than that of B_r.

Fig. 5. Top: Normalized critical field strength as a function of radius. We also show the maximum absolute value of the three components of the Prendergast field (see Sect. 2.2), which is scaled to Bcen = Bcrit, cen. Furthermore, we show the radial component of an improvised magnetic field that has a plateau and scales with the critical field strength at the hydrogen-burning shell Bcrit, H. The strength of the improvised field at the hydrogen-burning shell is set to BH = Bcrit, H. The hydrogen-burning shell is indicated by the vertical dotted line. Bottom: Trapped fraction as a function of the central field strength or the field strength at the hydrogen-burning shell for model 1A with ω = 10 ωdyn and l = 1. The black symbols correspond to the case where all components of the Prendergast field are taken into account. They are the same as the black symbols in the top panel of Fig. B.1. The blue symbols indicate that only the radial component of the Prendergast field has been used and the other components were set to zero. The red symbols correspond to magnetic field featuring the plateau that is scaled by Bcrit, H instead of Bcrit, cen. The dotted line marks either Bcrit, cen or Bcrit, H. The line connecting the symbols is the fit to the trapped fraction (Eq. (13)).

Interestingly, the behavior of the trapped fraction as a function of field strength when considering the improvised magnetic field does not differ significantly from the case where the Prendergast field is used. This applies under the condition that the improvised field is scaled by the field strength at the hydrogen-burning shell B_H/B_crit, H instead of B_cen/B_crit, cen. It indicates that it does not matter where in the core the field strength first approaches the critical field strength, and that the ratio between the actual field strength and the critical field strength at that location determines f_T. A prime location for this to happen is the hydrogen-burning shell, as the critical field strength is at a minimum there. In our models, the Prendergast field always reaches the critical field strength in the center of the star first, even if r_f lies outside the hydrogen-burning shell. Since the critical field strength at the center and at the hydrogen-burning shell differ by up to almost two orders of magnitude for some models (see Table 1), different field topologies that approach the critical field strength at either the center or the hydrogen-burning shell could cause a dichotomy in the trapped fraction.

When representing a set of rays as a unisphere parameterized by θ₀ and α (see Fig. A.2), we can see a clear structure in the distribution of reflected and trapped rays. The reflected rays are localized at two poles that are located opposite to each other. If the trapped fraction of the set of rays increases, the poles of reflected rays become smaller; if the trapped fraction decreases, the poles expand. The reflected rays are forming these poles because they share two properties: A small initial launch latitude θ₀ and a small latitudinal component of the wavevector k_θ. This means that the reflected rays are those rays that propagate in an equatorial band. The width of this band is being controlled by B_cen. This is related to the fact that the radial component of the magnetic field B_r is proportional to cos θ (Eq. (3)), which means that B_r vanishes at the equator. Since this is a characteristic of dipole fields, the same trend is likely to be observed for all dipole-type field topologies. However, it is not the case for quadrupole fields or other higher-order topologies. Therefore, we expect the rays to behave differently for these configurations⁴.

In conclusion, we expect Eq. (13) to be applicable for internal magnetic fields with |B_r|≳|B_θ|,|B_ϕ| and a dipole-like dependence on the latitude (i.e., B_r ∝ cos θ). Furthermore, Eq. (13) is also prone to the limitations discussed in Sect. 4.4. These include the use of the WKBJ approximation and the fact that we only considered internal magnetic fields with a dipole-like dependence on latitude.

3.4. Evolution of the trapped fraction

The critical field strength generally decreases as the star ascends the RGB (see Table 1). In addition, the strength of an internal magnetic field contained in the radiative core of an RGB star increases as the star evolves due to the conservation of magnetic flux and the contraction of the core region. Both of these effects cause the ratio B_cen/B_crit, cen to grow, which means that f_T is expected to become higher during evolution in the RGB. In Fig. 6, we show the evolution of the trapped fraction of both the dipole and quadrupole modes at ω = 2πν_max for the four evolutionary tracks presented in Fig. 1. For each track, we show the trapped fraction for three different central initial strengths of the magnetic field, which corresponds to the field strength at the end of the MS. These initial field strengths are scaled according to the conservation of magnetic flux during the evolution of the star.

Fig. 6. Trapped fraction evaluated with Eq. (13) at ω = 2πνmax, where νmax is the frequency of maximum oscillation power, as a function of the evolution of the star indicated by νmax (stars evolve from right to left). Each panel features one of the four evolutionary tracks presented in Fig. 1. Rows show tracks with the same stellar mass and columns show tracks with the same metallicity. The blue (red) lines correspond to the dipole (quadrupole) modes. Different linestyles (solid, dashed, dotted) represent different initial strengths of the radial component of the magnetic field at the stellar center, which corresponds to the field strength the end of the MS. Due to the conservation of magnetic flux, the magnetic field strength increases as the stars evolve. The small knees that occur at the points where fT becomes greater than 0, which are visible in some of the curves, are caused by the limited temporal resolution of the evolutionary tracks. We show the evolution of the trapped fraction considering BH/Bcrit, H instead of Bcen/Bcrit, cen in Fig. B.4.

As expected, f_T is initially 0% for all of our evolutionary tracks because the ratio B_cen/B_crit, cen is at first small. If the magnetic field strength has become strong enough, f_T increases and approaches 100%. The onset of the trapping occurs earlier in the evolution of the star when the initial field strength is higher. Similarly, the suppression of the quadrupole mode always starts at higher values of the frequency of maximum oscillation power ν_max compared to the dipole modes when considering the same initial magnetic field strength. This is a direct consequence of Eq. (6).

The stellar mass also affects the trapped fraction, with trapping starting at higher values of ν_max for lower mass stars. This effect becomes weaker as the star evolves. It is questionable whether this feature is observable, as the typical magnetic field strength at the end of the MS of stars with M_* = 1 M_⊙ and 2 M_⊙ may be different. In fact, the study by Bugnet et al. (2021) suggests that this may not be possible. They estimate the strength of the fossil magnetic field left by a potential convection-driven dynamo phase and show that a change in stellar mass and metallicity does not significantly alter the magnetic frequency shift of mixed modes of similar frequency on the RGB. This implies that the effect of the internal magnetic field may be comparable for different stellar masses and metallicities at a certain ν_max⁵. The influence of stellar metallicity on the value of f_T in Fig. 6 appears to be smaller than that of the initial field strength, the spherical degree, and the stellar mass.

The values selected for the initial magnetic field strength range from 10 to 400 kG. They are an order of magnitude larger than the MS field strengths inferred from the magnetic frequency shifts of the dipole modes by Deheuvels et al. (2023) for an axisymmetric dipolar field configuration. However, these field strengths cannot be directly compared to each other. This is because the central field strength of the Prendergast field is the maximum strength of the magnetic field, while the field strengths estimated by Deheuvels et al. (2023) are weighted averages dominated by the contribution of the hydrogen-burning shell (e.g., Li et al. 2022; Bhattacharya et al. 2024). In particular, the radial component of the magnetic field in the case of the Prendergast field decreases rapidly with increasing distance from the center, while a uniform distribution is assumed for the estimation of the MS field strengths by Deheuvels et al. (2023). Therefore, the central field strengths of the Prendergast field required to suppress the multipole modes are higher than the average values estimated from the observed frequency shifts.

The improvised magnetic field introduced in Sect. 3.3 is more similar to the field configuration assumed by Deheuvels et al. (2023), because the radial component of the field is uniform in the region at and interior to the hydrogen-burning shell. Since the improvised field is scaled by the field strength at the hydrogen-burning shell B_H, the MS field strengths estimated by Deheuvels et al. (2023) should be compared with B_crit, H instead of B_crit, cen. In Fig. B.4, we show that the evolution of the trapped fraction is comparable to that in Fig. 6 when using initial field strengths ranging from 1 to 40 kG and the ratio B_H/B_crit, H instead of B_cen/B_crit, cen. This indicates that the MS field strengths estimated by Deheuvels et al. (2023) might be sufficient to cause dipole mode suppression starting from ν_max ≈ 250 μHz, which is in good agreement with observations (e.g., Stello et al. 2016b; Mosser et al. 2017; Coppée et al. 2024). Outliers from the population of suppressed dipole mode stars with higher values of ν_max could be caused, for example, by slightly higher MS field strengths or different field configurations.

4. Discussion

4.1. Implications for observed power spectral densities

The trapped fraction obtained in this work is not a quantity that can be measured directly in observations. Instead, the suppressed dipole mode stars are classified based on the mode visibility, which is defined as the ratio of the amplitude of the multipole modes of a spherical degree l to the radial mode of the same pressure radial order (e.g., Fuller et al. 2015; Stello et al. 2016b; Cantiello et al. 2016). The mode visibility depends on the trapped fraction, but there is currently no general expression linking these two quantities. As of yet, this is only possible for the two boundary cases, f_T = 0% and f_T = 100% (e.g., Fuller et al. 2015; Takata 2016a; Mosser et al. 2017). If f_T = 0%, the mixed modes are not suppressed and have normal amplitudes. If, on the other hand, f_T = 100%, all the energy that leaks into the gravity mode cavity is dissipated by the internal magnetic field, which corresponds to the case discussed by Fuller et al. (2015) (see also Shibahashi 1979; Unno et al. 1989). For intermediate values of f_T, the observed visibility is expected to lie between these two limiting cases. Therefore, we can make qualitative comparisons between the expected behavior of f_T and the observed mode visibilities.

The results of our magneto-gravity ray tracing analysis show that the trapped fraction is a function of the frequency ν = ω/(2π). For a stellar model with an internal magnetic field of a certain strength, f_T is equal to 0% at higher frequencies and increases at frequencies below a threshold value. It is therefore expected that the suppression of the multipole modes occurs at lower frequencies. It depends on the respective model and magnetic field strength whether the suppression falls within the observable frequency range. The PSD of the RGB star KIC 6975038 (Fig. 7) features clear signs of dipole mode suppression at lower values of ν, while the amplitudes of the dipole modes appear to be normal at higher ν. This phenomenon has been associated with the onset of the dipole mode suppression at lower ν caused by an internal magnetic field (e.g., Mosser et al. 2017; Loi 2020a; Deheuvels et al. 2023). The case of KIC 6975038 is discussed in more detail in Sect. 4.3.

Fig. 7. Background-corrected PSD of KIC 6975038 obtained from Kepler data using the background model of Kallinger et al. (2014) (for details, see Coppée et al. 2024). Radial modes are indicated by the gray shaded area, dipole modes by the blue area, and quadrupole modes by the reddish area. In addition, we show the trapped fraction as a function of frequency computed for a stellar model with M* = 1.3 M⊙, Zinit = 0.02, νmax ≈ 128 μHz, and Bcen ≈ 6.23 MG (see Sect. 4.3) for both the dipole modes (in blue) and quadrupole modes (in red).

As soon as f_T is greater than 0%, it initially increases rapidly while ν decreases. Eventually, this sharp increase flattens out when f_T begins to approach 100%. The flattening of the frequency dependence of the trapped fraction is consistent with observations. So far, no trend in the observed visibilities with the frequency has been detected for the suppressed dipole mode stars, unless they exhibit features expected to be associated with the onset of dipole mode suppression (such as KIC 6975038).

The trapped fraction also depends on the spherical degree l. Our computations show that modes of higher l are suppressed at lower magnetic field strengths. For the PSD of a star, this means that the suppression of quadrupole modes starts at a higher frequency than the suppression of dipole modes (see Fig. 7 and Loi 2020a). In observations, however, the visibility of the quadrupole modes of the suppressed dipole mode stars is usually higher than that of the dipole modes. This is because the mode visibility depends not only on the trapped fraction, but also on the coupling strength of the mixed modes (e.g., Shibahashi 1979; Fuller et al. 2015; Takata 2016a,b; Mosser et al. 2017). The coupling strength is a measure of how strongly the oscillations in the pressure and gravity mode cavities interact. At higher spherical degrees, the interaction between the mode cavities is weaker (e.g., Hekker & Christensen-Dalsgaard 2017), such that the coupling strength is lower and the mode visibility is higher (e.g., Cantiello et al. 2016). Overall, our results suggest that the quadrupole modes of stars with suppressed dipole modes should also be suppressed, although their mode visibilities might be higher.

4.2. Trapped fraction along the red-giant branch

The observed dipole mode visibilities of RGB stars show a clear dichotomy. On one hand, there is a population of stars with normal dipole mode visibilities. One the other hand, there are the suppressed dipole mode stars, which have been observed for ν_max ≲ 250 μHz (Stello et al. 2016b; Mosser et al. 2017). While the two populations start to merge at ν_max ≲ 70 μHz (Stello et al. 2016b), there are only a few stars with ν_max > 70 μHz whose dipole mode visibility lies between the two branches and cannot be assigned to either of them.

In Sect. 3.4, we investigated the evolution of the trapped fraction under the assumption of conservation of magnetic flux for different magnetic field strengths at the end of the MS. We have shown that for stars with solar mass, the dipole mode suppression can start at ν_max ≲ 200 μHz if the central field strength at the end of the MS is ≲400 kG. If the strength of the field at the hydrogen-burning shell is considered instead, a field strength of ≲40 kG at the end of the MS is sufficient to cause dipole mode suppression at ν_max ≲ 250 μHz. These magnetic field strengths at the hydrogen-burning shell agree with those inferred by Deheuvels et al. (2023) from observed average frequency shifts of the dipole modes of RGB stars assuming an axisymmetric dipolar magnetic field configuration. Note however that these field strengths might be overestimated (Rui et al. 2024).

The small number of stars with observed dipole mode visibilities that lie between the visibilities of the population of normal RGB stars and the population of suppressed dipole mode stars could be related to the sharp increase in f_T immediately after the onset of dipole mode suppression (see Figs. 6 and B.4). It is not unreasonable to assume that within the uncertainties, there is a threshold value of f_T above which the dipole mode visibility is indistinguishable from that of a star experiencing complete dipole mode suppression (i.e., f_T = 100%). Once f_T > 0%, the trapped fraction could quickly reach this threshold as the star ascends the RGB since f_T rises sharply after the start of the dipole mode suppression. This would lead to a short time window for the detection of an RGB star with an intermediate visibility in terms of stellar evolution, which is in line with the current state of observations.

Contrary to Loi (2020a), we find that the trapped fraction always follows Eq. (13), even for more evolved RGB stars. However, the observed dipole mode visibilities of the suppressed dipole mode stars increase with decreasing ν_max (Stello et al. 2016b; Mosser et al. 2017). This trend can simply be explained by the decrease in coupling strength during the evolution of an RBG star (Cantiello et al. 2016). The fact that f_T always saturates at 100% in our simulations therefore does not contradict the observed behavior of dipole mode visibility.

Lastly, assuming that all RGB stars have an internal magnetic field of a certain strength at the end of the MS, we expect the fraction of stars with suppressed dipole modes to increase with decreasing ν_max. This is because the trapped fraction increases as the star evolves along the RGB (see Figs. 6 and B.4). However, this trend could be obscured in the observations due to the generally higher dipole mode visibility of the suppressed dipole mode stars at low ν_max (see Stello et al. 2016b; Mosser et al. 2017), which is caused by the decreasing coupling strength. Furthermore, the radiative damping in the core of RGB stars becomes stronger at lower ν_max (e.g., Dupret et al. 2009; Grosjean et al. 2014; Hekker & Christensen-Dalsgaard 2017; Bordadágua et al., in prep.), which makes the identification of effects caused by an internal magnetic field more challenging.

4.3. Frequency shift of KIC 6975038

The star KIC 6975038 is of particular interest for the hypothesis of an internal magnetic field, because the dipole modes are suppressed at lower frequencies, while their amplitudes appear to be normal at higher frequencies (see Fig. 7). Deheuvels et al. (2023) show that the dipole modes split into asymmetric multiplets at higher frequencies, of which one component per multiplet was identified. They measure the average frequency shift of the dipole modes at ν_max as δω₀/(2π) = 7.18 ± 0.03 μHz and attribute it to the presence of an internal magnetic field. Since only one component of each multiplet has been identified, the configuration of the internal magnetic field of KIC 6975038 is currently not constrained (Deheuvels et al. 2023).

Here, we used the first-order perturbation approach presented by Li et al. (2022) to calculate the expected average frequency shift for a selected stellar model with an internal magnetic field (see also Bugnet et al. 2021; Loi 2021; Bugnet 2022). The average frequency shift only depends on the radial component of the magnetic field. It is important to note that this approach assumes that the magnetic field only constitutes a small perturbation to the restoring forces of the oscillations. This assumption is broken when the field strength approaches the critical field strength. Therefore, our predictions of the average frequency shift are expected be underestimated (Rui et al. 2024). This means that our estimates of the observed average frequency shift must be lower than the observed value, but likely of the same order of magnitude (Rui et al. 2024), to be consistent with the current understanding of magnetic frequency shifts.

We set the internal field strength to 70% of the critical field strength of the dipole modes at ω = 2πν_max, so that the onset of dipole mode suppression is at about ν ≈ 135 μHz (see Fig. 7). Yu et al. (2018) estimate the mass of KIC 6975038 to be around 1.3 M_⊙. Therefore, we computed an additional evolutionary track with this mass and Z_init = 0.02 using the MESA inlist shown in Appendix C. We then selected a model with ν_max ≈ 128 μHz from this track, which is close to the observed value of ν_max (e.g., Deheuvels et al. 2023). In Fig. 7, we show f_T as a function of ν for this model.

We tested three topologies of the internal magnetic field: two Prendergast configurations and the improvised field introduced in Sect. 3.3. They were scaled such that B_cen = 0.7 B_crit, cen and $B_{\mathrm{H}}^{\mathrm{pla}}=0.7B_{\mathrm{crit},\mathrm{H}}$, respectively. When considering Prendergast fields, the average frequency shift depends on the maximum extent of the field r_f. We investigated two scenarios: A Prendergast field that is confined by the region at and within the hydrogen-burning shell and a Prendergast field that extends over the entire radiative core. Since the magnetic frequency shift is only sensitive to the region at and within the hydrogen-burning shell (e.g., Li et al. 2022; Bhattacharya et al. 2024), the average frequency shift induced by the improvised field does not change significantly when the maximum extent of the field is changed, as it behaves like a plateau in this region.

In Table 3 and Fig. 8, we show magnetic field strengths and the resulting average frequency shifts of the dipole modes for the three different field configurations for the azimuthal orders m = 0 and m = ±1 (see Li et al. 2022; Deheuvels et al. 2023)⁶. Note that due to the limitations of the Hamiltonian ray tracing approach described in Sect. 4.4 and the simplicity of the stellar model, the magnetic field strengths determined for KIC 6975038 should not be considered as measurements, but as an order of magnitude estimate (see Sect. 4.4). Since recent studies estimate the mass of KIC 6975038 to be less than 1.3 M_⊙ (e.g., Coppée et al., in prep.), we repeated the process described above for a model of the evolutionary track with M_* = 1 M_⊙ and Z_init = 0.02, as well as for a model of the one with M_* = 1 M_⊙ and Z_init = 0.005. The field strengths and average frequency shifts inferred for these models are also shown in Table 3 and Fig. 8.

Fig. 8. Average frequency shift of the dipole modes caused by various internal magnetic field configurations as a function of the magnetic field strength for three stellar models. The field strengths have been chosen such that Bcen = 0.7 Bcrit, cen and BH = 0.7 Bcrit, H. The critical field strengths were calculated at ω = 2πνmax with l = 1. Colors correspond to different stellar models and symbols indicate that a different magnetic field configuration has been used. Average frequency shifts corresponding to the azimuthal order m = 0 (m = ±1) are shown as solid symbols (outlines). The field configurations and the calculation of the average frequency shifts are described in Sect. 4.3 and in the notes under Table 3, which also summarizes the data shown in this figure. Note that the horizontal axis is broken.

All of our estimates for the average frequency shift are systematically lower than the observed value. However, they agree with it by less than a factor of 5.5. Taking into account the limitations of the first-order perturbation approach and the uncertainties in the estimation of the magnetic field strength based on the ray tracing method, all tested topologies of the internal magnetic field are in agreement with the observations. Based on these results, we cannot favor any of the field configurations over the others. Overall, our order-of-magnitude analysis suggests that both the suppression of the dipole modes at lower ν and the frequency shift of the modes at higher ν can be explained by the presence of an internal magnetic field. Our results thus support the hypothesis that the cause of these phenomena is indeed magnetic in nature. However, a more sophisticated treatment of the magnetic effects is required to firmly validate this interpretation of the features observed in the oscillations of KIC 6975038.

4.4. Limitations

As mentioned in Sect. 2.3, the Hamiltonian ray tracing approach used to estimate the trapped fractions relies on the WKBJ approximation of zeroth order. This means that higher order effects, such as the evanescence of the oscillation in certain regions of the star, mode excitation and damping, as well as phase dependencies, are not taken into account. However, none of these effects should have an influence on the trapped fraction (see also Loi 2020a). More important for our application is the assumption that the horizontal wavelength is much smaller than the horizontal variation of the internal magnetic field. Since the observable oscillation modes have a low spherical degree, this assumption is violated in our procedure and could lead to a bias in our estimates of the trapped fraction.

Nevertheless, there are a number of arguments that suggest that Eq. (13) describes the trapped fraction with reasonable accuracy. First, the trapped fraction is an increasing function of the strength of the internal magnetic field and transitions smoothly from 0% to 100%. This is exactly what is needed to explain the observation that some of the suppressed dipole mode stars have dipole modes with a mixed character (Mosser et al. 2017, Coppée et al. in prep.). Second, the transition from 0% to 100% takes place when the strength of the internal magnetic field approaches the critical field strength. This was predicted by Fuller et al. (2015) and later verified by Rui & Fuller (2023). Third, the average frequency shift estimated for KIC 6975038 corresponding to the field strengths obtained with Eq. (13) is systematically lower than the average frequency shift observed by Deheuvels et al. (2023). This is in line with the results of Rui et al. (2024). Fourth, Eq. (13) is a good fit for the trapped fraction as a function of the internal field strength also for rays with a high spherical degree of l = 10 and l = 100 (see Fig. B.3). These rays have a much smaller horizontal wavelength and therefore correspond more to the short wavelength assumption. By extrapolation, we conclude that Eq. (13) should also apply to the observable modes with smaller l, since it appears to be independent of the spherical degree. However, in order to verify these results, more sophisticated simulations of the interaction of the waves with the internal magnetic field are required. As long as this is not the case, Eq. (13) should be considered as an order of magnitude estimate that can be used to describe general qualitative trends of the multipole mode suppression.

A further limitation arises from the fact that we have only considered internal magnetic fields with a dipole-like dependence on the latitude. In reality, other field geometries are possible and have been observed (Li et al. 2022, 2023; Hatt et al. 2024). The applicability of Eq. (13) to different field configurations was discussed in Sect. 3.3 and we refer the reader to this section for more information. In short, we expect a different behavior of the trapped fraction for magnetic field topologies with higher order contributions such as a quadrupole component. However, it can be assumed that the transition of the trapped fraction from 0% to 100% takes place when the internal field strength approaches the critical field strength, which means that the resulting changes are most likely of a quantitative nature.

Finally, we have neglected rotation. The typical rotation frequencies of RGB stars are small compared to the frequencies of the oscillations. Therefore, the inclusion of rotation is not expected to significantly affect our conclusions (see also Loi 2020a,b), although this remains to be verified. For the dipole modes, Loi (2020b) speculate that the internal magnetic field could primarily damp the zonal oscillations (i.e., the component with m = 0 with respect to the magnetic axis). The reason for this is that the reflected rays are those that propagate in an equatorial band, which corresponds to sectoral oscillations. In the particular case where the magnetic axis and the rotation axis are misaligned, the oscillation modes are expected to separate into (2l + 1)² individual peaks (e.g., Gough & Thompson 1990; Loi 2021). However, following the argument of Loi (2020b), only the rotationally split peak corresponding to m = ±1 with respect to the magnetic axis could be visible. This would reduce the number of observable modes.

5. Summary

In this work, we built on the study of Loi (2020a) and performed a Hamiltonian ray tracing analysis of magneto-gravity waves using numerical models of RGB stars with different masses, metallicities and ages as background. From these simulations, we estimated the trapped fraction, which is a measure of the mode energy lost through interaction with an internal magnetic field. We confirm the results of Loi (2020a) that the trapped fraction transitions from 0% to 100% with increasing strength of the internal magnetic field. In particular, the trapped fraction depends on the magnetic field strength, the frequency of the oscillations, the spherical degree, and the background model (Sect. 3.1). However, we do not reproduce the result of Loi (2020a) that the trapped fraction saturates at values smaller than 100% for more evolved RGB stars. This could be a consequence of the changes we have made to the integration of the rays (Sect. 2.4) and the ray classification scheme (Appendix A.4).

Furthermore, we find that these dependencies are encompassed by the dependence of the critical magnetic field strength (Eq. (6)) on these parameters. We show that the trapped fraction depends only on the ratio of the central strength of the internal magnetic field to the central critical field strength. The behavior of the trapped fraction is described by Eq. (13) for dipole-like configurations of the internal magnetic field. In addition, we show in Sect. 3.3 that Eq. (13) is also applicable when the ratio of the magnetic field strength to the critical field strength is evaluated at the location of the hydrogen-burning shell where the critical field strength has a minimum, instead of at the center. Our results suggest that the ratio of the internal field strength to the critical field strength at the location where the field strength first approaches the critical field strength determines the trapped fraction.

Equation 13 provides a flexible tool that can be used to improve our understanding of the multipole mode suppression observed for RGB stars. It can, for example, be used to describe the trapped fraction as a star ascends the RGB, taking into account the evolution of the strength of the internal magnetic field (Sect. 3.4). In addition, Eq. (13) can be used to study the onset of mode suppression in the observed power spectra densities of RGB stars (Sect. 4.3). Note however that due to the limitations of the Hamiltonian ray tracing approach, Eq. (13) should be regarded as an order of magnitude estimate and needs to be validated by future studies using a more complete treatment of the magnetic effects (for details, see Sect. 4.4).

Based on our results and those of previous studies (Fuller et al. 2015; Mosser et al. 2017; Loi 2020a), we predict that the quadrupole modes of stars with suppressed dipole modes should always be suppressed (Sect. 4.1). In addition, we expect that the ratio of suppressed dipole modes to RGB stars with normal dipole mode amplitudes should increase along the RGB (i.e., with decreasing ν_max). However, this trend could be obscured by the generally higher dipole mode visibilities of the suppressed dipole mode stars at low ν_max (Sect. 4.2; Stello et al. 2016b; Cantiello et al. 2016).

In the next article in this series, we will investigate the relationship between intermediate values of the trapped fraction (i.e., 0%< f_T < 100%) and the observable visibility of the multipole modes. Furthermore, we will relate the different morphologies observed in the PSD of stars with suppressed dipole modes (Coppée et al. in prep.) to physical quantities.

Acknowledgments

We thank the anonymous referee for their valuable and constructive feedback that improved the manuscript. We acknowledge funding from the ERC Consolidator Grant DipolarSound (grant agreement # 101000296). In addition, we acknowledge support from the Klaus Tschira Foundation.

References

Appendix A: Additional information on the magneto-gravity ray tracing

A.1. Buoyancy frequency profiles

In Fig. A.1, we show the buoyancy frequency profiles of the models 2A, 2B, and 2C. As the models were selected to have a similar radius to models A, B and C of Loi (2020a), the overall qualitative trends remain unchanged compared to these models. However, there are small quantitative differences between our buoyancy frequency profiles and the profiles of Loi (2020a) due to the changes in the MESA version, the inlist, and the interpolation onto another grid performed by Loi (2020a).

Fig. A.1. Profiles of the buoyancy frequency as a function of radius for the models 2A, 2B, and 2C. The launch radius of the rays r0 is shown in red for each model, the maximum extend of the internal magnetic field rf is shown as the dashed black line.

A.2. Hamilton’s equations for magneto-gravity waves

Here, we show the set of first-order ordinary differential equations describing the propagation of a magneto-gravity wave packet in time in spherical polar coordinates (Loi 2020a):

$$\begin{array}{cc}& \frac{\mathrm{d}r}{\mathrm{d}t}=\frac{{\omega }_{\mathrm{A}}}{\omega }{V}_{\mathrm{A}r}-\frac{N^2{\kappa }_{\perp }^2{\kappa }_r}{\omega k},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& \frac{\mathrm{d}\theta }{\mathrm{d}t}=\frac{{\omega }_{\mathrm{A}}}{\omega }\frac{V_{\mathrm{A}\theta }}{r}+\frac{N^2{\kappa }_{\theta }{\kappa }_r^2}{\omega kr},\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& \frac{\mathrm{d}\varphi }{\mathrm{d}t}=\frac{{\omega }_{\mathrm{A}}}{\omega }\frac{V_{\mathrm{A}\varphi }}{rsin\theta }+\frac{N^2{\kappa }_{\varphi }{\kappa }_r^2}{\omega krsin\theta },\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& \frac{\mathrm{d}{k}_r}{\mathrm{d}t}=\frac{{\omega }_{\mathrm{A}}}{\omega }(\frac{k_{\theta }{V}_{\mathrm{A}\theta }+k_{\varphi }{V}_{\varphi }}{r}-\mathbf{k}·\frac{\partial {\mathbf{V}}_{\mathrm{A}}}{\partial r})\hfill \\ \hfill & +\frac{N{\kappa }_{\perp }^2}{\omega }(\frac{N{\kappa }_r^2}{r}-\frac{\mathrm{d}N}{\mathrm{d}r}),\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& \frac{\mathrm{d}{k}_{\theta }}{\mathrm{d}t}=-\frac{{\omega }_{\mathrm{A}}}{\omega }\frac{1}{r}(k_{\theta }{V}_{\mathrm{A}r}-k_{\varphi }{V}_{\mathrm{A}\varphi }cot\theta +\mathbf{k}·\frac{\partial {\mathbf{V}}_{\mathrm{A}}}{\partial \theta })\hfill \\ \hfill & +\frac{N^2{\kappa }_r}{r\omega }({\kappa }_{\theta }{\kappa }_{\perp }^2+{\kappa }_{\varphi }^2{\kappa }_{r}cot\theta ),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& \frac{\mathrm{d}{k}_{\varphi }}{\mathrm{d}t}=-\frac{{\omega }_{\mathrm{A}}}{\omega }\frac{k_{\varphi }}{r}(V_{\mathrm{A}r}+V_{\mathrm{A}\theta }cot\theta )+\frac{N^2{\kappa }_{\varphi }{\kappa }_r}{r\omega }({\kappa }_{\perp }^2-{\kappa }_{\theta }{\kappa }_{r}cot\theta ).\hfill \end{array}$$(A.6)

In these equations, κ = (κ_r,κ_θ,κ_ϕ) = k/k and ${\kappa }_{\perp }=\sqrt{{\kappa }_{\theta }^2+{\kappa }_{\varphi }^2}$.

A.3. Correction of the wavevector

As described in Sect. 2.4, we adjusted the radial component of the wavevector k_r at each timestep to conserve ω during the ray tracing simulation. This was achieved by finding the roots of the following equation which can be derived from the dispersion relation (Loi 2020a):

$$\begin{array}{c}\hfill 0=b_1{k}_r^4+b_2{k}_r^3+b_3{k}_r^2+b_4{k}_r+b_5,\end{array}$$(A.7)

with

$$\begin{array}{cc}\hfill b_1& =V_{\mathrm{A}r}^2,\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill b_2& =2V_{\mathrm{A}r}(k_{\theta }{V}_{\mathrm{A}\theta }+k_{\varphi }{V}_{\mathrm{A}\varphi }),\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill b_3& =V_{\mathrm{A}r}^2{k}_{\perp }^2+{(k_{\theta }{V}_{\mathrm{A}\theta }+k_{\varphi }{V}_{\mathrm{A}\varphi })}^2-{\omega }^2,\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill b_4& =2V_{\mathrm{A}r}{k}_{\perp }^2(k_{\theta }{V}_{\mathrm{A}\theta }+k_{\varphi }{V}_{\mathrm{A}\varphi }),\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill b_5& =k_{\perp }^2[{(k_{\theta }{V}_{\mathrm{A}\theta }+k_{\varphi }{V}_{\mathrm{A}\varphi })}^2+N^2-{\omega }^2].\hfill \end{array}$$(A.12)

Loi (2020a) used a Newton-Raphson iteration scheme to find the roots of Eq. A.7. They performed a maximum of five iterations after each timestep during the integration of each ray. If the scheme has not converged after five iterations, they set the value of k_r to the current estimate of the root and continued the integration of the ray. In this case, ω deviates from its initial value. Thus, the value of ω after each timestep can be used as a measure of accuracy.

Instead of using a Newton-Raphson scheme, we directly calculated the roots and selected the one with the real part closest to the current estimate of k_r as the new value of k_r (see Sect. 2.4 for details). Since we have discarded the imaginary part, ω can also be used as a measure of accuracy for our method. We show a comparison of a set of 1200 rays using the Newton-Raphson method to correct k_r with a set of rays where we used our method to correct k_r in Appendix A.5. Overall, we find that calculating the roots directly improves the stability of the integration and reduces the uncertainties.

A.4. Ray classification

Since it is not feasible to check each ray manually, we used a ray classification scheme to decide whether a ray is trapped, reflected, or unclassified. It consists of a series of conditions that are checked in sequence. If the ray fulfilled a condition, it was classified accordingly and the other steps were ignored. If the ray did not meet the condition, the next one was checked and so on. The ray classification scheme used in this work is an updated version of the one used by Loi (2020a). It works as follows:

1. The value of ω is different than the original value by more than 0.0001 % at any point during the integration. → Unclassified

2. The ray crosses r_f going outward.

   1. The final value of k_r is equal or larger than 1000 times the initial value of k_r⁷. → Unclassified

   2. The final value of k_r is smaller than 1000 times the initial value of k_r. → Reflected

3. The final value of k_r exceeds the initial value k_r by a factor of 10.

   1. The initial value of k_r is not always smaller than the value of k_r during the second half of the integration⁸. → Unclassified

   2. The initial value of k_r is always smaller than the values of k_r during the second half of the integration. → Trapped

4. None of the above conditions is fulfilled. → Unclassified

We have introduced additional quality checks to the ray classification scheme of Loi (2020a) before identifying a ray as reflected or trapped (see Fig. A.2). This was done because there were rays with clear numerical artifacts that were classified as either reflected or trapped by the classification scheme of Loi (2020a). The majority of these rays were identified as unclassified by our scheme. Due to the unpredictability of the numerical artifacts, we cannot exclude that there still is a small number of rays that exhibit unphysical behavior and were classified as either trapped or reflected. This may lead to a slight underestimation of the uncertainties of the trapped fraction, but should not have a significant impact on the results presented in this work. Any ray that was identified as unclassified according to the classification scheme of Loi (2020a) was also unclassified according to our scheme.

Fig. A.2. Sets of 1200 rays plotted as as unispheres parameterized by (θ, ϕ) = (θ0, α) for the model 2A with Bcen = Bcrit, cen, ω = 10 ωdyn, and l = 1. Trapped rays are shown as blue dots, reflected rays are shown as red diamonds, and unclassified rays are shown as gray squares. In the top row, we show rays for which we have used a Newton-Raphson iteration scheme for the correction of kr, which was used by Loi (2020a). In the bottom row, we show rays where we have used our method described in Sect. 2.4 to correct kr. In the left column, we employed the ray classification scheme of Loi (2020a) to decide whether a ray is trapped, reflected, or unclassified. In the right column, we used our ray classification scheme described in Appendix A.4. The trapped fraction fT of each set is indicated under its corresponding unisphere.

A.5. Comparison of the ray tracing setups

In Fig. A.2 we compare a set of 1200 rays using a Newton-Raphson iteration scheme to correct k_r with a set of rays where we used our method to correct k_r (i.e., calculating the roots of Eq. A.7 directly; see Sect. 2.4 and Appendix A.3). For both sets of rays, we show the ray classification according to the scheme of Loi (2020a), as well as according to our ray classification scheme (see Appendix A.4). Comparing the two classification schemes, the number of unclassified rays is higher when using our classification scheme than when using that of Loi (2020a). Furthermore, we find that a small number of rays that clearly exhibit numerical artifacts were not identified as unclassified when using the classification of Loi (2020a). Our classification scheme seems to be more reliable in finding these numerical artifacts. We therefore conclude that the uncertainties of the trapped fraction (Eq. 12) might be underestimated when using the classification scheme of Loi (2020a).

The comparison of the two methods for correcting k_r shows that there is more scatter in the behavior of the rays when using the Newton-Raphson iteration than when using our method to determine the roots of Eq. A.7. This can be recognized by the fact that in the sets of rays corresponding to the Newton-Raphson iteration, there is a small number of trapped rays that are completely surrounded by reflected rays and vice versa. There are no such rays when using our method to correct k_r. Nevertheless, the trapped fraction estimated using the Newton-Raphson iteration is not significantly different from the trapped fraction determined using our method when our ray classification scheme is utilized. In this case however, the uncertainties of the trapped fraction are an order of magnitude smaller when using our method (i.e., calculating the roots of Eq. A.7 directly), which is why we applied it in this work.

A.6. Comparison to Loi (2020a) and Loi (2020b)

A.6.1. Comparison to Loi (2020a)

The numerical integration of the magneto-gravity ray trajectories performed in this work is based on the study of Loi (2020a). We used the same set of differential equations (Eqs. A.1 to A.6), as well as the same ray initialization and a similar forward Runge-Kutta integration scheme. Therefore, we expect to obtain results consistent with Loi (2020a) for our models 2A, 2B, and 2C, which are reproductions of their models A, B, and C. While the trapped fraction qualitatively shows the same trends with B_cen, ω, and l as reported by Loi (2020a), we systematically obtain smaller trapped fractions for the same B_cen. If for example B_cen = B_crit, cen, we obtain a trapped fraction of roughly 50 %, while Loi (2020a) obtains trapped fractions ≳ 75 %. We find that neither the MESA profiles with increased spacial resolution (see Appendix A.1) nor the changes to the correction of k_r during the integration (see Appendix A.3) are responsible for this discrepancy. The difference in the critical field strengths of our model 2A and model A of Loi (2020a) (see Table 1) does not make up for this either.

Since the discrepancy in the trapped fraction is present for all models and combinations of B_cen, ω, and l, it must be caused by a systematic difference of the behavior of the individual rays in our implementation and that of Loi (2020a). In Figs. A.3 and A.4, we show reproductions of Figures 4 and 5 of Loi (2020a) using our implementation of the magneto-gravity ray tracing. In addition to the magnetic field strength used by Loi (2020a) (i.e., B_cen = 0.3 B_crit, cen), we also show the behavior of the ray for B_cen = 1.2 B_crit, cen, for which the ray in Fig. A.3 undergoes the same number of reflections as the ray shown in Figure 4 of Loi (2020a).

In Fig. A.3, we show that we need a magnetic field strength that is about four times greater than that of Loi (2020a) so that our ray undergoes the same number of reflections as the ray of Loi (2020a). Likewise, our ray shown in Fig. A.4 is not trapped when using the field strength of Loi (2020a), although they find that the ray is trapped. Using a field strength that is four times greater, the ray shown in Fig. A.4 is trapped, but its behavior is still significantly different to that of the ray shown in Figure 5 of Loi (2020a). This suggests that the discrepancy cannot be explained by an offset in B_crit, cen.

Figure 4 of Loi (2020a) shows a peculiarity that indicates that there are some systematic differences between their and our numerical integration of the rays. In the third panel where Loi (2020a) shows the evolution of the wavevector components in time, the maximum absolute value of k_ϕ is significantly larger than the maximum absolute value of k_r. During the propagation of our rays, the maximum absolute value of k_ϕ is always comparable or smaller than the maximum absolute value of k_r (see bottom panel of Figs. A.3 and A.4).

In Fig. A.5, we show a reproduction of Figure 6 of Loi (2020a). Here, the internal magnetic field has been turned off (i.e., B_cen = 0). In their Figure 6, the maximum absolute value of k_ϕ is larger than the maximum absolute value of k_r, similar as in their Figure 4. The same is true for the maximum absolute value of k_θ. This is never the case with our rays (see bottom panel of Fig. A.5). Since the ray tracing equations (Eqs. A.1 to A.6) depend on the current value of k_θ and k_ϕ, different values of k_θ and k_ϕ can influence the local behavior and thus the trajectory of the ray, as well as its eventual fate (i.e., trapped vs. reflected).

Fig. A.3. Reproduction of the lower two panels of Figure 4 of Loi (2020a) using our implementation of the ray tracing described in Sect. 2 for the model 2A. The limits of the axes are also similar to Figure 4 of Loi (2020a). In the left column, we show the behavior of the ray with Bcen = 0.3 Bcrit, cen, which is the value used by Loi (2020a). In the right column, we show how the ray behaves with Bcen = 1.2 Bcrit, cen. This value has been selected such that the ray undergoes the same number of reflections as the ray shown in Figure 4 of Loi (2020a). The remaining initial parameters are identical to Loi (2020a) (i.e., θ0 = 99°, α = 40.5°, ω = 10 ωdyn, and l = 1).

Fig. A.4. Reproduction of the lower two panels of Figure 5 of Loi (2020a) using the same setup as for Fig. A.3. Here, θ0 = 141° and α = 58.5°. The remaining initial parameters are the same as for Fig. A.3.

Fig. A.5. Reproduction of the lower two panels of Figure 6 of Loi (2020a) for the ray shown in Fig. A.4, but without an internal magnetic field (i.e., Bcen = 0).

A.6.2. Comparison to Loi (2020b)

The magneto-gravity ray tracing method presented by Loi (2020a) has also been used by Loi (2020b). It is a follow-up work of Loi (2020a) where the author investigates the effect of an internal magnetic field on the period spacing. The only difference to the ray tracing analysis presented by Loi (2020a) is that Loi (2020b) performed the integration of the rays on an analytic background instead of a finite grid. Since we want to verify our procedure, which is based on MESA profiles consisting of a finite radial grid, we projected the analytic expressions of Loi (2020b) for N and ρ on a radial logarithmic grid with 5000 grid cells. We used the resulting grids as mock MESA profiles to compare the behavior of our rays to Loi (2020b).

In Figs. A.6 and A.7, we show reproductions of Figures 4 and 5 of Loi (2020b) using our implementation of the magneto-gravity ray tracing. Since θ₀ and α are not stated by Loi (2020b), we chose them manually. Both Figs. A.6 and A.7 have a striking resemblance to Figures 4 and 5 of Loi (2020b) respectively. In particular, the maximum absolute values of k_θ and k_ϕ are always comparable or smaller than the maximum absolute value of k_r. This is in contrast to the figures shown by Loi (2020a).

Fig. A.6. Reproduction of Figure 4 of Loi (2020b) using our implementation of the ray tracing described in Sect. 2 but with synthetic N and ρ profiles (see Appendix A.6 for details). Since θ0 and α are not stated by Loi (2020b), they have been chosen manually (θ0 = 11.0° and α = 81.0°). The remaining initial parameters are identical to Loi (2020b) (i.e., r0 = 0.04 R*, rf = 0.01 R*, Bcen = Bcrit, cen, ω = 10 ωdyn, and l = 1).

Fig. A.7. Reproduction of Figure 5 of Loi (2020b) for the ray shown in Fig. A.6, but without an internal magnetic field (i.e., Bcen = 0).

In conclusion, the trajectories of our rays are consistent with the most recent implementation of the author of Loi (2020a). We hypothesize that the difference in the trapped fraction of Loi (2020a) and this work can be explained by the systematically different behavior of k_θ and k_ϕ. The reason for this discrepancy in the behavior of the wavevector components remains elusive.

Appendix B: Additional figures

B.1. Trapped fractions of the remaining models

In Figs. B.1 and B.2, we show the trapped fractions for the models 1A, 1B, and 1C, as well as 1a, 1b, 2a, and 2b. They follow the same trends as the trapped fractions shown in Figs. 3 and 4. In addition, we shown in Fig. B.3 the trapped fractions corresponding to a spherical degree of l = 10 and 100 for model 1A. They also follow the same trends.

Fig. B.1. Same as Fig. 4, now for the models 1A, 1B, and 1C.

Fig. B.2. Same as Fig. 4, now for the models 1a, 1b, 2a, and 2b.

Fig. B.3. Same as Fig. 4, now for the model 1A with ω = 10 ωdyn and l = 10 and 100.

B.2. Evolution of the trapped fraction

In Fig. B.4, we show the evolution of the trapped fraction for the four evolutionary tracks presented in Fig. 1 considering the ratio of the strength of the internal magnetic field and the critical field strength at the hydrogen-burning shell.

Fig. B.4. Same as Fig. 6, now considering the strength of the radial component of the internal magnetic field at the hydrogen-burning shell. Here, we consider the ratio BH / Bcrit, H instead of Bcen / Bcrit, cen when using Eq. 13.

Appendix C: MESA inlist

The MESA inlist used to generate the evolutionary tracks shown in Fig. 1 is as follows:

&star_job

create_pre_main_sequence_model=.false.

! .true. for PMS model

load_saved_model = .true.

! .false. for PMS model

load_model_filename = ’2M_ms.mod’

! load MS or PMS model

save_model_when_terminate = .true.

save_model_filename = ’2M_rgb.mod’

! chose model name

pgstar_flag=.false.

/

MYAMPeos

/

MYAMPkap

Zbase = 0.02 ! 0.005 for low Z star

/

MYAMPcontrols

initial_mass=2 ! 1 for 1 M_sun star

initial_z=0.02 ! 0.005 for low Z star

log_directory=’LOGS_RGB’

! chose output directory name

set_min_D_mix = .true.

min_D_mix = 5d-1

energy_eqn_option=’dedt’

use_gold_tolerances=.true.

calculate_Brunt_N2=.true.

terminal_interval = 10

do_history_file = .true.

history_interval = 1

write_profiles_flag = .true.

profile_interval = 50

mesh_delta_coeff=0.3d0

varcontrol_target=1d-5

min_allowed_varcontrol_target=1d-7

max_years_for_timestep=1.d7

/&

Note that the above inlist has no stopping condition. For the pre-MS models, we used the following stopping condition:

Lnuc_div_L_zams_limit=0.99d0

stop_near_zams=.true.

For the MS models, we used this stopping condition:

xa_central_lower_limit_species(1)=’h1’

xa_central_lower_limit(1)=1d-3

For the RGB models, we adjusted the run_star_extras.f90 file such that the computation stops when ν_max is smaller than 10 μHz. We made the MESA inlist and the run_star_extras.f90 file used in this work publicly available on Zenodo (doi: 10.5281/zenodo.14924586).

All Tables

All Figures

Fig. 1. Hertzsprung-Russell diagram of the four evolutionary tracks with different stellar masses and metallicities described in Sect. 2. The evolution of the MS is shown as a dotted line and the evolution after the MS as a solid line. We do not show the full evolution after the MS. The locations of the stellar models that served as the background for the ray tracing analysis are indicated by the various symbols. Star symbols indicate that the central field strength Bcen, the frequency ω, and the spherical degree l were varied during the ray tracing analysis of a particular model. Full dots indicate that Bcen and ω were varied and circles indicate that only Bcen was varied.

In the text

Fig. 2. Trajectory of a reflected and trapped ray (top row) and the corresponding wavevector components as a function of time (bottom row) for the model 2A with Bcen = Bcrit, cen, ω = 10 ωdyn, and l = 1. The reflected ray with θ0 = 128.9° and α = 3.7° is shown on the left, the trapped ray with θ0 = 98.7° and α = 90.9° is shown on the right. The dotted black circle in the top row represents the maximum extend of the magnetic field rf. Note that the limits of the axes do not coincide between left and right for the sake of clarity.

In the text

Fig. 3. Trapped fraction as a function of the central field strength of the internal magnetic field for the models 2A (top), 2B (middle), and 2C (bottom). Different colors indicate different combinations of the frequency ω and the spherical degree l (see legend). The dotted lines mark the corresponding central critical field strength. The error bars are determined by the uncertainty of fT (see Eq. (12)). The line connecting the symbols is a fit to the trapped fraction presented in Sect. 3.2.

In the text

	Fig. 4. Same as Fig. 3, now as a function of the central field strength of the internal magnetic field divided by the corresponding central critical field strength. The upper limit of the horizontal axis is now 10. Note that Bcrit, cen is different for each model, as well as for each combination of ω and l (see Eq. (6)).
In the text

Fig. 5. Top: Normalized critical field strength as a function of radius. We also show the maximum absolute value of the three components of the Prendergast field (see Sect. 2.2), which is scaled to Bcen = Bcrit, cen. Furthermore, we show the radial component of an improvised magnetic field that has a plateau and scales with the critical field strength at the hydrogen-burning shell Bcrit, H. The strength of the improvised field at the hydrogen-burning shell is set to BH = Bcrit, H. The hydrogen-burning shell is indicated by the vertical dotted line. Bottom: Trapped fraction as a function of the central field strength or the field strength at the hydrogen-burning shell for model 1A with ω = 10 ωdyn and l = 1. The black symbols correspond to the case where all components of the Prendergast field are taken into account. They are the same as the black symbols in the top panel of Fig. B.1. The blue symbols indicate that only the radial component of the Prendergast field has been used and the other components were set to zero. The red symbols correspond to magnetic field featuring the plateau that is scaled by Bcrit, H instead of Bcrit, cen. The dotted line marks either Bcrit, cen or Bcrit, H. The line connecting the symbols is the fit to the trapped fraction (Eq. (13)).

In the text

Fig. 6. Trapped fraction evaluated with Eq. (13) at ω = 2πνmax, where νmax is the frequency of maximum oscillation power, as a function of the evolution of the star indicated by νmax (stars evolve from right to left). Each panel features one of the four evolutionary tracks presented in Fig. 1. Rows show tracks with the same stellar mass and columns show tracks with the same metallicity. The blue (red) lines correspond to the dipole (quadrupole) modes. Different linestyles (solid, dashed, dotted) represent different initial strengths of the radial component of the magnetic field at the stellar center, which corresponds to the field strength the end of the MS. Due to the conservation of magnetic flux, the magnetic field strength increases as the stars evolve. The small knees that occur at the points where fT becomes greater than 0, which are visible in some of the curves, are caused by the limited temporal resolution of the evolutionary tracks. We show the evolution of the trapped fraction considering BH/Bcrit, H instead of Bcen/Bcrit, cen in Fig. B.4.

In the text

Fig. 7. Background-corrected PSD of KIC 6975038 obtained from Kepler data using the background model of Kallinger et al. (2014) (for details, see Coppée et al. 2024). Radial modes are indicated by the gray shaded area, dipole modes by the blue area, and quadrupole modes by the reddish area. In addition, we show the trapped fraction as a function of frequency computed for a stellar model with M* = 1.3 M⊙, Zinit = 0.02, νmax ≈ 128 μHz, and Bcen ≈ 6.23 MG (see Sect. 4.3) for both the dipole modes (in blue) and quadrupole modes (in red).

In the text

Fig. 8. Average frequency shift of the dipole modes caused by various internal magnetic field configurations as a function of the magnetic field strength for three stellar models. The field strengths have been chosen such that Bcen = 0.7 Bcrit, cen and BH = 0.7 Bcrit, H. The critical field strengths were calculated at ω = 2πνmax with l = 1. Colors correspond to different stellar models and symbols indicate that a different magnetic field configuration has been used. Average frequency shifts corresponding to the azimuthal order m = 0 (m = ±1) are shown as solid symbols (outlines). The field configurations and the calculation of the average frequency shifts are described in Sect. 4.3 and in the notes under Table 3, which also summarizes the data shown in this figure. Note that the horizontal axis is broken.

In the text

	Fig. A.1. Profiles of the buoyancy frequency as a function of radius for the models 2A, 2B, and 2C. The launch radius of the rays r0 is shown in red for each model, the maximum extend of the internal magnetic field rf is shown as the dashed black line.
In the text

Fig. A.2. Sets of 1200 rays plotted as as unispheres parameterized by (θ, ϕ) = (θ0, α) for the model 2A with Bcen = Bcrit, cen, ω = 10 ωdyn, and l = 1. Trapped rays are shown as blue dots, reflected rays are shown as red diamonds, and unclassified rays are shown as gray squares. In the top row, we show rays for which we have used a Newton-Raphson iteration scheme for the correction of kr, which was used by Loi (2020a). In the bottom row, we show rays where we have used our method described in Sect. 2.4 to correct kr. In the left column, we employed the ray classification scheme of Loi (2020a) to decide whether a ray is trapped, reflected, or unclassified. In the right column, we used our ray classification scheme described in Appendix A.4. The trapped fraction fT of each set is indicated under its corresponding unisphere.

In the text

Fig. A.3. Reproduction of the lower two panels of Figure 4 of Loi (2020a) using our implementation of the ray tracing described in Sect. 2 for the model 2A. The limits of the axes are also similar to Figure 4 of Loi (2020a). In the left column, we show the behavior of the ray with Bcen = 0.3 Bcrit, cen, which is the value used by Loi (2020a). In the right column, we show how the ray behaves with Bcen = 1.2 Bcrit, cen. This value has been selected such that the ray undergoes the same number of reflections as the ray shown in Figure 4 of Loi (2020a). The remaining initial parameters are identical to Loi (2020a) (i.e., θ0 = 99°, α = 40.5°, ω = 10 ωdyn, and l = 1).

In the text

	Fig. A.4. Reproduction of the lower two panels of Figure 5 of Loi (2020a) using the same setup as for Fig. A.3. Here, θ0 = 141° and α = 58.5°. The remaining initial parameters are the same as for Fig. A.3.
In the text

	Fig. A.5. Reproduction of the lower two panels of Figure 6 of Loi (2020a) for the ray shown in Fig. A.4, but without an internal magnetic field (i.e., Bcen = 0).
In the text

Fig. A.6. Reproduction of Figure 4 of Loi (2020b) using our implementation of the ray tracing described in Sect. 2 but with synthetic N and ρ profiles (see Appendix A.6 for details). Since θ0 and α are not stated by Loi (2020b), they have been chosen manually (θ0 = 11.0° and α = 81.0°). The remaining initial parameters are identical to Loi (2020b) (i.e., r0 = 0.04 R*, rf = 0.01 R*, Bcen = Bcrit, cen, ω = 10 ωdyn, and l = 1).

In the text

	Fig. A.7. Reproduction of Figure 5 of Loi (2020b) for the ray shown in Fig. A.6, but without an internal magnetic field (i.e., Bcen = 0).
In the text

	Fig. B.1. Same as Fig. 4, now for the models 1A, 1B, and 1C.
In the text

	Fig. B.2. Same as Fig. 4, now for the models 1a, 1b, 2a, and 2b.
In the text

	Fig. B.3. Same as Fig. 4, now for the model 1A with ω = 10 ωdyn and l = 10 and 100.
In the text

	Fig. B.4. Same as Fig. 6, now considering the strength of the radial component of the internal magnetic field at the hydrogen-burning shell. Here, we consider the ratio BH / Bcrit, H instead of Bcen / Bcrit, cen when using Eq. 13.
In the text