Structure and stability of the magnetic solar tachocline

Rather weak fossil magnetic fields in the radiative core can produce the solar tachocline if the field is poloidal and almost horizontal in the tachocline region, i.e. if the field is confined within the core. This particular field geometry is shown to result from a shallow (<~1Mm) penetration of the meridional flow existing in the convection zone into the radiative core. I.e., two conditions are crucial for a magnetic tachocline theory: (i) the presence of meridional flow of a few meters per second at the base of the convection zone, and (ii) a magnetic diffusivity inside the tachocline smaller than 10^8 cm^2s^-1. Numerical solutions for both confined poloidal fields and the resulting tachocline structures are presented. We find that the tachocline thickness runs asB_p^-1/2 with the poloidal field amplitude falling below 5% of the solar radius forB^p>5 mG. The resulting toroidal field amplitude inside the tachocline of about 100 G does not depend on theB^p. The hydromagnetic stability of the tachocline is only briefly discussed. For the hydrodynamic stability of latitudinal differential rotation we found that the critical 29% of the 2D theory of Watson (1981) are reduced to only 21% in 3D for marginal modes of about 6 Mm radial scale.


Introduction
The tachocline is a thin shell inside the Sun where the rotation pattern changes strongly. Beneath the tachocline, the solar rotation is rather uniform. Above the tachocline, the rotation rate varies with latitude to decrease from equator to pole by about 30%. The transition from differential to rigid rotation detected by helioseismology (Wilson et al 1997;Kosovichev et al 1997;Schou et al 1998) is shown in figure 1. The tachocline thickness is about 5% of the solar radius, its midpoint radius is (0.692 ± 0.005)R ⊙ , and it is slightly prolate in shape (Kosovichev 1996;Antia et al 1998;Charbonneau et al 1999b). The tachocline is located mostly if not totally beneath the base of convection zone at R in = 0.713R ⊙ (Christensen-Dalsgaard et al 1991; Basu & Antia 1997).
The given picture of the internal solar rotation suggests that some coupling must exist between the base of convection zone and the inner core. Otherwise, the rotational braking of the Sun on evolutionary time-scales would produce a rapidly rotating core (Dicke 1970). More than that, an even stronger link between low and high latitudes operates immediately below the base of convection zone in order to produce the rotation profile within the tachocline. The magnetic theory which explains the tachocline structure as a consequence of a weak internal fossil magnetic field simultaneously provides both the links (Rüdiger & Kitchatinov 1997;MacGregor & Charbonneau 1999). The present paper summarizes the basic ideas of the magnetic theory and describes new results of numerical MHD to model the tachocline and probing the internal field structure. Other approaches has been recently reviewed by Gilman (2005) and Garaud (2007).
Section 3 formulates the equations of the tachocline model and shows their solution for prescribed geometry of the fossil magnetic field within the radiative core. The role of meridional flow penetrating the uppermost region of the radiative core from the convection zone is discussed in section 4. In this section, the eigenmodes of internal poloidal field are computed with inclusion of the penetrating flow. Tachocline models with consistently defined internal field modes are discussed in section 5. In the final section 6 a first step is presented to understand the stability problem of the solar tachocline. In a hydrodynamical approach the 2D Watson theory for the stability of latitudinal differential rotation is reformulated in 3D with surprising results. The modes with small radial wavelengts are most unstable so that it is indeed essential for tachocline theories to take into account all the radial stratifications.

