# A Model of Magnetic Braking of Solar Rotation That Satisfies Observational Constraints

## Abstract

The model of magnetic braking of solar rotation considered by Charbonneau & MacGregor (1993) has been modified so that it is able to reproduce for the first time the rotational evolution of both the fastest and slowest rotators among solar-type stars in open clusters of different ages, without coming into conflict with other observational constraints, such as the time evolution of the atmospheric Li abundance in solar twins and the thinness of the solar tachocline. This new model assumes that rotation-driven turbulent diffusion, which is thought to amplify the viscosity and magnetic diffusivity in stellar radiative zones, is strongly anisotropic with the horizontal components of the transport coefficients strongly dominating over those in the vertical direction. Also taken into account is the poloidal field decay that helps to confine the width of the tachocline at the solar age. The model’s properties are investigated by numerically solving the azimuthal components of the coupled momentum and magnetic induction equations in two dimensions using a finite element method.

## 1 Introduction

Helioseismology has revealed that the Sun’s radiative core rotates nearly uniformly, at least down to the radius , with a rate that is close to the mean rotation rate of its convective envelope (e.g., Tomczyk, Schou, & Thompson 1995; Couvidat et al. 2003). It is also known that the Sun, like other solar-type stars, has been experiencing an effective braking of its surface rotation via a magnetized solar (stellar) wind (Schatzman 1962; Weber & Davis 1967; Kawaler 1988). Because the magnetic field lines that sling charged particles into space are rooted in the photosphere, it is the envelope rotation that should be decelerated by the wind torque, while the conservation of angular momentum in the radiative core should keep it rapidly rotating. This expectation appears to be in a stark contrast with the measurements of internal solar rotation, which indicate that the difference in the angular velocity between the Sun’s radiative core and convective envelope must have been erased by some angular momentum redistribution mechanism. A similar conclusion can be reached from comparisons of the rotation period distributions for solar-type stars in open clusters of different ages. It turns out that the core and envelope rotation have to be coupled with a characteristic timescale of angular momentum transfer between them increasing from a few Myr to nearly a hundred Myr between the fastest and slowest rotators (Denissenkov et al. 2009; cf. Irwin et al. 2007). A simple model of rotational evolution of solar-type stars with angular momentum transport described as a purely diffusive process shows that this coupling time interval roughly corresponds to a diffusion coefficient decreasing from cm s to cm s. However, such large diffusion coefficients would lead to rates of surface Li destruction strongly exceeding the trend observed in the Sun and its twins (Section 2.4). Therefore, we have to admit that the angular momentum redistribution in solar-type stars is unlikely to be accompanied by an equally rapid mass transfer. This restricts our search for possible mechanisms of angular momentum transport in solar-type stars to those employing waves and magnetic fields.

Presently, the most popular mechanisms of angular momentum redistribution in the Sun’s radiative core are the smoothing of its rotation profile by internal gravity waves (g-modes) generated by turbulent eddies in the convective envelope (Charbonnel & Talon 2005) and magnetic braking of differential rotation (Charbonneau & MacGregor 1992, 1993; Spruit 1999, 2002). Using a simplified analytical prescription for the spectrum of internal gravity waves and other approximations, Charbonnel & Talon (2005) have succeeded in obtaining a nearly flat rotation profile in the present-day Sun simultaneously with a time evolution of its internal rotation consistent with the observed solar Li abundance. However, because of the complexity of physical processes involved in the generation, propagation and dissipation of internal gravity waves in the Sun, it is difficult to assess the validity of this solution. For instance, it would not be valid if a flatter spectrum, like that computed by Rogers & Glatzmaier (2005), were used. Non-linear wave-wave interactions, as well as the interaction of gravity waves with a magnetic field, could also change the solution drastically (Rogers, MacGregor, & Glatzmaier 2008; Rogers & MacGregor 2009). Even in a case very similar to that considered by Charbonnel & Talon (2005) the question arises as to why gravity waves do not disturb the uniform solar rotation? Indeed, the power of gravity waves generated in the present-day Sun should approximately be the same as it was in the young Sun; consequently, given their peculiar anti-diffusive nature in shear flows (e.g., Plumb 1977; Ringot 1998), gravity waves should have forced the solar rotation profile away from the uniform one (Denissenkov, Pinsonneault, & MacGregor 2008). However, this has not occurred.

For a magnetic field configuration to influence differential rotation of the solar radiative core, it must have a non-vanishing azimuthal (toroidal) component. To achieve this, it is usually assumed that there is an axially symmetric poloidal magnetic field frozen into the plasma inside the radiative core which can be stretched around the rotation axis by the differential rotation itself to form the necessary toroidal field. As far as the origin of the poloidal field is concerned, it is widely believed that the field was either inherited from a protostellar cloud or generated by a convective dynamo during the Sun’s pre-MS evolution when its convective envelope occupied a much larger volume. Spruit (1999) proposed the original alternative explanation that the poloidal field might be constantly replenished by radially stretching the toroidal field in a dynamo-like cyclic process, the necessary radial plasma displacements being caused by the non-axisymmetric () kink instability of the toroidal field (the Tayler-Spruit dynamo).

Magnetic braking is implemented by the azimuthal component of the magnetic tension term of the Lorentz force. Given that its corresponding engagement time is of order several Alfvén timescales for the poloidal field, magnetic breaking appears to be very efficient even for a small field amplitude (Mestel & Weiss 1987). The only problem associated with this mechanism is that it leaves behind a torus-shaped “dead” zone inside the radiative core with the angular velocity increasing toward the center of the torus cross-section. This is a direct consequence of Ferraro’s law of isorotation (rotation that keeps the angular velocity constant along each field line) (Ferraro 1937) and the presence of a constant dipole poloidal field. To solve this problem, Charbonneau & MacGregor (1993) (hereafter, referred to as CHMG93) and other researchers (e.g., Ruediger & Kitchatinov 1996) artificially amplified the viscous transport with the justification that rotation and associated hydrodynamic instabilities should lead to turbulence, and hence that the microscopic viscosity should be augmented by a much stronger turbulent viscosity. However, the minimum value of the enhanced viscosity ( cm s) with which it is possible to get a nearly flat rotation profile in the present-day Sun results in too much Li being transported below the bottom of convective envelope and thus destroyed, inconsistent with the observed solar Li abundance. Therefore, the proposed isotropic viscosity enhancement cannot be considered a satisfactory solution.

Only recently has the physics-based modeling of angular momentum transport in solar-type stars begun to use the extensive data sets of rotation periods for solar analogs in open clusters as an observational constraint. The main goal of previous studies was to demonstrate that a particular model could (or could not) produce the solid-body rotation of the Sun. In particular, Eggenberger et al. (2005) showed that the Tayler-Spruit dynamo could account for the flat rotation profile of the Sun. However, the correct model must also explain the cause of the transition from the solid-body rotational evolution of the fastest rotators to the evolution wherein a radial differential rotation is sustained for a hundred Myr in the slowest rotators. It turns out that the Tayler-Spruit dynamo cannot pass this observational test because it always enforces a nearly uniform rotation in solar-type stars, no matter how slowly they rotate (Denissenkov et al. 2009).

