Modeling of Magneto-Rotational Stellar Evolution I. Method and first applications

While magnetic fields have long been considered to be important for the evolution of magnetic non-degenerate stars and compact stars, it has become clear in recent years that actually all of the stars are deeply affected. This is particularly true regarding their internal angular momentum distribution, but magnetic fields may also influence internal mixing processes and even the fate of the star. We propose a new framework for stellar evolution simulations, in which the interplay between magnetic field, rotation, mass loss, and changes in the stellar density and temperature distributions are treated self-consistently. For average large-scale stellar magnetic fields which are symmetric to the axis of rotation of the star, we derive 1D evolution equations for the toroidal and poloidal components from the mean-field MHD equation by applying Alfven's theorem, and a conservative form of the angular momentum transfer due to the Lorentz force is formulated. We implement our formalism into a numerical stellar evolution code and simulate the magneto-rotational evolution of 1.5 M$_\odot$ stars. The Lorentz force aided by the $\Omega$ effect imposes torsional Alfven waves propagating through the magnetized medium, leading to near-rigid rotation within the Alfven timescale. Our models with different initial spins and B-fields can reproduce the main observed properties of Ap/Bp stars. Calculations continued to the red-giant regime show a pronounced core-envelope coupling, which reproduces the core and surface rotation periods determined by asteroseismic observations.