An heuristic approach
The basic parameters of the tachocline theory can be estimated as follows. The estimates are made with a simplified shear flow model in Cartesian geometry sketched in figure 2. The plane of z = 0 mimics the bottom of the convection zone where a shear flow u = (U(y, z), 0, 0) imitating the differential rotation is prescribed as (1) The Cartesian coordinates x, y and z correspond to azimuth, latitude and depth beneath the convection zone. An uniform (poloidal) field B 0 along the y-axis is prescribed in the radiative zone. The shear flow penetrates the z > 0 region and produces an x-component B(y, z) of the field so that can be written. The equations for the stationary (toroidal) x-components of both magnetic and velocity fields then read where ν is the viscosity and η the magnetic diffusivity. The boundary condition for the flow is (1). For the toroidal field we impose the vacuum condition B = 0 at z = 0 motivated by the very large turbulent magnetic diffusivity inside the convection zone compared with microscopic diffusivity of the radiative interior. The remaining two conditions require U and B to vanish for z → ∞. The solution of the equations (3) simply reads U(y, z) = U 0 exp (−λ 1 kz) cos (λ 2 kz) sin ky, where Pm = ν/η is the magnetic Prandtl number and the parameters λ 1 = 1 + Ha 2 1/4 cos 1 2 arctanHa , λ 2 = 1 + Ha 2 1/4 sin 1 2 arctanHa (5) depend on the Hartmann number which is the basic parameter of the whole theory. With Ha = 0, the equation (5) gives λ 1 = 1, λ 2 = 0, and (4) converts into the nonmagnetic solution, B = 0, U = U 0 e −kz sin ky. The resulting shear flow penetrates deep inside the z > 0 region. No slender tachocline can be formed without magnetic fields. But with solar parameters one finds Ha ≃ 10 7 B 0 (B 0 in G) just below the convection zone. Very large Ha can thus be expected. For this case, λ 1 = λ 2 = Ha/2, and the solution (4) provides two important consequences: (i) The tachocline thickness, D tach , is strongly reduced compared to its nonmagnetic value, D 0 , The tachocline is so thin because its extension in z-direction is only due to viscous stress while the smoothing in y-direction is performed by the much stronger Maxwell stress.
(ii) The amplitude of the toroidal magnetic field in the steady tachocline does not depend on the poloidal field strength, i.e.
The ratio of magnetic to kinetic energy in the tachocline equals Pm. In terms of the Alfvén velocity, predicting B ∼ 1000 G for the Sun. After (7) even a weak poloidal field of only B 0 ∼ 10 −3 G can reduce D tach below 5% of the solar radius. Equation (7) does also apply when B 0 is not an internal field of the radiative core. In case the poloidal field of the solar cycle diffuses into the core (Forgács-Dajka & Petrovay 2002) ν and η must be replaced by the turbulent diffusivity values. As their magnetic Prandtl number is not much smaller than unity the resulting toroidal magnetic field after (9) becomes much stronger than 1000 G.

Tachocline model in spherical geometry
The tachocline equations will be formulated for axial symmetry. Then the magnetic field can be expressed in terms of A and B in accord to where r, θ and φ are spherical coordinates and e φ is the azimuthal unit vector. In the next section we show that meridional flow is only significant for the structure of the poloidal field. It can be neglected in the tachocline equations below. Here the only mean flow is the differential rotational with Ω = Ω(r, θ).
If the poloidal field is given, the tachocline is described by the equations for the toroidal field B and the angular velocity Ω η r ∂ ∂θ  (14) and viscosity (15) calculated with the solar structure model by Stix & Skaley (1990).
The rotation profile at the top boundary is produced by the convection zone for which the expression Ω = 2.9 1 − 0.15 cos 2 θ µrad s −1 at r = R in , is used derived by Charbonneau et al (1999a). The remaining boundary conditions are the vacuum condition for toroidal field on the top and the conditions for the regularity of the fields for r = 0: The used diffusivity profiles are shown in figure 3. We do not expect turbulence in the radiative core and apply, therefore, the microscopic magnetic diffusivity and the viscosity (cf. Kippenhahn & Weigert 1994) including molecular (ν micro ) and radiative (ν rad ) parts; κ is opacity. The equations (11) with their boundary conditions (12), (13) provide tachocline solutions if the poloidal field is known. For a very first view the field model is used to probe whether slender tachoclines can indeed be found. B p is the free amplitude of the poloidal field. Figure 4 compares the rotation pattern in a radiative zone computed with either zero magnetic field and with B P = 0.01 G. B P is the maximum strength of poloidal field inside the core. The maximum is attained in the core center. The field in the tachocline region is much weaker. Even this rather weak field suffices to produce a very strict and slender tachocline.

