© K. Takahashi and N. Langer 2021

1. 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. 2015, 2016; Folsom et al. 2016, 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 correlations 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; GK-giants (Aurière et al. 2015) as well as asymptotic-giant-branch (AGB) and post-AGB stars (Vlemmings 2014, 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.

2. Methods

To compute the time evolution of stellar models, we use the 1D stellar evolution code HOSHI (Takahashi et al. 2016, 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 reaction–diffusion equation as

(1)

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¹ 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 ¹²C(α,γ)¹⁶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 , where v_cv is the velocity of convective eddies, l_cv = α_MLTmin(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

(2)

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 D_therm ≡ (1/C_Pρ)(4acT³/3κρ), ∇_rad ≡ (3κ/16πacG)(PL/T⁴M), ∇_ad ≡ (∂lnT/∂lnP)|_{s = const.}, ∇_μ ≡ d log μ/d log P, δ ≡ −(∂lnρ/∂lnT)_P, μ, ϕ ≡ (∂lnρ/∂lnμ)_PT, 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⁻¹² M_⊙ yr⁻¹ and ∼10⁻¹¹ − 10⁻⁸ M_⊙ yr⁻¹ during the main-sequence and red-giant phases, respectively.

2.1. 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 . 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

(3)

where P₂ is the second-degree Legendre polynomial. Here, indicates the degree of rotation compared with the local gravity, and the length scale, a, satisfies the relation (Denissenkov & VandenBerg 2003),

(4)

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)

(6)

where ∇_MLT is the convective temperature gradient determined by the mixing-length theory. The parameters f_P and f_T are calculated as

(7)

(8)

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 .

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. (2010, 2012) as

(9)

where |Ṁ(v_rot = 0)| is the mass loss rate of a nonrotating counterpart having the same luminosity and effective temperature, v_rot and 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

(10)

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.

2.2. Stellar magnetism

As the most fundamental assumption, we assume that a large-scale stable magnetic field is embedded inside a star. Here, “stable” means that the magnetic field is in a magneto-hydrostatic 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; 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 long-time 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.

2.2.1. 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:

(11)

(12)

(13)

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

(14)

To express the poloidal magnetic field, we utilize the toroidal component of the vector potential, A_tor = A_ϕe_ϕ, as

(15)

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

(16)

(17)

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

(18)

which yields

(19)

(20)

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 sin2 θ. Thus, we write

(21)

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).

2.2.2. 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 Φ ≡ ∫_SP ⋅ ndS = ∫_SP ⋅ 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

(22)

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:

(23)

where v 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

(24)

in which Stokes’ theorem is applied, and Φ_B ≡ ∫_SB ⋅ dS and U ≡ v − V.

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

(25)

Similarly, the right hand side of Eq. (24) becomes

(26)

We take the Lagrangian expansion velocity, v_m, as the velocity of the polar cap, thus V = v_me_r. The vector field, U ≡ v − 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

(27)

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

(28)

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 divergence-free assuming MHD. Using the evolution equation of the electric current, the evolution equation of the electric current flux,

(29)

is obtained. Owing to the axial symmetry, the electric current flux can be written as

(30)

(31)

The right-hand side of Eq. (29) is reduced to

(32)

in which the latitudinal dependencies of the magnetic field described above are used, and a relation , which is valid for a spherically symmetric flow, is assumed. By taking the terms with sin 2θ dependence, the evolution equation of B(r) is obtained as

(33)

The remaining terms with sin θ dependence may be used to describe the evolution of B_stb as

(34)

However, considering the function of stabilizing the poloidal component, a simpler relation,

(35)

with f_tp ≳ 𝒪(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,

(36)

(37)

It is noteworthy that both Ar and B/ρr scale as |B|r² if we assume ρ ∝ r⁻³. 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.

2.2.3. 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.

2.3. Evolution equation of the angular velocity

The equation of fluid momentum conservation can be written as

(38)

where g, Π, and are the gravity, the Reynolds stress tensor, and the Maxwell stress tensor, respectively. Here, we use a relation 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

(39)

where axial symmetry is assumed. By averaging over a sphere, we obtain

(40)

where S ≡ ∫2πr² sin θdθ is the surface area of the sphere and i ≡ ∫2πr⁴ sin³θdθ/S ∼ 2r²/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;

(41)

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 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⁻¹, which can be checked by applying ν_cv = ν_eff to Eq. (41).

2.4. Input physics

2.4.1. 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:

(42)

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,

(43)

is used for rotating models. Similarly, the effective magnetic viscosity is estimated as

(44)

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.

2.4.2. 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 , where 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.

2.5. 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.

2.6. 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).

3. 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²/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 and . The corresponding invariant becomes constant along the characteristic dr/dt = ±c. Here, 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.

3.1. 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⁻¹, and correspondingly, the estimated wave-crossing time from the center to the surface is ∼1.3 × 10⁹ s. A step-function distribution of Ω = 10⁻⁴ rad s⁻¹ for M ≤ 1 M_⊙, and Ω = −10⁻⁴ rad s⁻¹ 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⁸ s, and is then reflected to follow the down-going wavefront. The down-going front reaches the stellar center at ∼1.4 × 10⁹ s and then is also reflected. At ∼1.8 × 10⁹ 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

(45)

The wave-crossing time in the simulation is estimated to be 1.71 × 10⁹ 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⁹ 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

(46)

This is proportional to the angular velocity difference at the wavefront, ΔΩ, but is independent of the strength of the poloidal 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⁷ G at M_r = 1 M_⊙ with ΔΩ = 2 × 10⁻⁴ rad s⁻¹ 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 former 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¹⁵ erg g⁻¹ for the gravitational and thermal energies, ∼10¹² erg g⁻¹ for the rotational kinetic energy, but only ∼10⁴ erg g⁻¹ 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.

3.2. 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 magneto-hydrodynamic dissipation, the latter works effectively only during the first one or two wave-crossing times. First, shear rotation propagates inside the star with the velocity of . 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

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Ω/∂lnr, the growth time of the Pitts–Tayler instability can be estimated as . Hence, the condition above can be expressed as . Our simulation with a strong poloidal field of B_r ∼ 1 kG has τ_wc ∼ 10⁹ s. Therefore, the condition can be satisfied with a canonical value of Ω ∼ 10⁻⁵ s⁻¹, 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, , 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⁻² s⁻¹ 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.

4. 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

(47)

(48)

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 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 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 ∂Ω/∂lnr ∼ 10⁻¹⁰ s⁻¹ 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 hydrogen-burning 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, , 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-momentum-loss 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 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 magnetic 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.

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⁻³ M_⊙ during its main-sequence phase. Comparably, the models with B_p, ini = 10 000 G lose only 2.1 × 10⁻⁴ to 3.4 × 10⁻⁴ 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⁻² 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⁻³ M_⊙ for the model with B_p, ini = 100 G and P_rot, ini = 1 d and 1.8 × 10⁻³ 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.

5. 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 (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⁴ 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⁴³ g cm² s⁻¹. 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⁻² 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⁶ G is induced.

6. Comparison with previous results

6.1. 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 (1999, 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. 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_rB_ϕ/4π. The magnetic stress is further rewritten in terms of the viscosity, ν_TS, by equating ρν_TSr∂_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, 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.

6.2. 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⁻¹ 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 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.

7. Comparison with observations

7.1. 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 erg cm⁻³, 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⁻¹ close to the surface. It is ∼10⁻⁴ 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⁻⁵ 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 non-magnetic 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⁻³ cm s⁻¹. Considering that a model rotating with a period of P_rot ∼ 1 d forms meridional flow with v_ES ∼ 1 cm s⁻¹, 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, turbulence 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. 2007, 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⁻¹ (Netopil et al. 2017) is considerably slower than the major peak at ∼200 km s⁻¹ of the wide distribution of rotation velocities of normal AB-type stars. Some of the Ap/Bp stars even show super-long rotational periods of P_rot >  1000 d (Kochukhov & Bagnulo 2006; Netopil et al. 2017; Mathys et al. 2019, 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⁻¹, while normal Herbig Ae/Be stars obey a wide distribution with a typical velocity of 50−250 km s⁻¹ (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⁻¹ 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. (2007, 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 root-mean-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. (2007, 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, 1992), in order to 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².

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 main-sequence phase, it is expected that the surface angular velocity of a magnetic star exponentially decreases with time as

where Ω_init is the initial angular velocity and 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⁻¹² M_⊙ yr⁻¹, 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, 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.

7.2. 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 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. 2009, 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 magneto-hydrostatic 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 Stȩ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.

7.3. 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⁻² 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 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⁻¹ in the convective envelope, while Ω is constant in the inner radiative region (Fig. 13). The distribution, Ω ∝ r⁻¹, 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⁻² R_⊙, much deeper within than ours. To reveal the differential rotation inside the convective envelope in a red-giant 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⁶ 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.

7.4. 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, and the radius of R_WD = 1.20 × 10⁻² 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 , 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,

is obtained. Our fiducial model has Ω_base ∼ 10⁻⁸ rad s⁻¹ and ΔM_acc = 0.28 M_⊙, hence a white dwarf angular velocity of Ω_WD ∼ 1.3 × 10⁻⁴ rad s⁻¹, and correspondingly a rotation velocity of the white dwarf of v_WD ∼ 1.1 km s⁻¹ 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⁻⁷ rad s⁻¹, and similar ΔM_acc = 0.25 M_⊙, and thus Ω_WD ∼ 1.2 × 10⁻³ rad s⁻¹ and v_WD ∼ 9.9 km s⁻¹ are expected. Both estimates are compatible with the spectroscopic upper limit of ∼10 km s⁻¹ (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⁻³ − 1.03 × 10⁻⁵ rad s⁻¹, which are measured by asteroseismic observations, and magnetic white dwarfs of P ∼ 0.2−429 h corresponding to Ω ∼ 8.7 × 10⁻³ − 4.1 × 10⁻⁶ rad s⁻¹, 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.

8. 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 solar-type 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, 2019), will be interesting targets in our future explorations.

---

Acknowledgments

K.T. appreciate invaluable discussions with John D. Landstreet, Matthew E. Shultz, Alexandre David-Uraz, Jiri Krtička, Jim Fuller, Götz Gräfener, Youhei Masada, Munehito Shoda, and Tomoya Takiwaki. K.T. is grateful to Zahra Mirzaiyan and Antonio Ferriz-Mas for discussions of mathematical expression.

References

Appendix A: Alfvén’s theorem

Theorem: We have ℳ as a 2d compact manifold with a boundary 𝒞. ℳ is embedded in a 3d manifold and moves with time with a velocity field v(x, t). Consider a flux Φ of a time-dependent vector field P(x, t) on ℳ. If and only if P is divergence-free then

(A.1)

Proof: We express the position of the element of ℳ at a given time t in a parametric form by x = α(ξ, η; t), where ξ ∈ [ξ_i, ξ_f] and η ∈ [η_i, η_f]. Then the velocity field is equated as v = ∂α/∂t. The parameters (ξ, η) are concentrically defined: the center of ℳ at time t is specified as x = α(ξ_i, η; t) = x_c(t), likewise the boundary 𝒞 as x = α(ξ_f, η; t). The geometry defined here is illustrated in Fig. A.1.

In the parametric form, a flux Φ of P on ℳ is defined as a function of time as

(A.2)

(A.3)

So now we equate the total time derivative of Φ,

(A.4)

Let V be the volume swept by ℳ for t ∈ [t_i, t]. The surface of V consists of three manifolds;

• S₁ = {x ∈ α(ξ, η, t)|ξ ∈ [ξ_i, ξ_f],η ∈ [η_i, η_f],t = t_i}

• S₂ = {x ∈ α(ξ, η, t)|ξ ∈ [ξ_i, ξ_f],η ∈ [η_i, η_f],t = t}

• S₃ = {x ∈ α(ξ, η, t)|ξ = ξ_f, η ∈ [η_i, η_f],t ∈ [t_i, t]}.

S₁ and S₂ are ℳ at t_i and t, and S₃ composes the side of V.

In the parametric form, volume integral of div P on V is done as

(A.5)

Here Gauss’s theorem is used. By taking the total time derivative,

(A.6)

is obtained. 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

(A.7)

and

(A.8)

where, Stokes’s theorem is used. Notice that . Therefore, Eq. (A.6) equates with

(A.9)

Here Gauss’s theorem is again applied to the divergence-free vector field ∂P/∂t.

So we obtain

(A.10)

Appendix B: Viscosity of the Pitts–Tayler instability

The vertical and horizontal length scales of the fluctuation driven by the Pitts–Tayler instability are respectively defined as l_v and l_h. The horizontal scale will be replaced with r. The vertical sale is limited to be

(B.1)

where is the Alfvén frequency of the toroidal magnetic field and N is the Brunt–Väisälä frequency, since the vertical replacement should work against the restoring force of buoyancy. This gives the maximum scale for the vertical fluctuation. On the other hand, the fluctuation will be damped by the dissipation, if the vertical fluctuation is too small. Therefore, the dissipation time should be longer than the growth time of the Pitts–Tayler instability, τ_PT. This gives

(B.2)

The two inequalities yield a relation

(B.3)

has been used in previous works, but in such a case τ_PT becomes zero when Ω is zero. In order to avoid too rapid growth, we apply instead so that the growth time approaches 1/ω_A in the limit of Ω → 0.

Following the original discussion by Spruit (2002), we assume that the inequality (B.3) reaches the equipoise situation when the instability saturates. This is because a turbulent viscosity induced by the Pitts–Tayler instability starts to act as the effective viscosity in the right-hand side of the inequality. This results in

(B.4)

Besides, according to Maeder & Meynet (2004), we assume the relation between the Brunt–Väisälä frequency and the effective viscosity as

(B.5)

where K is the thermal diffusivity and N_T and N_μ are oscillation frequencies associated with thermal and chemical gradients, respectively.

Finally, the two equations yield a quadratic equation for ν_PT/K,

Since it has a negative zeroth-degree coefficient, it has one positive solution when is positive. We solve this equation for the estimate of the turbulent viscosity of the Pitts–Tayler instability.

Appendix C: Wind-magnetic field interaction

According to 2D axisymmetric MHD simulations with an aligned dipole magnetic field by ud-Doula & Owocki (2002, 2008, 2009), we estimate the effects of the magnetic wind confinement and the magnetic braking as follows.

First, the magnetic confinement parameter η_* is calculated as

(C.1)

where B_eq = B_θ(r = R, θ = π/2) is the surface magnetic field strength at the equator and and v_{∞, B = 0} are the wind mass-loss rate and the terminal wind velocity for a nonmagnetic model. In the present work we approximately estimate the terminal wind velocity as . Next, the Alfvén radius, R_A, where the radial components of the field and the matter flow have an equal energy density, is estimated as

(C.2)

Efficiencies of the magnetic confinement and the magnetic braking are estimated as

(C.3)

and

(C.4)

where R_c = R + 0.7(R_A − R) is a maximum closure radius of magnetic loops. We note that an additional term of , in which R_K = R(v_K/v_rot)^2/3 is the Kepler corotation radius, is included in the original formula of the magnetic confinement in Ud-Doula et al. (2008). It is discussed that this term accounts for the breakout of the gas from the closed loops due to the fast rotation of the star. However, since the enhancement happens even for a non-magnetic (η_* = 0) model, this additional term likely partly accounts for the Ω effect, which is already taken into account in our simulation by Eq. (9). To avoid double-counting of the Ω effect, this term is omitted from our simulation.

Appendix D: Code test

D.1. Magnetic flux conservation

As the simplest test case, we have calculated a 1.5 M_⊙ stellar evolution with the magnetic field but switching off the magnetic dissipation and the Ω effect to test whether the magnetic field satisfies the flux conservation. The initial magnetic field is arbitrarily set to have a r⁻³ radial dependence for both the poloidal and toroidal components. The evolution is followed from the ZAMS phase through the TAMS and the red-giant phase until the star experiences helium flash and starts core helium burning, entering into the red clump in the Hertzsprung-Russell Diagram. During the evolution, the star experiences significant contraction in the central core and expansion in the outer envelope.

Figure D.1 shows the resulting evolution of the internal magnetic field. We note that B_r in the figure shows the polar value of the radial magnetic component, thus B_r = 2A(r)/r, and, B_ϕ in the figure is the toroidal component at θ = π/4, thus B_ϕ = B(r). 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.4 M_⊙. 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⁻³. The star at this phase has a convective hydrogen-burning core of ∼0.2 M_⊙. The convective turbulence explains the magnetic diffusivity of ∼10¹² cm² s⁻¹, 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⁹ 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⁰ 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¹³ s for the further extension of the steady region. Similarly to the poloidal component, the toroidal component also reaches a steady state at ∼10⁹ s, and the steady region extends further with a longer timescale of ∼10¹³ 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² 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 be determined by solving the equation

which yields 0 = m² − (n − 1)m − (n + 6), where m ≡ ∂lnB/∂lnr and n ≡ −∂lnη/∂lnr. The dissipative region can be divided into the fully convective region of ≲0.17 M_⊙ and the surrounding overshooting region at ∼0.17−0.21 M_⊙. 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⁹ s, the magnetic field finds a steady state of B_r ∝ B_ϕ ∝ r⁻³, 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.

Movie

Movie of Fig. 3 Access here

All Tables

All Figures