OPEN ACCESS

<figure><img src="image_1.png" /></figure>

CrossMark

# Effect of Rotation on Wave Mixing in Intermediate-mass Stars

A. Varghese<fnref n="1" />,<fnref n="2" /> ID, R. P. Ratnasingam<fnref n="1" /> ID, R. Vanon<fnref n="3" /> ID, P. V. F. Edelmann<fnref n="4" /> ID, S. Mathis<fnref n="2" />, and T. M. Rogers<fnref n="1" />,<fnref n="5" /> ID

Received 2023 November 9; revised 2024 May 24; accepted 2024 June 4; published 2024 July 23

## Abstract

Internal gravity waves are likely to cause mixing in stellar interiors. Studies show that the mixing by these waves changes drastically across age and mass. Here, we study the effect of rotation on this wave mixing by considering a 7 $M_\odot$ model at zero-age main sequence and mid-main sequence. We compare the mixing profiles at a range of rotation rates ($1 \times 10^{-5}$, $2 \times 10^{-5}$, $3 \times 10^{-5}$, $4 \times 10^{-5}$, and $1 \times 10^{-4}$ rad s$^{-1}$) and observe that the mixing decreases with decreasing Rossby number. This can be attributed to the effect of rotation on convection, which influences the amplitude with which the waves are excited near the convective–radiative interface.

Unified Astronomy Thesaurus concepts: Stellar evolution (1599); Hydrodynamics (1963)

## 1. Introduction

Internal gravity waves (IGWs) are waves propagating in stably stratified fluids with gravity as the restoring force. These are the waves responsible for the quasibennial oscillation (Baldwin et al. 2001) in Earth's atmosphere and for the turbulent mixing in oceans (Munk & Wunsch 1998). In stars, these waves are stochastically excited either by convective plumes penetrating into the stable stratified region (Hurlburt et al. 1986; Montalbán & Schatzman 2000; Pinçon et al. 2016) or from the Reynolds stress in the bulk of the convection zone (Kumar et al. 1999; Belkacem et al. 2009; Samadi et al. 2010; Lecoanet & Quataert 2013). The amplitudes of these waves are affected by a number of factors such as the density stratification, Brunt–Väisälä frequency, thermal damping, and geometric effects (Ratnasingam et al. 2018).

Studies suggest that, along with other instabilities, IGWs can transport angular momentum and cause mixing in the radiative regions of stars. IGWs were proposed to solve the solar neutrino discrepancy (Press 1981) and explain the Li gap in F-type stars and the observed Li abundance in the Sun (Garcia Lopez & Spruit 1991; Montalban 1994; Montalban & Schatzman 1996). Rogers & McElwaine (2017) showed that the mixing by IGWs could be treated as a diffusive process with an amplitude that is proportional to the square of wave amplitudes in the radiation zone. Following their work, Varghese et al. (2023) extended the analysis to stars of different masses and ages. They showed that the mixing is stronger in more massive stars and decreases as stars age.

Recent studies demonstrated that the excitation of IGWs is strongly influenced by rotation. The action of the Coriolis acceleration on these waves modifies the damping and the spatial structure of these waves, resulting in the appearance of new types of waves (Mathis et al. 2008, 2014). One among these are the gravito-inertial waves (GIWs) where the lower frequencies are significantly affected under the restoring action

of buoyancy and Coriolis force (Berthomieu et al. 1978; Lee & Saio 1997). These waves are excited near the convective–radiative interface, similar to IGWs, either through small-scale eddies or by plumes with the coupling between the waves and the turbulence strongly influenced by the Coriolis acceleration (Mathis et al. 2014). GIWs are strongly coupled with the turbulence when they are in the subinertial regime ($\omega < 2\Omega$), remaining as propagative inertial waves in the convection zone, and are weakly coupled in the super-inertial regime ($\omega > 2\Omega$) when they become evanescent (Mathis et al. 2008; Mathis 2009) in the convection zone. Rotation, therefore, modifies the coupling between the waves and the turbulent convective flows, subsequently influencing the excitation of gravity waves and GIWs in the stellar radiative interiors (Mathis et al. 2014).

With asteroseismology probing the internal dynamics of stars (Van Reeth et al. 2015, 2016; Pápics et al. 2017; Szewczuk et al. 2021; Breton et al. 2022), gravito-inertial modes (Mathis 2009; Dintrans & Rieutord 2000; Ballot et al. 2010) have been detected (Neiner et al. 2012, 2020), leading to a new path for exploring stellar interiors. Excitation of these modes is considered to be a possible explanation of the variability of certain rapidly rotating stars (Balona et al. 1996; Dintrans et al. 1999) as the pure-gravity modes and rotation alone could not explain the observations (Brunsden et al. 2018; Antoci et al. 2019; Ouazzani et al. 2019). Detection of these modes illustrates the need to include Coriolis acceleration in stellar modeling (Aerts 2021; Mombarg et al. 2021; Pedersen et al. 2021) and provides asteroseismic calibrations of various parameters such as the overshoot parameter. Gravito-inertial asteroseismology can give us information related to convection, buoyancy, and rotation (Aerts 2021) and therefore, is expected to provide better constraints for numerical simulations.

The influence of rotation, along with that of IGWs, on mixing has been studied previously by Talon & Charbonnel (2005) and Charbonnel & Talon (2007) where they explained the surface Li abundance in low-mass stars. They computed 1D stellar models by incorporating the angular momentum transport by meridional circulation and shear turbulence following the diffusive/advective formalism (Chaboyer &

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Newcastle University, Newcastle Upon Tyne, UK</footnote>

<footnote n="2">²Université Paris-Saclay, Université de Paris, Sorbonne Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³University of Leeds, Leeds, UK</footnote>

<footnote n="4">⁴Los Alamos National Laboratory, Los Alamos, NM, USA</footnote>

<footnote n="5">⁵Planetary Science Institute, Tucson, AZ, USA</footnote>
Zahn 1992; Zahn 1992; Maeder & Zahn 1998) and IGWs. Their model considered the wave excitation by stochastic eddies similar to the work of Kumar et al. (1999). They showed that the shear layer developed near the convection zone due to the deposition of angular momentum acts as a wave filter (Gough & McIntyre 1998) and replaced the excited wave spectrum together with the action of shear-layer oscillation by a filtered spectrum and diffusion coefficient. Inclusion of the angular momentum transport by waves along with the meridional circulations and rotational instabilities in the 1D models by these studies suggested a decrease in the mixing compared with that of rotation alone. These studies have taken into account the indirect influence of the angular momentum transport on the chemical mixing, meaning that IGWs modify the shear that modifies the chemical mixing. Our work focuses solely on the direct chemical mixing induced by IGWs in the stellar interiors, considering stars with a convective core and radiative envelope without taking into account the role of angular momentum transport by these waves and the subsequent indirect chemical mixing induced.