Another important problem that the correct mechanism of angular momentum transport in the Sun has to explain is the thinness of the solar tachocline. The tachocline is a layer in which the differential rotation (both radial and latitudinal) of the convective envelope sharply changes to the nearly solid-body rotation of the radiative core (Spiegel & Zahn 1992). Its thickness, as measured by helioseismology, is only about 4% of the solar radius (Basu & Antia 2003). It has been suggested that no purely hydrodynamic mechanism can explain its existence and that a large-scale magnetic field in the radiative core is therefore needed to confine the tachocline structure (Gough & McIntyre 1998; Ruediger & Kitchatinov 2007).

To summarize: because large-scale magnetic fields have certainly been playing crucial roles during different phases of the Sun’s rotational evolution — such as the synchronizing of rotation of the young Sun and its surrounding protoplanetary disk (the disk locking) (e.g., Shu et al. 1994; Matt & Pudritz 2005), providing the leverage for the solar wind, confining the tachocline, being generated through a convective dynamo and then driving the solar activity — we believe that it is worthwhile to develop a mechanism of magnetic breaking of differential rotation of the Sun’s radiative core.

In this paper, we elaborate upon and extend the magnetic braking mechanism that was originally proposed by Mestel & Weiss (1987) and then studied in detail numerically by CHMG93 followed by Ruediger & Kitchatinov (1996). Our only modification of the mechanism is to assume that the turbulence induced by rotation-driven hydrodynamic instabilities in the Sun’s radiative core is highly anisotropic with its corresponding horizontal component of turbulent viscosity strongly dominating over that in the vertical direction . The hypothesis of anisotropic turbulent diffusion in stellar radiative zones was first advanced by Zahn (1992) and it had since been productively used by many other researchers to study rotational mixing in both MS and red giant branch stars (Talon & Zahn 1997; Talon et al. 1997; Maeder 2003; Palacios et al. 2003, 2006). Initially, it was even considered to explain the thinness of the solar tachocline by Spiegel & Zahn (1992). We will show that this hypothesis alone permits a solution of almost all of the aforementioned problems associated with this particular mechanism of magnetic braking, such as the Li problem (a sufficiently small value of the diffusion coefficient can be chosen so that it does not lead to excessive Li destruction), the dead-zone problem, and the difference in spin-down of the fastest and slowest rotators among solar-type stars in open clusters.

The paper is organized as follows. In Section 2, we briefly discuss simple models of rotational evolution of solar-type stars and review the main results that were obtained by comparison of the model predictions with the observed period distributions for solar counterparts in open clusters. In Section 3.1, we reproduce the 2D MHD computations of CHMG93 and show that they cannot account for the spin-down of the fastest rotators in open clusters, especially when one reduces the vertical viscosity to a value more or less compatible with the surface Li abundance evolution in solar twins. In Section 3.2, we present our new dynamic model of magnetic braking with the anisotropic turbulent diffusion () and discuss its advantages compared to the original model. Section 3.3 describes a number of stationary solutions that address the problem of the penetration of the latitudinal differential rotation from the bottom of convective envelope into the radiative interior in the present-day Sun. Finally, we summarize our conclusions in Section 4.

## 2 Simple Models of Rotational Evolution of Solar-Type Stars

The simple models of rotational evolution of solar-type stars — i.e., the double-zone model and the purely diffusive one — have been discussed in detail and compared with each other by Denissenkov et al. (2009). Here, we will only briefly describe the models and summarize the main results obtained by comparing their predictions with distributions of angular velocity of the (assumed) rigidly rotating convective envelope for solar analogs in open clusters, , which are derived from the corresponding observed rotation period () distributions. We will also show that the purely diffusive model fails to simultaneously comply with the observational constraints imposed by the spin-down of solar analogs in open clusters and the time evolution of the surface Li abundance in solar twins. Although the simple models do not identify the physics behind the internal transport of angular momentum in solar-type stars, they nevertheless provide us with useful estimates of characteristic timescales of relevant processes and their dependence on the rotation rate. Therefore, they can serve as a good starting point for our analysis.

### 2.1 The Pre-Main Sequence Disk Locking

In Fig. 1, crosses with vertical error bars represent the upper 90th and lower 10th percentiles of the distributions for stars having masses in the interval sampled from open clusters with ages between 2 Myr and 600 Myr (for details, see Denissenkov et al. 2009). The data for the youngest clusters set up the initial distributions for rotational evolution computations. Of course, these are not the true initial distributions inherited by the stars during their formation epoch because they have already been modified (reduced) by a disk-locking process compared to what they would have been if the stars’ pre-MS evolution had taken place in complete isolation. The fact is that after its birth a protostar continues to contract before it settles down on the zero-age MS (ZAMS) at an age of about 40 Myr (for ). During this period of time, the star’s radius and total moment of inertia are decreasing, hence its should be growing as a consequence of angular momentum conservation. Contrary to this, it is the surface angular velocity rather than the angular momentum that is observed to be nearly preserved in contracting pre-MS stars (e.g., Rebull, Wolff, & Strom 2004). This is usually explained by a complex magnetic interaction that pumps out angular momentum from the protostar to its surrounding protoplanetary dust disk (e.g., Shu et al. 1994; Matt & Pudritz 2005). As a result, the magnetic disk locking keeps nearly constant (this is modeled by horizontal fragments of dotted curves in Fig. 1) for a period of time up to 10 – 20 Myr at most before the disk disappears. After that, the star’s residual contraction and conservation of angular momentum cause an increase of its that continues until the star reaches the ZAMS.

### 2.2 The Surface Loss of Angular Momentum

From the moment when the young Sun arrives at the ZAMS ( Myr) until the Sun’s present-day age of Gyr, its internal structure does not change significantly. In particular, the dimensionless moment of inertia of the Sun’s convective envelope remains constant to within 6%. Note that amounts only to 14% of the Sun’s total moment of inertia . If there were no angular momentum transfer from the Sun’s radiative core to its convective envelope then the only process involved in the Sun’s rotational evolution during its life on the MS () would be the braking of envelope rotation by the magnetized solar wind. The corresponding rate of angular momentum loss is often approximated by the following equation:

(1) |

where is a constant calibrated to give rad s (for days), and is a mass-dependent velocity above which the wind gets saturated (Chaboyer, Demarque, & Pinsonneault 1995; Krishnamurthi et al. 1997; Andronov, Pinsonneault, & Sills 2003). The value of has been adjusted so that the upper solid curve in Fig. 1 simulating the solid-body rotational evolution of rapidly rotating open-cluster solar-type stars could fit the upper 90th percentiles of their distributions as closely as possible (Denissenkov et al. 2009).

The differential equation (1) can easily be solved. If the initial angular velocity exceeds the saturation threshold then the solution consists of two different parts. The first one, which is valid for , is

(2) |

where the wind braking timescale is , while is the time when has decreased to a value of . To shorten the equations, we have used the notations and . The exponential decay of is replaced by the Skumanich relation (Skumanich 1972)

(3) |

for the time interval . Finally, the solar-calibrated wind constant has to be equal to

(4) |

Alternatively, if the star arrives at the ZAMS with then its surface angular velocity will be declining according to the Skumanich law from the very beginning. The above equations can still be used provided that one substitutes in all of them.