Introduction
Magnetic fields are visible for many different types of stars.As the most evident example, the Sun exhibits magnetic activity in the form of spots, prominences, flares, and mass ejections (e.g., Solanki et al. 2006).The Sun is considered to be prototypical of FGK-type main-sequence stars that are known to have convective envelopes and most of these cool stars are known to host magnetic fields (Landstreet 1992;Donati & Landstreet 2009).These magnetic fields are thought to have a "dynamo origin", meaning that the field is continuously amplified through hydrodynamic induction and would otherwise decay within the Alfvén timescale.This understanding is supported by observed correlations of the field strengths detected in these main-sequence stars with their fundamental parameters, such as their mass, age, and rotation periods (Vidotto et al. 2014;See et al. 2015See et al. , 2016;;Folsom et al. 2016Folsom et al. , 2018)).
In the upper main sequence, where stars have radiative envelopes, about 10% of the stars are magnetic (Landstreet 1992;Wade et al. 2014).Their strong (typically ∼1 kG) fields are characterized by a large scale (∼dipole) structure and since neither convection nor rotation could currently produce these fields, it is believed that they are stable over a significant part of the stellar lifetime (Wade et al. 2000;Silvester et al. 2014).In contrast to stars with convective envelopes, clear correla-Movie associated to Fig. 3 is available at https://www.aanda.orgtions between the field strength and major stellar parameters have not been found thus far.Such properties indicate the "fossil" origin of the field, namely: the strong magnetic field was somehow amplified in the past and ultimately reached a stable configuration through magneto-hydrodynamical relaxation (Braithwaite & Spruit 2017).
Magnetic fields are also found among evolved stars; GKgiants (Aurière et al. 2015) as well as asymptotic-giant-branch (AGB) and post-AGB stars (Vlemmings 2014(Vlemmings , 2019, and references therein), AFGK-yellow supergiants (Grunhut et al. 2010), and M-supergiants including α Ori (Aurière et al. 2010;Tessore et al. 2017) -despite the fact that field strengths in evolved stars are often small because of the large radii and slow rotation rates involved (Aurière et al. 2008).Furthermore, ∼10% of the white dwarfs show magnetic fields and while all neutron stars appear to have magnetic fields, some 10% of them (known as magnetars) possess extremely strong surface magnetic fields (Chanmugam 1992;Ferrario et al. 2015).
Magnetic fields can have significant effects on stellar evolution.It has been known for a long time that the spin evolution of solar-like stars is essentially coupled with the magnetic field evolution.The slow rotation rate of the Sun has been understood as the result of the magnetic braking (Weber & Davis 1967) and the solar wind itself is largely driven by the surface magnetic activities (e.g., Parker 1958;Ofman 2010).The magnetic stress of the amplified field inside the convective envelope contributes to the angular momentum transfer together with the turbulent viscosity and the mean meridional flow to determine the rotation profile in the envelope (Brandenburg 2018).Moreover, the convective dynamo is the result of the symmetry breaking due to stellar rotation (Brun & Browning 2017).
The evolution of main-sequence stars with radiative envelopes can also be significantly affected by magnetic fields.Because of the intrinsically strong stellar wind, magnetic braking for massive OB-type stars can be strong enough to become directly observable (e.g., σ Ori E; Townsend et al. 2010).A strong surface field can trap the wind material into a corotating magnetosphere, which provides an extensive explanation for the X-ray emission of, for instance, IQ Aur (Babel & Montmerle 1997) and the phase variations of Balmer-line emissions of σ Ori E (Townsend & Owocki 2005).
Theoretical works have also revealed the importance of the magnetic field in stellar evolution.For instance, magnetic stress furnished by internal fields may account for the efficient angular momentum transfer in radiatively stratified regions (Spruit 1999).Many stellar evolution simulations take magnetic angular momentum transfer into account, following a model originally proposed by Spruit (2002) and reassessed by other authors (Maeder & Meynet 2003;Heger et al. 2005;Suijs et al. 2008;Denissenkov & Pinsonneault 2007;Fuller et al. 2019;Ma & Fuller 2019).The effect of magnetic braking, as well as magnetic wind confinement, have been evaluated in the context of massive star evolution (Meynet et al. 2011;Petit et al. 2017;Georgy et al. 2017;Keszthelyi et al. 2019).Sufficiently strong internal fields can modify the stellar structure via magnetic pressure and tension and, more importantly, by affecting the adiabatic indices, which, in turn, modifies the efficiency of convective energy transport.Such effects have been considered in Feiden & Chaboyer (2012, 2013, 2014) for low-mass star evolution, following prescriptions developed by Lydon & Sofia (1995).
The majority of theoretical works have modeled magnetic stars by considering the specific magnetic effects individually.However, all effects, in reality, must relate to each other because they are governed fundamentally by precisely the same field.To consider the integrated effect, the global structure of the magnetic field must be modeled.In addition, since the magnetic field ought to evolve over time along with other physical quantities, the time-dependent treatment is most desirable.Such a treatment was recently derived by Potter et al. (2012).In their framework, the global field structure is significantly simplified to have axial symmetry and dipole-like structure.Furthermore, they formulated the evolution equation of the global magnetic field based on a mean-field MHD equation.As a consequence, important processes are effectively incorporated, such as the αand η-effects, which express the interaction between turbulence and magnetic field, as well as the Ω effect, in which differential rotation winds up poloidal magnetic component to enhance the toroidal component.
Meanwhile, the prescription in Potter et al. (2012) still displays some incompleteness.For example, it is likely that their evolution equations for magnetic fields do not reproduce the magnetic flux conservation, which is one of the most fundamental outcomes of the ideal magnetohydrodynamic (MHD) assumptions.Similarly, their expression of the angular momentum transport by the Lorentz force does not reproduce the angular momentum conservation.These problems arise because of the ambiguity that exists in the averaging process for formulating 1D evolution equations starting from the more general 3D MHD equation.
In this work, we present a new framework in which a mutual interaction of the magnetic field and the rotation during the evolution of the star is treated in a physically consistent manner.In the next section, equations that describe the evolution of the stellar magnetic field and stellar rotation, assumptions and approximations made for the formulation, and a brief description of the numerical construction are provided.In Sect.3, we show that our formulation, including the Ω effect and the Lorentz force, naturally leads to a torsional Alfvén wave, in which differential rotation and toroidal magnetic field propagate together inside the star.We further discuss that the Alfvén wave accounts for a highly efficient mechanism of angular momentum redistribution when realistic dissipation is taken into account.While we plan to apply this new formulation to general-purpose stellar evolution simulations in the future, here we provide results of magneto-rotational evolution calculations for stars of 1.5 M to demonstrate the capabilities and limitations of our formulation.The corresponding main-sequence evolution is analyzed in Sect.4, and the red-giant phase is presented in Sect. 5.In Sect.5, we discuss comparisons between our simulation results and other theoretical models (Sect.6), and relevant observations (Sect.7).We present our conclusions in Sect.8.

Methods
To compute the time evolution of stellar models, we use the 1D stellar evolution code HOSHI (Takahashi et al. 2016(Takahashi et al. , 2018)).The code iteratively solves the four structure equations of mass conservation, the momentum balance equation in hydrostatic or hydrodynamic form, the equation of energy conservation in the form of an entropy equation, and the energy transport equation by the so-called Henyey method.The equation of state in the code consists of a mixture of ideal gases of photon, averaged nuclei, electron, and positron.An analytical treatment of Blinnikov et al. (1996) is applied for the electron-positron gas.The free energy of the Coulomb interaction for degenerate states is included (Salpeter & van Horn 1969;Slattery et al. 1982).The ionization of hydrogen, helium, carbon, nitrogen, and oxygen is also treated by solving the Saha equation.For the opacity, the Rosseland mean opacity of the OPAL project (Iglesias & Rogers 1996) is used together with the conductive opacity by Potekhin et al. (2006) and the molecular opacity by Ferguson et al. (2005).
In addition to the structure equations, the evolution of the abundances of the chemical species is solved through a reactiondiffusion equation as where Y i , Ẏi,reac , and D eff are the mole fraction of ith isotope, the rate of change of Y i due to nuclear reactions, and the effective chemical diffusivity, respectively.49 isotopes 1 are considered in this work, and the reaction rates are taken from the current version of JINA REACLIB (Cyburt et al. 2010), except for the 12 C(α,γ) 16 O-rate, for which we use the rate from Caughlan & Fowler (1988) multiplied by a factor of 1.2.The Ledoux criterion is used to evaluate convective instability.The standard mixing-length-theory (Böhm-Vitense 1958) is applied to compute the amount of energy transported by convection.
The mixing-length theory also provides the diffusion coefficient for chemical mixing in convection zones as D cv = 1 3 v cv l cv , where v cv is the velocity of convective eddies, l cv = α MLT min(H P , r) is the length scale of the convective flow, and α MLT and H P are the mixing-length parameter and the pressure scale height, respectively.To consider the effect of convective overshooting, the eddy velocity for regions surrounding the convective region is calculated as where f ov is an adjustable parameter determining the e-folding length scale, v cv,0 and h P,0 are the eddy velocity and the pressure scale height at the edge of the convective region, and ∆r is the distance from the edge.This treatment yields a similar diffusion coefficient distribution to the exponential diffusive overshoot described in Herwig (2000).
In semiconvective layers, we apply a diffusion coefficient of the form and thermohaline convection is treated with (Kippenhahn et al. 1980;Wellstein et al. 2001;Siess 2009).Here are the thermal diffusivity, the radiative temperature gradient, the adiabatic temperature gradient, the µ-gradient, and relevant thermodynamic derivatives, respectively, with κ being the Rosseland mean opacity.For the two control parameters, we use f sc = 0.3 and f thh = 1.0.
The HOSHI-code includes stellar wind induced mass and angular momentum loss.For models presented in this work with low effective temperatures of log T eff [K] < 3.9, an empirical mass-loss formula by de Jager et al. (1988) is applied.Consequently, 1.5 M nonrotating model experiences wind mass loss of ∼3 × 10 −12 M yr −1 and ∼10 −11 −10 −8 M yr −1 during the main-sequence and red-giant phases, respectively.

Stellar rotation
The effects of stellar rotation are similarly taken into account as described in Takahashi et al. (2014).To describe a rotating star in a 1D formulation we assume shellular rotation, where all material on an isobaric surface shares the same angular velocity (Zahn 1992).We define the volume and mass enclosed by an isobar as either V P and M P .Accordingly, the mean radius is defined as 4πr 3 P /3 = V P .Thus the angular velocity Ω is defined as a function of the mass coordinate, M P , in our simulation.The derivation of the evolution equation of Ω is shown in Sect.2.3, following a description of our treatment of the stellar magnetic field.
To consider the effect of deformation due to the centrifugal force, the isobar is assumed to have a shape described by where P 2 is the second-degree Legendre polynomial.Here, ≡ (Ω 2 r 3 P /2GM P )(a/r P ) 3 indicates the degree of rotation compared with the local gravity, and the length scale, a, satisfies the relation (Denissenkov & VandenBerg 2003), The centrifugal force not only deforms the isobars but also affects the pressure balance and temperature gradient in the star.Following Endal & Sofia (1976), these effects are taken into account in the structure equations by introducing the parameters of f P and f T as (5) where ∇ MLT is the convective temperature gradient determined by the mixing-length theory.The parameters f P and f T are calculated as where g is the local effective gravity, in which the centrifugal force is taken into account.The pointy brackets imply an average over the isobaric surface S P as q ≡ 1 S P qdσ.The wind mass-loss rate is enhanced when the stellar rotation is considered and the surface condition approaches the ΩΓ limit (Langer 1998;Maeder & Meynet 2000).We apply an enhancement according to Yoon et al. (2010Yoon et al. ( , 2012) ) as where | Ṁ(v rot = 0)| is the mass loss rate of a nonrotating counterpart having the same luminosity and effective temperature, are the rotation velocity and the critical rotation velocity at the stellar equator, τ KH is the Kelvin-Helmholtz timescale, and L Edd ≡ 4πcGM/κ is the Eddington luminosity.
Even without the help of a magnetic field, the angular momentum loss by stellar winds may be one order of magnitude more efficient than the accompanying mass loss when efficient angular momentum transport takes place at the subsurface region of the star (Langer 1998).The rate of angular momentum loss for a non-magnetic star is calculated as where j surf is the specific angular momentum at the surface of the star.In our current models, j surf is treated as a constant during the mass change, which implies the angular momentum is quickly redistributed throughout the subsurface region.

Stellar magnetism
As the most fundamental assumption, we assume that a largescale stable magnetic field is embedded inside a star.Here, "stable" means that the magnetic field is in a magnetohydrostatic quasi-equilibrium and maintains its structure, particularly the geometry, for a considerably longer timescale than the Alfvén time.From a theoretical point of view, the structure, or even the existence, of such a stable magnetic field is far from trivial, which has been shown both, analytically and numerically (e.g., Tayler 1973;Wright 1973; A19, page 3 of 28 Markey & Tayler 1973, 1974;Braithwaite & Spruit 2004, 2017;Duez & Mathis 2010;Akgün et al. 2013 and references therein).Nonetheless, our assumption may be reasonable at least for the radiative stellar envelopes close to the surface because observations have shown that the surface magnetic field of radiative main-sequence stars displays a large scale structure and longtime stability (more than decades, i.e., longer than the Alfvén timescale of ∼10 yr; e.g., Landstreet 1992).Indeed, 3D MHD simulations have shown that initially random fields can relax into a stable hydrostatic equilibrium in radiative stellar envelopes, following a magneto-hydrodynamical adjustment with an Alfvén time.The equilibrium states obtained in this way often have axial symmetry and a comparable or stronger toroidal component as compared to the poloidal component (Braithwaite & Spruit 2004;Braithwaite & Nordlund 2006;Braithwaite 2009).Justifications are more scarce for the convective layers since the fields are expected to vary on the convective timescale and to show a small-scale structure.However, even in this case, our method may still allow us to follow the time-averaged evolution of the mean magnetic field strength in the convective layers.
Although the large-scale stellar fields are assumed to be in a magneto-hydrostatic equilibrium, they can evolve over time responding to environmental changes.The most evident example is the density evolution due to the background thermal or elemental change, which leads to the field-strength evolution coming about as a consequence of the magnetic flux conservation.The next such mechanism is the Ω effect, in which differential rotation inside the star winds up the poloidal magnetic component to induce the additional toroidal component.We note that differential rotation can exist in a hydrostatic object since rotation does not require any driving force.Finally, we also treat the effects of turbulence on the magnetic fields as such a background effect.More specifically, any turbulence in our model is regarded as a perturbation that does not affect the hydrostatic pressure balance and its effects on thermal, elemental, rotational, and magnetic transports are expressed through simple effective theories.For the turbulent effects on the magnetic field, we take into account the α effect, which expresses the turbulent induction of the magnetic field, and the η effect, which shows the turbulent magnetic diffusivity, according to the mean-field MHD-dynamo equation (Brandenburg 2018).
Although deformation effects due to rotation are included in the stellar structure equations as described above, those are small unless the star rotates close to critical.Therefore, spherical symmetry is assumed when deriving the governing equations of the magnetic field.This treatment makes the derivation much easier.Generalization of the derivations to the case of a star with a deformed structure is beyond the scope of the current work.

Simplification of the magnetic field
The stellar magnetic field is assumed to be axially symmetric with respect to the rotation axis of the star.The field is divided into a poloidal and a toroidal component: where B r , B θ , and B φ are the r, θ, and φ components of the magnetic fields, respectively, and they are functions of the radius, r, and the colatitude, θ.
Because the magnetic field satisfies the solenoidal (the divergence-free) condition, we can find a vector potential A that satisfies To express the poloidal magnetic field, we utilize the toroidal component of the vector potential, A tor = A φ e φ , as so that not only B tor but also B pol naturally satisfies the solenoidal condition.Because of the axial symmetry, the poloidal components of the magnetic field can be related to A φ as To express the magnetic field evolution by a 1D method, the latitudinal dependence of the magnetic field has to be somehow determined.As for the simplest case, we approximate that the poloidal field has the same latitudinal dependence as a dipolar field, thus which yields To achieve a magneto-hydrostatic state, this poloidal field will be accompanied by a toroidal component, with a comparable or stronger field strength than the poloidal field (Braithwaite & Nordlund 2006;Braithwaite 2009); otherwise, the poloidal field is unstable (Wright 1973;Markey & Tayler 1973, 1974).We may assume that this toroidal component has the polarity of sin θ.Besides, we take a toroidal field into account that is additionally induced by the Ω effect.Such a toroidal component should have the same polarity as the original poloidal field and we take one of the simplest latitudinal dependencies as sin 2θ.Thus, we write However, we show below (Sect.2.2.2) that only the toroidal component induced by the Ω effect is capable of transferring angular momentum.Therefore, we only discuss the evolution of B(r) in the subsequent sections (Sects.3-5).

Evolution equations for the magnetic field
We start by considering a surface, S , that is embedded in a star and moves with time with velocity field, V(r, t).We let P be the time-dependent vector field, which will be later regarded as the magnetic field or the electric current field.Then the flux of P on S is defined as Φ ≡ S P • ndS = S P • dS, where n is the normal vector of S .In this situation, Alfvén's theorem tells that the total time derivative of Φ is written as if and only if ∇ • P = 0 is satisfied.For interested readers, an elementary proof is given in Appendix A.
The macroscopic evolution of the magnetic field may be described by the mean-field MHD-dynamo equation: where u is the fluid velocity, η is the magnetic diffusivity, and α and η t are the pseudo-tensors indicating the α-effect and the turbulent magnetic diffusivity (Brandenburg 2018).Because of the large magnetic Reynolds number in a stellar environment, η is almost dominated by the turbulent magnetic diffusivity.Hence, going forward, we rename the total magnetic diffusivity as η + η t → η.
Since ∇ • B = 0, Eq. ( 22) can be applied to the magnetic field.Substituting Eq. ( 23) for the ∂P/∂t term in Eq. ( 22), the evolution equation of the magnetic flux is obtained as in which Stokes' theorem is applied, and For practical purposes, the surface, S , together with the boundary, C, have to be defined.Here, we take S as a polar cap with radius, r, and opening angle θ, and accordingly, the boundary, C, is taken to be a conic section that is located at the angle θ (Fig. 1).Considering the axial symmetry of the magnetic field, magnetic flux penetrating the surface is written as Similarly, the right hand side of Eq. (24) becomes We take the Lagrangian expansion velocity, v m , as the velocity of the polar cap, thus V = v m e r .The vector field, U ≡ u − V, now indicates a flow other than that which is due to stellar expansion or contraction.For the sake of simplicity, we consider that only the rotational flow contributes to U. In other words, we tentatively neglect the advection of the magnetic field by the meridional flow.Ultimately, the evolution equation of A φ is obtained as Here and hereafter, a time derivative at constant radius is shown by ∂/∂t, while a time derivative at constant mass coordinate is shown by d/dt.As we assume A φ = A(r) sin θ, we obtain the evolution equation of A(r) as Equation ( 22) can also be applied to the electric current field j(r, t) since j is proportional to the curl of B and thus divergencefree assuming MHD.Using the evolution equation of the electric current, the evolution equation of the electric current flux, is obtained.Owing to the axial symmetry, the electric current flux can be written as The right-hand side of Eq. ( 29) is reduced to However, considering the function of stabilizing the poloidal component, a simpler relation, with f tp O(1) as the ratio between poloidal and toroidal field strengths, might also provide a reasonable estimate.
We have not yet determined the functional form of the electromotive force induced by the α effect, namely, (α • B)/c.In this work, we keep this issue open and thus, our results omit the α-effect.This is because α is a pseudo-tensor, which can be a function of the magnetic field strength, the rotation frequency and the local thermodynamic quantities as well, so that it will be complex in the general case.Nevertheless, the α effect accounts for the induction of the poloidal field from the toroidal field, which is not provided by other means considered in this work.Such a term is therefore indispensable for obtaining magnetic amplification especially in convective regions (Featherstone et al. 2009;Augustson et al. 2016;Hotta et al. 2016).Thus, we will include this term in a future work.
In summary, the evolution of the stellar magnetic field is described by the two equations, It is noteworthy that both Ar and B/ρr scale as |B|r 2 if we assume ρ ∝ r −3 .Therefore, the terms on the left-hand side of both Eqs.((36) and (37)) show the magnetic flux conservation.For Eq. ( 36), the rest term describes the magnetic diffusion.The first, second, and third terms on the right-hand side of Eq. (37) account for the Ω effect, the magnetic diffusion, and the magnetic advection caused by the gradient of the magnetic diffusivity, respectively.

Boundary conditions
We solve two diffusion-advection equations for the magnetic field evolution.Hence, in total, four boundary conditions are needed for closures.At the center of the star, we set A = 0 and B = 0 such that the magnetic field does not diverge.At the surface of the star, we set 1 + (1/A)(∂Ar/∂r) = 0 and B = 0.The first condition is obtained by approximating that the poloidal magnetic field outside the star coincides with the dipole field.The second condition is obtained by assuming that there is no radial electric current that penetrates the stellar surface to outer space.

Evolution equation of the angular velocity
The equation of fluid momentum conservation can be written as where g, Π, and M ≡ 1 4π BB − 1 8π |B| 2 I are the gravity, the Reynolds stress tensor, and the Maxwell stress tensor, respectively.Here, we use a relation 1 c j × B = ∇ • M describing the balance between the Lorentz force and the magnetic stress, which is satisfied for MHD.We note that the viscous stress owing to the fluid viscosity is neglected because of the extremely large Reynolds number of the stellar system.
Starting from the basic equation, one can write the evolution equation of the specific angular momentum as where axial symmetry is assumed.By averaging over a sphere, we obtain where S ≡ 2πr 2 sin θdθ is the surface area of the sphere and i ≡ 2πr 4 sin 3 θdθ/S ∼ 2r 2 /3 is the specific moment of inertia of the sphere.Finally, by assuming the latitudinal dependencies of the magnetic field and by applying an empirical formula of the viscous force due to shear, the evolution equation of the angular velocity is obtained; where ν cv and ν eff are the effective viscosities due to convective turbulence and turbulence induced by other instabilities.Here, it is noteworthy that B φ in the magnetic tensor should omit B stb , the toroidal component that would exist for stabilization of the poloidal field, but only account for the toroidal component that is induced by the Ω effect.This is done for the purposes of consistency with the assumption on the dynamical equilibrium of the original field; otherwise nonzero Lorentz forces would disrupt the original configuration.This requirement is naturally satisfied if we assume that B stb has ∼ sin θ polarity, such that the effect on the Lorentz force is canceled out after averaging over the sphere.
The first term on the right-hand side of Eq. ( 41) describes the angular momentum transfer owing to the Reynolds stress due to convective turbulence.In general, angular momentum distribution in a convective region is affected by the interplay between the convective advection and turbulence and rotation.However, the theoretical treatment is uncertain (e.g., Tassoul 2000).Thus, our model incorporates the parameter, n cv , which indicates which kind of angular velocity structure forms in a convective region in equilibrium.Two extreme cases are either n cv = 0 (favoring rigid body rotation) or n cv = 2 (favoring isotropic specific angular momentum).We apply n cv = 0 for our calculations unless otherwise noted.For the second term, we assume that a region reaches rigid rotation in equilibrium because shear rotation is the source of the energy that drives instability in most of the considered cases.The third term accounts for the angular momentum transfer due to the Maxwell stress.
It is noteworthy that all of these stress terms have a conservative form, that is, only the surface term appears after integrating over the whole star.For the magnetic stress, this property originates from the fact that the electromagnetic field, assuming A19, page 6 of 28 MHD, only possesses a negligible amount of momentum compared to the matter.We also note that even in the case of n cv = 2, the rotation law does not necessarily reach the case of isotropic specific angular momentum because the second term counterbalances towards a rigid rotation.In our code, the most efficient rotation induced viscosity in a convective region is the dynamical shear instability.Considering the huge uncertainty involved in the efficiency estimates, we assume that ν DS = ν MLT in case of n cv = 2.This yields Ω(r) ∝ r −1 , which can be checked by applying ν cv = ν eff to Eq. (41).

Diffusion coefficients
We assume that several hydrodynamical instabilities may develop in the rotating stellar interior and that turbulence driven by such instabilities accounts for the Reynolds stress.The convective viscosity is estimated as ν cv = D cv and ν eff consists of six terms, namely: where each term stands for the viscosity corresponding to the Eddington-Sweet circulation (ν ES ), the dynamical shear instability (ν DS ), the secular shear instability (ν SS ), the Solberg-Høiland instability (ν SH ), the Goldreich-Schubert-Fricke (GSF) instability (ν GSF ), and the Pitts-Tayler instability (ν PT ), respectively.The viscosity coefficients ν ES , ν DS , ν SH , and ν GSF are calculated according to Pinsonneault et al. (1989) and Heger et al. (2000), with the modification that we use the minimum of the pressure scale height, H P , and the radius to estimate the typical length scale for each instability.In the original works, they use the velocity scale height of the respective flow, which is further limited by the radius or the width of the unstable region, as the typical length instead.We find that this modification has only a limited effect on the overall stellar evolution.To compute ν SS , we follow the prescription by Maeder (1997).An m = 1 instability is assumed to grow in a region with a strong toroidal magnetic field.The effective viscosity owing to this Pitts-Tayler instability, ν PT , is estimated according to Spruit (2002) and Maeder & Meynet (2004).For clarity, we give the corresponding equations in Appendix B. We note that these prescriptions involve a control parameter, f µ , which is multiplied to the µ-gradient.This is an influential parameter of the stellar simulation, as it affects the stability conditions.
We assume that turbulence driven by (magneto-)hydrodynamical instabilities accounts for the chemical mixing as well.Another control parameter f c is set, which indicates the ratio between the chemical diffusivity and the viscosity.Thus, is used for rotating models.Similarly, the effective magnetic viscosity is estimated as however, f m = 1 is set in the current work.We note that the Eddington-Sweet circulation accounts for the most efficient "turbulent" diffusivity in the radiative envelope in the present models.The Eddington-Sweet circulation, which is also referred to as the Eddington-Vogt circulation, was first postulated as a laminar meridional flow driven by a thermal imbalance that results from the difference of temperature gradient in a latitudinal direction in a rotating star (von Zeipel 1924a,b;Eddington 1925;Vogt 1925;Sweet 1950).As the effect on the angular momentum transport could be modeled as an advection (cf.Maeder & Zahn 1998), magnetic advection due to the Eddington-Sweet circulation would be formulated by more stringent consideration of the U × B term in Eq. ( 24).However, considering the existence of the baroclinic instability, which will operate with a dynamical timescale in a radiative zone even with a very small differential rotation (Fujimoto 1988;Kitchatinov 2014), it will also be natural to assume that the Eddington-Sweet circulation in an actual star would typically be accompanied by turbulence.Bearing the uncertainties involved in the theoretical modeling in mind, we leave this discussion open in the current work.

Wind-magnetic field interaction
With a strong surface magnetic field, part of the stellar wind blowing from a closed field region will be trapped to form a magnetosphere surrounding the star (Donati et al. 2002).Indeed, the time variation and the Balmer-line emission profiles observed in a well-known magnetic Bp star, σ Ori E, have been reproduced by considering such a rigidly rotating magnetosphere (Townsend & Owocki 2005;Townsend et al. 2005).Because the net mass loss rate can be significantly reduced due to the magnetic confinement, we take this effect into account in our simulations according to the work of Ud-Doula et al. (2008).The magnetic confinement parameter, η * , is estimated first and then used to derive the confinement efficiency f conf (η * ) ≡ Ṁ/ Ṁ(B = 0), where Ṁ(B = 0) is the mass loss rate of a non-magnetic star.The detailed procedure is explained in Appendix C.
While the strong surface field reduces the mass-loss rate, it can enhance the rate of the angular momentum loss in contrast.This is because the stellar angular momentum is not only reduced by the material flow but also by the Maxwell stress (Weber & Davis 1967;Ud-Doula et al. 2009).The braking efficiency f break (η * ) ≡ J/ J(B = 0), where J(B = 0) is the angular momentum loss rate of a non-magnetic star, is estimated in this case.We account for the effect of the magnetic braking according to Ud-Doula et al. (2009).The interaction of wind material with a dipole magnetic field is assumed in their analysis, while a simple monopole geometry is assumed in the classic analysis by Weber & Davis (1967), which is often used to model the evolution of solar-like stars.We give further details of our treatment in Appendix C.

Numerical settings and code test
Equations ( 36), (37), and (41), and the diffusion part of Eq. ( 1) are numerically solved with a finite-difference method.We use a first-order backward difference for the time derivative and a second-order central difference for the space derivative.Together with the boundary conditions, these difference equations are iteratively solved at the same time and while they are decoupled from the equations of stellar structure.One time step for the structure equations is further divided into numerous (typically ∼1000) time steps for the evolution equations of the magnetic field and the angular momentum.The latter time step is controlled such that the relative differences in A, B, and Ω are restricted to be smaller than λ in each step, where λ is an arbitrarily chosen control parameter of about 10%.The basic features of the numerical code such as magnetic flux conservation and magnetic dissipation are confirmed, and details of the code tests are described in Appendix D.

Other possible magnetic effects
Several other magnetic effects are not considered in the present work.Even though they do not significantly affect the interplay between the evolution of the magnetic field and the stellar rotation, they can have a significant effect on the stellar structure in some cases.Here, we briefly review these effects, which we plan to implement on top of the present formulation in forthcoming papers.
For strong magnetic fields, the stellar structure can be modified by the magnetic pressure and the magnetic tension.Feiden & Chaboyer (2012, 2013, 2014) take these effects into account for low-mass stellar evolution calculations using a geometry parameter introduced by Lydon & Sofia (1995).Duez et al. (2010) develop a more rigorous and general treatment and apply their formulation to model the young Sun.
Lydon & Sofia (1995) also showed how large-scale magnetic fields can affect the equation of state, especially the adiabatic index, and the equations in the mixing length theory for convection.When the local magnetic field is strong enough, the pressure change during an adiabatic motion can differ from the non-magnetic case since part of the work goes into the form of magnetic energy.The modification of the adiabatic index further affects the criterion for convective instability.In addition, it changes the specific heat and thus the efficiency of the convective energy transport.
The internal energy equation is in principle also modified when considering the magnetic effects of Jule heating and the Poynting flux.Because our formulation includes magnetic dissipation due to turbulent magnetic diffusivity (the η effect), a corresponding Jule heating term may be taken into account in the internal energy equation for consistency in the future.
Strong and stable magnetic fields inside the star might suppress hydrodynamical flows such as convection and meridional circulations.The stabilization effect may be included in stellar evolution models by modifying the convective criterion (e.g., Lydon & Sofia 1995;Petermann et al. 2015).

Angular momentum transport via dissipating torsional Alfvén wave
Differential rotation winds up the poloidal magnetic field to enhance the toroidal component, which is the Ω effect.As the toroidal component gets strong, the magnetic stress increases as well, which counteracts to reduce the differential rotation.Here, we consider what would happen in a stellar model when these two effects are incorporated.We can simplify our set of equations such that the Ω effect and the magnetic stress are resolved by two linear differential equations as where the effects of magnetic diffusion and viscous angular momentum transport are neglected, the rates of change of Ar, radius, density, and specific moment of inertia i are assumed to be small, and the relations i ∼ 2r 2 /3 and A ∼ rB r /2 are used.This hyperbolic system may be analyzed with methods used in fluid dynamics.By diagonalizing the matrix we are able to obtain a set of eigenvalues and eigenvectors as ±c ≡ 1 √ 5 v A and r ± = (1 1/ 5πρr 4 ) t .The corresponding invariant dw ± = dΩ ∓ d(Br 3 )/ 5πρr 4 becomes constant along the characteristic dr/dt = ±c.Here, v A ≡ B r / 4πρ is the Alfvén velocity of the radial magnetic field.Therefore, we expect that a wave that propagates with the Alfvén velocity forms in this system.
In this section, we analyze how this wave propagation manifests itself in our numerical models.Furthermore, we demonstrate that with the help of viscosity and magnetic dissipation, this torsional Alfvén wave serves as a highly efficient mechanism for the redistribution of angular momentum.

Formation and propagation of torsional Alfvén wave
We explore the wave propagation using a 1.5 M main-sequence model which has a radius of ∼1.5 R .We calculated the coupled evolution of the toroidal field and the stellar rotation, including the Ω effect and the magnetic stress, but we set the viscosity and magnetic diffusivity to zero.A uniform radial field as B r (=2A(r)/r) = 1 kG is set in the beginning, so that the wave velocity becomes c ∼ 80 cm s −1 , and correspondingly, the estimated wave-crossing time from the center to the surface is ∼1.3 × 10 9 s.A step-function distribution of Ω = 10 −4 rad s −1 for M ≤ 1 M , and Ω = −10 −4 rad s −1 otherwise, is imposed as the initial angular velocity distribution.The toroidal magnetic field B φ (=B(r)) is set to be zero everywhere.
Figure 2 shows how the rotation and the magnetic field evolve with time.In the beginning, two wavefronts launch from the discontinuity and start propagating towards the stellar surface and the center.A wave, which is similar to a rarefaction wave in fluid dynamics, is formed in between the two wavefronts.The up-going front reaches the surface at ∼3 × 10 8 s, and is then reflected to follow the down-going wavefront.The downgoing front reaches the stellar center at ∼1.4 × 10 9 s and then is also reflected.At ∼1.8 × 10 9 s, the two wavefronts cross and penetrate each other.
Although the wavefronts have a diffuse structure, especially in the central region, Fig. 1 shows that their propagation agrees well with the prediction, which is shown as black dashed and dotted lines.The positions of the two lines are directly calculated by the time integral of The wave-crossing time in the simulation is estimated to be 1.71 × 10 9 s, as the two wavefronts meet again at this time after traveling either through the stellar surface or the center.This agrees well with the simple estimate of ∼1.3 × 10 9 s given above.The timescale of shear rotation does not appear in the propagation timescale because the efficiency of the winding-up of the poloidal magnetic field is proportional to the required torque to affect the angular momentum of the material, and therefore they cancel each other out.However, the timescale of the Ω effect relates to the time that the wave passes the width of the wavefront, λ/c, and therefore, the strength of the toroidal magnetic component at the wavefront can be estimated as This is proportional to the angular velocity difference at the wavefront, ∆Ω, but is independent of the strength of the poloidal A19, page 8 of 28 magnetic component.This is because the stronger the seed poloidal magnetic field is, the shorter is the wave-crossing time of the width of the wavefront, and therefore they cancel each other out.This relation yields B φ ∼ 7 × 10 7 G at M r = 1 M with ∆Ω = 2 × 10 −4 rad s −1 and explains the simulation result.
In another test calculation with a ten-times stronger initial poloidal magnetic field, the wave velocity increases by a factor of 10, but the toroidal magnetic component does not change.On the other hand, in a calculation with a ten-times smaller initial angular velocity, the toroidal magnetic component decreases by a factor of 10, but the wave velocity stays constant.
The standard Alfvén wave and the Ω-B wave discussed here are essentially identical, as they share the same driving force and a similar propagation velocity.Hence, we refer to the wave solution in our simulation as a torsional Alfvén wave.While a small fluctuation propagates along the magnetic field in the for-mer case, a large number of windings are required to launch the wave in the latter.This difference arises from the weak magnetic field considered in the present case.We compute ∼10 15 erg g −1 for the gravitational and thermal energies, ∼10 12 erg g −1 for the rotational kinetic energy, but only ∼10 4 erg g −1 for the magnetic energy of the poloidal component.With such a weak magnetic field, a strong magnetic amplification due to the Ω effect is required for the magnetic stress to affect the dynamics of the rotating flow.The toroidal component induced in the Alfvén wave always has comparable specific energy to the rotational kinetic energy; see Eq. ( 46).The required number of windings decreases for higher initial poloidal field strengths.The torsional Alfvén wave will converge to the standard Alfvén wave if the poloidal component is so strong that a small number of windings is sufficient to drive the wave.A19, page 9 of 28 A&A 646, A19 (2021 9.9-10.3×10 9sec with hydro inst.For this stellar model the wave-crossing time is t wc = 1.71 × 10 9 s.The top seven plots of each panel cover about two wave-crossing times (3.4 × 10 9 s) with five snapshots per plot, one every 10 8 s, for the period given in the top right corner of each plot.In the bottom row of plots, results after about six wave-crossing times are shown, i.e., from 9.9 × 10 9 s to 10.3 × 10 9 s.The temporal evolution is available as an online movie.

Torsional Alfvén wave with dissipation
Here, we discuss the wave propagation for calculations where dissipation effects are included, using the same stellar model and initial conditions as in the previous section.Figure 3 shows the resulting evolution of angular velocity, toroidal magnetic, and magnetic viscosity.The cases of neglecting dissipation effects (left) can be compared with results including dissipation due to hydrodynamic instabilities only (middle), and considering both hydrodynamic and magneto-hydrodynamic dissipation (right).
As shown above, without dissipation, the waves travel freely through the star, and the wave-fronts meet each other once per wave-crossing time.Due to numerical diffusion, the angular velocity distribution after two wave-crossing times (seventh plot in the first column), which should coincide with the initial distribution, has become somewhat more diffuse.Nevertheless, we can follow the back-and-forth sloshing of the wave for more than 30 crossing times.
The middle column of Fig. 3 shows that turbulent dissipation due to hydrodynamic instabilities affects wave propagation.Whereas turbulence due to secular shear and GSF instabilities accounts for some viscosity in the radiative envelope, this is too small to matter here.However, our 1.5 M model has a hydrogen-burning convective core, in which rigid rotation is established by the large convective viscosity (we note that we apply n cv = 0 in this calculation).This has a big impact on the angular momentum redistribution.Since waves with shorter wavelengths have shorter dissipation time, the convective region effectively filters out waves with shorter wavelengths, which originally compose the step function used for the initial condition.As a result, a standing wave is quickly formed in the model, with a wavelength of twice the stellar radius.Hence, the standing wave corresponds to the n = 1 fundamental mode oscillation, which has one node at ∼1.1 M (see the online movie of Fig. 3).
In the calculation including hydrodynamic and magnetohydrodynamic dissipation, the latter works effectively only A19, page 10 of 28 no dissp.w/ hydro w/ hydro+PT Fig. 4. Surface angular velocity as a function of time in test calculations of our 1.5 M magneto-rotational main-sequence model, in which constant radial magnetic field of B r = 1 kG and a step function with an arbitrary amplitude for the initial angular momentum distribution are imposed as initial conditions.The star starts to oscillate on the Alfvén time scale.Lines correspond to the case without dissipation (red, solid), including only hydrodynamic instabilities (green, dashed), and the case with hydrodynamic and magneto-hydrodynamic instabilities (blue, dash-dotted), respectively.The wave-crossing time of t wc = 1.71 × 10 9 s is referenced by the thick black bar at the bottom left corner.
during the first one or two wave-crossing times.First, shear rotation propagates inside the star with the velocity of ∼v A / √ 5.In regions passed by this shear wave, a toroidal field is induced by the Ω effect, which soon triggers the Pitts-Tayler instability.As a result, the large diffusivity due to the Pitts-Tayler instability covers the entire star at t ∼ t wc , which diffuses the angular velocity very effectively, together with the convective diffusivity.At later times, however, the diffusivity due to the Pitts-Tayler instability becomes weak since there is no more strong shear rotation amplifying the toroidal field component.Thus, the overall evolution becomes comparable to the case that only considers the hydrodynamic instabilities.
Figure 4 shows the first 500 yr evolution of the surface angular velocity of our model.The damped, quasi-sinusoidal variation for the models including dissipation is due to the torsional oscillation.The period is comparable to the wave-crossing time t wc , which is given by dr with the stellar radius, R. If such waves were excited, such oscillations could exist in real stars.Because the oscillation timescale relates to the internal density and poloidal magnetic field distributions, observations of the changing surface rotation frequency may allow deriving the internal magnetic field strength by observing the change in the surface rotation frequency.We discuss the comparison between our model and relevant observations in Sect.7.2.The torsional Alfvén wave oscillation gradually decreases its amplitude.The model including only hydrodynamic instabilities approaches a trivial stationary state of ∂Ω/∂r = B φ = 0 with a decay timescale of t ∼ 10 t wc .This means that the angular momentum in the star is effectively redistributed to achieve rigid rotation, as torsional Alfvén wave propagates and dissipates throughout the star.
It is noteworthy that an integrated evolution of the toroidal magnetic field and the rotational flow in a radiative region of the Sun has been modeled by 2D axisymmetric simulations in Charbonneau & MacGregor (1992, 1993).In their simulations, the phase shift across poloidal field lines can be followed.This leads to efficient wave dissipation because large gradients in the toroidal field are formed (Charbonneau & MacGregor 1992).Although there is a difference in the detailed mechanism of the dissipation, their results show that the radiative region in the sun approaches rigid rotation as the dissipative Alfvén wave propagates, which is consistent with our model.
The onset of Pitts-Tayler instability could affect the wave propagation if the turbulent dissipation modifies the wavefront structure.Therefore, we consider a differentially rotating region in a star, where considerable strength of the poloidal field exists but initially zero toroidal component.The condition for the Pitts-Tayler instability to grow within a propagation time of the Alfvén wave is τ PT < τ wc , where τ PT is the growth time of the Pitts-Tayler instability and τ wc is the wave-crossing time.When the toroidal component is induced by the Ω effect with a differential rotation parameter q ≡ ∂ ln Ω/∂ ln r, the growth time of the Pitts-Tayler instability can be estimated as τ PT = √ 3τ 2 wc /(Ωq).Hence, the condition above can be expressed as √ 3Ωτ wc q > 1.Our simulation with a strong poloidal field of B r ∼ 1 kG has τ wc ∼ 10 9 s.Therefore, the condition can be satisfied with a canonical value of Ω ∼ 10 −5 s −1 , although q can have a variety of value of 1.
Nonetheless, the condition for the wavefront to be disturbed by the growing turbulence is normally not satisfied in a radiative stellar envelope.For the turbulence driven by the Pitts-Tayler instability to affect the wavefront, l v > v A τ PT may be required.The vertical length scale of the turbulence, l v , is estimated as r(ω A/N ), using the toroidal Alfvén angular frequency, ω A ≡ B tor / 4πρr, and the Brunt-Väisälä frequency, N (Spruit 2002).Thus, the condition becomes Ωq/N > 1.For our 1.5 M model, we find N ∼ 10 −2 s −1 in the radiative envelope, which implies that in most of the case this condition will not be satisfied.Also, this condition does not depend on the poloidal field strength.Therefore, in our simplified 1D picture, although the Pitts-Tayler instability is expected to grow at the wavefront, the wave propagation will not be affected by the turbulence induced by the instability because the unstable region will be too thin.

Main-sequence evolution of 1.5 M stars
To explore the capabilities of our new modeling approach, we calculate the evolution of solar metallicity 1.5 M models with the full framework as described in Sect. 2. We set α MLT = 1.8, f ov = 0.01, f µ = 0.1, and f c = 0.125, for the mixing-length parameter, the overshoot parameter, the µ-barrier parameter, and the chemical diffusion/viscosity ratio parameter, respectively.Initially, rigid rotation is applied, and the initial rotation period is chosen from P rot,ini = 1, 10, 100, and 1000 d.For the magnetic field, we apply the simplest possible functions as This type of vector potential yields a uniform poloidal magnetic field along the rotation axis inside the star.The strength of the surface magnetic field at the pole, B p,ini , is chosen as 10, 100, 1000, or 10 000 G. The model characteristics are summarized in Table 1.
In the left panel of Fig. 5, evolution in the HR diagram are compared for rapidly rotating (P rot,ini = 1 d) 1.5 M models with Notes.τ MS is the main-sequence lifetime, ∆M MS is the total mass lost during the MS phase, Ṁ MS = ∆M MS /τ MS is the averaged mass loss rate for the MS phase, τ break,ZAMS = −J ZAMS / JZAMS is the braking timescale measured at ZAMS, η * ,ZAMS is the magnetic confinement parameter at ZAMS, and P rot,TAMS , B p,TAMS , and J TAMS are the rotation period, the field strength at the pole, and the total angular momentum at TAMS.B p,ini = 10 G and 10 kG and without magnetic effects.In the right panels, internal profiles of angular velocity, hydrogen mass fraction, radial and toroidal magnetic field strength, and the ratio between the magnetic pressure and the pressure for the same models at central hydrogen mass fractions of X c = 0.1 are shown.In Fig. 6, the reconstructed 2D magnetic field structures for the A19, page 12 of 28 weakly magnetic (B p,ini = 10 G) model at X c = 0.7 and 0.1 are shown as well.
The fast rotation with P rot,ini = 1 d, which accounts for ∼20% of the Keplerian rotation at the surface, supports the stellar surface, reducing the effective temperature.However, the surface velocity of the strongly magnetic (B p,ini = 10 kG) model quickly decreases because of the efficient magnetic braking.This explains the offset in the HR diagram in the early main sequence, in which the strongly magnetic model shows a slightly higher effective temperature than the others.
During the main-sequence evolution, the stellar envelope expands, whereas the convective core shrinks.As a result, in the model without magnetic fields, significant differential rotation develops in the radiative envelope.The hydrogen-burning core is slightly extended by rotation induced mixing due to the secular shear instability.On the other hand, magnetic models evolve remaining close to rigidly rotating due to the highly efficient angular momentum redistribution by the dissipating torsional Alfvén wave.Even in the weakly magnetic model with B p,ini = 10 G, which develops a toroidal field of only ∼−10 G in its radiative envelope, only a tiny amount of differential rotation of ∂Ω/∂ ln r ∼ 10 −10 s −1 develops at the core-envelope boundary.This behavior appears consistent with the results from asteroseismology which finds only a small deviation from rigid rotation in main-sequence stars (e.g.Aerts et al. 2017).Consequently, no rotation induced mixing develops close to the hydrogenburning core in the magnetic model.This explains the ∼5% fainter terminal-age main-sequence (TAMS) luminosity of the weakly magnetic model.The difference in the TAMS luminosities between B p,ini = 10 G and 10 kG models will be explained by the slight difference of the total TAMS masses.
The poloidal magnetic field strength in the magnetic model decreases over time, partly due to the envelope expansion, but more importantly, as a result of magnetic dissipation due to rotation induced turbulence.Here, turbulence triggered by the Eddington-Sweet circulation, the efficiency of which is roughly proportional to the square of the rotation frequency but not to the shear rotation, accounts for the magnetic dissipation.The ratio between the magnetic pressure, B mag = (B 2 r + B 2 θ + B 2 φ )/8π, and the pressure (or the inverse of the plasma-β) is shown in the right-bottom panel of Fig. 5.This ratio will also indicate the significance of the influence of the magnetic pressure and tension and magnetic modification of the adiabaticity to the stellar structure, which are not currently treated in our magnetic models (Sect.2.6).This figure shows this ratio is fairly small in particular in the inner region of the star, justifying the current treatment.Meanwhile, a strong magnetic field may influence the surface structure, which will indirectly affect the stellar evolution by changing the mass-loss and angular-momentumloss histories.Also, it is noteworthy that material mixing due to the magnetic Pitts-Tayler instability does not take place during the main-sequence evolution in magnetic models because it requires differential rotation.Therefore, in our current prescription, magnetic fields with different strength affect the evolution only in that they change the rotation velocity and, thus, the rotation induced mixing.
The evolution of the surface rotation periods for our models with different initial spins and magnetic field strengths is shown in Fig. 7. Models without wind mass loss (black solid lines) conserve their initial angular momenta.This figure shows that the stronger the surface magnetic field, the stronger the magnetic braking takes place.With the adopted mass loss rate, a B-field of ∼100 G is strong enough to spin down the stars by an order of magnitude within their main-sequence lifetime, except for our fastest rotating model.The efficiency of the magnetic braking scales with the strength of the surface magnetic field.Thus, the A19, page 13 of 28 spin-down of the models with B p,ini = 1000, and 10 000 G is ∼10 and 100 times faster than that of the models with B p,ini = 100 G.However, the magnetic braking gets weaker in our fastest rotating models (P rot,ini = 1 d).This is because the surface magnetic fields are dissipated due to the efficient η effect.
Figure 8 shows the time evolution of the polar surface magnetic strength B p .For comparison, the results of models without the η effect, which have P rot,ini = 1000 d, are also shown (black solid lines).In these, the surface magnetic field decreases due to mass loss.When a mass is lost, the layers below the surface expand, and magnetic flux conservation leads to weaker fields.A much more significant drop in magnetic field strength takes place in models with the η effect.We see that the field strengths decrease faster for more swiftly rotating models, which argues for dissipation driven by rotation induced turbulence.For example, for B p,ini = 10 G (top-left panel), the strongest dissipation takes place in the model with P rot,ini = 1 d, but for models with P rot,ini = 1000 and 10 000 d, the surface magnetic field evolves as it would in a model without an η effect.A similar behavior is obtained for models with B p,ini = 100 G.
Our most highly magnetic models (B p,ini = 1000 and 10 000 G) follow a more complicated evolution.The two slowly rotating models maintain stronger surface fields than the model without the η effect.This is so because the strong surface mag-netic fields suppress stellar wind mass loss by the magnetic confinement.This also takes place in the model with P rot,ini = 10 d, but not in the fastest rotating one (P rot,ini = 1 d).Here, the surface magnetic field is quickly dissipated by rotation induced turbulence.We conclude that the surface magnetic field strength is affected by wind mass loss and by rotation induced dissipation.However, at the same time, these two effects are influenced by a strong surface magnetic field, that is, the wind mass loss is suppressed by magnetic confinement and the rotation is also slowed down by magnetic braking.
The evolutionary paths of B p and P rot are shown together in Fig. 9.As dots are placed on the tracks after a constant elapsed time, the density of dots in this figure is somewhat representative of the likelihood of observing a star at a given location.The very low density in the lower right corner of Fig. 9 implies that the likelihood of finding a strongly magnetized star (B p > 1000 G) with a rotation period below ∼1000 d is very small.This is the case given that for the considered field strength, the magnetic spin-down timescale becomes shorter than ∼1% of the main-sequence lifetime.At the same time, the likelihood of finding a rapidly rotating (P rot < 1 d), strongly magnetized star (B p > 100 G) is also low because fast rotation results in efficient magnetic dissipation.
A19, page 14 of 28 A19, page 15 of 28 Figure 9 also indicates the evolution of the total amount of mass lost by the models due to stellar winds.The model with B p,ini = 10 G and P rot,ini = 1000 d is best understood, as it has the most modest magnetic confinement and negligible rotational mass-loss enhancement.This model loses 5.9 × 10 −3 M during its main-sequence phase.Comparably, the models with B p,ini = 10 000 G lose only 2.1 × 10 −4 to 3.4 × 10 −4 M thanks to efficient magnetic confinement.On the other hand, the model with B p,ini = 10 G and P rot,ini = 1 d experiences magnetic dissipation and, therefore, the magnetic confinement on the model is minimal.Its wind mass loss is enhanced by the fast rotation and it loses 1.33 × 10 −2 M .The rotational enhancement also takes place for other models with P rot,ini = 1 d, but with the stronger initial surface fields of B p,ini = 100 and 1000 G, the lost amounts are reduced because either the magnetic confinement or the reduction of the rotational enhancement due to the magnetic braking also taking place.The lost masses are 8.3 × 10 −3 M for the model with B p,ini = 100 G and P rot,ini = 1 d and 1.8×10 −3 M for the model with B p,ini = 1000 G and P rot,ini = 1 d.We conclude that by considering a surface magnetic field, the total wind mass loss displays a complex behavior because of the interplay between wind mass-loss rate, rotation, magnetic confinement, and magnetic spin-down, even for a fixed stellar mass and metallicity.

Red-giant branch evolution of 1.5 M stars
Here, we describe the results of our model calculation of a 1.5 M star of solar metallicity from the zero-age main sequence (ZAMS) up to core helium ignition at the tip of the red-giant branch, where the model experiences a violent helium flash and the calculation is concluded.We include the same physics as that used in the models of the previous section.Here, we use an initial rotation period of P rot = 1.4 d, a zero toroidal field, and a uniform poloidal field with B p,ini = 10 G.These values may be representative of normal A-type stars.Our fiducial model is computed using n cv = 0, and later we will discuss results obtained with n cv = 2 for an otherwise identical model.
The evolution of the fiducial model in the HR diagram is shown in the left panel of Fig. 10.We divide the evolution into five phases, which are indicated by lines with different colors (red, orange, yellow, green, and blue, from the beginning) in the figure.Representatives, including the model at helium ignition are selected from each phase and shown as dots on the HR diagram.Corresponding internal profiles of physical key properties are shown in the right panels using the same line color (dark-blue for the model at helium ignition).
The top-left figure shows the increasing density contrast between the core and the envelope.The angular velocity profiles demonstrate that rigid rotation is maintained even after the core contraction (yellow and green) until the core growth phase (blue) due to the efficient magnetic angular momentum transfer.This efficiency can be estimated through the crossing time of the torsional Alfvén wave, which is defined as τ wc = r 20πρ/B r (middle-right panel).Profiles of the radial magnetic field strength are shown in the top-right panel.The radial component keeps 0.1 G until envelope convection develops and, consequently, the local wave-crossing time remains to be shorter than 10 4 yr.
The magnetic diffusivity is significantly enhanced in the convective envelope.Hence, it decreases the radial component after the development of envelope convection, increasing the local wave-crossing time (green and blue).Finally, it becomes longer than the evolutionary timescale of ∼10 Myr.Consequently, the model starts to develop differential rotation during the core growth phase (blue).
The rotation of the red-giant core significantly accelerates during this phase.The acceleration (of more than two orders of magnitude) can be evaluated as follows: at first, the matter accreted by the helium core has a large specific angular momentum of where R base and Ω base are the radius and the angular velocity of the base of the convective envelope, respectively.Then, the total accreted angular momentum can be estimated as where M acc ∼ 0.14 M is the accreted amount of mass.This accreted amount of angular momentum largely exceeds the original core angular momentum of ∼10 43 g cm 2 s −1 .Since the magnetic field maintains rigid rotation in the helium core, the saccreted angular momentum is quickly redistributed throughout the core.Assuming that the accreted amount of angular momentum is redistributed within a core of M core = 0.43 M and R core = 1.6 × 10 −2 R , we expect a core angular velocity of which corresponds well to the simulation result.As a consequence of the accelerated core rotation, strong differential rotation occurs at the core-envelope boundary.Consequently, a strong toroidal field of ∼10 6 G is induced.

Models including the Tayler-Spruit dynamo
As the significant impact of magnetic fields on the internal angular momentum distribution of stars has become more and more evident during the last two decades (cf., Sects. 1 and 7 below), many stellar evolution models have attempted to account for this by incorporating magnetic angular momentum diffusion as proposed by Spruit (1999Spruit ( , 2002) ) in various versions (Maeder & Meynet 2003, 2004;Heger et al. 2005;Yoon & Langer 2005;Denissenkov & Pinsonneault 2007;Suijs et al. 2008;Brott et al. 2011a;Yoon et al. 2012;Fuller et al. 2019).This model relies on a dynamo picture, the so-called Tayler-Spruit dynamo, which is assumed to operate in differentially rotating, radiative layers in stars.
The original picture of the Tayler-Spruit dynamo consists of three steps.First, due to the Ω effect, a strong toroidal magnetic field develops in a radiative layer due to differential rotation.Second, the Pitts-Tayler instability induces turbulence in this region.Finally, a radial magnetic field with considerable strength is induced by the turbulent stretching of the toroidal field.The induced radial component plays a role as the next seed field of the Ω effect such that it closes the dynamo loop.
A19, page 16 of 28 Following this picture, we may estimate the strength of local magnetic components of B r and B φ as well as the magnetic stress S TS = B r B φ /4π.The magnetic stress is further rewritten in terms of the viscosity, ν TS , by equating ρν TS r∂ r Ω = S TS .This magnetic viscosity is then used in evolutionary calculations in the form of a diffusion coefficient for angular momentum transport.The Tayler-Spruit picture is incompatible with our modeling in two aspects.First, the time evolution of the toroidal component B φ is significantly simplified such that B φ is directly proportional to ∂ r Ω.This simplification leads to a diffusion approximation of the magnetic stress, so that angular momentum redistribution takes place locally as a form of diffusion, while the most natural consequence of the Lorentz force and the Ω effect is the formation of the torsional Alfvén wave (Sect.3).Second, the Tayler-Spruit dynamo relies on the α effect of the Pitts-Tayler instability, such that the "radial" component that contributes to the Lorentz force is dominated by the secondary generated field.In contrast, the poloidal component that contributes to the wave propagation in our case is the original field and a hydrodynamic induction to reproduce the poloidal component is not assumed.
It is evident that when a star has a structured poloidal magnetic field, the dominant phenomenon occurring after the Ω effect will be wave propagation, at least for the simplified 1D geometry adopted here.For a structured but weak poloidal field, the field strength only affects the timescale, and the torsional Alfvén wave will still form.This is because the vertical length scale of the Pitts-Tayler instability is always smaller than the wave-crossing length scale (see Sect. 3.2).
Therefore, the important question is bound to consider what happens in a star that has a very weak and unstructured initial poloidal field (as a strong and unstructured poloidal field would be unstable and does not exist).Also, in this case, a strong toroidal component may develop due to the Ω effect in a region with a differential rotation.Because the initial poloidal is weak and unstructured, the torsional Alfvén wave launched by the Ω effect will be locally trapped.Then the region will be predominantly affected by the Pitts-Tayler instability.Turbulence may affect the magnetic field through both the α and the η effect.
If the η effect wins, then the initially weak poloidal component will soon dissipate, so that the region will have a pure toroidal magnetic field.Because purely toroidal fields are unstable (leading to the development of such instabilities as the Pitts-Tayler instability; e.g., Tayler 1973), eventually, all the magnetic energy will dissipate into heat.Thus, the magnetic field does not affect the differential rotation in the region.If the α effect wins, A19, page 17 of 28 the poloidal field is amplified to have a complicated 3D structure that is embedded within the instability region.Characterized by a short crossing time, torsional Alfvén waves will now propagate along the amplified poloidal field, and non-linear interaction (e.g., phase mixing) may efficiently dissipate the waves.Since there will be no preferred direction for the induced poloidal field, the dissipation will take place nearly chaotically.This might result in a similar outcome to what is discussed in the Tayler-Spruit dynamo.
Therefore, the balance between the η and α effects is of particular importance.Whether the α effect works in this situation as assumed in the Spruit-Taylor picture is debatable since contradictory results based on multi-dimensional MHD simulations have been reported (Braithwaite 2006;Zahn et al. 2007).Unfortunately, multi-dimensional simulations with realistic thermal and magnetic diffusivities, which are desirable for testing the dynamo picture, remain a challenge given the current computational resources (see discussion in Braithwaite & Spruit 2017).The question of whether the Tayler-Spruit dynamo works in a star remains unclear at this time.We compare both the results obtained with the Taylor-Spruit dynamo and our results with the relevant observations in Sect.7.

Other approaches
Several works have incorporated other magnetic effects than the magnetic viscosity based on the Tayler-Spruit dynamo into stellar evolution simulations.The interaction between a surface magnetic field and the stellar wind has been considered in the context of modeling the evolution of massive stars (Meynet et al. 2011;Petit et al. 2017;Georgy et al. 2017;Keszthelyi et al. 2019); in particular, for studying the expected spin-down due to their intrinsically large mass loss rate.While these works have revealed the significance of the effects of magnetic braking or magnetic confinement, they do not account for the evolution of the surface magnetic field but, rather, they assume either a constant magnetic field strength or constant magnetic flux during the evolution.Our simulation implies that neither assumption can be realistic because the surface magnetic field may change due to other mechanisms in addition to the flux conservation.In particular, the mass loss will affect the surface magnetic field for massive stars, as it replaces the surface material with a matter that originally stayed below the surface, which has a different magnetic flux compared to the original surface.The Ohmic decay of the field can also be significant.Therefore, a self-consistent global simulation for the stellar magnetism may improve the current estimate of the field interaction with the stellar wind in massive stars.
Furthermore, in some of the above-quoted calculations, the assumptions for the internal magnetic fields are unrelated to the assumptions for the surface field.For example, Meynet et al. (2011) present massive star models with a strong surface field and spin-down but while assuming that no B-field is present inside the star, with the consequence of introducing strong internal differential rotation.Such inconsistencies are avoided with the present approach.
The models of Feiden & Chaboyer (2012, 2013, 2014), and Feiden (2016) take into account several magnetic field effects, such as the magnetic pressure and the magnetic tension, and especially the efficiency change of the convective energy transport.The stellar structure of low mass stars with convective envelopes is shown to be sensitive to these effects.However, these models also lack a detailed theory of the evolution of the stellar magnetic field.Instead, they assume a constant surface magnetic field with a simple radial profile for the internal field strength.The substantial effects on the stellar structure shown by their models argue in favor of the importance of incorporating an appropriate evolution theory for the stellar magnetic field in stellar evolution calculations, which satisfies the essential MHD conditions, such as flux conservation and divergence-free magnetic field configurations.
Such global simulations were performed for the first time by Potter et al. (2012) and their formalism has also been used in their later works (Quentin & Tout 2018).Although the physics included in the modeling is similar to ours, their formulation has two fundamental shortcomings.Their evolution equations for the magnetic field do not reproduce the magnetic flux conservation and, similarly, the angular momentum conservation is not guaranteed with their expression of the Lorentz force.This could be because their choice of the surface-averaging of the original 3D expressions is too simple: it is likely that simple weighted surface-averaging is applied for both the mean-field MHD-dynamo equation and the 3D Lorentz force.
Interestingly, Potter's models successively reproduce the observed population of slowly rotating but nitrogen-enhanced massive stars (Hunter et al. 2008).These stars cannot be obtained by standard rotating single stars models (Brott et al. 2011b), whereas binary evolution can produce such stars (Langer et al. 2008;Marchant 2016).Their magnetically braked models not only slows down the surface rotation but also allows for strong differential rotation to develop inside the star.Furthermore, the Pitts-Tayler instability aided by the Ω effect develops, which allows for efficient chemical diffusion to account for the surface nitrogen enhancement.The assumptions in Potter's work are comparable to those in the Tayler-Spruit dynamo, thus drawing a common picture with other evolutionary models (e.g., Meynet et al. 2011).In contrast, our model predicts a nearly rigid rotation inside a magnetic star and there is no efficient matter-mixing resulting from any instabilities powered by differential rotation.Considering the large impact of chemical mixing on stellar evolution, we will investigate the relationship between nitrogen enhancement, rotation, and the magnetic field in a future work.
Another effect that is omitted from our present study is the magnetic suppression of hydrodynamic flows.Its relevance in accounting for the stability of atmospheres of Ap/Bp stars, which has been required for atmospheric diffusion processes to take place, was discussed by Michaud (1970).Similarly, a small macroturbulence of only ∼ a few km s −1 was measured in the O-type star NGC 1624-2, which has a strong dipolar surface field of ∼20 kG.Usually, the macroturbulent velocities, which are thought to be caused by pressure waves emitted by sub-photospheric convection zones (Grassitelli et al. 2015), are at least one order of magnitude larger in such stars (Simón-Díaz et al. 2017), however, Sundqvist et al. (2013) show quantitatively that the magnetic pressure in this star is strong enough to suppress the sub-surface convection.The impact of the magnetic inhibition of the core convection on blue supergiant evolution has been investigated by Petermann et al. (2015), where convectively unstable regions in a star are artificially reduced by modifying the convective criterion.They have indeed shown that this modification significantly affects the stellar structure to reproduce the enigmatic surface temperature of the progenitor of supernova 1987A, which otherwise will require a stellar merger during the evolution (Menon & Heger 2017;Urushibata et al. 2018).
In the convection criterion proposed by Lydon & Sofia (1995), a magnetic field can both stabilize and destabilize A19, page 18 of 28 the region, depending on the radial gradient of the field strength.In a 3D radiation magneto-hydrodynamic simulation by Tremblay et al. (2015), which has a cuboid-shaped computational domain embedded in an atmosphere of a white dwarf ("box-in-a-star"), it has been shown that convective transport is significantly impeded when the plasma-β (the ratio between the thermal and the magnetic pressure) is less than the unity.However, it is unclear whether such a strong field can remain within a convectively unstable layer because theoretical works strongly indicate that there is no stable magnetic structure in a marginally convective unstable barotropic region (e.g., Reisenegger 2009;Mitchell et al. 2015).A subsequent investigation of convective instability under a strong magnetic field will be required to understand.

ApBp stars
Ap/Bp stars are main-sequence A-and B-type stars that show enhancements in surface chemical abundances of elements such as Si, Cr, Fe, and Eu.Observationally, there is a strong coincidence of the peculiar surface abundance pattern and the strong surface magnetic field (Babcock 1958).For example, a surface magnetic field stronger than 100 G has been detected for 41 out of 97 Ap/Bp stars, while no magnetic field was found in 138 normal AB-type stars (Bagnulo et al. 2006).The contemporary understanding of this finding is that the strong surface magnetic field stabilizes the stellar subsurface layers such that peculiar chemical abundance at the surface can result from long-term gravitational settling and radiative levitation (Michaud 1970).
Possibly, there are two different ways to improve the stability of stellar subsurface layers.The first possibility is that strong enough magnetic fields will regulate the fluid flow such that the flow that erases the subtle chemical imprints in Ap/Bp stars is prevented.In particular, the quantitative assessment in Sundqvist et al. (2013) indicates that a strong magnetic field can indeed suppress sub-photospheric convection.In our 1.5 M models, this subsurface convective turbulence is estimated to have the energy density of 1 2 ρv 2 cv ∼ 6.3 × 10 3 −1.6 × 10 5 erg cm −3 , which yields an equipartition field strength of B eq ∼ 300−1400 G.Such strong fields are maintained during the greater part of the main-sequence evolution in our magnetic models with B p,ini ≥ 1 kG and P rot,ini ≥ 10 d (Fig. 8).Similarly, meridional flows that would exist in a rotating star have been postulated to disrupt the chemical inhomogeneity.By applying an estimate of Kippenhahn (1974), our nonmagnetic model with rapid rotation of P rot,ini = 1 d is estimated to have v ES ∼ 1 cm s −1 close to the surface.It is ∼10 −4 times slower than the turbulence velocity of the subsurface convection and, hence, it could be significantly affected by very weak fields with strengths of 10 −5 G.
The second possibility is more indirect: because magnetic stars are also known to be slowly rotating stars, their subsurface layers will be less affected by rotation induced flows than nonmagnetic stars.For instance, the flow velocity of the Eddington-Sweet circulation is assumed to be proportional to the square of the rotation frequency.Michaud (1970) has discussed that the subsurface layers of Ap/Bp type stars need to lack flows faster than ∼10 −3 cm s −1 .Considering that a model rotating with a period of P rot ∼ 1 d forms meridional flow with v ES ∼ 1 cm s −1 , long rotation periods with P rot 30 d would be required to achieve such slow flow velocities.With the help of magnetic braking, this condition is again satisfied in our models with B p,ini ≥ 1 kG and P rot,ini ≥ 10 d (Fig. 7).Besides this, turbu-lence powered by instabilities due to differential rotation, such as dynamical and secular shear instabilities and the Pitts-Tayler instability, may disrupt the chemical inhomogeneity.Our magnetic model is also compatible with that since our stellar models keep nearly rigidly rotating during the whole main-sequence phase (Fig. 5).
Correlations among the stellar age (τ), mass (M), rotation period (P rot ), and surface magnetic field strength (B z ) of Ap/Bp stars have been studied by several authors (Mathys et al. 1997;Hubrig et al. 2000;Bagnulo et al. 2006;Kochukhov & Bagnulo 2006;Landstreet et al. 2007Landstreet et al. , 2008;;Mathys 2017;Netopil et al. 2017).This showed that magnetic Ap/Bp stars are, at the same time, slow rotators.The peak rotation velocity of ∼40 km s −1 (Netopil et al. 2017) is considerably slower than the major peak at ∼200 km s −1 of the wide distribution of rotation velocities of normal AB-type stars.Some of the Ap/Bp stars even show superlong rotational periods of P rot > 1000 d (Kochukhov & Bagnulo 2006;Netopil et al. 2017;Mathys et al. 2019Mathys et al. , 2020)).Interestingly, a similar tendency is also observed for pre-MS stars known as Herbig Ae/Be stars: magnetic Herbig Ae/Be stars are concentrated to have slow rotation velocities of 50 km s −1 , while normal Herbig Ae/Be stars obey a wide distribution with a typical velocity of 50−250 km s −1 (Alecian et al. 2013).
Our results on the evolution of the surface rotation period with magnetic braking and magnetic dissipation (Fig. 9) is qualitatively consistent with these observations.For strong-enough initial surface fields of 100 G, the magnetic braking efficiently reduces the stellar angular momentum to increase the rotation period by a factor of 10 or even more, whereas the rotation period stays nearly constant within a factor of ∼2 difference for the whole main-sequence phase for models with weaker fields.Moreover, models with the stronger initial field of 1000 G can account for the super-slow rotators with P rot > 1000 d if they have initial rotation periods of 10 d.Furthermore, since magnetic dissipation becomes more efficient with faster rotation, surface fields of models with fast initial rotation velocities of ∼100 km s −1 with P rot,ini = 1 d rapidly decrease below ∼100 G, thus contributing to the lack of fast-rotating magnetic stars.
Observational tests of the evolution of the surface magnetic field are more complicated.An absence of young magnetic stars has been reported by Hubrig et al. (2000).However, later studies have not confirmed this conclusion but some did identify some young magnetic stars in their sample (Landstreet et al. 2007).The authors have attributed the inconsistency to the different ways of bolometric and effective temperature corrections, which implies a significant uncertainty in the age determination by fitting the position on the HR diagram.Moreover, Landstreet et al. (2007Landstreet et al. ( , 2008) ) have collected their stellar sample from open star clusters, which allows an age determination by isochrone fitting with a considerably improved accuracy.Landstreet et al. (2008) suggested that a strongly magnetized star with rootmean-square (rms) fields larger than 1 kG only appears close to the ZAMS in the HR diagram.This is in line with the finding by Fossati et al. (2016) for O-type stars that the fraction of magnetic stars decreases for larger evolutionary age.
Our current simulations are not intended to explain any of the observations.Also, our 1.5 M models are not completely consistent with the mass range of the Ap stars of 2−3 M observed by Landstreet et al. (2007Landstreet et al. ( , 2008)).Nevertheless, as a demonstration of the capability of our magnetic stellar evolution code, our results are compared with observational results in Fig. 11.We note that here we multiply a factor of 3.3 to the observed field strength of B rms , which is the median rms of the so-called mean longitudinal field, B l (Landstreet 1988(Landstreet , 1992)), in order to Fractional Age B p,ini =10 kG 1 kG 100 G Fig. 11.Comparison of the surface magnetic field evolution of our 1.5 M models and observational results for 2−3 M Ap stars analysed by Landstreet et al. (2007Landstreet et al. ( , 2008)).Our models with B p,ini = 10 kG, 1 kG, and 100 G are shown by red, blue, and green lines, respectively (see legend).Similar as in Fig. 8, solid, dashed, dotted, and dash-dotted lines indicate initial periods of P rot,ini of 1, 10, 100, and 1000 d, respectively.Dots with error bars correspond to observed stars, where black filled symbols show stars for which a field is detected, while gray open symbols are probable magnetic stars.Numbers associated with several stars denote V sin i in units of km s −1 .The threshold magnetic field of B th ∼ 300 G, below which essentially no Ap/Bp stars are found (Aurière et al. 2007), is shown by the black dashed line.A numerical factor of 3.3 is multiplied to the field strengths reported in the literature to convert B rms into B p .compare it with our simulation results of the polar field strength B p .We refer to Aurière et al. (2007) for the factor of 3.3 2 .
The majority of Ap stars appear to have a polar field strength of ∼1 kG during the entirety of the main-sequence phase.This may be compatible with our models assuming an initial field strength of B p,ini ∼ 1 kG.On the other hand, although it is not (yet) statistically significant, the dearth of strongly magnetized (B p 10 kG), evolved (fractional age 0.5) Ap stars is not reproduced by our models.This implies that the magnetic flux conservation is insufficient with regard to explaining the observation.However, since efficient magnetic braking stops rotations for those strongly magnetized models, no efficient rotational turbulence is expected from the current modeling (see models with 1 kG in Fig. 9).Furthermore, several Ap stars, including evolved ones (fractional age 0.5), still show certain surface rotation velocities.This is inconsistent with our current models, since our models with B p,ini 1 kG (except for the model with the fastest initial rotation) essentially stop rotation during the early main-sequence phase (see also Fig. 9 and Table 1).This might imply an inaccuracy of our treatment of the wind mass loss or the magnetic braking.
Due to its intrinsic weakness, the wind mass-loss rates of AB-type stars are highly uncertain (see discussion in Krtička 2014).Using the finite rotational periods of magnetic stars, it may be possible to constrain the unknown mass-loss rates in the less massive stars.For example, by assuming the moment of inertia and the mass loss rate are constant during the mainsequence phase, it is expected that the surface angular velocity of a magnetic star exponentially decreases with time as Ω(t) = Ω init e −t/τ break , 2 This simple conversion will still provide a reasonable comparison with the theoretical results, although it is not B rms but, rather, B max l that is used in the original estimate, and this conversion will actually yield a lower limit of the polar strength.
where Ω init is the initial angular velocity and τ break = J/ J is the braking timescale.The angular momentum loss rate linearly correlates with the wind mass-loss rate of a nonmagnetic star and the magnetic braking efficiency, which also linearly correlates with the surface field strength in a strong field limit (see Appendix C).Under this simplified case, ∼97% reduction of the currently applied value of ∼10 −12 M yr −1 , which is estimated according to de Jager et al. (1988), will be required to stay within the rotation period of HD 50169 of 10 600 d, the slowest rotating Ap star with accurate magnetic field determinations so far discovered, for a magnetic stellar model with P rot,ini = 100 d and the field strength of HD 50169 of B p ∼ 4300 G (Mathys et al. 2019(Mathys et al. , 2020)).Where a magnetic star of B p ∼ 100 G is assumed to acquire 40 times slower rotation period at TAMS phase (Table 1) with the mass loss rate from de Jager et al. (1988).We note that such a reduction of the mass-loss rate would not significantly change our present prediction for weakly magnetic models (B p,ini 10 G; Fig. 9).
Another important issue is to understand the threshold of B thr ∼ 300 G below which no Ap or Bp stars with definite field detection have been found.According to Aurière et al. (2007), the threshold does not result from the observational bias of limiting the sample to Ap or Bp stars that satisfy the condition for stabilizing the surface layers.In this case, non-Ap stars with large scale fields below 300 G should exist, which are, however, not observed (e.g., Bagnulo et al. 2006).Aurière et al. (2007) discuss the notion that differential rotation can form in the subsurface region of weakly magnetized stars and, accordingly, the field can decay through the Pitts-Tayler instability.However, our model is incompatible with this hypothesis since a magnetic model develops essentially no differential rotation even for a case with a weak initial field strength of B p,ini = 10 G (Fig. 5).Jermyn & Cantiello (2020) investigate the relation between the surface magnetic field and the subsurface convection and propose a possible explanation of a bimodal distribution of the field strength, which would be observed in O-type stars (Grunhut et al. 2017).They consider the idea that strong fields suppress subsurface convection; otherwise, the emergence of the subsurface convection erases the magnetic field.Their critical field strength to suppress the subsurface convection might correspond to the 300 G threshold.As discussed earlier, our simple estimate also shows that magnetic fields with 300 G may be sufficiently strong to affect the subsurface convection.
On the other hand, the origin of this threshold could relate to the formation of magnetic stars.For example, magnetic stars might be predominantly formed via stellar mergers (Schneider et al. 2019) that yield magnetic fields stronger than this threshold; otherwise, star formation could form only non-magnetic stars.Indeed, for O-type stars, several magnetic stars below the 300 G threshold have been found (Fossati et al. 2015), despite the detection limit of 600 G for this type of star (Fossati et al. 2016).The threshold may be present and may have lower values for massive stars (see also Jermyn & Cantiello 2020).However, it is also possible that a lower threshold is the result of magnetic field decay during stellar evolution, which could have different efficiencies depending on the stellar mass, and does not relate to the stability of the stellar fields (Fossati et al. 2016).If this is the case, it is not inconsistent with our models when suitable birth probability distributions for magnetic field and rotation were adopted.

Rotation period changes in Bp stars
Some magnetic stars show the variability of their rotation frequency on timescales of ∼10−100 yr, which is found A19, page 20 of 28 spectroscopically, photometrically, or through magnetic field measurements (Shultz et al. 2018;Pyper & Adelman 2020, and references therein).The rotation period of the magnetic star can be determined by analyzing the variations (e.g., Krticka et al. 2009Krticka et al. , 2015)), and in some cases, even a time evolution of the rotational period can be measured.For instance, comparing the long base-line photometric data, an increasing rotation period, and thus a decreasing rotational rate, has been detected for a magnetic B type star, σ Ori E, which has been explained by magnetic braking (Townsend et al. 2010).
For other B type stars of CU Vir, HD 37776, HD 142990 (Pyper et al. 1998;Mikulášek et al. 2008;Mikulášek 2016;Shultz et al. 2019), and, more recently, 13 And and V913 Sco (Pyper & Adelman 2020) -decreasing rotational periods have even been found, indicating that these stars accelerate their surface rotation.Obviously, this is inconsistent with the prediction of the magnetic braking theory.Considering the high occurrence rate of spin-up, it is possible that the period derivative in σ Ori E may also change its sign in the future (Shultz et al. 2019).
Although the mechanism for the observed spin-up has not yet been clarified, an interesting explanation has been proposed by Krtička et al. (2017).They considered a star in magnetohydrostatic equilibrium and derived an incompressible wave equation describing the propagation of a torsional magnetohydrodynamic wave.By assuming axial symmetry and applying a simple poloidal field structure, they found that the periodic cycle of 67.6(5) yr estimated for CU Vir can be reproduced well by the basic resonant frequencies of 51 yr in their model.Hence, they have shown that the torsional surface oscillation, which results from the propagation of the torsional Alfvén wave, can provide a possible explanation for the observed period decrease (see also Ste ¸pień 1998).For other possible mechanisms, see discussions in Shultz et al. (2019) and references therein.
Assuming similar physics parameters as done by Krticka et al. leads us to quite similar magneto-hydrodynamical waves in our model.As shown in Sect.3.2, our model starts to oscillate torsionally when an initial perturbation is added to the background equilibrium state (Fig. 4).In both models -Krticka's and our own -the propagation of the torsional Alfvén wave accounts for the oscillation.The magnetic strength distribution assumed in Krticka et al. is simpler than ours: they considered a field configuration of (B R , B φ , B z ) = (B, 0, −zB/R) with constant B in cylindrical coordinates R, φ, and z.On the other hand, the poloidal component in our model has a general radial dependency.In addition, we consider the long term evolution of the magnetic field.In fact, our model predicts that the field close to the surface can evolve to have a different field strength than the internal region.Although the models presented by Krtička et al. (2017) cannot reproduce the period change of HD 37776, it might be possible to explain the rapid rotational period by considering more realistic radial distributions of the stellar magnetic field.If this hypothesis is correct, then the period of the variation of the rotation rate will roughly correlate with the strength of the surface magnetic field.
This oscillation period, which has a similar timescale to the wave-crossing time of the torsional Alfvén wave, depends on the strength of the internal poloidal field component.Therefore, it will be fundamentally possible to determine the internal field strength distribution by observing the high-order rotational period change at the surface of the star.The idea of determining the internal magnetic field by observing the surface oscillations is reminiscent of asteroseismology.However, instead of gravity or pressure waves, here, it is the propagation of Alfvén waves that produces the observable signal.
Estimates for the presence of an internal magnetic field are provided indirectly, namely, by considering that a strong enough magnetic field can reduce the amplitude of non-radial dipolar (l = 1) mode oscillations (Fuller et al. 2015;Stello et al. 2016;Cantiello et al. 2016).However, so far this method has only provided the information of fields close to red-giant cores that result from long and complicated evolutionary histories.Perhaps, torsional oscillation, if observed in reality, can reveal the field distribution in the radiative envelope of a main-sequence star, which can be directly compared with corresponding stellar evolution models.
Applying asteroseismology to a magnetic star still requires further developments of the theory (Kiefer et al. 2017;Loi & Papaloizou 2017, 2018, 2020;Kiefer & Roth 2018).Also, 2D modeling of the background stellar structure will be important in considering a realistic poloidal field configuration (Rincon & Rieutord 2003;Reese et al. 2004;Prat et al. 2019).Nonetheless, MHD-oscillations may play a crucial role in detecting and analyzing internal stellar magnetic fields in the future.