Following the work of Rogers & McElwaine (2017) and Varghese et al. (2023), we study the effect of rotation on wave mixing by considering a $7 M_\odot$ model at zero-age main sequence (ZAMS) and mid-main sequence (midMS). We do not include TAMS models because, as described in Ratnasingam et al. (2020) and Varghese et al. (2023), these waves are highly attenuated. We achieve this by running 2D simulations with different rotation rates ($\Omega$) using a background reference model from Modules of Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019). Section 2 introduces the numerical techniques used for obtaining the mixing profiles. Section 3 presents our findings, and Section 4 discuss our conclusions.

## 2. Background and Numerical Setup

### 2.1. Two-dimensional Hydrodynamical and Tracer-particle Simulations

We conducted 2D simulations to study IGWs in stellar interiors and coupled them with tracer-particle simulations to determine the mixing by these waves. We generated the background reference models from MESA for a nonrotating $7 M_\odot$ star at ZAMS (core hydrogen mass fraction, $X_c = 0.70$) and at midMS ($X_c = 0.35$) with the metallicity set to $Z = 0.02$. We set all stellar parameters similar to those in Ratnasingam et al. (2023) and Varghese et al. (2023). The inlists used to generate the models are available in Zenodo doi:10.5281/zenodo.2596370. The simulations solve the Navier-Stokes equations in the anelastic approximation by considering an equatorial slice of the star with stress-free, isothermal, and impermeable boundary conditions similar to that of Rogers et al. (2013) and Ratnasingam et al. (2023).

The equations are given by

$$ \nabla \cdot (\bar{\rho} \mathbf{v}) = 0, \qquad (1) $$

$$ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\nabla\left(\frac{P}{\bar{\rho}}\right) - C\bar{g}\hat{r} + 2(\mathbf{v} \times \mathbf{\Omega}) \\ + \bar{\nu}\left(\nabla^2\mathbf{v} + \frac{1}{3}\nabla(\nabla \cdot \mathbf{v})\right), \qquad (2) $$

$$ \frac{\partial T}{\partial t} + (\mathbf{v} \cdot \nabla)T = -v_r \left( \frac{\partial \bar{T}}{\partial r} - (\gamma - 1) \bar{T} h_{\rho} \right) \\ + (\gamma - 1) T h_{\rho} v_r + \frac{1}{c_v \bar{\rho}} \nabla \cdot (c_p \bar{k} \bar{\rho} \nabla T) \\ + \frac{1}{c_v \bar{\rho}} \nabla \cdot (c_p \bar{k} \bar{\rho} \nabla \bar{T}), \qquad (3) $$

where $\bar{\rho}$ and $\bar{k}$ are the reference state density and thermal diffusivity, respectively.

Equation (1) is the mass conservation equation in the anelastic approximation with $\mathbf{v}$ as the fluid velocity, while Equation (2) represents the momentum conservation equation with $\mathbf{\Omega}$ as the rotation rate and $\bar{\nu}$ as the kinematic viscosity. $P$ is the reduced pressure defined by Braginsky & Roberts (1995) as

$$ P = \frac{p}{\bar{\rho}} + U, \qquad (4) $$

where $U$ is the gravitational potential perturbation, which is neglected in our simulations (Cowling 1941), and $p$ is pressure perturbation. The codensity perturbation $C$ as defined in Rogers & Glatzmaier (2005) is

$$ C = -\frac{1}{\bar{T}} \left( T + \frac{1}{g\bar{\rho}} \frac{d\bar{T}}{dz} p \right). \qquad (5) $$

Equation (3) is the energy-conservation equation written as a temperature equation, where $v_r$ is the radial velocity and $c_v$ is the specific heat capacity at constant volume. The first term in the right-hand side of Equation (3) represents the super- or subadiabaticity (Rogers et al. 2013), which drives the convection in our simulation. It is set to a constant positive value in the convection zone and is calculated from the 1D MESA model in the radiation zone.

Equations (1)–(3) are solved using a Fourier decomposition method in the azimuthal direction ($\theta$) and a finite difference scheme on a nonuniform grid in the radial direction ($r$). The variables are updated using the Adams-Bashforth explicit method for the nonlinear terms and Crank-Nicolson implicit time-stepping method for the linear terms.

To maintain numerical stability, the simulation domain is cut off at 90% of the total stellar radius, where the density drops beyond ~6 orders of magnitude from the stellar center as shown in Figure 1. We also maintain constant thermal ($\kappa$) and viscous ($\nu$) diffusivities at $5 \times 10^{12} \text{ cm}^2 \text{ s}^{-1}$ in all our simulations and thus the Prandtl number $Pr$,

$$ Pr = \frac{\nu}{\kappa}, \qquad (6) $$

is equal to 1 throughout the simulation domain. We have considered this value, which is much higher than the actual value in the stellar interior, to ensure numerical stability. Recent work by Vanon et al. (2023) could achieve a mild improvement on the Pr (Pr ~ 5–15 in the convection zone and ~0.02–0.4 in the radiative envelope) in their three-dimensional simulations (see also Breton et al. 2022).

We then set the initial solid body rotation rate as $1 \times 10^{-5}$, $2 \times 10^{-5}$, $3 \times 10^{-5}$, $4 \times 10^{-5}$, and $1 \times 10^{-4} \text{ rad s}^{-1}$, which in terms of critical rotation rate, $v_c = \sqrt{GM/R^3}$, are 0.035$v_c$, 0.07$v_c$, 0.1$v_c$, 0.14$v_c$, and 0.35$v_c$ for the ZAMS models and 0.069$v_c$, 0.13$v_c$, 0.20$v_c$, 0.27$v_c$ and 0.69$v_c$ for the midMS
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Density and Brunt–Väisälä frequency as a function of fractional radius for the 7 $M_\odot$ at ZAMS (purple) and midMS (cyan). The vertical dashed line in the density plot indicates the cut off radius of the simulation domain.</figcaption></figure>