Equations (2) – (3) suffice to describe the Sun’s spin-down for the two limiting cases in which its core and envelope rotation are either completely decoupled or, alternatively, rigidly coupled. For the first case, which is only of academic interest because it is in conflict with uniform rotation in the solar interior, equation (4) gives cm g s for , and cm g s for (the upper and lower solid curves in Fig. 1). In the second case, which appears to adequately describe the fastest rotators, we have to substitute the total moment of inertia instead of into (4), resulting in cm g s and cm g s for the same initial angular velocities. It is interesting that the solutions (2) – (3) do not depend on the moment of inertia; therefore, the curves representing them for the two limiting cases completely coincide in Fig. 1. The entire difference between the cases is contained in the values of the wind constant . These values have to be larger, in proportion to the ratio , for the solid-body rotational evolution because in this case the wind has to be faster to remove the correspondingly increased amount of angular momentum, so that finally we could get . Note that the last constraint comes directly from the observed rotation period distributions for solar-type stars in open clusters that show a convergence of to by the solar age. It is also clear that the correct mechanism of angular momentum redistribution in solar-type stars should not tolerate a large variation of the wind constant between its applications to the fastest and slowest rotators as long as we believe that the dependence of on has correctly been factored out in equation (1).

The illustrated degeneracy of the solutions (2) – (3) with respect to the two limiting cases of coupling between the core and envelope rotation is lifted when we begin to consider more realistic cases in which the coupling is allowed to increase progressively in time. To demonstrate this via a simple example, let us assume that the moment of inertia of the star’s rigidly rotating envelope that is decoupled from the rest of it obeys the following rule:

(5) |

where . Although this case does not have a physical justification, it admits an analytical solution of equation (1) and it also guarantees that (no core/envelope coupling at the beginning) and (full coupling at the end). For , the solution of equation (1) for the dependence (5) is plotted in Fig. 1 with a dot-dashed curve. It needs cm g s. Thus, in spite of the fact that the evolution begins with angular momentum being lost only from the convective envelope, the wind constant has a value much closer to the previously considered case of solid-body rotational evolution. This is explained by the fact that the star will later have to get rid of more angular momentum than in the old decoupled case because the moment of inertia of its decoupled and uniformly rotating envelope (not just the convective envelope but also an outer part of the radiative core adjoint to it) increases well above . Therefore, it is necessary to lose additional angular momentum from the surface of the convective envelope in the future, which will be supplied by the radiative core as its rotation gets more and more coupled to that of the envelope, demanding a larger wind constant which, in turn, leads to a steep initial decline of .

### 2.3 The Double-Zone Model

The double-zone model (MacGregor 1991) assumes that the Sun (or a solar-type star) consists of two uniformly rotating zones, the core and envelope, and that an excess of angular momentum in the core, compared to the case in which the whole star rotates rigidly, is transferred to the envelope on a specified constant timescale . It is evident that, for this core-envelope rotational coupling to replenish the angular momentum content of the envelope faster than it is drained by the wind, one needs to have Myr (for ), and vice versa. Computations show (e.g., Denissenkov et al. 2009) that the solid-body rotational evolution of a solar-type star, in the case that its core and envelope rotate nearly synchronously all the time, can only be achieved with Myr.

It turns out that the upper 90th percentiles of distributions (the fastest rotators) for solar analogs in open clusters are very well reproduced by the double-zone model with Myr (upper dotted curves in panels A and B in Fig. 2). On the contrary, the slowest rotators (the lower 10th percentiles) at ages of order 100 Myr are found to be located below the curves representing the solid-body rotational evolution ( Myr) that starts from the 10th percentiles of distributions for the youngest clusters, even when the safe upper limit of Myr is chosen for the disk-locking time (the lower dotted curve in Fig. 2A). A rigorous statistical analysis that uses Monte Carlo simulations and compares the full distributions (modeled versus observed ones) rather than just their percentiles confirms this conclusion. It also gives the following estimates of the coupling time for the slowest rotators: Myr for stars with and Myr for stars with (Denissenkov et al. 2009).

### 2.4 The Purely Diffusive Model

The purely diffusive model uses appropriate initial and boundary conditions, including (1), to solve the following equation:

(6) |

where is a constant viscosity whose physical nature is not specified. The time-dependent density distribution in the radiative core and other necessary stellar structure parameters are taken from full stellar evolution computations. Denissenkov et al. (2009) have demonstrated that the solid-body rotational evolution of solar-type stars can only be simulated with cm s (solid curves in Fig. 2). They have also established an approximate correspondence between the diffusion coefficient and the coupling time from the double-zone model that produce a similar rotational evolution. In particular, this relation says that the spin-down of a slowly rotating star computed using the double-zone model with Myr looks almost identical to that obtained with the purely diffusive angular momentum redistribution for cm s (dot-dashed curves in Fig. 2). This solution gives an example similar to (5) when the moment of inertia of the rigidly rotating envelope effectively increases with time as the diffusive transport of angular momentum gradually penetrates into the radiative core. It should be noted that even for a coupling time as long as 90 Myr the corresponding diffusion coefficient still remains as large as cm s.

The angular momentum redistribution in solar-type stars is known to be accompanied by a chemical element transport, the most pronounced manifestation of which is the strongly reduced abundance of Li in the solar atmosphere. Meléndez et al. (2009) have recently reported that the Sun does not seem to be unique in this respect. It turns out that solar twins (stars with mass and metallicity not different from the solar values to within a few percent) show a dependence of the surface Li abundance on age to which the Sun also belongs (Fig. 3). A natural interpretation of this relation is that a rate of Li destruction in the atmospheres of solar twins is approximately the same, the Sun simply being a relatively old star.

We have added diffusion terms to a network of nuclear kinetics equations relevant for the Sun’s evolution (the pp chains and CNO cycle reactions). This network has been solved in a post-processing way using, as a background, previously stored files with all necessary information about the Sun’s internal structure from its fully convective pre-MS evolution to its present age. The resulting time evolution of the surface Li abundance is plotted in Fig. 3 for several values of the diffusion coefficient. It is seen that only a value of - cm s satisfies the observed anti-correlation. On the other hand, a value of cm s that was used by CHMG93 to model the magnetic braking of solar rotation leads to an excessive Li depletion inconsistent with the observations. In our modified model of magnetic braking of solar rotation with anisotropic turbulent diffusion, we will employ a value of cm s which is close to the values constrained by the Li abundance. Note that, without being assisted by another angular momentum transport mechanism, the purely diffusive mixing with this small coefficient can reproduce neither the Sun’s nearly uniform rotation nor the spin-down of the fastest rotators in open clusters (dashed curves in Fig. 2).

Finally, it is important to note that Be, another fragile element, does not seem to share the fate of Li in evolved MS stars with close to solar effective temperatures, including the Sun (Randich 2008). This means that the chemical mixing that causes the Li destruction does not penetrate below in the Sun (eg., see Fig. 1 of Barnes, Charbonneau & MacGregor 1999), being the radius of its core-envelope interface, where the temperature is high enough for protons to begin destroying Be. Alternatively, it is possible that the diffusion coefficient rapidly declines with depth (e.g., Richer, Michaud, & Turcotte 2000). If this is true then the Li mixing may have nothing to do with the global angular momentum redistribution in the Sun. On the other hand, to be in accord with observations the latter process cannot be assisted or accompanied by a mass transfer whose rate, when described in terms of a diffusion coefficient, greatly exceeds cm s.

## 3 Dynamic and Stationary Models of Solar Rotation with Large-Scale Magnetic Fields

### 3.1 Magnetic Braking of Solar Rotation: the Old Solution