Magnetic field confinement by meridional flow
Our tachocline model is not yet consistent. It works with a prescribed poloidal field. We have shown that the model is not very sensitive to the poloidal field amplitude. Even a weak field can produce tachoclinic structures. The model is, however, very sensitive to the field geometry. The tachocline can be found only if the field lines are almost horizontal near the top of radiative zone. Steady rotation and magnetic fields can be realized in highly conducting fluids only with constant angular velocity along the field lines (Ferraro 1937). Therefore, the open field geometry with field lines crossing the boundary between convection and radiative zones is not appropriate for tachocline formation. The poloidal field prescribed by (16) is of a closed type. The field can, however, change to an open structure due to magnetic diffusion (Brun & Zahn 2006). We shall show in the following that already a small penetration of meridional flow from the convection zone into the radiative core produces the confined field geometry required for the tachocline formation (see Kitchatinov & Rüdiger 2006). If the electric conductivity in the core is high enough compared to that of the convection zone then the flow has massive consequences and the poloidal field lines in the transition zone become parallel to the meridional flow (cf. Mestel 1999).

Penetration of meridional circulation
A global poleward flow is observed at the solar surface (Komm et al 1993) persisting to a depth of at least 12 Mm (Zhao & Kosovichev 2004). There must be a return flow towards the equator somewhere deeper. Theoretical models indeed predict an equatorward flow of ∼ 10 m s −1 at the base of convection zone (Kitchatinov & Rüdiger 1999;Miesch et al 2000;Rempel 2005). This flow can penetrate beneath the bottom of the convection zone into the radiative core (figure 5). This penetration has been discussed recently in relation to dynamo models for solar activity (Nandy & Choudhuri 2002;Gilman & Miesch 2004;Rüdiger et al 2005). The penetration results from the viscous drag imposed by the meridional flow at the base of convection zone on the fluid beneath, and this is opposed by Coriolis force. The thickness of such an Ekman layer results as D pen ∼ ν/Ω (see Gilman & Miesch 2004). The plot in the left of figure 5 is indeed approximated by D pen ≃ 2.3 · ν/Ω or, equivalently, by for solar parameters ‡. The penetration distance is about 100 m for microscopic viscosity and remains shorter than 1 Mm for any reasonable value of eddy viscosity. The distance is small even compared to the tachocline depth, i.e. D pen ≪ D tach . This shallow penetration is, nevertheless, highly important for the geometry of the internal poloidal field. The ratio of the time τ diff = D 2 pen /η for the magnetic field diffusion across the penetration layer to the characteristic time of induction of a meridional magnetic field from a radial one, τ s = D pen /u m , gives the magnetic Reynolds ‡ D pen of figure 5 was defined as distance from the convection zone bottom to the position where the meridional flow falls to zero (with u m ≃ 10 m s −1 for the meridional velocity). This Reynolds number is large, Rm∼ 10 3 , for microscopic diffusivities but it remains above this limit with eddy diffusivities up to 10 6 cm 2 s −1 . The large ratio of (18) means that the latitudinal field inside the penetration layer is large compared to the radial field component. With other words, the field has just the confined geometry required for the tachocline formation.
Only with eddy diffusivities of 10 12 cm 2 s −1 the Reynolds number (18) sinks to unity, so that only in this case the influence of penetration on the internal field geometry is too weak.
The ratio of τ diff to the advection time is small independently of whether microscopic or eddy diffusivity is used. If the tachocline is stable in the hydrodynamic regime (see section 6) the microscopic diffusivities should be used. The small ratio (19) justifies the neglect of the meridional flow in the tachocline equations (11). The diffusion time, τ diff , is too short as could the penetration layer be dynamo-relevant. The configuration of the internal poloidal field is probably the only process for which the penetration of the meridional circulation into the stable radiative zone is important.