<table><caption>**Table 1**
The Total Time Interval the Simulations Are Run in Terms of Convective Turnover Times for ZAMS and midMS Models</caption><thead><tr><th rowspan="2">Ω<br>(rad s<sup>-1</sup>)</th><th colspan="2">Convective Turnovers</th></tr><tr><th>ZAMS</th><th>midMS</th></tr></thead><tbody><tr><td>1 × 10<sup>-5</sup></td><td>97</td><td>133</td></tr><tr><td>2 × 10<sup>-5</sup></td><td>107</td><td>119</td></tr><tr><td>3 × 10<sup>-5</sup></td><td>93</td><td>110</td></tr><tr><td>4 × 10<sup>-5</sup></td><td>103</td><td>74</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>115</td><td>44<sup>a</sup></td></tr></tbody></table>

**Note.**

models to study the effect of rotation on waves.<fnref n="6" /> We run all the simulations to a total of $4.5 \times 10^7$ s, which is approximately 225 wave crossing times of a 5 µHz wave (see Table 1 for time interval in terms of convective turnover times). Studies by Rogers et al. (2013) determined that there was no significant variations between the models with an initial differential rotation and the ones with an initial solid body rotation. Hence, we also expect the diffusion coefficients calculated in this work to have a very small dependence on whether the model has an initial differential rotation or solid body rotation. Figure 2 shows the root-mean-squared velocity averaged over the convection zone as a function of time for ZAMS and midMS models. We note that the $v_{rms}$ attained a steady state from $t = 1 \times 10^7$ s (indicated as vertical dashed lines in Figure 2), and therefore we chose the velocity data saved at regular time intervals from this value for our further analysis.

To determine the mixing by IGWs, we introduced $N$ tracer particles into our simulations and tracked their trajectories over a time $T$. We calculated the diffusion coefficient, $D$, based on the equations given by Rogers & McElwaine (2017),

$$D(r, \tau) = \frac{Q(r, \tau)}{2\tau n(r, \tau)} - \frac{P(r, \tau)^2}{2\tau n(r, \tau)^2}, \quad (7)$$

where $n(r, \tau)$ is the number of subtrajectories starting at $r$ with a duration of $\tau$; $P(r, \tau)$ is the sum of the lengths of the subtrajectories; and $Q(r, \tau)$ is the sum of the square of these subtrajectories. More details on the calculation of $n$, $P$, $Q$, and $D(r, \tau)$ can be found in Rogers & McElwaine (2017).

We plotted the diffusion coefficients as a function of radius at different time differences,<fnref n="7" /> $\tau$, and chose a time difference ($\tau = 2 \times 10^7$ s) such that it has contributions from a sufficient number of time steps and the amplitude of the profiles are converged (Figure 3). We then used this profile for all our further analysis.

# 3. Results

## 3.1. Mixing Profiles at Different Rotations

Figure 3 shows the radial diffusion profiles for the 7 $M_\odot$ model at ZAMS (left) and midMS (right) for rotation rates $\Omega = 1 \times 10^{-5}$ (orange), $2 \times 10^{-5}$ (blue), $3 \times 10^{-5}$ (pink), $4 \times 10^{-5}$ (green), and $1 \times 10^{-4}$ (yellow) rad s$^{-1}$. We note that all the profiles follow the same trend, which is increasing toward the surface as we move from the convective–radiative interface as expected from Rogers & McElwaine (2017). This increase can be explained by the change in the wave amplitude due to the balance between decreasing density stratification and the increasing thermal damping toward the stellar surface. However, there are three notable features in the plot. (i) There is a particular drop in the amplitude of mixing profile (almost an order of magnitude) as we move from $3 \times 10^{-5}$ to $4 \times 10^{-5}$ in the midMS model. (ii) The difference between the amount of mixing at the slowest rotation and the fastest rotating model is different in ZAMS and midMS. The zoomed-in plot in Figure 3(a) gives a clearer picture of the trend in the ZAMS model. In ZAMS, we see a modest change of a factor of approximately 3, while at midMS we see a substantial change of more than an order of magnitude. (iii) We note that the overall mixing decreases with increasing rotation rate for both ages.

Before looking in detail at these features, we compared the diffusion coefficients at a variety of radii as a function of the

<footnote>a The simulation is run until the $v_{rms}$ attained steady state as shown in Figure 2.</footnote>

<footnote n="6">⁶ This kind of rotation likely distorts the star such that our cylindrical geometry is not appropriate.</footnote>

<footnote n="7">⁷ The smaller $\tau$ has contributions from more time steps compared to a larger value of $\tau$. As an example, consider $T = 100$ with time step = 1. Then, for $\tau = 5$, it has contributions from 5–0, 6–1, 7–2, and so on, whereas $\tau = 97$ can result only from 100–3, 99–2, and 98–1.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2. $v_{rms}$ averaged over the convection zone for (a) ZAMS and (b) midMS at different angular velocity values. The vertical dashed line indicates the time from which we chose the data for tracer-particle simulations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Diffusion coefficient as a function of fractional radius for $7 M_{\odot}$ at (a) ZAMS and (b) midMS models at $\Omega = 1 \times 10^{-5}$ (orange), $2 \times 10^{-5}$ (blue), $3 \times 10^{-5}$ (pink), $4 \times 10^{-5}$ (light green), and $1 \times 10^{-4}$ (yellow) rad s$^{-1}$. The vertical dashed line represents the convective–radiative interface for each model.</figcaption></figure>

convective Rossby number, $R_o$ defined as

$$R_o = \frac{v_{rms}}{2\Omega L}, \qquad (8)$$

where $v_{rms}$ is the averaged root-mean-squared velocity in the convection zone, $L$ is the extent of the convection zone, and $\Omega$ is the rotational frequency. Figure 4 shows the diffusion coefficients as a function of $R_o$ at different radii for all the models studied. We note that there is a clear break in the midMS models at $R_o \sim 0.04$, below which where we see a steep decline in the diffusion coefficient. Regardless of the age of the star or the radius chosen, we see that the trend remains the same, that is, the diffusion coefficient increases with increasing Rossby number. To understand this observed trend, i.e., the wave-induced mixing diminishing with rotation, we will study the variation of the damping rate and of the