Because our main goal is to modify the well-known model of magnetic braking of solar internal differential rotation for it to enable to explain the recent observational data on the distributions for solar-type stars in open clusters, without causing a conflict with the Li abundance data for solar twins, we will describe the model employed by CHMG93 that we want to develop only briefly. For more details about the model, the interested reader is referred to the original paper for an excellent analysis of numerical results obtained with it and also a discussion of its shortcomings. To save space, we introduce, at the outset, the general forms of relevant equations that encompass both the original case considered by CHMG93 and our modified case with anisotropic turbulent diffusion. In fact, we have chosen an equivalent system of equations similar to that used in the follow-up study by Ruediger & Kitchatinov (1996), which allows the poloidal field potential to fade with time.

In the spherical polar coordinates , the magnetic breaking of solar rotation in the presence of anisotropic turbulence in the radiative core is described by the following equations

(7) | |||||

(8) | |||||

(9) | |||||

where is a normalized angular momentum loss rate, with and being, as before, the moment of inertia and angular velocity of the convective envelope. For the rate of angular momentum loss from the surface , we continue to use equation (1), whereas CHMG93 used the formulation of Weber & Davis (1967).

Equations (7) – (8) represent the - (azimuthal) components of the momentum and induction equations. Following CHMG93, we neglect any meridional motion of either rotational or magnetic origin. We also omit the Coriolis force, while solving the problem in an inertial (non-rotating) frame of reference in which the polar axis is directed along the rotation axis. The right-hand side of the momentum equation (7) contains the viscous term, which has been split into its vertical and horizontal components, the magnetic tension part of the Lorentz force , where the total (poloidal plus toroidal) field is in which the poloidal field has been expressed via its potential as follows

(10) |

and, finally, the term simulating the surface angular momentum loss by smearing over the entire convective envelope, which therefore must vanish in the radiative core.

Equation (9) is decoupled from the first two equations, therefore it can be solved separately. Following Ruediger & Kitchatinov (1996), we expand the poloidal field potential into a series of the associated Legendre polynomials and keep only the longest living dipole term

(11) |

which we substitute into (9). As a result, we have to solve the eigenvalue problem

(12) |

for the boundary conditions , where , and the potential will thus be determined.

CHMG93 considered the isotropic case in which and . They artificially amplified the microscopic kinematic viscosity and magnetic diffusivity by factors of order and , respectively, making them as large as cm s and cm s and keeping their values constant throughout the radiative core. It is convenient to transform equation (12) into its dimensionless form

(13) |

where , , and . Taking the values of and cm s used by CHMG93, we solve the eigenvalue problem (12 – 13) and find that, in this case, the poloidal field potential can be approximated by the following relation:

(14) |

where Myr for . The same, or very similar, form of the potential was used by both CHMG93 and Ruediger & Kitchatinov (1997). However, in both of those investigations the time dependence of was ignored, which might not be a good assumption, given the short decay time for obtained with that for the enhanced magnetic diffusivity.

The decay time can be as high as 5.4 Gyr, which is comparable to the solar age, only when the maximum value of the microscopic magnetic diffusivity in the Sun’s radiative core cm s is substituted into the above expression for instead of (cf. Ruediger & Kitchatinov 1996). But, in this case, magnetic braking will hardly work, especially, after our having correspondingly decreased the amplified macroscopic viscosity either to the maximum microscopic value of cm s or to its minimum turbulent value of (assuming that is still of turbulent origin but has also reached its minimum possible value). The latter substitution brings the turbulent magnetic Prandtl number close to its most plausible value of unity.

Keeping this inconsistency in mind, we proceed with our revision of the original model of CHMG93. Like them, we apply equations (7) – (8) and (14) to the whole star assuming that and jump to their large convective turbulent values of order cm s at the core-envelope interface. Solutions for the angular velocity and toroidal magnetic field are sought in the space domain on the time interval . Following CHMG93, it is convenient to use as the second independent variable instead of , in which case the relevant space domain turns into the rectangular area , where .

We start our computations with the following initial conditions: , and . The model’s internal structure is considered to be fixed and described by the density distribution for the present-day solar model taken from our full stellar evolution computations (Denissenkov et al. 2009). For the dimensionless moment of inertia of the convective envelope, we use the time-averaged value of . Given that the surface loss of angular momentum is already included in equation (7) and that the equations are solved for the whole star, the boundary conditions become simple: , at each point on the rotation axis, and on all the boundaries, including the core-envelope interface, the latter condition being imposed because of the extremely large turbulent magnetic diffusivity in the convective envelope that makes it feel like a vacuum to the internal magnetic field (Ruediger & Kitchatinov 1996).

We consider two poloidal field configurations of the four cases investigated by CHMG93. They are specified by the parameter that defines the radial extent of non-vanishing potential with respect to (equation 14). We have chosen the qualitatively different configurations D3 and D2 in the notations of CHMG93. In the most promising, according to CHMG93 and MacGregor & Charbonneau (1999), configuration D3 the potential occupies the whole radiative core up to the bottom of convective envelope, i.e. we have (Fig. 5A). The configuration D2 is set up by using the value of , which means that can now penetrate into the lower part of the convective envelope (Fig. 5B) but, what is more important, the poloidal field has now a non-vanishing radial component at for some . Given that in the configuration D2 is thus directly coupled to the envelope rotation, it can more effectively transfer into the radiative core the torque applied to the convective envelope by the magnetized wind. In the D3 case, there is a time delay before envelope braking will be felt in the core. This time is taken for the viscous transport of angular momentum to build up a sufficiently strong rotational shear just beneath the core-envelope interface which will then begin to interact with . For a comparison, we also introduce a purely diffusive configuration D0 with (Fig. 7A). Finally, our modified magnetic braking models with anisotropic turbulent diffusion that have the D3 and D2 poloidal fields and strong horizontal components of turbulent viscosity in the radiative core will be referred to as D3H and D2H.

To numerically solve the partial differential equations (7) – (8) in the 2D space domain with the potential given by equation (14) we make use of the COMSOL Multiphysics software package. Like the custom-made computer code employed by CHMG93, the COMSOL Multiphysics solves PDEs using the finite element method. In our computations, we use quadrilateral mesh elements to partition the space domain into a total number of elements that slightly varies from one solution to another but that is approximately equal to and comparable to the mesh sizes used by CHMG93. Note that we seek a solution in a half meridian, whereas CHMG93 did their computations for a quadrant. This reduces the resolution of our mesh by a factor of two compared to theirs. Otherwise, apart from the different prescription for the surface loss of angular momentum and the option of taking into account the exponential decay of , our analysis of magnetic braking of solar rotation with the isotropic viscosity and magnetic diffusivity is very similar to that by CHMG93.

Analyzing the results of their extensive computations, CHMG93 have come to the conclusion that the rotational evolution of the Sun weakly depends on the poloidal field amplitude (equations 10 and 14). On the contrary, the average strength of the toroidal field that is generated by differential rotation in the radiative core changes significantly with , a weaker poloidal field resulting in a stronger toroidal field. This anti-correlation emerges from the necessity of the magnetic torque to compensate for the one exerted by the magnetized solar wind, the former being proportional to the products and in which and because the total magnetic field in the radiative core is .

