Magnetohydrodynamic Simulations of the Tayler Instability in Rotating Stellar Interiors

The Tayler instability is an important but poorly studied magnetohydrodynamic instability that likely operates in stellar interiors. The nonlinear saturation of the Tayler instability is poorly understood and has crucial consequences for dynamo action and angular momentum transport in radiative regions of stars. We perform three-dimensional MHD simulations of the Tayler instability in a cylindrical geometry, including strong buoyancy and Coriolis forces as appropriate for its operation in realistic rotating stars. The linear growth of the instability is characterized by a predominantly $m=1$ oscillation with growth rates roughly following analytical expectations. The non-linear saturation of the instability appears to be caused by secondary shear instabilities and is also accompanied by a morphological change of the flow. We argue, however, that non-linear saturation likely occurs via other mechanisms in real stars where the separation of scales is larger than those reached by our simulations. We also observe dynamo action via the amplification of the axisymmetric poloidal magnetic field, suggesting that Tayler instability could be important for magnetic field generation and angular momentum transport in the radiative regions of evolving stars.


INTRODUCTION
The interplay between rotation and magnetism is crucial for understanding the evolution of stars and the compact objects they produce. Differential rotation generated by contracting stellar cores may source various magnetohydrodynamic (MHD) instabilities that can amplify magnetic fields and/or transport AM outwards to slow the rotation of the stellar core. However, the instabilities at work and their saturation mechanisms remain poorly understood.
Asteroseismic observations have helped by providing internal rotation rate measurements for stars on the main sequence (Kurtz et al. 2014;Saio et al. 2015;Benomar et al. 2015;Van Reeth et al. 2018), sub-giant/red giant branch (RGB) (Beck et al. 2012;Mosser et al. 2012;Deheuvels et al. 2014;Triana et al. 2017;Gehan et al. 2018), red clump (Mosser et al. 2012;Deheuvels et al. 2015), and finally in WD remnants (Hermes et al. 2017). The conclusion drawn from these measurements is unambiguous: core rotation rates are relatively slow, and the vast majority of AM is extracted from stellar cores as they evolve. The spin rates of red giant cores and WDs are slower than theoretically predicted by nearly all hydrodynamic AM transport mechanisms (Fuller et al. 2014;Belkacem et al. 2015;Spada et al. 2016;Eggenberger et al. 2017;Ouazzani et al. 2019).
The non-axisymmetric MHD Tayler instability (Tayler 1973;Spruit 1999;Goldstein et al. 2019) is likely the most important MHD instability in radiative regions of stars. Tayler instability is a kinktype instability of toroidal (azimuthal) fields that have been created by winding up a radial seed field through differential rotation. Above a critical field strength, field loops slip sideways relative to one another with a predominantly non-axisymmetric = 1 wavenumber. While the dynamics of the linear instability are well understood, the nonlinear (and likely turbulent) saturation of the instability, and the resulting AM transport are poorly understood and controversial (e.g., Braithwaite 2006;Zahn et al. 2007).
The Tayler-Spruit (TS) dynamo (Spruit 2002) is one possible saturation mechanism of the Tayler instability. In this theory, toroidal magnetic field energy is turbulently dissipated by the fluid motions, ★ suoqing@shao.ac.cn and the instability saturates when the turbulent dissipation rate is equal to the energy input via winding of the radial field. In the presence of a composition gradient, the resulting torque density due to Maxwell stresses is where = ln Ω/ ln is the dimensionless radial shear, and eff is the effective stratification, which is usually nearly equal to the compositional stratification in post-MS stars. The TS dynamo has been implemented into many stellar evolution codes, but it predicts much faster core rotation than observed in post-MS stars ).
However, Fuller et al. (2019) argued that Spruit (2002) overestimated the energy damping rate of the instability, because only magnetic energy in the disordered (perturbed) field can be turbulently damped. By calculating an energy damping rate due to weak magnetic turbulence, Fuller et al. (2019) argued the instability can grow to larger amplitudes, producing a larger Maxwell stress in its saturated state. Fuller et al. (2019) find the Tayler torques are where is a saturation parameter of order unity. The different scaling is very important because Ω/ eff ∼ 10 −4 in RGB stars, so the prescription of Fuller et al. (2019) allows for significantly more AM transport.
Because the saturation of the Tayler instability is a complex and nonlinear process, it is important to examine this process via numerical simulations. The Tayler instability has been seen in a few simulations, but only in limited configurations not including both rotation and realistic stratification. Weber et al. (2015) and Gellert et al. (2008) used = 0 (i.e., no stratification) and Guerrero et al. (2019) used Ω = 0 (i.e., no rotation). The first simulation of the Tayler instability with shear and buoyancy (Braithwaite 2006) was compressible, limiting the dynamic range and time scale over which simulations could be performed. Those simulations used Ω/ = 1, in stark contrast to the values of Ω/ ∼ 10 −4 1 expected in real stars. The analytic predictions of Spruit (2002) and Fuller et al. (2019) also assume Ω , so it is important to simulate that parameter regime. Recently, Petitdemange et al. (2023) presented a suite of simulations of Tayler instability, finding apparent agreement with the prediction of Eq. (1), which we discuss further in Section 4.
In this work, we perform three-dimensional MHD simulations of the Tayler instability, including both stratification and rotation. We also vary dimensionless parameters over a small range in an attempt to determine scaling relations and extrapolate the nature of the saturated state to parameters characteristic of real stars. Our paper is organized as follows. The numerical methods and selected parameters are described in §2. In §3, we discuss the results regarding the linear and non-linear evolution of the Tayler instability. We finally conclude in §4.