excitation source as a function of rotation in the following sections.

## 3.2. Linear Theory

In this section, we discuss the theoretical prescription from Rogers & McElwaine (2017), where they determined the diffusion coefficient $D$,

$$D = A v_{wave}^2(\omega, l, r), \qquad (9)$$

with the coefficient$^8 A \sim 1$ s (Rogers & McElwaine 2017; Varghese et al. 2023) and the wave amplitude calculated using the linear theory without considering rotation from

<footnote>$^8$ The constant $A$ is unknown and likely depends on simulation parameters such as the total time domain considered. We expect that this parameter can be constrained using asteroseismology (Pedersen et al. 2021).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Diffusion coefficient as a function of Rossby number for a 7 M☉ model at (a) ZAMS and (b) midMS for $\Omega = 1 \times 10^{-5}$, $2 \times 10^{-5}$, $3 \times 10^{-5}$, $4 \times 10^{-5}$, and $1 \times 10^{-4}$ rad s⁻¹ at three different radii.</figcaption></figure>

Ratnasingam et al. (2018),

$$v_{\text{wave}}(\omega, l, r) = v_0(\omega, l, r) \left(\frac{\rho}{\rho_0}\right)^{-\frac{1}{2}} \left(\frac{r}{r_0}\right)^{-1} \times \left(\frac{N^2 - \omega^2}{N_0^2 - \omega^2}\right)^{-\frac{1}{4}} e^{-\frac{\tau}{2}}, \quad (10)$$

where $\rho_0$, $r_0$, and $N_0$ are the density, radius, and the Brunt–Väisälä Frequency at the initial reference point. $v_0(\omega, l, r)$ is the initial wave amplitude for a given frequency and wavenumber. The damping coefficient $\mathcal{T}$ as given by Kumar et al. (1999) and Zahn et al. (1997) is expressed as

$$\mathcal{T}(\omega, l, r) = \int_{r_0}^{r} \frac{16\sigma T^3}{3\rho^2\kappa c_p} \left( \frac{l(l+1)^{\frac{3}{2}}N^3}{r^3\omega^4} \right) \left( 1 - \frac{\omega^2}{N^2} \right)^{\frac{1}{2}} dr, \quad (11)$$

where $\sigma$, $T$, $c_p$, $\kappa$, and $l$ are the Stefan–Boltzmann constant, temperature, specific heat capacity at constant pressure, opacity, and wavenumber.

Using these equations Varghese et al. (2023) found that the dominant waves contributing to the mixing profiles are low-frequency waves within the range of 4–9 µHz, with 5 µHz being the dominant wave for the 7 M☉ ZAMS and midMS nonrotating models. As discussed in Ratnasingam et al. (2018), these are the waves generated with amplitudes large enough to escape the thermal damping but low enough that they are efficiently generated by convection. While many frequencies contribute to the mixing profile, the aim of identifying the dominant waves contributing the mixing profiles was to make the inclusion of wave mixing in 1D stellar evolution models easier.

### 3.3. Linear Theory Including Rotation

To theoretically determine the effect of rotation on waves, we followed the approach of Press (1981) and Ratnasingam et al. (2020) in solving the linearized hydrodynamic equations in the anelastic approximation but with the inclusion of rotation term in the momentum equation (Equation (13)). The linearized equations neglecting thermal and viscous diffusivities are given

below:

$$\nabla \cdot (\bar{\rho}v) = 0, \quad (12)$$

$$\frac{\partial v}{\partial t} = -\nabla \left(\frac{P}{\bar{\rho}}\right) - C\bar{g}\hat{r} + 2(v \times \Omega), \quad (13)$$

and

$$\frac{\partial T}{\partial t} = -v_r \left( \frac{\partial \bar{T}}{\partial r} - (\gamma - 1) \bar{T} h_{\rho} \right). \quad (14)$$

We reduce the above to get the following second-order differential equation for wave propagation considering a wave ansatz, $v_r(r, \theta, z) \propto v_r(r)e^{im\theta}e^{-i\omega t}$,

$$0 = \frac{\partial^2 \alpha}{\partial r^2} + \left[ \frac{N^2}{\omega^2} - 1 \right] \frac{m^2}{r^2} \alpha + \frac{1}{2} \left[ \frac{\partial h_{\rho}}{\partial r} - h_{\rho}^2 + \frac{h_{\rho}}{r} \right] \alpha + \frac{1}{4r^2} \alpha - \frac{2\Omega h_{\rho} m}{\omega r} \alpha. \quad (15)$$

Here, $\alpha = v_r \bar{\rho}^{\frac{1}{2}} r^{\frac{3}{2}}; m$ is the 2D Fourier basis wavenumber; $\omega$ is the angular frequency; and $h_\rho$ is the inverse-density scale height. In this work, we introduce the term including the rotational effects, $2\Omega h_\rho m/\omega r$. Here, we look at the combined effect of radiative damping, density stratification, and rotation.

### 3.4. Radiative Damping and Rotation

Earlier studies of Press (1981) and Kumar et al. (1999) define the wave-damping opacity as

$$\mathcal{T}(\omega, l, r) = \int_{r_0}^{r} \frac{\gamma[\omega, l, r]}{|v_g[\omega, l, r]|} dr, \quad (16)$$

where the damping rate $\gamma$ is given as

$$\gamma(\omega, l, r) = K k_r^2, \quad (17)$$

with the radial wavenumber $k_r^2$,

$$k_r^2 = \left( \frac{N^2}{\omega^2} - 1 \right) k_h^2. \quad (18)$$
The vertical group velocity $v_g$ and the thermal diffusivity $K$ are given as

$$v_g(r) = \frac{\partial \omega}{\partial k_r} = - \frac{(N^2 - \omega^2)^{\frac{1}{2}} \omega^2}{k_h N^2}, \qquad (19)$$

$$K(r) = \frac{16\sigma T^3}{3\rho^2\kappa c_p}. \qquad (20)$$

In the absence of rotation, the horizontal wavenumber $k_h$ is given as

$$k_h^2 = \frac{l(l+1)}{r^2}, \qquad (21)$$

and the Equation (16) simplifies to Equation (11).