CHMG93 have distinguished the following three important epochs in the magnetic braking of solar rotation that differ from each other by the relative strength of the magnetic and wind torques. During the first epoch, the magnetic torque is much weaker than the wind torque and, as a consequence, the rotational shear produced by the spin-down of the convective envelope at its interface with the radiative core is spread by the viscosity almost without hindrance to the core surface. So, the shear in the outer core grows almost linearly with time and, as a result, the toroidal field is amplified there almost quadratically in time. This short epoch, lasting between a few thousand to a few million years, is followed by the second epoch during which time the differential rotation profile begins to interact with the growing toroidal field. For the poloidal field configurations in which the core and envelope are magnetically coupled (e.g., the D2 case), this epoch is characterized by vigorous phase mixing throughout the radiative core. This process is initiated by large scale toroidal field oscillations that originate close to the bottom of convective envelope and propagate through the radiative core along the poloidal field lines using them as elastic rubber strings. Because the strings are rooted in the convective envelope (for the D2 configuration) at different colatitudes and therefore have different lengths, the Alfvén waves carrying the perturbations back and forward along the strings get out of phase as time progresses. Eventually, this causes large gradients of toroidal magnetic field to build up locally at different places and, as a result, the large-scale toroidal field oscillations are diffused away either via ohmic dissipation or via turbulent magnetic reconnection, depending on the nature of (e.g., Spruit 1999). Finally, during the third epoch, a state of dynamical balance that is achieved at the end of the second epoch between the magnetic torque in the radiative core and the total stresses (magnetic plus viscous) at the bottom of convective envelope is maintained as the surface wind torque decreases with time.

First, we have repeated the computations of CHMG93 for one of their most favorite models, namely the one that starts with and has the D3 configuration of time-independent poloidal field with the amplitude G. Our results are presented in Fig. 4 (thick solid curves in panels A and C as well as panels B and D). The geometry and magnitude of the residual toroidal field in the present-day Sun (panel D) look very similar to those obtained by CHMG93 (upper quadrant of the lower right panel in their Fig. 6). As they reported, the model has a nearly flat rotation profile by the solar age (panels B and C). However, it seems unlikely that it will be able to fit the upper 90th percentiles of the distributions for open-cluster solar analogs after its initial rotation rate is increased to a value of corresponding to the observed maximum (the peak on the dotted curve representing the solid-body rotational evolution simulated using the double-zone model). This supposition is confirmed by our additional test computations that started with (thin solid curve in panel A).

It is also not clear what could possibly force the rotational evolution of the CHMG93 model to change its morphology in the vs. Age plane from the convex one for the fastest rotators to the concave one for the slowest rotators, as is indicated by both direct observations (e.g., compare the solid and dot-dashed curves in panels A and B in Fig. 2, respectively) and their rigorous statistical analysis (Denissenkov et al. 2009). Physically, the change implies a transition from solid-body rotational evolution to evolution with a differential rotation sustained for a few tens to a hundred Myr. As has been mentioned, this cannot be produced by a variation of . We can see now that an increase of the initial rotation rate does not make the rotational evolution of the CHMG93 model approach the dotted curve either. In other words, this model does not have enough internal degrees of freedom to explain the whole spectrum of observed rotation rates. Besides, its high viscosity ( cm s) causes the model to conflict with the observed Li abundances in solar twins (Section 2.4). When we reduce the viscosity to the value of cm s, which is more or less compatible with the Li data, then the model’s rotational evolution becomes absolutely incapable of reproducing the upper 90th percentiles of the distributions (dashed curve in Fig. 4A).

As far as the importance of taking into account the exponential decay of the poloidal field is concerned, it depends on how big the values of and are in the model in question. For the above considered case of , and cm s, a more consistent solution with decaying on the timescale of Myr (equation 14) is plotted in Fig. 4A as a dot-dashed curve. Although it does not appear to strongly deviate from the old solution (the thin solid curve), the difference between the two solutions will probably be significant for Monte Carlo simulations and a rigorous statistical analysis of rotation period distributions for solar analogs in open clusters like those performed by Denissenkov et al. (2009).

Finally, it is worth noting that the purely diffusive model with the same viscosity cm s as the one used by CHMG93 seems to produce both the Sun’s rotational evolution and its final internal rotation profile (the dot-dashed curves in panels A and C in Fig. 2) that are not very different from those computed with the model of CHMG93, especially when the poloidal field potential is allowed to decrease with time (the dot-dashed curves in panels A and C in Fig. 4). The secondary role of magnetic braking compared to the viscous transport of angular momentum is emphasized when we take a much smaller value of cm s, consistent with the solar Li data (compare dashed curves in the same four panels). In this case, the only noticeable effect produced by the presence of magnetic fields is a slightly slower rotation of the radiative core at the solar age, especially close to its center, which still strongly deviates from the helioseismology data.

### 3.2 Magnetic Braking With Anisotropic Turbulent Diffusion

We start with the aforementioned failed D3 model of magnetic braking of solar rotation that has an isotropic viscosity cm s and magnetic diffusivity cm s and which begins its rotational MS evolution with (thin dashed curves in Fig. 6). Note that its viscosity is already increased by nearly two orders of magnitude compared to the maximum value of the microscopic (molecular for the Sun) viscosity in the radiative core. This is not a problem though, given that vertical diffusion with a coefficient of a similar or slightly stronger magnitude is indeed predicted to occur in radiative zones of solar-metallicity MS stars with , especially close to the ZAMS (Palacios et al. 2003), where internal gradients of the mean molecular weight are too weak to effectively hinder it. This diffusion is associated with either rotation-driven meridional circulation or small-scale turbulence caused by rotational shear instabilities. It has been successfully used to explain chemical composition anomalies in the atmospheres of massive MS stars and their descendants (e.g., Meynet & Maeder 2000; Denissenkov 2005).

Whereas the model of magnetic angular momentum transport based on the hypothetical Tayler-Spruit dynamo has a problem reproducing differential rotation in the slowest rotators (Denissenkov et al. 2009), the D3 model of CHMG93, on the contrary, cannot provide a rapidly rotating solar-type star with a means to spin-down its radiative core fast enough for the star’s rotational evolution to fit the upper 90th percentiles of the distributions for solar analogs in open clusters. The question is whether we can make a reasonable modification of the CHMG93 model such that it will become capable of reproducing the rotational evolution of both the fastest and slowest rotators while being consistent with the other observational constraints? We have shown that a variation of the isotropic viscosity does not help to solve the problem self-consistently with or without a magnetic field. The major obstacles are the observational limit of cm s imposed by the rate of the evolutionary Li depletion in solar twins (Fig. 3) and the rapid decay of the poloidal field caused by high viscosity and magnetic diffusivity. On the theoretical side, stellar hydrodynamics does not predict a build-up of vertical turbulent viscosities as large as cm s in the radiative cores of rapidly spinning MS solar-type stars, as required for the stars to evolve like solid-body rotators. Besides, such large diffusion coefficients would solve the problem without being assisted by magnetic fields.

It should be noted that the poloidal magnetic field does help the viscous transport to reduce, even a little, the degree of differential rotation in the Sun’s radiative core, as claimed by CHMG93. Indeed, panels C in our Figs. 7 and 5 indicate that, whereas the contours of constant angular velocity in the present-day Sun’s D0 model are concentric spheres, they become consecutive layers in meridian cross sections of an axisymmetric torus in the D3 model. This torus, called the “dead zone” by Mestel & Weiss (1987), is shaped by the poloidal field, whose lines form a similar toroidal structure (Fig. 5A), according to Ferraro’s law . Thus, in the D3 configuration, the viscous transport has to redistribute angular momentum throughout a smaller volume than in the absence of magnetic fields (the D0 case). However, when cm s, it does not have enough time to finish its work by the solar age. As CHMG93 reported, the situation is slightly improved when the poloidal field configuration is switched to D2. In this case, there is no waiting time before the spin-down of the convective envelope gets communicated to the radiative core because poloidal field lines are now rooted in the envelope, hence they can start transferring the wind torque to the core from the very beginning. Besides, because not every poloidal field line is now parallel to the bottom of convective envelope (some of them emanate from it, being directed inside the core), a larger cylindrical region around the rotation axis (not just the axis and its immediate adjacent layers as it was in the D3 case) acquires a rotation rate close to that of the envelope. As a result, the dead zone becomes thinner (Fig. 5D). However, even in this case the rotational evolution is still far from being close to that of a solid-body rotator (thick dashed curves in Fig. 6).