Simulation code
We use the spectral MHD code Dedalus 1 (Burns et al. 2020) for our simulations. We use version 2.2006 of the code with commit hash 9bf7eb1. Because of its spectral nature, Dedalus can achieve comparable accuracy with relatively lower resolutions compared with extremely high-resolution simulations using finite-volume codes, and it parallelizes efficiently using MPI. This feature is particularly useful for our 3D simulations, since only in three dimensions can the Tayler instability develop. Dedalus has already demonstrated its ability to handle different types of MHD problems including effects of stratification and nearly incompressible dynamics including convection, waves, and magnetic fields (e.g., Lecoanet et al. 2015Lecoanet et al. , 2017Couston et al. 2018).

Initial conditions
For convenience, we non-dimensionalize the initial conditions by setting the characteristic scales (width box , averaged gas densitȳ and gravity ) to unity 1. To mimic a latitudinal band of a star, our simulations are performed in 3D cylindrical coordinates ( , , ) with the domain size of where in = box , out = 2 box and = box /2 2 . We set up initial conditions as a magnetized, gravitationally stratified medium with density gradient 0 / and gravitational acceleration along the -axis: Here, 0 denotes the density of the unperturbed state which is a function of the scale height . The averaged density¯, and the density gradient 0 / are constants.
1 http://dedalus-project.org 2 The simulation domain has an aspect ratio of 1 in -plane. We choose this aspect ratio for convenience, and the wave numbers along and directions can also be sufficiently resolved with this aspect ratio, given that ∼ 2 / A (where the Alfvén frequency A ≡ / √¯) and ∼ 2 / box are at similar orders of magnitude for parameters used in our simulations (see the following §2.6).
We initialize the simulation with a toroidal magnetic field along -direction, with a power-law profile in radius: where 0 and are constant. Here we use = 2 for the initial magnetic field profile which is expected to be Tayler unstable (Tayler 1973). The system is initially in magnetostatic equilibrium with the unperturbed pressure 0 satisfying: such that pressure gradients are balanced by magnetic forces and gravity respectively in the and directions.

Initial perturbations
Since the Tayler instability is a non-axisymmetric instability, initially non-axisymmetric perturbations are needed, otherwise, perfect axisymmetry will be maintained throughout the simulations. To maintain ∇ · = 0, we effectively evolve the magnetic vector potential in our equation sets with ≡ ∇ × . We initialize white-noise perturbations to the magnetic fields by setting the magnetic potential vector as: where [0, 1] is a random number uniformly distributed between 0 and 1. By taking = ∇ × , we obtain divergence-free magnetic fields with in desired form in Eq. (8) and white noise perturbations with a magnitude of ∼ 10 −10 0 in and . Because the white noise distribution does not introduce any characteristic length scales, these initial conditions do not add to the initial magnitude of any particular modes.