Core-envelope decoupling in red-giant star
Without angular momentum transport, the spin periods of the cores of red-giant stars would be expected to be very short, ∼10 −2 d.Longer periods would result if an effective angular momentum transport is assumed, but the naive expectation has been that red-giant core rotation periods would be much shorter than those of the surfaces.However, over the last decade, asteroseismology has revealed rotation periods of red-giant cores of ∼10 d (Beck et al. 2012;Mosser et al. 2012;Deheuvels et al. 2014).This implies that even the most efficient angular momentum transport mechanism proposed at that time (the Tayler-Spruit dynamo) left a discrepancy of predicted and observed core rotation periods of more than one order of magnitude (Cantiello et al. 2014;Spada et al. 2016;Eggenberger et al. 2017).Recently, Fuller et al. (2019) revised the diffusion coefficient of the TS dynamo, suggesting that the saturation of the dynamo cycle takes place much later than assumed in Spruit (2002), such that a stronger toroidal component and magnetic torque are obtained.The red-giant core periods have been well reproduced in their simulation, however, the basic picture of the TS dynamo theory still remains uncertain (Sect.6.1).Although other mechanisms such as angular momentum transfer by internal gravity waves (Fuller et al. 2014;Pinçon et al. 2017) have been proposed, these are generally not efficient enough to overcome the problem.
The rotation periods at the center and the surface of our magnetic models are compared with observations in Fig. 12. Central spin periods of a model with n cv = 0 is shown by the red solid line, and that with n cv = 2 is by the red dotted line.The surface period of the n cv = 0 model is shown by the black dashed line, while that of the n cv = 2 model is omitted since they are nearly identical.In addition to the magnetic models, central spin periods of a model with no internal angular momentum transfer (labeled as "j cons.",green) and a model with the Tayler-Spruit dynamo ("TS dynamo", blue) are shown as well.
The central period of our fiducial magnetic model with n cv = 0 coincides with its surface period up to the point of log R/R = 1.5.In other words, the model sustains near-rigid rotation even after the formation of the red-giant envelope.While this model is incompatible with the observations, we stress the importance of this result: a stellar evolution calculation that self-consistently accounts for the interaction between differential rotation and magnetic field obtains the enigmatically large A19, page 21 of 28 efficiency of the angular momentum transfer needed to understand the red-giant core periods.Interestingly, in the magnetic model with n cv = 2, the core develops a shorter rotation period than the surface after the surface rotation period increases to ∼20 d and keeps its period roughly constant up to the point when log R/R = 1.4.Therefore, not only the core rotation period but also the surface period agrees with the observations.This is because this model develops an angular momentum distribution as Ω ∝ r −1 in the convective envelope, while Ω is constant in the inner radiative region (Fig. 13).The distribution, Ω ∝ r −1 , results from the approximation, ν DS = ν MLT (Sect.1).On the other hand, the rigid rotation in the inner radiative region is solely due to the magnetic effect.
Our n cv = 2 result predicts that strong shear exists at the base of the convective layer, at ∼1 R , in a red-giant star.This is in contrast to the result of Fuller et al. (2019) (Fig. 5), in which the shear rotation mainly forms at the hydrogen-burning shell, at ∼3 × 10 −2 R , much deeper within than ours.To reveal the differential rotation inside the convective envelope in a redgiant star by asteroseismology, the detection of quadrupole and octopole pulsation modes may be required (Ahlborn et al. 2020).We expect that future asteroseismic observations will discriminate among these two cases.
We have also noticed an interesting correspondence between the very strong toroidal field of ∼10 6 G obtained in our simulation at the core-envelope boundary (B φ at M r ∼ 0.4 M , Fig. 10) and possible observational constraint for the internal magnetic field for red-giant stars of 10 6.5 G (Fuller et al. 2015;Stello et al. 2016;Cantiello et al. 2016).More stringent simulations and discussions for the internal magnetic field of red-giant stars will be done in the future.
We note that more than two orders of magnitude acceleration takes place in the cores of magnetic models at log R/R = 1.4, the mechanism of which has been discussed in Sect.5, whereas it does not take place in the model with the Tayler-Spruit dynamo.This is because the core has already a larger specific angular momentum than the accreted material and, additionally, the timescale of the angular momentum transfer inside the core at this point is already too long to affect the rotation period at the center of the core.