When analyzing the above results, especially contour plots C and D in Fig. 5, we have noticed that a strong horizontal turbulence, a concept with different approximate descriptions that have already been used in other stellar evolution computations, might solve the problem of magnetic braking of solar rotation as well. In fact, the same hydrodynamic models that give an estimate of – cm s for the vertical component of turbulent diffusion also predict a much stronger value, of order – cm s, for its horizontal component in the radiative core of the solar-metallicity MS star with and (Palacios et al. 2003). In a solar-mass star the corresponding turbulent viscosities are expected to be lower by a factor of 10, which is roughly proportional to the ratio of luminosities for MS stars with and , because the baroclinic and secular Kelvin-Helmholtz instabilities driving this anisotropic turbulence have growth rates inversely proportional to the thermal timescale. On the other hand, the coefficient of horizontal turbulent diffusion has been suggested to take on much greater values (by more than four orders of magnitude) than previously thought (cf. Zahn 1992; Maeder 2003; Mathis, Palacios, & Zahn 2007).

Zahn (1992) was the first to bring anisotropic turbulent diffusion to the attention of the stellar astrophysics community. Anisotropic turbulence has both quite a natural physical justification and an experimental basis in geophysics and oceanography. Tassoul (2000) notes that “in the Earth’s lower atmosphere one has , whereas this ratio may be as large as in the surface layer of the ocean where large-scale currents are observed”. It is even more likely to occur in stellar radiative zones with their extremely large Reynolds numbers and strong thermal stratifications in the vertical direction.

Unfortunately, neither nor can be derived from first principles yet. Therefore, crude estimates are usually adopted and then tuned up by making an order of magnitude evaluation of the relevant energy balance, using appropriate results from laboratory or numerical experiments, and, last but not least, trying to comply with available observational constraints concerning mixing in stars. Given this uncertainty, while pursuing our modest goal of producing a modified CHMG93 model that should be as simple as the original one, we have chosen the following prescriptions for the anisotropic diffusion coefficients in solar-type MS stars:

(15) |

The first expression represents the approximate upper limit for the diffusion rate associated with the radial mass transport that is more or less consistent with Li data in solar twins. The more realistic value would be closer to cm s (dot-dashed curve in Fig. 3). However, the observational data do not exclude the possibility that could be as large as cm s during the first 300 Myr and then reduced to cm s (dotted curve in Fig. 3; however, our rotation evolution computations show that this variable viscosity does not solve the problem, as is illustrated by Fig. 8). In other words, our assumed constant value for should instead be considered a time average, therefore it can be a bit larger, say like cm s (short-dashed curve in Fig. 3). But it definitely cannot be as large as cm s (solid curve in Fig. 3) and remain constant, as CHMG93 assumed. As for our choice of , its amplitude seems to have a correct order of magnitude (in agreement with other most recent applications) for a ZAMS solar model rotating with the angular velocity , taking into account that Mathis, Palacios, & Zahn (2007) have estimated – cm s for a young (252 Myr) solar-metallicity model star that started its MS evolution with approximately one half the rotation rate. On the other hand, the assumed dependence of on does not seem unnatural and it provides us with the necessary additional degree of freedom, such that a variation of its value can probably be used to model the transition from the fastest to slowest rotators, without causing a conflict with other observational constraints.

Thick solid curves in Fig. 6 have been computed using the anisotropic diffusion coefficients (15) for the D2 poloidal field configuration. We have assumed that cm s and . The poloidal field decay has been taken into account (equation 14) with an appropriate value of the decay time obtained as a solution of equation (13). We can see that now the surface rotation evolution curve in panel A fits the upper 90th percentiles of the distributions for solar analogs in open clusters, while its corresponding internal rotation profile at the solar age appears to be nearly uniform (panel B). A comparison of four similar models in panel A demonstrates that it is the combination of the D2 poloidal field geometry and strong horizontal turbulence (the D2H model) that results in the rotational evolution resembling that of the fastest rotators because neither the D2 nor the D3H models are able to produce a convex evolutionary curve. The contours of constant angular velocity and toroidal magnetic field for the D2H model of the solar age are shown in panels B and D in Fig. 7. It is important to note that, in the absence of internal large-scale magnetic fields, strong horizontal turbulent diffusion could not erase the differential rotation in the radiative core because it is the dipole poloidal field and its interaction with the rotational shear that change the geometry of differential rotation from the spherically symmetric one (Fig. 7C) to the toroidal one (Fig. 5D). On the other hand, the assumption of anisotropic turbulence with allows us to choose a sufficiently small value of , thus avoiding a conflict with the observed surface Li abundances in solar twins.

Finally, when we apply the D2H configuration to model the rotational evolution of solar-type stars with initial angular velocities ranging from to , we indeed obtain a nice transition from the convex evolutionary curves for the fastest rotators to the concave ones for the slowest rotators (Fig. 9), as hinted by observations. Note that all of these models end up with nearly flat internal rotation profiles, and their solar-calibrated wind constants do not differ by very much: whereas the fastest rotating model has cm g s, the slowest rotating model needs cm g s.

### 3.3 Stationary Solutions: the Width of the Solar Tachocline

The conclusion made by MacGregor & Charbonneau (1999) that it is the D3 configuration, rather than the D2 one, that seems to be the most appropriate one for the model of magnetic braking of solar rotation is based on their consideration of the thinness of the solar tachocline. Because of the direct magnetic coupling of the convective envelope with the radiative core through the D2 poloidal field, the latitudinal differential rotation at the core-envelope interface penetrates into the bulk of the radiative core more easily in the D2 case than in the D3 case. The differential rotation of the convective envelope in the present-day Sun is maintained by the inhomogeneous turbulence of a rotating fluid (e.g., Kitchatinov & Ruediger 1993; Ruediger & Kitchatinov 1996). Following Charbonneau, Dikpati, & Gilman (1999) and Ruediger & Kitchatinov (2007), we use the simple approximate relation

(16) |

Its values at , , , and are plotted in panels C and D in Fig. 2 as horizontal line segments at . It is seen that the variation of the interface angular velocity with the colatitude in the present-day Sun is much smaller than the decrease of during the Sun’s possible previous MS evolution (Fig. 2A). Therefore, we have neglected this variation in our rotational evolution computations. However, now, when we consider the final model of the solar age, we have to be sure that our assumptions and, in particular, our preference for the D2 poloidal field configuration, do not lead to the formation of a thick tachocline.

A stationary solution that has to be tested for the spread of the interfacial differential rotation can easily be obtained with the same COMSOL Multiphysics code that we have used to compute the rotational evolution of solar-type stars. We have only to set to zero the left-hand sides of equations (7) – (8), omit the term proportional to , use the relation (16) as a new outer boundary condition, and shrink the computational space domain from to . We have used this approach to obtain a number of stationary solutions, the results of which are discussed below.