Governing equations
We express fluid quantities as the sum of the unperturbed fields (denoted by the subscript 0) and the variations (denoted by the prime symbol), e.g., = 0 + , = 0 + , etc., and solve the following fundamental governing equations of incompressible magnetohydrodynamics: where / denotes / + · ∇, is the fluid velocity, and is the Brunt-Väisälä frequency defined as: We use a Boussinesq approximation that | −¯| |¯|, which is acceptable because of the incompressible nature of the Tayler instability and its short radial length scale. Vertical stratification appears through the buoyancy term, which appears in spite of the Boussinesq approximation. This mimics a simulation of a star over a radial length scale much less than the density scale height. We transform our simulations into the rotating frame by adding the Coriolis term 2 0 × , with bulk angular velocity 0 = Ω 0ˆ. We include explicit diffusivity, viscosity and magnetic resistivity as , and respectively, where the diffusivity mimics the compositional diffusivity in a real star. Temperature perturbations are not included because the instability operates in an isothermal regime in post-MS stars. Since the magnetic diffusivity is usually larger than microscopic viscosity in real stars, we adopt a relatively large magnetic diffusivity with > ∼ .

Boundary conditions
We apply periodic boundary conditions along the and -directions. Note that although a density profile in the direction is implied as described by Eq. (6), what actually solved in the governing equations (13) -(17) are variations of fluid quantities (e.g., , , , , etc.), therefore periodic boundary conditions can be applied to the -axis. We apply = 0, = 0 and = 0 at the inner and outer boundaries. We apply the electric scalar potential = 0 and the magnetic potential = 0 and on both inner and outer boundaries, and out − in at the inner and outer boundaries respectively, in order to maintain continuity in the magnetic potential described by Eq. (12). These boundary conditions are consistent with the initial field profile but allow the magnetic field to evolve. The boundary conditions enforce the toroidal magnetic flux to be conserved, but the magnetic energy can decrease (although it cannot go to zero).

Simulation parameters
The non-dimensionalized parameters used in our simulations are summarized in Tab. 1. Since the simulations span a range of parameters (mainly the rotation frequency Ω 0 , initial Alfvén frequency A , and Brunt-Väisälä frequency ), a combination of these name elements in Tab. 1 is used to refer to one simulation where a certain combination of parameters are adopted, e.g., the notation "Om.5_OmA.25_N1" refers to the simulation with Ω 0 = 0.5, A = 0.25 and = 1.
We note that like most other numerical work, our simulations depart from the actual parameters due to limited computational power. The Reynolds number Re ≡ 1/ and the magnetic Reynolds number Re m ≡ 1/ used here are much smaller than those in real stars. In addition, real RGB stars likely have A /Ω ∼ 10 −1 1, Froude number Ω/ ∼ 10 −4 1 and magnetic Prandtl number Pr m ≡ / 1; however, our simulations can only reach much smaller scale separations with A /Ω ∼ Ω/ ∼ 0.5 and Pr m ∼ 0.5.
The scale separation in our simulations is constrained due to the following reasons: (1) since the vertical ( direction in code setup) length scales of the Tayler instability must satisfy (Ω/ ), Ω/ thus cannot be too small otherwise cannot be resolved; (2) since the growth rate of the Tayler instability roughly scales as 2 A /Ω, the initial A (or 0 ) cannot be too small, otherwise the growth of the Tayler instability will be too slow for the simulations to follow; 3) since the Tayler instability occurs when (Spruit 2002;Zahn et al. 2007), given A /Ω < 1 and Ω/ < 1, the magnetic resistivity needs to be small enough to allow the Taylor instability to develop, but not too small to be unresolved on the grid scale. Similarly, although Pr m ≡ / 1 would be ideal, here we use Pr m = 0.5 so that will not be too small to be resolved either.
We shall see that since the scale separation in our simulations ( A Ω and Pr m 1) is much smaller than that in real stars ( A Ω and Pr m 1), we do not expect the scaling relations measured from the simulations to perfectly replicate theoretical predictions made under the limit of large-scale separations (e.g., Fuller et al. 2019). Nevertheless, our simulations probe more realistic parameter spaces with the correct ordering of scales A < Ω < that occur in real stars. This parameter space has not been fully explored yet in previous studies, such as Braithwaite (2006) who used Ω/ = 1 and compressible fluid equations, Weber et al. (2015) and Gellert et al. (2008), who used = 0, and Guerrero et al. (2019) who used Ω = 0. As will be seen in the following sections, this setup enables a clean and detailed numerical study of the Tayler instability and its saturation while retaining much of the key physics for stellar interiors.