From Equation (15), the modified radial wavenumber $k_{rr}$ in the presence of rotation is given by

$$k_{rr}^2 = \left( \frac{N^2}{\omega^2} - 1 \right) k_h^2 - \frac{2\Omega h_\rho k_h}{\omega}, \qquad (22)$$

where $k_h = m/r$. The equation including the effect of rotation on the vertical group velocity with rotation $v_{grr}$ is given as

$$v_{grr} = \frac{\partial \omega}{\partial k_{rr}} = - \frac{k_{rr}}{\frac{N^2 k_h^2}{\omega^3} - \frac{\Omega h_\rho k_h}{\omega^2}}. \qquad (23)$$

Therefore, the damping rate is modified by rotation and is given as $\gamma_{rr}$:

$$\gamma_{rr}(\omega, l, r) = K k_{rr}^2, \qquad (24)$$

thereby influencing the spatial damping opacity in rotating models, $T_{rr}$:

$$T_{rr}(\omega, l, r) = \int_{r_0}^{r} \frac{\gamma_{rr}[\omega, l, r]}{|v_{grr}[\omega, l, r]|} dr, \qquad (25)$$

with

$$T_{rr}(\omega, l, r) = \int_{r_0}^{r} \frac{16\sigma T^3}{3\rho^2\kappa c_p} \left( \left[ 1 - \frac{\omega^2}{N^2} \right] k_h^2 - \frac{2\Omega h_\rho \omega k_h}{N^2} \right)^{1/2} \times \left( \frac{N^3 k_h^2}{\omega^4} - \frac{\Omega h_\rho N k_h}{\omega^3} \right). \qquad (26)$$

The wave amplitude given by Equation (10) can be rewritten as

$$v_{\text{wave}_{rr}}(\omega, l, r) = v_0(\omega, l, r) \left(\frac{\rho}{\rho_0}\right)^{-\frac{1}{2}} \left(\frac{r}{r_0}\right)^{-1} \times \left(\frac{(N^2 - \omega^2)m^2 - 2\Omega h_\rho m\omega r}{(N_0^2 - \omega^2)m^2 - 2\Omega h_\rho m\omega r_0}\right)^{-\frac{1}{4}} e^{-\frac{\tau_{rr}}{2}}. \qquad (27)$$

We now compare the damping experienced by the waves at different rotation rates for the ZAMS and midMS models. This is carried out by plotting the attenuation coefficient ($\exp-(T_{rr}/2)$) calculated with Equation (26) as a function of the fractional radius for a 5 µHz wave at $l=1$ shown in Figure 5.<fnref n="9" /> We observe an overall decrease in the attenuation factor with rotation. However, we find that this decrease in the attenuation factor with rotation is not sufficient enough to

explain the difference seen in the mixing profiles (Figure 3). Hence, we will now evaluate the influence of rotation on convection and on the subsequent wave excitation in the following sections. We begin by comparing the theoretical diffusion coefficient calculated using Equations (26), (27), and (9) with the simulation profiles.

### 3.5. Theoretical Mixing Profiles

We calculated the theoretical diffusion profiles using the Equations (26), (27), and (9) for each frequency and wavenumber. Considering the following relation between the tangential and radial velocity based on the assumption that the frequencies are much smaller than the Brunt–Väisälä frequency

$$\frac{v_0}{v_r} = \frac{N}{\omega}, \qquad (28)$$

we set $v_0(\omega, l, r)$ to the root-mean-squared velocity at a radius just outside the convection zone averaged over the total time domain multiplied by a factor of $\omega/N$. The wave launch point is chosen such that the Brunt–Väisälä frequency is greater than 10 µHz since the dominant waves found from Varghese et al. (2023) were in the range of 4–9 µHz for ZAMS and midMS models at $l=1$. Figure 6 shows the theoretical profiles calculated for a 5 µHz wave for the ZAMS and midMS models. We observe that the theoretical profiles do not follow a trend similar to that of the simulation in the case of ZAMS models. The profiles from the theory follows a similar trend as that of the simulation in midMS models, with a particular drop in the amplitude of the mixing profiles as we move from $3 \times 10^{-5}$ to $4 \times 10^{-5}$ rad s$^{-1}$. However, the amplitude difference is only ~2.6 compared to an order of magnitude difference in the simulation. We know that the initial velocity $v_0$ is set by the convection, but it is still unclear how to define this value. Hence, to get a deeper understanding on the choice of the initial velocity and of the difference observed between the different rotation rates, we look in detail at the influence of rotation on convection and the theory from Augustson et al. (2020), which also includes further effects of rotation on wave generation and propagation.

### 3.6. Influence of Rotation on Convection

Figure 2 shows an overall decrease in the rms velocity averaged over the convection zone with rotation. Chandrasekhar (1961) demonstrated that the rotation influences the onset of turbulent convection. In addition, Stevenson (1979) predicted that the mode of convection that carries more heat is inhibited by the action of rotation building a mixing-length theory for rotating convection. Barker et al. (2014) verified this theoretical prediction by conducting high-resolution nonlinear numerical simulations in a Cartesian box assuming the Boussinesq approximation. Later, Augustson & Mathis (2019) generalized the initial study by Stevenson (1979) by taking into account viscous and heat diffusions and confirmed the results obtained by Stevenson (1979). Hence, the decrease in the convective velocities noted in Figure 2 can be attributed to the influence of rotation on convection.

A notable feature in Figure 3 is the large decrease in the amplitude of the mixing profile in the radiation zone from $\Omega = 3 \times 10^{-5}$ to $4 \times 10^{-5}$ rad s$^{-1}$ in the case of the midMS

<footnote n="9">⁹ We chose $l=1$ in accordance with Varghese et al. (2023), where they had shown that the higher wavenumbers experience higher thermal damping as expected from Equation (11).</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 5. The attenuation coefficient ($\exp(-(T_{rr}/2))$) as a function of fractional radius for the 7 $M_\odot$ at (a) ZAMS and (b) midMS at different angular velocity values considering a 5 µHz wave at $l=1$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Theoretical diffusion coefficient calculated using Equations (9), (26), and (27) for (a) ZAMS and (b) midMS for different $\Omega$ as a function of fractional radius for a 5 µHz wave. The vertical dashed line indicates the convective–radiative interface.</figcaption></figure>