Rotation periods of white dwarfs
Finally, we discuss the rotation velocities of white dwarfs that are expected from our simulations.The ∼2 orders of magnitude acceleration of the core shown in Fig. 12 is due to the matter accretion of the helium core, as explained in Sect. 5. Assuming that the angular momentum of a white dwarf chiefly originates from the accreting matter, the angular velocity of a white dwarf will be estimated as where J acc and I WD are the angular momentum accreted onto the central core after rotational decoupling and the moment of inertia of the white dwarf, respectively.Here, I WD = 0.205 × R 2 WD M WD and the radius of R WD = 1.20 × 10 −2 R and the mass of M WD = 0.575 M are taken from a 1.5 M model of Suijs et al. (2008).Furthermore, we have J acc ∼ (2R 2 base Ω base /3)∆M acc , where R base ∼ 1 R and Ω base are the radius and the rotation rate at the base of the convective envelope, and ∆M acc is the mass that is accreted onto the core after the decoupling.Therefore, A19, page 22 of 28 is obtained.Our fiducial model has Ω base ∼ 10 −8 rad s −1 and ∆M acc = 0.28 M , hence a white dwarf angular velocity of Ω WD ∼ 1.3 × 10 −4 rad s −1 , and correspondingly a rotation velocity of the white dwarf of v WD ∼ 1.1 km s −1 is expected.The n cv = 2 model has a rotation rate at the base of the envelope which is one order of magnitude faster, Ω base ∼ 10 −7 rad s −1 , and similar ∆M acc = 0.25 M , and thus Ω WD ∼ 1.2 × 10 −3 rad s −1 and v WD ∼ 9.9 km s −1 are expected.Both estimates are compatible with the spectroscopic upper limit of ∼10 km s −1 (Berger et al. 2005).
For the present models, the α effect has not yet been taken into account.A red-giant model with the α effect will show stronger magnetic fields than the present models, and the magnetic braking will reduce the angular momentum of accreted material and thus reduce the rotation rate of the white dwarf.This reduction will improve the comparison to observations, since rotation periods of both groups, non-magnetic white dwarfs with P ∼ 1−169 h corresponding to Ω ∼ 1.7 × 10 −3 −1.03 × 10 −5 rad s −1 , which are measured by asteroseismic observations, and magnetic white dwarfs of P ∼ 0.2−429 h corresponding to Ω ∼ 8.7 × 10 −3 −4.1 × 10 −6 rad s −1 , which are measured by spectroscopic modulation, include slow rotators (Kawaler 2015;Córsico et al. 2019, and references therein).
Another relevant question here is whether a large-scale poloidal field linking the whole radiative mantle region, which covers the central core and is covered by the convective envelope, can form in a red-giant star.Because the magnetic fields induced by the α effect will show a time variability or would probably be dominated by small scale fields, the magnetic fields could quickly cancel each other out after the material migrates from the convective region into the radiative region.If the magnetic link is lost, the rotational decoupling happens no matter how strong the field is formed in the convective envelope.The rotation velocity of the white dwarf depends on when the rotational decoupling takes place because the accreted mass on the helium core after the decoupling determines the angular momentum of the white dwarf.