Morphologies
We first study the run Om.5_OmA.25_N1 as a fiducial case. 10 -11 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 u rms u r u z u φ 0 20 40 60 80 100 t 10 -11 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 (right two columns of Fig. 1), the amplitudes have saturated, and the = 1 structure has mixed together with higher modes. Fig. 2 shows the time evolution of the same set of quantities as Fig. 1, but viewed in the − plane. The structures on theplane are mostly horizontally aligned, with a short wavelength in the −direction. Because the dominant perturbations have = 1, the values of , , etc., oscillate in time as the flow pattern propagates. In the non-linear and saturation phases, the banded structure is mostly maintained, but with a slightly longer wavelength in the -direction and distorted structure. However, the -components of the flow, and , become highly turbulent during the non-linear saturation phase, losing the banded structure. This appears to be the result of secondary shear instabilities that develop as the instability saturates. The unstable eigenmodes are initially confined to 1.5. Simultaneously with saturation, the flow pattern migrates inward, reaching nearly the inner boundary. We will discuss the transition to the nonlinear stage and the shear instability in the following sections. Fig. 3 shows the time evolution of the root-mean-squared velocities (top) and magnetic fields (bottom) in the fiducial simulation Om.5_OmA.25_N1. The evolution of both velocities and magnetic fields goes through a linear growth stage for ∼ 80 before reaching saturation. As expected, the magnitude of is smaller than and by a factor of several, due to the buoyancy force that restricts motion in the −direction. Similarly, remains several times smaller than or . At saturation, the perturbed field components remain more than a factor of ten weaker than the background field , which weakens only slightly by the end of the simulation.

Linear growth rate
We further examine the growth of magnetic fields by decompos- . Time evolution of volume-averaged root-mean-squared azimuthal modes of magnetic fields | | in the fiducial simulation Om.5_OmA.25_N1. The = 1 mode is dominant in the linear growth phase, and higher modes start to grow at later times, e.g., the = 0 and = 2 mode grows from ∼ 40 at twice the rate of the = 1 mode.
ing them into different azimuthal modes with the following equation with as azimuthal mode numbers, and ... , denoting averaging over and under cylindrical coordinates. We plot the time evolution of the amplitudes of different modes in Fig. 4. The = 1 mode is dominant over higher modes in the linear growth phase, consistent with the apparent = 1 mode structure in Fig. 1. Modes of ≥ 2 start to grow at ∼ 40, and all modes ultimately reach saturation at ∼ 80. We will discuss the nonlinear coupling and saturation in §3.3 and §3.4. We measure the linear growth rates of the Tayler instability from simulations with varying initial-conditions parameters, including the angular frequency Ω 0 , Aflvén frequency A and Brunt-Väisälä frequency . The growth rates measured from simulations (orange stars) are plotted against the eigenmodes (blue dots) of the Tayler instability setup calculated with the eigenvalue problem solver in Dedalus. The simulated linear growth rates are well-predicted by the linear eigenvalue calculations, following a scaling relation of ∝ 4.5 A Ω −2 0 0 (Fig. 5). This scaling is different than expected from Spruit (1999)  ∼ 0.04 -0.05, Ω 0 ∼ 0.1 -0.2, = 1 and = 5 × 10 −6 ) and consequently much smaller growth rates ( ∼ a few ×10 −3 ) than the fiducial simulations. With larger scale separations between A , Ω 0 and , the best power-law fits (dashed line) is roughly consistent with the theoretical prediction of ∝ 2 A Ω −1 0 (Spruit 1999), and the growth rates are much smaller than those in Fig. 5. in the limit that A Ω 0 . This occurs because our simulations are not actually in the asymptotic limit of A Ω used for the analytic estimates in those works: in Fig. 6, we further carry out linear eigenvalue calculations with larger scale separations with ∼ 0.04 -0.05, Ω 0 ∼ 0.1 -0.2, = 1 and = 5 × 10 −6 , and find that the obtained scaling relations of the growth rates are quite consistent with the predictions by Spruit (1999). However, the resulting growth rates are as low as a few 10 −3 , which are much smaller than those with fiducial parameters ( ∼ a few ×10 −1 in Fig. 5) and are prohibitively small for numerical simulations to follow the growth of the Tayler instability. Therefore, we stick to the fiducial parameters for the simulations even though they have limited scale separations, and bear it in mind when comparing our results with analytic estimates that the simulations are not fully in the limit of A Ω 0 as used in many analytic works.