A model of the poloidal field
For both axisymmetric flows and fields the poloidal magnetic field equation decouples from the tachocline equations (11) and reads where u is the velocity field. The vacuum boundary condition for the poloidal field can be applied on the top of the radiative zone because of the very large turbulent diffusivity in the convection zone, The condition is usually formulated in terms of a Legendre polynomials expansion nA n (r) P 1 n (cos θ) .
The other condition is A = 0 at r = 0.
The meridional flow in the bulk of the radiative core is very slow. The characteristic time of the Eddington-Sweet circulation even exceeds the solar age (Tassoul 2000). Only the flow penetrating from convection zone is significant for the tachocline. The stream function, ψ(r, θ), of the penetrating flow drops to zero at a small depth of order D pen inside the core. The function is, however, finite at the top boundary of the core where it scales as whereψ is dimensionless function of order unity. u m is the meridional velocity amplitude at the boundary. The penetration depth is so small even compared to tachocline thickness that we are motivated to integrate equation (20) across the penetration layer instead of resolving this layer explicitely. Such an integration reformulates the top boundary condition which now reads where Rm is the magnetic Reynolds number (18) and the RHS is defined in (22). For Rm = 0 (25) reduces to the vacuum condition (21). The true Reynolds number is Rm ∼ 10 3 . Penetrating meridional flow is now included via the condition (25) and poloidal field structure can be found by solving eigenvalue problem for diffusion equation, where eigenvalue τ is the time of resistive decay of normal modes. The results discussed below were obtained in computations with the simplest stream function, whereP 1 2 is the normalized Legendre polynomial. The penetrating flow is expected to confine the internal poloidal field within the radiative core. The parameter is used to estimate the field confinement. The δφ-parameter (28) estimates the ratio of the magnetic flux through the surface of the core to the characteristic value of the flux within the core. The linear equations (26) and (25) define the field structure but not its amplitude. The amplitude B p of the poloidal field is a free parameter of the model.

Normal modes of the internal field
How the structure of the internal field changes towards a confined geometry as Rm increases is shown in figure 6. The internal field computed with Rm = 0 has an open structure. Even a moderate flow with Rm = 10 changes the field considerably towards the confined geometry. The internal field for the solar value Rm = 1000 shown in figure 6 is almost totally confined. Less than 1% of magnetic flux belongs to 'open' field lines in this case (figure 7). The figures 6 and 7 only display the most slowly  decaying dipolar modes. Other normal modes have shorter scales in radius or latitude and also shorter decay times. The three modes following the longest-living dipole ordered for decreasing lifetimes are shown in figure 8. The higher-order modes also have a confined structure. The decay times of the given normal modes are long enough to allow tachocline computations with steady poloidal fields.
For sufficiently high Rm the direction of the penetrating meridional flow does not play a role. The flow from equator to poles also makes efficient confinements of the internal field. The effect can be understood as a magnetic field expulsion from the region  of circulating motion (Weiss 1966). The transition from convection zone to radiative core was treated as a sharp boundary in our model. A smooth decrease of eddy diffusivity in convection zone towards its base may also contribute to the field confinement via the diamagnetic effect of inhomogeneous turbulence (cf. Krause & Rädler 1980). Figure 9 shows the angular velocity distributions inside the core computed with dipolar eigenmodes of internal field of various strengths. Dependence of the tachocline thickness on poloidal field amplitude is shown in figure 10 (D tach was defined as depth of exponential decrease of equator-to-pole difference in angular velocity). Even a very weak field, B p ∼ 10 −3 G, can produce tachocline. The depth D tach drops below 5% of the solar radius for B p > 5 mG.