Conclusion
In this work, we develop a new formalism for stellar evolution calculations, which includes the interaction of stellar rotation and stellar magnetic field as self-consistently as possible, in a 1D model.With this method, we computed evolutionary models of 1.5 M stars, adopting various initial magnetic field strengths and rotation rates.We compared the theoretical models with relevant observations of (1) ages, rotation rates, and magnetic field strengths of Ap stars; (2) surface rotation variations observed in Bp stars; (3) core and surface rotation periods of red-giant stars; and (4) rotation periods of white dwarfs.Even though we have not manipulated the model to explain any of these observations, we found a generally good agreement between our modeling and observations.This work demonstrates the first results in a series in which we intend to develop and apply a new scheme to simulate the evolution of magneto-rotating stars.In forthcoming papers, additional magnetic effects on the stellar evolution, such as the α effect, magnetic pressure, and magnetic modification of convective criterion will be included and analyzed.There are plenty of possible applications of the new scheme: the evolution of solartype stars, for which abundant works have been done; the evolution of massive stars is an interesting target as well because strong surface fields will significantly affect the evolution as they have a relatively fast rotation and strong wind mass loss and, of course, because they are the progenitors of (non-)magnetized compact remnants; the evolution of binary systems, such as a merger remnants (Beloborodov 2014;Schneider et al. 2019), as well as tidally interacting binaries (Vidal et al. 2018(Vidal et al. , 2019))  The top two panels showing radial and toroidal magnetic field components exhibit the effect of core contraction and envelope expansion: the magnetic field in the central core of 0.4 M is amplified by about two orders of magnitude, while that in the outer envelope is reduced by about four orders of magnitude.Nevertheless, the two conserved quantities, which are shown in the bottom panels, are entirely conserved during the whole evolutionary phases.Only a small fluctuation is seen for B φ /ρr at ∼0.3−0.4M .This results from automated mesh refinement, which is done to capture a thin structure of the hydrogen-burning shell that surrounds the helium core.