Nonlinear coupling
From Fig. 4, we can see that during the linear growth phase, the amplitude of the = 1 mode grows exponentially as expected. The = 0 and ≥ 2 modes initially decay because they are stable, but eventually they also start growing exponentially. This is a consequence of non-linear power transfer from the large amplitude = 1 mode to other values of . From the induction equation (16), we can see that the magnetic field grows as (neglecting diffusive effects which are small in this analysis) During the linear growth stage, fluctuations in and are dominated by the = 1 component, which has time and spatial dependence ∝ sin( − ) , and similar for . Here is the real part of the frequency of the fastest growing mode, and is its growth rate. Hence, to lowest order, the ≠ 1 components grow as for = 0 and ≥ 2 modes. Multiplying each side of Eq. (21) by sin( ) and integrating over volume gives the contribution to the non-linear growth of ì at a desired wavenumber . We see that to the lowest order, only the = 0 and = 2 modes have a non-vanishing integral, arising from the first and second terms inside the parentheses in Eq. (21), respectively. Hence, we see that the = 0 and = 2 modes grow at exactly twice the rate as the = 1 mode, as long as the = 1 mode has much larger amplitude. At a given time , the amplitude of the non-linearly excited modes scales as ì ∝ ∫ ì ∝ ( , , = 1, = 0) 2 2 . Hence at a given time , the ratio of the non-linearly excited mode to the linearly excited mode is where is a constant of proportionality that captures the non-linear coupling between modes. From Fig. 4, we estimate ∼ 0.1/ , where is the background magnetic field strength. Hence, the non-linearly excited modes have much smaller amplitudes in the linear regime where ( = 1) . Extending this calculation to > 2 requires higher chains of non-linear interaction that results in a non-linear growth rate for modes with ≥ 2. This scaling is verified by the growth rates and amplitudes shown in Fig. 4 in the linear regime, with 70. Hence, in the linear growth phase, we clearly see a non-linear transfer of power to smaller scales. By the time the instability saturates, the high-modes have reached amplitudes comparable to (but smaller than) the = 1 mode, and the weakly non-linear analysis presented above begins to break down.

Non-linear saturation
The end of linear growth and apparent saturation of the instability at around = 80 in our simulations is accompanied by a remarkable change in the motions and morphology of the simulation domain (see Fig. 2). This includes the appearance of turbulent non-layered structure in and (Fig. 2), as well as the inward motion of the perturbed flow.
The turbulent structure appears to stem from secondary shear instabilities. We believe the evolution is similar to the magnetized Rayleigh-Taylor instability, where the primary instability also drives opposing flows which then exhibit shear instabilities (e.g., Stone & Gardiner 2007). For the Tayler instability as we consider here, buoyancy provides a restoring force in the direction. The horizontal flows driven by the instability will become unstable to secondary Kelvin-Helmholtz instabilities when where ⊥ = √︃ 2 + 2 is the fluid velocity perpendicular toˆ.
Near the end of the simulation, the dominant modes have roughly seven wavelengths in the −direction and hence ∼ 40. Therefore in our simulations with = 1, we expect shear instabilities to occur when ⊥ ∼ 0.05. Indeed, Fig. 3 shows that the system reaches its saturated state very near this scale, and that this saturation is accompanied by turbulence shown in Fig. 2. Despite the shear instabilities, , , , and all maintain a banded structure in the -direction, and they also maintain a predominantly = 1 azimuthal structure.
At the same time that secondary instabilities develop, the unstable region also "migrates inward". During linear growth, the instability is restricted to large radii with 1.5, because the Tayler instability only occurs in the outer part of the domain where the background magnetic field is stronger. As the instability saturates around = 80, however, the unstable motions and magnetic field perturbations move inwards to the inner boundary at = 1. The inward-moving fingers have a predominantly = 1 azimuthal structure and feature inward radial motion at alternating heights, accompanied by an outward return flow in between. Their structures resemble convective plumes/fingers seen in the early phases of convective/thermohaline instability.
We posit that the secondary instabilities cause saturation both by dissipating energy from the growing modes, and by allowing for the inward motion of the instability. Once the fluid elements obtain velocities large enough to overcome the Richardson criterion of Eq. (23), inertial forces become comparable to buoyancy forces, and fluid elements can flow with fewer constrictions. This may allow for circulation to smaller radii that was previously prevented by buoyancy/Coriolis forces. It is also possible that non-linear inertial terms effectively change the linear dispersion relation, allowing unstable motions to grow at smaller radii. Future work should investigate this effect in detail, and determine whether latitudinal migration of the Tayler instability could occur in real stars.