Tachocline models
The tachocline thickness in relation to the amplitude to the prescribed fossil field is given in figure 10. It is very close to the rough estimate (7). The simulations also confirm the above finding that the toroidal field amplitude, B t = 100...200 G, hardly varies while the poloidal field amplitude changes by several orders of magnitude. Also the higher order modes given in figure 8 also produce tachoclinic structures provided that the fields have the confined structure which is always the case for Rm= 100 or larger. Figures 9 and 10 are valid for Rm = 1000. Already Rm = 100 suffices for tachocline formation, but for smaller Rm the poloidal field remains 'too open' for the formation of a tachocline. With other words, a meridional flow with minimum 1 m s −1 amplitude at the base of convection zone is necessary for the success of the magnetic tachocline theory. The rather shallow penetration of the meridional flow of the convection zone into the radiative core influences the internal field geometry strongly enough and in such a way that the field becomes appropriate for the tachocline formation. We suggest that a layer thiner than 1 Mm beneath the convection zone where a meridional flow of 1-10 m s −1 enters the radiative core is responsible for the existence of the solar tachocline.

Hydrodynamic stability: the Watson approach in 3D
Latitudinal differential rotation can be unstable even without magnetic field if the shear ∂Ω/∂θ is positive and sufficiently strong (Watson 1981). The critical value of 29% latitudinal shear found by Watson for the nonaxisymmetric mode with m = 1 resulted from a theory with strong stratification but without any radial velocity and radial wavelengths. The value has also appeared in a 3D numerical examination of marginal stability of a shell rotating fast enough with the rotation law but of fully incompressible material (Arlt et al 2007). The critical shear increases to higher values, however, if the real rotation law (including its radial variations) of the solar tachocline is adopted. In this case the superrotation observed in the equatorial region of the rotation law strongly stabilizes the shear instability. In this section the Watson approach with strong density stratification (for m = 1) is extended to 3D, i.e. to the inclusion of finite radial velocities and finite radial Figure 11. The buoyancy frequency (30) in the upper radiative core after the solar model of Stix & Skaley (1990). wavelengths. The stabilizing effect of a subadiabatic stratification on the mean flow is characterized by the buoyancy frequency, where S = C v ln (P/ρ γ ) is the specific entropy of ideal gas. The frequency is very large in the solar radiative core (see figure 11). The larger the N the more the radial fluctuations are suppressed by the buoyancy force. Radial velocities should therefore be small. Our stability analysis is local in the radial dimension, i.e. we use Fourier modes exp ikr with kr ≫ 1. The equation system remains, however, global in horizontal dimensions. Instabilities of rotating shear flows are fast. Their growth rates are of order Ω or smaller. In the short-wave approximation the velocity field can be assumed as divergence-free, divu ′ = 0. The next assumption concerns the pressure. Local thermal disturbances occur at constant pressure so that ρ ′ /ρ = −T ′ /T or S ′ = −C p ρ ′ /ρ where perturbations have been marked by dashes. This assumption is again justified by the incompressible nature of the fluid. We start from the linearized equations of velocity, i.e.
and entropy The basic flow is the rotation law (29). Perturbations of velocity and pressure are now expressed in terms of scalar potentials like with the operator The identities are used to reformulate the relations in terms of the potentials. The perturbations are considered as Fourier modes for time, azimuth and radius in the form exp (i(−ωt + mφ + kr)). For an instability the eigenvalue ω must possess a positive imaginary part. Only the highest order terms in kr for the same variable have been considered. The pressure term in (31) can be formulated as In order to denormalize the variables time is measured in units of Ω −1 0 and velocities are scaled with rΩ 0 . The remaining variables are Then the equation for the poloidal flow readŝ with µ = cos θ and the normalized radial wavelengtĥ The first term on the RHS describes the stabilizing effect of the stratification. It vanishes for smallλ. Apart from this stabilizing buoyancy term, the wavelength only exists in diffusive terms. The second term of the RHS includes the action of viscosity. Here we have used the quantities for the diffusive parameters ν and χ. Now the limit N → 0 can no longer be used for the transition to nonstratified but viscous media. The numerical values for the solar tachocline are The second and the following lines in (38) describe the influences of the basic rotation. Note that only latitudinal derivatives of Ω appear.
The complete system of equations also includes the equation for toroidal flow, and the entropy equation Note that for the Reynolds number Re = Ωr 2 /ν follows which with (41) is very large. For positive a in (29) the modes A1 and S2 become unstable for sufficiently large a = a crit . Figure 12 shows the dependence of the a crit on the normalized wavelengthλ. Of course, for large enough radial wavelength the 29%-value of the Watson theory with 2D fluctuations is reproduced. It is reduced, however, to a = 0.21 in the full 3D model calculations. We see that in contrast to the original Watson approach short rather than long radial scales are preferred. The minimum a crit appears forλ ≃ 0.6, so that the characteristic wavelength of λ ≃ 6 Mm results for the solar tachocline. For shorter scales the instability for m = 1 disappears due to the action of the dissipation processes but modes with higher m remain still unstable.