D.2. Magnetic dissipation
Here, the effect of magnetic dissipation is tested.Firstly, magnetic dissipation due to core convection is calculated taking the initial condition from the above test calculation.It is in the main-sequence phase with the central hydrogen mass fraction of 0.3.Both field components initially have a radial dependence of ∼r −3 .The star at this phase has a convective hydrogen-burning core of ∼0.2 M .The convective turbulence explains the magnetic diffusivity of ∼10 12 cm 2 s −1 , which is large enough to establish a steady state for the magnetic field within a much shorter timescale than the evolutionary time.Figure D.2 shows the result of the evolution of the radial (top) and toroidal (middle and bottom) magnetic field components.As a result of the magnetic diffusion, the radial component reaches a steady state inside the convective region within a timescale of ∼10 9 s, where the time derivative of the poloidal field becomes uniformly nearly zero.This steady state is achieved as the magnetic diffusive flux becomes uniform in the convective region: with B r ∝ r 0 and A ∝ r.The convective region is surrounded by an overshoot region, in which magnetic diffusivity exponentially decreases with radius.The small magnetic diffusivity limits the magnetic diffusive flux, explaining the longer timescale of ∼10 13 s for the further extension of the steady region.Similarly to the poloidal component, the toroidal component also reaches a steady state at ∼10 9 s, and the steady region extends further with a longer timescale of ∼10 13 s.The middle panel shows the calculation result in which the magnetic diffusion term in Eq. ( 37) is considered but the magnetic advection term is not; because of the B φ ∝ r 2 distribution, the diffusive flux of the toroidal field becomes uniform as well.
In the bottom panel, the toroidal field evolution of a calculation in which both the magnetic diffusion and the magnetic advection are taken into account is shown.We might suspect that the toroidal component in this case does not reach the steady state because B φ at late times does not converge.However, the rate of the change of B φ is significantly smaller than the fluxes of the diffusion and the advection.In fact, both fluxes cancel each other out to achieve ∂B/∂t ∼ 0. In this meaning, again, we consider that the toroidal field evolves keeping a steady state.In this steady state, the radial gradient of the toroidal magnetic field can The dissipative region can be divided into the fully convective region of 0.17 M and the surrounding overshooting region at ∼0.17−0.21M .In the former region, n is so small that m ∼ 2 is achieved.Meanwhile, n is as large as ∼200 in the latter region, resulting in m ∼ 200 to compensate for the advection flux by the diffusion flux.
Secondly, magnetic dissipation due to the envelope convection is tested in a similar way taking an initial condition with a red-giant envelope of a radius of R = 10 R .The uniform initial distributions of B r = B φ = 1 kG are applied for the magnetic components in this case.The results are shown in Fig. D.3.Within a diffusion timescale of ∼10 9 s, the magnetic field finds a steady state of B r ∝ B φ ∝ r −3 , which corresponds to having a zero magnetic diffusive flux.The advection term in this case has only a minor effect, probably due to the much narrower convective overshooting region located at the base of the convective envelope.
In conclusion, the magnetic field is destined to reach a steady state under the efficient magnetic dissipation effects of convection.In a convective core, the magnetic field distributes such that the magnetic flux uniformly distributes.In contrast, the distribution of the magnetic field is such that the magnetic flux becomes zero in a convective envelope.A19, page 28 of 28