Dynamo action
A major question regarding the Tayler-Spruit dynamo is whether a dynamo actually occurs. In the picture advanced by Spruit (2002) to components from the simulation Om.5_OmA.25_N1, demonstrating the non-linear induction process needed to amplify the axisymmetric poloidal field to close a dynamo loop. differential rotation amplifies by winding up a weak poloidal field. The dynamo loop is supposed to be closed by the amplification of the poloidal field via the unstable motions associated with Tayler instability. However, as pointed out by Zahn et al. (2007) (see also Fuller et al. 2019), the unstable motions are = 1 and do not directly produce an axisymmetric component of the poloidal field that is needed for amplification of via differential rotation. Hence, a non-linear coupling process is needed to amplify the axisymmetric poloidal field, and this step of the dynamo is poorly understood. Fig. 7 shows the growth of the axisymmetric component of , , =0 , which is the poloidal component of the field needed to be generated by the Tayler instability in order for a dynamo to occur. Our simulations do not include differential rotation and thus cannot capture the winding of , =0 needed to close the dynamo loop. However, our simulations clearly demonstrate nonlinear amplification of , =0 , as also discussed in Section 3.3. We have verified that the proportionality of the scaling / ∼ A / predicted by Spruit (2002) and Fuller et al. (2019) roughly holds, but with a much smaller normalization factor of ∼ 10 −2 for the saturated state of our simulations, rather than the predicted ∼ ( A / ) ∼ 1/4 , i.e., the axisymmetric component of in our simulations saturates at much lower values than predicted by the scaling relation. However, both of those models are based on dynamos sustained by energy input by shear, which cannot occur in our simulations, so the disagreement is not surprising.
Future work will be necessary to understand the saturated state of the Tayler-Spruit dynamo, but our simulations do indicate that a dynamo based on non-linear induction can greatly amplify the axisymmetric component of the poloidal field.

Scaling Relations
A key result of our simulations is how the properties of the saturated state depend on input parameters ( A , Ω 0 , etc.). Fig. 8 shows how the mean amplitude of the magnetic field perturbations in the saturated state, ,sat , scale with Ω 0 , A , and . The saturated values of ,sat scale strongly with rotation rate and initial magnetic field, with ,sat ∝ Ω −1.5 0 and ,sat ∝ 3 A , and almost no dependence on . These scalings are similar to (but slightly weaker than) the linear growth rate scaling presented in Fig. 5.
Based on the discussion above, we believe the instability saturates in our simulations due to secondary shear instabilities. For the large-scale growing modes with ∼ 1, the non-linear inertial term ( · ∇) produces an effective damping rate of order NL ∼ 2 / . In the rapidly rotating limit, we expect ∝ ( A /Ω 0 ) ,sat (Spruit 2002;Fuller et al. 2019), so the non-linear damping rate 3 × 10 -1 4 × 10 -1 6 × 10 -1 will scale as NL ∝ ,sat ( A /Ω). Setting this equal to the linear growth rate, which scales as ∝ 4.5 A /Ω 2 0 (Fig. 5) yields an expected scaling of the non-linear saturation amplitude of ,sat ∝ 3.5 A /Ω 0 , similar to the scaling Fig. 8. The slight mismatch between these scaling laws could stem from the small dynamic range covered by our simulations, or the fact that they are not quite in the rapidly rotating limit. Fuller et al. (2019) argued that Tayler instability saturates via weak magnetic turbulence with a non-linear damping rate of NL ∼ A ∝ ,sat . Setting this expectation equal to the linear growth rates would entail ,sat ∝ 4.5 A /Ω 2 0 , steeper than the scalings shown in Fig. 8. Hence, it appears that our simulations do not saturate via weak magnetic turbulence as predicted by Fuller et al. (2019). However, in the next section, we discuss why the saturation mechanism is likely to be different when the Tayler instability operates in stellar interiors.