Concluding remarks
The presence of a meridional flow with an amplitude of a few meters per second at the bottom of the convection zone is shown as necessary for the magnetic tachocline model. The flow provides the confined geometry (with field lines parallel to the outer spherical boundary) of the internal magnetic field which is necessary for the tachocline formation. The bottom flow is also a key ingredient of advection-dominated dynamo models for the solar cycle (see Rüdiger & Hollerbach 2004 for detailed references). The flow was predicted theoretically but its existence is not yet confirmed by observations. We hope, of course, that helioseismology will probe the deep meridional flow soon.
The polar cusps in the tachoclines of figure 9 are consequences of the assumed axisymmetry. This assumption may not be realistic. The penetrating meridional flow is, however, expected to confine even nonaxisymmetric fields to the core so that the fields should also be able to produce the tachoclinic structures.
The tachocline computations result in toroidal fields of 100 to 200 G, much stronger than the original poloidal fields. These toroidal fields can be subject to currentdriven pinch-type instabilities (Tayler 1973). Hence, the question is whether the Figure 12. Stability map for latitudinal differential rotation (3D hydrodynamics). Most unstable are the perturbation modes with the vertical scaleλ ≃ 0.6. The critical magnitude of latitudinal shear is reduced to 0.21 compared to the 0.29 value for large wavelengths. For short wavelengths the instability for m = 1 rapidly disappears by dissipation effects. steady axisymmetric solutions of our model are stable against all the nonaxisymmetric disturbances. The threshold field strength for the instability were estimated to be of order 10 2...3 G (Spruit 1999;Arlt et al 2007;Kitchatinov & Rüdiger 2007) what is (only) somewhat larger than the field amplitudes in our model.
It is important for the magnetic tachocline model that the magnetic Reynolds number (18) is large. This is true for microscopic diffusivity or for not too large eddy diffusivities up to about 10 8 cm 2 s −1 . The tachocline should, therefore, be stable or only mildly turbulent to allow poloidal field confinement by meridional flow. Some low level of turbulence may, however, be even necessary. The magnetic field is so efficient in producing tachocline structures that poloidal fields of only 1 G lead to radial scales smaller than 1% of the solar radius. It is, however, D tach ≃ 0.04R ⊙ . Hence, the internal poloidal field must be (much) smaller than 1 G (Kitchatinov et al 2001) or, if it is not, some kind of instability may prevent the tachocline thickness to reduce below the helioseismologically detected level.
The stability issue is a clear perspective for future tachocline studies. As a corresponding application we have developed in the last section the 2D theory of Watson (1981) for the hydrodynamic instability of latitudinal differential rotation to a fully 3D theory. Now the critical wavelength of the unstable mode with m = 1 is only 6 Mm while the critical latitudinal shear (see (29)) is reduced from 29% to 21%. This (linear) theory works with a Prandtl number of order 10 −6 and a very high Reynolds number of order 10 12 .