model. We see a decrease of approximately an order of magnitude compared to the other models (approximately a factor of 3 from $\Omega = 1 \times 10^{-5}$ to $2 \times 10^{-5}$ rad s$^{-1}$ and 2 from $\Omega = 4 \times 10^{-5}$ to $1 \times 10^{-4}$ rad s$^{-1}$). This can be attributed to the variation of rms velocity in the convection zone. Even though we note an overall decrease in the rms velocity averaged over the convection zone with rotation, this decrease is very evident for $\Omega = 4 \times 10^{-5}$ at midMS as seen from Figure 2. Under the action of rotation, the critical Rayleigh number ($Ra_c$) is given as

$$Ra_c \propto Ek^{-\frac{4}{3}} \qquad (29)$$

where the Ekman number, Ek is given as

$$Ek = \frac{\nu}{2\Omega L}. \qquad (30)$$

Hence, $Ra_c$ increases with rotation, implying that at higher rotation, stronger forcing is necessary for driving the convection to be in the supercritical regime compared to that of a

slowly rotating model (Chandrasekhar 1961; Plumley & Julien 2019). This can also be noted from the values of the Reynolds number Re given in Table 2. Here Re is defined as

$$Re = \frac{v_{rms}L}{\nu}, \qquad (31)$$

where we recall that $L$ is the radial extent of the convection zone. We see that Re remains within 750–950 for the ZAMS model, whereas there is a larger variation, a 30% decrease, between $\Omega = 3 \times 10^{-5}$ and $4 \times 10^{-5}$ rad s$^{-1}$, compared to 7% decrease between $\Omega = 2 \times 10^{-5}$ and $3 \times 10^{-5}$ rad s$^{-1}$, and a 10% decrease between $\Omega = 1 \times 10^{-5}$ and $2 \times 10^{-5}$ rad s$^{-1}$ as we move from the slow-rotating to fast-rotating model in midMS.

We therefore infer from the above that the difference in the rms velocity observed is due to the convection being inhibited at higher rotation. Consequently, higher rotation rates influence the convection significantly, such that the waves are generated
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Interfacial flux calculated using Equation (32) for (a) ZAMS and (b) midMS for different Ω as a function of frequency.</figcaption></figure>