DISCUSSION AND CONCLUSION
Although secondary shear instabilities appear to saturate the Tayler instability in our simulations, it is likely that a different mechanism will operate in real stars. Fuller et al. (2019) predict that non-linear Alfvén wave dissipation produces damping at the rate NL ∼ A / , which in our units equates to NL ∼ / ∼ 0.04 in the saturated state. The dimensionless growth rate (Fig. 5) is larger than this, roughly 0.4 for the fiducial simulation Om.5_OmA.25_N1. This indicates that the instability saturates via shear instabilities before it reaches an amplitude high enough for Alfvén wave dissipation to dominate. In a real star, however, Alfvén wave dissipation would likely occur before shear instability. As described above, the dissipation rate from shear instabilities is NL ∼ 2 / , and we expect ∼ ( A /Ω 0 ) A in the limit A /Ω 0 1 applicable to real stars. Because our simulations only reach A /Ω 0 ∼ 1/2, dissipation via shear instability and Alfvénic dissipation are comparable, and evidently, shear instability is more important by a factor of a few. But in a real star with A /Ω 0 1, shear instabilities become less important, and Alfvénic dissipation is expected to dominate. Hence, our simulations cannot grow to an amplitude high enough to test the prediction of Fuller et al. (2019), because of the small ratio of Ω to A . Future work should push towards larger scale separations to better test these predictions.
To determine whether Alfvén wave damping or shear instability will saturate Tayler instability, we can determine which mechanism requires a lower threshold to operate. As described above, shear instability requires ⊥ ∼ 2 / . The fastest growing modes have 2 ∼ (orange) and total + (green) energies from the fiducial simulation Om.5_OmA.25_N1. Together with the order-of-magnitude growth of the kinetic energy driven by the Tayler instability, the magnetic energy slightly dissipates at a rate roughly consistent with the prediction of ∼ | |, where is the growth rate of the Tayler instability.
Another argument against shear instability being important in real stars is as follows. According to the saturated state of Fuller et al. (2019), the value of ⊥ reaches ⊥ ∼ Ω 2 / , where = ln Ω/ ln is the dimensionless shear. This is much smaller than unless Ω 2 2 , which cannot occur in real stars, since a background state with Ω 2 2 would already satisfy the Richardson criterion for overturning the stratification. Therefore, at the saturated state due to Alfvén wave dissipation, the value of ⊥ would be much less than required to overturn the stratification and cause shear instability. Simulations including shear and pushing to smaller ratios of Ω 0 / will be needed to test this prediction.
Our simulations appear to disagree with the model of Spruit (2002), in which turbulent dissipation in the saturated state dissipates magnetic energy at the rate ∼ 2 , where is the linear growth rate (equal to the non-linear damping rate) of the Tayler instability. Even though our saturated state is somewhat turbulent, it does not dissipate magnetic energy at this rate, as can be seen from the bottom panel of Fig. 9. Since the growth rate of the instability is ∼ 0.4, this would imply that all of the magnetic energy would dissipate within a time span of Δ ∼ 2.5. In contrast, Fig. 9 shows that the magnetic energy only decreases by ∼10% in the final Δ = 20 of our simulation, corresponding to a magnetic energy dissipation rate of ∼ 2 × 10 −2 as shown in Fig. 10. This corresponds to magnetic energy dissipation ∼2 orders of magnitude slower than ∼ 2 . Fig. 10 also shows that kinetic energy dissipation is negligible relative to magnetic energy dissipation. However, we note that our imposed background field is ordered (uniform in the −direction), and a field built up by a dynamo could be disordered, allowing for more magnetic energy dissipation. Also, our boundary conditions prevent from being totally erased, which may limit the turbulent dissipation of . Future work with different initial conditions and boundary conditions will shed light on this issue.
Shortly before this paper was submitted, Petitdemange et al. (2023) presented a suite of dynamo simulations driven by differential rotation and exhibiting the Tayler instability in a spherical shell. Since their methods were very different from ours, we do not attempt a direct comparison with our results. In their setup, shear was created by enforcing the outer boundaries to rotate at different rates, though this also created Ekman boundary layer effects and a related "weak dynamo" that complicated the analysis. The Tayler instability was triggered once the magnetic field generated by the weak dynamo exceeded the threshold derived from linear theory. The saturation mechanism of the Tayler instability was not clear from that work, but the magnetic torques in the saturated state appeared to scale according to the predictions of (Spruit 2002). However, most of those simulations appeared to have Ω/ ∼ 1 3 , where the predictions of (Spruit 2002) and Fuller et al. (2019) are similar. While the work of Petitdemange et al. (2023) is a great leap forward in modeling the Tayler instability, our understanding of the saturation mechanism and magnetic torques in real stellar interiors remains incomplete.