Fig. 1 .
Fig. 1.Schematic illustration of the geometry used to define fluxes Φ B and Φ j .S , which is the surface of a polar cap with a radius r and the opening angle θ, is shown by the red shaded area.Fluxes are defined as Φ B ≡ S B • ndS and Φ j ≡ S j • ndS .

Fig. 2 .
Fig.2.Angular velocity, Ω (red lines), and strength of the φ-component of the magnetic field, B φ (blue lines), for ten different times during the evolution of a 1.5 M main-sequence stellar model, illustrating the propagation of torsional Alfvén waves in the stellar interior.Here, the Ω effect and the Maxwell stress are taken into account, but the η effect is not.The x-axis for the left panels, showing times up to 8 × 10 8 s, is the Lagrangian mass coordinate, the right panels, depicting later times, use the radius coordinate.On the y-axis, The logarithm of Ω divided by Ω min (or log B φ,min /B φ,min ) is plotted with the normalization values of Ω min = 1 × 10 −6 rad s −1 and B φ,min = 1 × 10 3 G.Prograde rotation and positive Bfield polarity are plotted in the upper half of the panels, retrograde rotation and negative polarity are plotted in the lower half.Small structure with −Ω min < Ω < Ω min and −B φ,min < B φ < B φ,min is omitted from this plot.The black dashed and black dotted lines show the positions of the wavefronts, r f , which are estimated by Eq. (45).Arrows indicate the direction of the wave propagation.