Let us start with non-magnetic solar models. Fig. 10A shows a stationary solution (contours of constant angular velocity) for one such model in which the kinematic viscosity is isotropic and has its near maximum microscopic value of cm s that is assumed (as before) to be constant throughout the radiative core. This solution is evidently in conflict with helioseismology data indicating that both the latitudinal and radial differential rotation of the solar convective envelope are smoothed out when one crosses the thin () tachocline, so that the radiative core rotates like a solid body at least down to (e.g., Ruediger & Kitchatinov 2007, and references therein). This result is not surprising because the isotropic viscosity is known to spread a horizontal velocity inhomogeneity of a given length scale to a radial velocity inhomogeneity of the same length scale (Ruediger & Kitchatinov 1997). Fig. 10B presents a stationary solution for the same non-magnetic model but with the anisotropic viscosity whose horizontal component is cm s cm s. In this case, the strong horizontal turbulence smooths out the interfacial differential rotation in a sufficiently thin layer located immediately beneath the core-envelope interface (cf. Spiegel & Zahn 1992).

Now, we return to the isotropic viscosity but add a poloidal magnetic field. Fig. 10C shows a stationary solution for the case similar to that considered by CHMG93, namely for the D3 poloidal field configuration of amplitude G interacting with the rotational shear via the (turbulent) viscosity cm s, the induced toroidal field being diffused away at the rate cm s. In this case, our contours of constant angular velocity indeed resemble very closely those plotted in the lower right panel in Fig. 1 by MacGregor & Charbonneau (1999). After returning to the isotropic viscosity, we would expect the interfacial differential rotation to again be spread radially deep into the core but this has not happened. The reason is that now it is the latitudinal transport of angular momentum by the magnetic stress, proportional to the product , that suppresses the differential rotation and restricts its penetration into the radiative interior (Ruediger & Kitchatinov 1997). However, the original CHMG93 model did not take into account the decay of the poloidal field in spite of the fact that its assumed extremely high value of the magnetic diffusivity implies the very short decay time of 35 Myr. With this e-folding time, the poloidal field will already have the amplitude of merely G by the age of 725 Myr. Fig. 10D shows contours of constant angular velocity in the D3 model of the solar age with this weak field. Because there is actually no poloidal field left, the solution looks almost identical to that for the pure D0 case.

We have chosen the D2 configuration because it more effectively transfers the wind torque and rotational shear into the radiative core than the D3 configuration. However, this same property plays a negative role when we try to prevent the interfacial differential rotation in the present-day Sun from penetrating the radiative interior. Fortunately, our D2H model combines the D2 geometry with the assumption of strong horizontal turbulent diffusion which, as we have seen in Fig. 10B, can potentially counteract the negative effect produced by the non-vanishing radial component of the D2 poloidal field at the core-envelope interface. To test this possibility, we have computed four stationary solutions for our favorite D2H model. Fig. 11A shows contours of constant angular velocity at the solar age for a D2H model in which neither the anisotropic diffusion coefficients nor the poloidal field amplitude were allowed to change with time. This is the worst possible solution in terms of the tachocline width. However, like the original CHMG93 solution, our first D2H test solution is not fully consistent because it does not take into account the poloidal field decay. For our assumed values of cm s and cm s (for the fastest solar-type rotator with ) equation (13) gives Gyr. Hence, by the solar age, the initial amplitude G will be reduced to G. Fig. 11B demonstrates that for a D2 poloidal field of such an order of magnitude the tachocline can be very thin, due to the continuing strong horizontal turbulent diffusion. Now, given our assumption about the dependence of on the angular velocity (equation 15), it would be fair to construct a solution for the case of having been reduced to cm s, proportionally to the decrease of from to . Its respective contours of constant are presented in Fig. 11C. We see that even in this case the obtained tachocline width is not more discrepant than in the original CHMG93 solution (Fig. 10C). Finally, Fig. 11D illustrates that the previous solution is greatly improved when the amplitude of the residual poloidal field is reduced by another order of magnitude. So, the anisotropic turbulence with the horizontal components of turbulent viscosity strongly dominating over those in the vertical direction plus the account of the poloidal field decay incorporated into our D2H model make it possible to obtain a sufficiently thin tachocline in the present-day Sun’s model even for the D2 poloidal field configuration.

## 4 Conclusion

The main goal of this work has been to modify the model of magnetic braking of solar rotation presented and discussed in detail by CHMG93 so that it could reproduce the most recent observational data concerning the rotation of solar analogs in open clusters of different ages, without coming into conflict with other available observational constraints. Although we have used a different prescription for the surface angular momentum loss rate (equation 1), its corresponding dependence of the envelope braking rate (yr) on the rotation rate does not deviate significantly from the dependence predicted by the Weber-Davis MHD wind model that was employed by CHMG93 (for comparison, we used Fig. 1 from Charbonneau 1992). It turns out that even in its original formulation the CHMG93 model experiences difficulties in fitting the upper 90th percentiles of the distributions for cluster solar-type stars. The discrepancy with the observationally constrained rotational evolution of the fastest rotators becomes much worse when we reduce the kinematic viscosity by one order of magnitude and/or take into account the poloidal field decay caused by the high magnetic diffusivity in the CHMG93 model. The viscosity reduction makes the model more consistent with the data for the atmospheric Li depletion in solar twins. The value of cm s used by CHMG93 leads to excessive Li destruction unless it describes, in terms of an effective diffusion coefficient, another angular momentum transport mechanism implemented as a process without a radial mass transfer, e.g. angular momentum redistribution by internal gravity waves. However, in this case the question would arise as to why that additional process could not do the entire work alone? Poloidal field decay also seems to be necessary, given that the model already follows the decay of the induced toroidal field. After these reasonable modifications are made, the CHMG93 model produces a rotational evolution of a rapidly rotating solar-type star that differs greatly from that of the fastest rotators among solar analogs in open clusters. Besides, at the solar age, the model still possesses strong differential rotation in the radiative core, in conflict with helioseismology measurements. This differential rotation has such a geometry that its corresponding contours of constant angular velocity form consecutive layers in cross sections of an axisymmetric torus, a region called the “dead zone” by Mestel & Weiss (1987). Its shape is a manifestation of Ferraro’s law of isorotation, , with the poloidal field having a dipole structure.

To save the CHMG93 model, we have assumed that the rotation-driven turbulence that is thought to amplify the viscosity and magnetic diffusivity compared to their microscopic values in stellar radiative zones is actually anisotropic with the horizontal components of the transport coefficients strongly dominating over those in the vertical direction. This is not a new hypothesis. It was first introduced to the stellar astrophysics community by Zahn (1992) and, since then, it has widely been used to model both chemical mixing and angular momentum transport in stars. The strong horizontal turbulent diffusion serves a two-fold goal in our modified model: it erases the dead zone along isobaric surfaces (a horizontal erosion of latitudinal differential rotation), while allowing us to choose a sufficiently small value for the coefficient of vertical diffusion that does not lead to a conflict with the Li data in solar twins. We have also found it necessary to switch from the D3 poloidal field configuration, that was considered to be the most appropriate one in the original CHMG93 model, to the D2 configuration. In the D2 geometry, the dipole poloidal field partially penetrates the convective envelope, which allows it to more effectively transfer the surface wind torque to the radiative interior. As a result, a thinner dead zone develops, which makes it easier for the horizontal turbulence to erase it. The horizontal component of the turbulent viscosity gives us an additional degree of freedom whose assumed linear dependence on the initial (ZAMS) angular velocity enables the models to properly reproduce the transition from the convex morphology of the surface rotation evolution of the fastest rotators to the concave curves for the slowest rotators, as suggested by observations.