Conclusion
In this study, we perform three-dimensional MHD simulations of the Tayler instability in a cylindrical annulus, with strong buoyancy and Coriolis forces incorporated to mimic realistic stellar environments. The simulations are initialized with a strong toroidal field which is in magnetostatic equilibrium. We explore a range of parameter space by varying the angular velocities Ω 0 , Alfvén frequency A and Brunt-Brunt-Väisälä frequency , with the the correct order of > Ω 0 > A that typically exists in realistic stars. We find that as theoretically expected, the initial conditions adopted are unstable to the Taylor instability. The = 1 mode clearly dominates the linear growth of the instability, and the linear growth rates are well-predicted by the linear eigenvalue calculations. The linear growth phase is later accompanied by a non-linear coupling between the = 1 and other azimuthal modes, leading to the growth of = 0 and ≥ 2 modes. The = 0 component of the poloidal field is amplified by the instability, signaling a dynamo that can regenerate the poloidal field as necessary for the Tayler-Spruit dynamo to occur.
Both the linear growth rates and the saturated magnetic field measured in the simulations scale strongly with angular velocities Ω 0 and the Alfvén frequencies A parameters, and are nearly independent of . While this has been predicted from linear theory and non-linear saturation models, the scaling in our simulations is steeper, due to the fact they are not in the asymptotic limit of A Ω 0 assumed in analytic work. With greater scale separations, the linear eigenvalue calculations are able to replicate the scaling relations of ∝ 2 A Ω −1 0 predicted by analytic work, even though the resulted growth rates are prohibitively too small to simulate numerically.
The linear growth is ultimately followed by a non-linear saturation of the instability, which appears to be caused by secondary shear instabilities. The saturation is also accompanied by an inward migration of unstable motions, whose cause is not clear. We argued that saturation via secondary shear instability is unlikely to operate in real stars where the stratification is greater (preventing shear instabilities), and where Alfvén wave damping becomes more important due to the larger separation of scales than we could achieve in our simulations. During the saturated phase, energy is dissipated primarily through magnetic diffusion. However, it is dissipated much more slowly than predicted by the model of Spruit (2002), likely entailing the instability can grow to larger amplitudes as expected in the model of Fuller et al. (2019).
Certain caveats apply to this study. First, due to limited com-putational capability, the scale separations between the parameters Ω 0 , A and are much less than that in real stars, the scaling relations obtained from simulations might differ from what occurs in realistic stellar environments. Second, although amplification of the axisymmetric poloidal magnetic field is observed in the simulations, differential rotation is not included in our simulations, which is necessary to close the loop of the Tayler-Spruit dynamo. We thus cannot directly test theoretical predictions of the angular momentum transport caused by this dynamo, leaving room for improvement in future work.