<table><caption>**Table 2**
Reynolds Number, Re, and the Averaged $v_{rms}$ in the Convection Zone for All the Models Studied</caption><thead><tr><th>Model</th><th>Ω<br>(rad s<sup>-1</sup>)</th><th>Re</th><th>Average v<sub>rms</sub><br>(cm s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="5">7 M<sub>☉</sub> ZAMS</td><td>1 × 10<sup>-5</sup></td><td>813</td><td>93,988</td></tr><tr><td>2 × 10<sup>-5</sup></td><td>891</td><td>103,048</td></tr><tr><td>3 × 10<sup>-5</sup></td><td>779</td><td>90,049</td></tr><tr><td>4 × 10<sup>-5</sup></td><td>862</td><td>99,642</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>956</td><td>110,548</td></tr><tr><td rowspan="5">7 M<sub>☉</sub> midMS</td><td>1 × 10<sup>-5</sup></td><td>921</td><td>116,810</td></tr><tr><td>2 × 10<sup>-5</sup></td><td>822</td><td>104,189</td></tr><tr><td>3 × 10<sup>-5</sup></td><td>764</td><td>96,811</td></tr><tr><td>4 × 10<sup>-5</sup></td><td>518</td><td>65,671</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>306</td><td>38,792</td></tr></tbody></table>

with a lower amplitude (Takehiro et al. 2020), which in turn explains a decrease in the mixing by these waves.

The significant difference observed in the Re values in the midMS models compared to that of the ZAMS can thus also explain the considerable variation in the mixing profiles in midMS as we move from the slow to fast-rotating models. We see a larger variation across midMS but not in ZAMS, suggesting that the convection needs to be forced even stronger in midMS to be at an equivalent Re. The evident decrease observed in the rms velocities of midMS models (Table 2) must influence the amplitude with which the waves are generated at the interface and the mixing they trigger. To confirm this, we now consider how the kinetic energy of the turbulent convection is transmitted to IGWs at the convective-core boundary as a function of rotation.

### 3.7. Influence of Rotation on Wave Excitation

In addition to the effect on convection, rotation affects the wave excitation as shown in Augustson et al. (2020). To evaluate this quantitatively for these simulations, we calculated the flux of the GIWs for all the models studied in this work. We adopted the approach of Augustson et al. (2020) adapted to our equatorial geometry. The flux was calculated by considering

Equation (40) in Augustson et al. (2020) for wave excitation by pressure fluctuations at the convection–radiation interface and applying it to the peculiar case of the equator. This leads to

$$F_r = \frac{1}{2} \rho v_{\text{int}}^3 \frac{\omega}{(N^2 + 4\Omega^2 - \omega^2)^{1/2}} \quad (32)$$

The above equation relates the wave energy flux in the radial direction to the convective velocity at the convective–radiative interface measured in our simulations, $v_{\text{int}}$, which is the root-mean-squared velocity averaged over the simulation time domain; the wave frequency, $\omega$; and the rotation rate, $\Omega$. Figure 7 shows the wave flux as a function of frequencies at each rotation rate for both the ZAMS and midMS models. We observe that the interfacial wave flux decreases with rotation, agreeing with the results of Augustson et al. (2020). We chose the wave launch point such that the Brunt–Väisälä frequency is greater than 10 µHz since the dominant waves found from Varghese et al. (2023) were in the range of 4–9 µHz for ZAMS and midMS models as in Section 3.5. The decrease in the wave energy flux seen in Figure 7 implies that the transport of chemicals induced by the waves in the stellar interior should decrease with rotation.

Our aim now is thus to compare the diffusion coefficients ($D_{\text{flux}}$) calculated from the wave flux with that of the simulation. Considering Equation (35) from Augustson et al. (2020) adapted to our equatorial geometry, we obtain the general relation between the radial flux of energy and the radial components of GIW velocity:

$$F_r = \frac{\rho}{2k_h} (N^2 + 4\Omega^2 - \omega^2)^{1/2} v_w^2 \quad (33)$$

We thus obtain the wave velocity at the interface, $v_w$, as

$$v_w^2 = \frac{2mF_r}{r\rho(N^2 + 4\Omega^2 - \omega^2)^{1/2}} \quad (34)$$

We now substitute $v_w$ as $v_0$ in Equation (27) and calculate the diffusion coefficient using Equation (9) for each frequency and wavenumber. Figure 8 shows the calculated diffusion coefficient for a 5 µHz wave. We see that the difference between
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Diffusion coefficient calculated from the interfacial flux for a 5 µHz wave for (a) ZAMS and (b) midMS for different Ω as a function of fractional radii. The vertical dashed line indicates the convective–radiative interface.</figcaption></figure>

these theoretical profiles at different rotation rates agrees reasonably well with that of the difference between the simulation mixing profiles shown in Figure 3 for models up to rotation rate Ω = 4 × 10⁻⁵ rad s⁻¹. The difference is approximately a factor of 1.6 from Ω = 1 × 10⁻⁵ to 2 × 10⁻⁵ rad s⁻¹, 1.3 from Ω = 2 × 10⁻⁵ to 3 × 10⁻⁵ rad s⁻¹, and 6.3 from Ω = 3 × 10⁻⁵ to 4 × 10⁻⁵ rad s⁻¹ in midMS. The theoretical profiles do not match with that of the simulation at the surface. This is because the radial velocity is forced to zero at the surface in our numerical simulation, causing a drop in the amplitude of the simulation profile near the surface.

On the other hand, we see a difference of ~10 in ZAMS and ~12 in midMS between the theoretical mixing profiles of 4 × 10⁻⁵ rad s⁻¹ and 1 × 10⁻⁴ rad s⁻¹ compared to a difference of ~2 seen in the simulation profiles (Figure 3) for both the models. We propose it could be due to the fact that we are comparing the theoretical profiles of a single frequency with the simulation profiles, which could have contributions from a number of frequencies. At higher rotation, the frequencies contributing to the mixing profiles could be significantly different from that of the slowly rotating model, causing the observed difference. Moreover, earlier studies suggest that the rotation reduces the overshoot depth in stellar models (Augustson & Mathis 2019; Korre & Featherstone 2021). This could influence the radius at which the waves are excited, particularly at high rotation. Furthermore, the considered theoretical model of GIW excitation by pressure fluctuations at the convection–radiation interface is a simplified one when compared to those of IGW excitation by turbulent plumes (e.g., Schatzman 1993; Pinçon et al. 2016). We attribute these factors to the difference noted between the theoretical and simulation profiles of the fastest rotating model. Overall, we can conclude that the decreased mixing we see in our simulations is due to the effect of rotation on reducing the convective motions and hence, reducing wave excitation and propagation.

**4. Conclusions**

We studied the influence of rotation on mixing induced by IGWs by considering a 7 M☉ model at ZAMS and midMS.

Overall, we noted that the mixing decreases with increase in rotation irrespective of the age considered. We attribute this to the influence of rotation on convection and the subsequent wave excitation at the convective–radiative interface and the wave propagation in the radiation zone.

The averaged rms velocity in the convection zone decreases with increasing rotation, agreeing with the studies of Augustson & Mathis (2019) and Augustson et al. (2020), particularly for midMS models. We argue that this modification of the convective velocity with rotation is the reason for a less efficient excitation of GIWs and as a consequence the reduced wave-induced mixing. We noted that the decrease in velocity is higher as we move from Ω = 3 × 10⁻⁵ to 4 × 10⁻⁵ rad s⁻¹ in the midMS model. We explain this decrease by the increasing critical Rayleigh number with rotation, resulting in a dramatically reduced Re. This is coherent with the reduction of the rms velocity of the convective mode that carries more heat predicted by the rotating mixing-length theory derived by Stevenson (1979) and Augustson & Mathis (2019), which has been confirmed in local Cartesian high-resolution nonlinear simulations computed by Barker et al. (2014).

While damping in the radiative regions is increased for increasing rotation, in real stars, this may be less important than the role rotation plays on reducing convective velocities and hence wave driving. However, we expect the damping to play a significant role near the surface of stars, particularly at later stages of evolution, where the thermal diffusivities are found to be much higher. Hence, considering the influence of rotation on convection in stars could provide a better constraint to the amplitude with which the IGWs are generated at the interface and the subsequent mixing caused by these waves in the radiation zone. This maybe of great importance for predicting the internal mixing in rapidly rotating intermediate-mass and massive stars (Pedersen et al. 2021) and in pre-main-sequence late-type stars (Charbonnel et al. 2013), which are known to rotate faster than our Sun (Gallet & Bouvier 2015).

The implementation of these mixing profiles in a 1D stellar evolution code can result in significant changes across ages. Computations of 1D models including time-dependent IGW mixing for different models presented in Varghese et al. (2023)
are currently ongoing. Preliminary results show that inclusion of the mixing by IGWs influences the composition gradient, resulting in a more smoothed Brunt-Väisälä frequency peak compared to the model with a constant envelope mixing that was used (J. Mombarg 2024, private communication)

## Acknowledgments

We thank the anonymous reviewer for the helpful comments. We acknowledge support from STFC grant ST/L005549/1 and NASA grant NNX17AB92G. Computing was carried out on Rocket High Performance Computing service at Newcastle University and DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk), funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. P.V.F.E. was supported by the U.S. Department of Energy through the Los Alamos National Laboratory (LANL). LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001). This work has been assigned a document release number LA-23-32363. S.M. acknowledges support from the European Research Council through HORIZON ERC SyG grant 4D-STAR101071505, from the CNES SOHO-GOLF and PLATO grants at CEA-DAP, and from PNPS (CNRS/INSU). R.V. was funded by STFC grant ST/W000873/1.

## ORCID iDs

A. Varghese https://orcid.org/0000-0002-0761-4545

R. P. Ratnasingam https://orcid.org/0000-0002-7250-6524

R. Vanon https://orcid.org/0000-0001-5747-8476

P. V. F. Edelmann https://orcid.org/0000-0001-7019-9578

T. M. Rogers https://orcid.org/0000-0002-2306-1362

## References

Aerts, C. 2021, *RvMP*, 93, 015001

Antoci, V., Cunha, M. S., Bowman, D. M., et al. 2019, *MNRAS*, 490, 4040

Augustson, K. C., & Mathis, S. 2019, *ApJ*, 874, 83