MacGregor & Charbonneau (1999) have chosen the D3 geometry as the most suitable one for the modeling of magnetic breaking of solar rotation because its corresponding poloidal field configuration does not lead to a penetration of the differential rotation of the Sun’s convective envelope into its radiative interior, consistent with the observations. The helioseismology measurements show that the envelope differential rotation gets smoothed out in a thin () transition layer (the solar tachocline) immediately beneath the bottom of the Sun’s convective envelope and that the radiative core below the tachocline rotates like a solid body at least down to . Contrary to this, it has been shown that the non-vanishing radial component of the D2 poloidal field at the core-envelope interface causes the interfacial latitudinal differential rotation to move into the bulk of the radiative core. However, after our modifications of the original CHMG93 model, the new D2H model appears to be free of this problem. Both the strong horizontal turbulent diffusion and the poloidal field decay work towards confining the width of the solar tachocline.

It is important to note that our D2H model is very simplistic and it is based on a number of assumptions whose validity has not been proven. For instance, it is not clear how the rotation-driven turbulence should interact with the strong large-scale magnetic fields. Also, the formation of the solar tachocline is a much more complex process than described by our model. In particular, it must take into account a penetration of the meridional circulation from the convective zone into the radiative core (e.g., Gough & McIntyre 1998; Ruediger & Kitchatinov 2007; Gough 2009), which our model completely ignores. We are aware of these shortcomings. Therefore, we consider our model only as a legitimate extension of the CHMG93 model in the sense that we have not made any modifications of the latter that are very different from the original assumptions and approximations used by CHMG93. On the other hand, we have shown that just one additional assumption, that of strong horizontal viscosity, helps to solve several problems with the old model.

### References

- Andronov, N., Pinsonneault, M. H., & Sills, A. 2003, ApJ, 582, 358
- Barnes, G., Charbonneau, P., & MacGregor, K. B. 1999, ApJ, 511, 466
- Basu, A., & Antia, H. M. 2003, ApJ, 585, 553
- Chaboyer, B., Demarque, P., & Pinsonneault, M. H. 1995, ApJ, 441, 865
- Charbonneau, P. 1992, in Cool stars, stellar systems, and the sun, Proceedings of the 7th Cambridge Workshop, ASP Conference Series (ASP: San Francisco), 26, 416
- Charbonneau, P., & MacGregor, K. B. 1992, ApJ, 387, 639
- Charbonneau, P., & MacGregor, K. B. 1993, ApJ, 417, 762
- Charbonneau, P., Dikpati, M., & Gilman, P. A. 1999, ApJ, 526, 523
- Charbonnel, C., & Talon, S. 2005, Science, 309, 2189
- Couvidat, S., García, R. A., Turck-Chièze, Corbard, T., Henney, C. J., & Jiménez-Reyes, S. 2003, ApJ, 597, L77
- Denissenkov, P. A. 2005, ApJ, 622, 1058
- Denissenkov, P. A., Pinsonneault, M., & MacGregor, K. B. 2008, ApJ, 684, 757
- Denissenkov, P. A., Pinsonneault, M., Terndrup, D. M., & Newsham, G. 2009, arXiv:0911.1121v1 [astro-ph.SR]
- Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9
- Ferraro, V. C. A. 1937, MNRAS, 97, 458
- Gough, D. O., & McIntyre, M. E. 1998, Nature, 394, 755
- Gough, D. O. 2009, arXiv:0905.4924v1 [astro-ph.SR]
- Irwin, J., Hodgkin, S., Aigrain, S., Hebb, L., Bouvier, J., Clarke, C., Moraux, E., & Bramich, D. M., 2007, MNRAS, 377, 741
- Kawaler, S. D. 1988, ApJ, 333, 236
- Kitchatinov, L. L., & Ruediger, G. 1993, A&A, 276, 96
- Krishamurthi, A., Pinsonneault, M. H., Barnes, S., & Sofia, S. 1997, ApJ, 480, 303
- MacGregor, K. B. 1991, in Angular Momentum Evolution of Young Stars, ed. S. Catalano, & J. R. Stauffer, (Dordrecht: Kluwer Academic Publishers), 315
- MacGregor, K. B., & Charbonneau, P. 1999, ApJ, 519, 911
- Maeder, A. 2003, A&A, 399, 263
- Mathis, S., Palacios, A., & Zahn, J.-P. 2007, A&A, 462, 1063
- Matt, S., & Pudritz, R. E. 2005, ApJ, 632, L135
- Meléndez, J., Ramírez, I., Casagrande, L., Asplund, M., Gustafsson, B., Yong, D., do Nascimento, Jr., J. D., Castro, M., & Bazot, M. 2009, arXiv:0910.5845v1 [astro-ph.SR]
- Mestel, L., & Weiss, N. O. 1987, MNRAS, 226, 123
- Meynet, G., & Maeder, A. 2000, A&A, 361, 101
- Palacios, A., Charbonnel, C., Talon, S., & Forestini, M. 2003, A&A, 399, 603
- Palacios, A., Charbonnel, C., Talon, S., & Siess, L. 2006, A&A, 453, 261
- Plumb, R. A. 1977, J. Atmos. Sci., 34, 1847
- Randich, S. 2008, Mem. S. A. It., 79, 516
- Rebull, L. M., Wolff, S. C., & Strom, S. E. 2004, AJ, 127, 1029
- Richer, J., Michaud, G., & Turcotte, S. 2000, ApJ, 529, 338
- Ringot, O. 1998, A&A, 335, L89
- Rogers, T. M., & Glatzmaier, G. A. 2005, MNRAS, 364, 1135
- Rogers, T. M., & MacGregor, K. B. 2009, 2009MNRAS.tmp.1534R
- Rogers, T. M., MacGregor, K. B., & Glatzmaier, G. A. 2008, MNRAS, 387, 616
- Ruediger, G., & Kitchatinov, L. L. 1996, ApJ, 466, 1078
- Ruediger, G., & Kitchatinov, L. L. 1997, AN, 318, 273
- Ruediger, G., & Kitchatinov, L. L. 2007, NJPh, 9, 302
- Schatzman, E. 1962, AnAp, 25, 18
- Shu, F., Najita, J., Ostriker, E., Wilkin, F., Ruden, S., & Lizano, S. 1994, ApJ, 429, 781
- Skumanich, A. 1972, ApJ, 171, 565
- Spiegel, E. A., & Zahn, J.-P. 1992, A&A, 265, 106
- Spruit, H. C. 1999, A&A, 349, 189
- Spruit, H. C. 2002, A&A, 381, 923
- Talon, S., & Zahn, J.-P. 1997, A&A, 317, 749
- Talon, S., Zahn, J.-P., Maeder, A., & Meynet, G. 1997, A&A, 322, 209
- Tassoul, J.-L. 2000, Stellar Rotation/Jean-Louis Tassoul. Cambridge ; New York : Cambridge University Press, 2000. (Cambridge astrophysics series ; 36)
- Tomczyk, S., Schou, J., & Thompson, M. J. 1995, ApJ, 448, L57
- Weber, E. J., & Davis, L., Jr. 1967, ApJ, 148, 217
- Zahn, J.-P. 1992, A&A, 256, 115