Fig. 3 .
Fig.3.Angular velocity Ω (red), toroidal magnetic component B φ (blue), and magnetic viscosity η (green) as a function of the Lagrangian mass coordinate at eight different moments in time in test calculations of a 1.5 M main-sequence stellar model, using different assumptions for magnetic and viscous dissipation.Left panel: results for no dissipation, results with dissipation only due to hydrodynamic instabilities are shown in the middle, and results with hydrodynamic dissipation and Pitts-Tayler instabilities are shown in the right panel.For this stellar model the wave-crossing time is t wc = 1.71 × 10 9 s.The top seven plots of each panel cover about two wave-crossing times (3.4 × 10 9 s) with five snapshots per plot, one every 10 8 s, for the period given in the top right corner of each plot.In the bottom row of plots, results after about six wave-crossing times are shown, i.e., from 9.9 × 10 9 s to 10.3 × 10 9 s.The temporal evolution is available as an online movie.

Fig. 5 .
Fig. 5. Left panel: main-sequence evolution of 1.5 M models in the HR diagram, computed with n cv = 0 and P rot,ini = 1 d.The result of magnetic models with an initial magnetic field strength of B p,ini = 10 G and 10 kG are shown by the red solid and blue dashed lines, while the result of the nonmagnetic model is shown by the green dash-dotted line.Epochs for which the central hydrogen become X c = 0.1 (and 0.7 for the B p,ini = 10 G model) are indicated by dots.Right panels: profiles of internal angular velocity (top left), radial and toroidal magnetic field strength (top right), hydrogen mass fraction (bottom left), and the ratio between the magnetic pressure and the pressure (bottom right) at X c = 0.1 as a function of Lagrangian mass coordinate.Again, profiles of the B p,ini = 10 G and 10 kG models are shown by the red solid and the blue dashed lines, while that of the nonmagnetic model are by the green dash-dotted lines.Top right panel: radial and toroidal field strength are shown by the thick and thin lines, respectively, and the radial component for the B p,ini = 10 kG model is multiplied with a factor of 10 −3 .The surface values are highlighted by dots in the right-top and right-bottom panels.

Fig. 6 .
Fig.6.Reconstructed 2D magnetic field structures at core hydrogen mass fractions of X c = 0.7 (left) and 0.1 (right), for our 1.5 M model with an initial magnetic field strength of B p,ini = 10 G and an initial rotation period of P rot,ini = 1 d.The x-axis goes through the equatorial plane, while the z-axis corresponds to the stellar rotation axis.Thick purple lines show the stellar surface, while thin purple lines designate the convective core boundary.Green lines indicate magnetic field lines, the interval of which is taken such that 1 line in 1 R corresponds to the field strength of 0.2 G.The field lines outside of the star are constructed assuming a dipole structure.The color indicates the strength of the toroidal field.Note that deformation due to centrifugal forces and structure change of the circumstellar magnetic field due to wind interaction are not taken into account here.

Fig. 7 .
Fig. 7. Evolution of the surface rotation periods as a function of time during core hydrogen burning for our 1.5 M models.Models with B p,ini = 10, 100, 1000, and 10 000 G are shown by red solid, dashed, dotted, and dash-dotted lines, respectively.Models with P rot,ini = 1000, 100, 10, and 1 d are shown in the top-left, top-right, bottom-left, and bottom-right panels, respectively.The solid black lines correspond to models without stellar wind mass loss and with B p,ini = 10 G for each initial rotation period.

Fig. 8 .Fig. 9 .
Fig. 8. Evolution of the polar surface magnetic field strength as a function of time for the same models as those shown in Fig. 7. Models with P rot,ini = 1, 10, 100, and 1000 d correspond to red solid, dashed, dotted, and dash-dotted lines, respectively.Models with B p,ini = 10 G are shown in the top left panel, those with B p,ini = 100, 1000, and 10 000 G are, respectively, shown in the top-right, bottom-left, and bottom-right panels.The black solid lines depict models without the η effect with P rot,ini = 1000 d for each initial magnetic field strength.

Fig. 10 .
Fig. 10.Left panel: evolution of a 1.5 M magneto-rotating model in the HR diagram, computed with n cv = 0.The line colors indicate the different evolutionary phases.Grey lines are iso-radius lines for log R/R = 0.0, 0.5, 1.0, 1.5, and 2.0.Right panels: profiles of internal density (top left), radial field strength (top right), rotation frequency (middle left), local wave-crossing time (middle right), specific angular momentum (bottom left), and toroidal field strength (bottom right) as a function of the Langrangian mass coordinate.Left panel: six epochs for which profiles are shown in the right panels are indicated by colors and dots that have the same color as the profiles.

Fig. 12 .Fig. 13 .
Fig. 12. Rotation period at the surface (black dashed line) and near the center (red line) of our fiducial 1.5 M model (n cv = 0) as function of its radius.The central rotation period of a similar model with n cv = 2 is shown by the red dotted line.For comparison, the central rotation periods of models with no angular momentum transfer (green line) and with transfer due to the Tayler-Spruit dynamo (blue line) are also shown.Theoretical models are compared with the results of asteroseismic observations.Core rotation periods obtained by asteroseismic observations(Mosser et al. 2012;Gehan et al. 2018) are shown by the magenta dashed line.Surface and core periods obtained byDeheuvels et al. (2014) are shown by black and magenta pluses.

Fig. D. 1 .
Fig. D.1.Internal magnetic field distributions of 1.5 M magnetorotational star with 5 evolutionary stages, at ZAMS (red), at the middle of the main-sequence phase (MS; green), at TAMS (blue), at the middle of the red-giant phase (RG; magenta), and at the red clump phase (Clump; cyan) are shown.In this test case, ηand Ω-effects are switched off to confirm magnetic flux conservation, and the Maxwell stress is also neglected.Top two panels: evolution of the radial component (B r ; left) and the toroidal component (B φ ; right).Bottom two panels: evolution of corresponding conserved quantities, B r r 2 by the left, and B φ /ρr by the right.
Fig. D.2.Evolution of magnetic field distributions during the core hydrogen-burning phase for the inner 0.34 M region of a 1.5 M magneto-rotational star, for five different epochs, which are indicated by the legends.In this test case, magnetic dissipation by the η effect is included and accounts for the short timescale evolution.Top panel: radial component of the field, B r .Middle panel: the toroidal component, B φ , is shown for a case in which only the diffusion term of the magnetic dissipation is included.Meanwhile, B φ evolution for a case with both diffusion and advection terms of the magnetic dissipation is shown in the bottom panel.
Fig. D.3.Same as Fig. D.2 but showing magnetic dissipation in a convective envelope during the red-giant phase of a 1.5 M magnetorotational star.
, will be interesting targets in our future explorations.The left hand side becomes zero because P is divergence-free.The first term of the right hand side is dΦ/dt.The second and third terms are equated as Fig. A.1.Illustration of the geometry.The 2d manifold M is shown by blue shaded areas.Below: M at t i .Above: M at t.An element of M, which is parametrically specified as α(ξ, η; t), is shown by a red point.It moves with velocity u(ξ, η; t).The divergence-free vector field, P, is shown as green lines.The thicker ones are P at t i and the thinner ones are P at t.
Here Gauss's theorem is again applied to the divergence-free vector field ∂P/∂t.A19, page 25 of 28