Augustson, K. C., Mathis, S., & Astoul, A. 2020, *ApJ*, 903, 90

Baldwin, M. P., Gray, L. J., Dunkerton, T. J., et al. 2001, *RvGeo*, 39, 179

Ballot, J., Lignières, F., Reese, D. R., & Rieutord, M. 2010, *A&A*, 518, A30

Balona, L. A., Böhm, T., Foing, B. H., et al. 1996, *MNRAS*, 281, 1315

Barker, A. J., Dempsey, A. M., & Lithwick, Y. 2014, *ApJ*, 791, 13

Belkacem, K., Samadi, R., Goupil, M. J., et al. 2009, *A&A*, 494, 191

Berthomieu, G., Gonczi, G., Graff, P., Provost, J., & Rocca, A. 1978, *A&A*, 70, 597

Braginsky, S. I., & Roberts, P. H. 1995, *GApFD*, 79, 1

Breton, S. N., Brun, A. S., & García, R. A. 2022, *A&A*, 667, A43

Brunsden, E., Pollard, K. R., Wright, D. J., De Cat, P., & Cottrell, P. L. 2018, *MNRAS*, 475, 3813

Chaboyer, B., & Zahn, J. P. 1992, *A&A*, 253, 173

Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability, International Series of Monographs on Physics (Oxford: Clarendon)

Charbonnel, C., Decressin, T., Amard, L., Palacios, A., & Talon, S. 2013, *A&A*, 554, A40

Charbonnel, C., & Talon, S. 2007, in AIP Conf. Ser. 948, Unsolved Problems in Stellar Physics: A Conference in Honor of Douglas Gough, ed. R. J. Stancliffe et al. (Melville, NY: AIP), 15

Cowling, T. G. 1941, *MNRAS*, 101, 367

Dintrans, B., & Rieutord, M. 2000, *A&A*, 354, 86

Dintrans, B., Rieutord, M., & Valdettaro, L. 1999, *JFM*, 398, 271

Gallet, F., & Bouvier, J. 2015, *A&A*, 577, A98

Garcia Lopez, R. J., & Spruit, H. C. 1991, *ApJ*, 377, 268

Gough, D. O., & McIntyre, M. E. 1998, *Natur*, 394, 755

Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, *ApJ*, 311, 563

Korre, L., & Featherstone, N. A. 2021, *ApJ*, 923, 52

Kumar, P., Talon, S., & Zahn, J.-P. 1999, *ApJ*, 520, 859

Lecoanet, D., & Quataert, E. 2013, *MNRAS*, 430, 2363

Lee, U., & Saio, H. 1997, *ApJ*, 491, 839

Maeder, A., & Zahn, J.-P. 1998, *A&A*, 334, 1000

Mathis, S. 2009, *A&A*, 506, 811

Mathis, S., Neiner, C., & Tran Minh, N. 2014, *A&A*, 565, A47

Mathis, S., Talon, S., Pantillon, F. P., & Zahn, J. P. 2008, *SoPh*, 251, 101

Mombarg, J. S. G., Van Reeth, T., & Aerts, C. 2021, *A&A*, 650, A58

Montalban, J. 1994, *A&A*, 281, 421

Montalban, J., & Schatzman, E. 1996, *A&A*, 305, 513

Montalbán, J., & Schatzman, E. 2000, *A&A*, 354, 943

Munk, W., & Wunsch, C. 1998, Deep Sea Research Part I: Oceanographic Research Papers, 45, 1977

Neiner, C., Floquet, M., Samadi, R., et al. 2012, *A&A*, 546, A47

Neiner, C., Lee, U., Mathis, S., et al. 2020, *A&A*, 644, A9

Ouazzani, R. M., Marques, J. P., Goupil, M. J., et al. 2019, *A&A*, 626, A121

Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, *A&A*, 598, A74

Paxton, B., Bildsten, L., Dotter, A., et al. 2011, *ApJS*, 192, 3

Paxton, B., Cantiello, M., Arras, P., et al. 2013, *ApJS*, 208, 4

Paxton, B., Marchant, P., Schwab, J., et al. 2015, *ApJS*, 220, 15

Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, *ApJS*, 234, 34

Paxton, B., Smolec, R., Schwab, J., et al. 2019, *ApJS*, 243, 10

Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, *NatAs*, 5, 715

Pinçon, C., Belkacem, K., & Goupil, M. J. 2016, *A&A*, 588, A122

Plumley, M., & Julien, K. 2019, *E&SS*, 6, 1580

Press, W. H. 1981, *ApJ*, 245, 286

Ratnasingam, R. P., Edelmann, P. V. F., & Rogers, T. M. 2018, *MNRAS*, 482, 5500

Ratnasingam, R. P., Edelmann, P. V. F., & Rogers, T. M. 2020, *MNRAS*, 497, 4231

Ratnasingam, R. P., Rogers, T. M., Chowdhury, S., et al. 2023, *A&A*, 674, A134

Rogers, T. M., & Glatzmaier, G. A. 2005, *MNRAS*, 364, 1135

Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, *ApJ*, 772, 21

Rogers, T. M., & McElwaine, J. N. 2017, *ApJL*, 848, L1

Samadi, R., Belkacem, K., Goupil, M. J., et al. 2010, *Ap&SS*, 328, 253

Schatzman, E. 1993, *A&A*, 279, 431

Stevenson, D. J. 1979, *GApFD*, 12, 139

Szewczuk, W., Walczak, P., & Daszyńska-Daszkiewicz, J. 2021, *MNRAS*, 503, 5894

Takehiro, S.-i., Brun, A. S., & Yamada, M. 2020, *ApJ*, 893, 83

Talon, S., & Charbonnel, C. 2005, *A&A*, 440, 981

Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, *A&A*, 593, A120

Van Reeth, T., Tkachenko, A., Aerts, C., et al. 2015, *ApJS*, 218, 27

Vanon, R., Edelmann, P. V. F., Ratnasingam, R. P., Varghese, A., & Rogers, T. M. 2023, *ApJ*, 954, 171

Varghese, A., Ratnasingam, R. P., Vanon, R., Edelmann, P. V. F., & Rogers, T. M. 2023, *ApJ*, 942, 53

Zahn, J. P. 1992, *A&A*, 265, 115

Zahn, J. P., Talon, S., & Matias, J. 1997, *A&A*, 322, 320