The Solar Tachocline:
A Self-Consistent Model
of Magnetic Confinement
Toby Wood
Queens’ College, Cambridge
A Dissertation submitted to
the University of Cambridge for
the degree of Doctor of Philosophy
April 2010
“There is no easy way from the earth to the stars”
Lucius Annaeus Seneca
Acknowledgements
First and foremost, I would like to thank my supervisor, Prof. M. E. McIn-
tyre, for his continuing guidance and infectious enthusiasm for science. I have
also benefited enormously from discussions with members of my research group,
including Prof. D. O. Gough, Dr. G. I. Ogilvie, Prof. J. C. B. Papaloizou,
Prof. M. R. E. Proctor and Prof. N. O. Weiss. I am indebted to Prof. J. Christensen-
Dalsgaard for providing the one-dimensional profiles of the solar structure used
in chapter 5.
I never would have got this far without the support of my family, friends and
colleagues, whose advice at every stage has been invaluable. I am particularly in-
debted to Mark Rosin, for suggesting the possibility of a life beyond mathematics,
and to Ce´line Guervilly, for proving him right.
Declaration
This dissertation is based on research done at the Department of Applied Math-
ematics and Theoretical Physics from October 2006 to February 2010. The ma-
terial presented is original except where reference is made to the work of others
in the text. Much of Parts III and IV are based upon research work performed in
collaboration with my supervisor Professor Michael E. McIntyre, and published
as Wood & McIntyre (2007) and Wood & McIntyre (2010).
TOBY WOOD
Cambridge,
November 15, 2010
iv
Abstract
In this dissertation we consider the dynamics of the solar interior, with particular
focus on angular momentum balance and magnetic field confinement within the
tachocline.
In Part I we review current knowledge of the Sun’s rotation. We summarise
the main mechanisms by which angular momentum is transported within the
Sun, and discuss the difficulties in reconciling the observed uniform rotation of
the radiative interior with purely hydrodynamical theories. Following Gough &
McIntyre (1998) we conclude that a global-scale interior magnetic field provides
the most plausible explanation for the observed uniform rotation, provided that
it is confined within the tachocline.
We discuss potential mechanisms for magnetic field confinement, assuming that
the field has a roughly axial-dipolar structure. In particular, we argue that the
field is confined, in high latitudes, by a laminar downwelling flow driven by tur-
bulence in the tachocline and convection zone above.
In Part II we describe how the magnetic confinement picture is affected by the
presence of compositional stratification in the “helium settling layer” below the
convection zone. We use scaling arguments to estimate the rate at which the
settling layer forms, and verify our predictions with a simple numerical model.
We discuss the implications for lithium depletion in the convection zone.
In Part III we present numerical results showing how the Sun’s interior magnetic
field can be confined, in the polar regions, while maintaining uniform rotation
within the radiative envelope. These results come from solving the full, nonlinear
equations numerically. We also show how these results can be understood in
terms of a reduced, analytical model that is asymptotically valid in the parameter
regime of relevance to the solar tachocline.
In Part IV we discuss how our high-latitude model can be extended to a global
model of magnetic confinement within the tachocline.
vi
Contents
I Solar Dynamics 1
1 Basic Solar Physics 3
1.1 Solar structure and observations . . . . . . . . . . . . . . . . . . . 3
1.2 Rotating, stratified fluids . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Magnetised fluids . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Solar spin-down . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 The Sun’s Differential Rotation 15
2.1 Standard solar models . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Helioseismology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3 Differential rotation . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 The tachocline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.5 Dynamo theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3 The Magnetic Interior 35
3.1 The Sun’s fossil field . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 The Ferraro constraint . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 The magnetic confinement problem . . . . . . . . . . . . . . . . . 38
vii
viii Contents
4 Magnetic Field Confinement 43
4.1 Magnetic pumping . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Mean meridional circulations (MMCs) . . . . . . . . . . . . . . . 46
4.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
II Solar Composition 63
5 Helium Settling 65
5.1 The µ-choke . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.2 The settling equation . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.3 The formation of the helium settling layer . . . . . . . . . . . . . 71
5.4 The compositional gradient . . . . . . . . . . . . . . . . . . . . . 75
5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6 Depletion of Lithium and Beryllium 79
III The Magnetic Confinement Layer 83
7 The High-Latitude Tachocline 85
7.1 The structure of the magnetic field . . . . . . . . . . . . . . . . . 85
7.2 The structure of the tachocline . . . . . . . . . . . . . . . . . . . 86
7.3 Angular momentum balance . . . . . . . . . . . . . . . . . . . . . 87
7.4 The equations and parameters . . . . . . . . . . . . . . . . . . . . 89
7.5 The limit of infinite stratification . . . . . . . . . . . . . . . . . . 92
7.6 Finite stratification . . . . . . . . . . . . . . . . . . . . . . . . . . 96
CONTENTS ix
7.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
8 The Magnetic Confinement Layer 99
9 The Helium Sublayer 105
9.1 Sublayer scalings . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
9.2 The subtail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
10 The Numerical Model 111
10.1 Numerical complexity . . . . . . . . . . . . . . . . . . . . . . . . . 111
10.2 The numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . 113
10.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . 114
10.4 The numerical solutions . . . . . . . . . . . . . . . . . . . . . . . 118
10.5 Upper-domain “slipperiness” . . . . . . . . . . . . . . . . . . . . . 122
IV Conclusions 125
11 Discussion 127
12 Future Work 131
12.1 Gyroscopic pumping . . . . . . . . . . . . . . . . . . . . . . . . . 131
12.2 Lithium burning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
V Appendices 133
A Parameter Values 135
B Meridional Circulation 137
C MAC Waves 141
D The Magnetic Margules Front 145
E Magnetic Pumping 149
F The Tayler Instability 153
G The Tayler–Spruit Dynamo 157
H The Numerical Scheme 161
H.1 Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
H.2 Finite-differences . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
H.3 Magnetostrophic balance . . . . . . . . . . . . . . . . . . . . . . . 166
I Numerical Test Cases 169
I.1 Transient adjustment . . . . . . . . . . . . . . . . . . . . . . . . . 169
I.2 Magnetostrophic balance . . . . . . . . . . . . . . . . . . . . . . . 170
I.3 Numerical boundaries . . . . . . . . . . . . . . . . . . . . . . . . . 171
I.4 Horizontal resolution . . . . . . . . . . . . . . . . . . . . . . . . . 173
I.5 Horizontal diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . 174
x
Part I
Solar Dynamics
1
Chapter 1
Basic Solar Physics
1.1 Solar structure and observations
1.1.1 Thermodynamics
The Sun formed from the gravitational collapse of a molecular cloud approxi-
mately 5 billion years ago. This collapse was halted once the radial (i.e. outward)
pressure gradient became large enough to balance the gravitational attraction,
leading to a state of hydrostatic balance that has persisted throughout the Sun’s
main-sequence lifetime. Departures from hydrostatic balance are restored on the
dynamical timescale (≈ 20 minutes).
Within the Sun’s core, which extends to around 30% of the Sun’s radius1 R,
high pressure and temperature cause hydrogen nuclei to fuse into helium, re-
leasing the energy that powers the Sun (and all life on Earth). The thermal
energy released in the core is carried outward by photon radiation. Throughout
the solar interior, pressure, density and temperature all decrease outward with
radius (see figure 1.1), but the radial entropy gradient changes sign at around
0.7R. The Sun’s “super-adiabatic” outer layers are convectively unstable, and
heat is transported outward through these layers by fluid motions as well as by
photon radiation. The boundary between the convective and radiative zones is
1Throughout this dissertation, we will use R for spherical radius and r for cylindrical radius.
The value of R and other relevant parameters are listed in appendix A.
3
4 1. Basic Solar Physics
blurred slightly by the presence of dynamical instabilities and penetration by
“overshooting” convective plumes, but these processes are limited by the strong
stable stratification of the radiative envelope, as will be described in chapter 2.
Figure 1.1: The interior structure of the Sun, shaded according to temperature. The dashed
line indicates the boundary between the convection zone and the radiative envelope.
1.1.2 Rotation
The one-dimensional solar model just described ignores the effects of the Sun’s
rotation. From an energetic perspective this is a valid approximation, because
the total rotational energy of the Sun is only a tiny fraction of its thermal energy.
But the Sun’s rotation significantly affects the dynamics of its interior, for reasons
described in §1.2.
The Sun rotates at a rate of about once a month. The total angular momentum
of the Sun is in fact only a small fraction of the angular momentum that was
present in the molecular cloud from which the Sun formed. This cloud would have
rotated relative to its centre at a rate given roughly by “Oort’s second constant”2
— i.e. half the vorticity of the local galactic rotation. If angular momentum had
been perfectly conserved during the gravitational collapse of this gas cloud, the
“ballerina effect” would have caused the Sun to form with a rotation period of
2Approximately one rotation every 500 million years.
1.2 Rotating, stratified fluids 5
just a few seconds. In reality most of this angular momentum was lost during
the early stages of the Sun’s formation, but nevertheless the Sun probably had a
rotation period of less than one day when it arrived on the main sequence (e.g.
Mestel & Weiss, 1987, & refs). The Sun must therefore have continued to shed
angular momentum even during its main-sequence lifetime. We return to the
issue of “solar spin-down” in §1.4.
1.1.3 Magnetism
The Sun has a magnetic field, which can be observed where it extends out from
the solar surface. The field at the surface has an average strength of around one
gauss (slightly stronger than the Earth’s surface field), although field strengths
of a few thousand gauss exist in localised patches called sunspots. The field is
generated within the Sun’s convection zone, where the fluid motions act as a
dynamo (§2.5), and decays in strength algebraically at large distances from the
solar surface. The spin-down of the Sun is closely connected to the presence
of this external magnetic field, as will be described in §1.4. First, however, we
introduce some relevant fluid dynamics associated with rotating, stratified and
magnetised fluids.
1.2 Rotating, stratified fluids
In a rotating system, small perturbations to the lines of absolute vorticity3
propagate as inertial waves (also known as epicyclic or Coriolis waves), with
a wavespeed proportional to the rate of rotation. The dynamical importance of
rotation is generally quantified by the Rossby number Ro, which is a measure
of inertial forces, in the rotating frame, against Coriolis forces (e.g. Greenspan,
1968). If motions in the fluid have a characteristic timescale τ , say, and the
overall rotation rate of the system is Ω, then the Rossby number is defined as
Ro = 1/(2Ωτ). For a rapidly rotating (i.e. low Rossby number) system, inertial
wave propagation inhibits fluid motions that deform the lines of absolute vortic-
ity — this is sometimes known as “rotational stiffness”. This physical principle
3Here, “absolute” vorticity refers to the total vorticity in an inertial frame.
6 1. Basic Solar Physics
is crystallised within the Taylor–Proudman theorem: for a barotropic fluid in a
balance of Coriolis, pressure and gravitational forces, any motions perpendicular
to the rotation axis must be invariant in the direction parallel to the rotation axis.
Flow of this kind is called “geostrophic”. In a geostrophic fluid, any differential
rotation must be “constant on cylinders”.
In a stratified, non-barotropic fluid, surfaces of constant density can be tilted
with respect to surfaces of constant pressure. This “baroclinicity” is a source of
vorticity in the direction perpendicular to the gradients of pressure and density.
The presence of baroclinicity can overcome the Taylor–Proudman constraint. In
particular, baroclinic production of azimuthal vorticity can balance the azimuthal
vorticity produced by axial variations in the centripetal acceleration; this is known
as “thermal-wind balance” (e.g. Pedlosky, 1979). Adopting cylindrical polar co-
ordinates (r, φ, z) centred on the axis of rotation, this balance can be expressed
as
∂
∂z (Ω
2r) = 1ρ2 (∇p×∇ρ) · eφ , (1.1)
where Ω is the local rotation rate, p and ρ are the pressure and density fields,
and eφ is the unit-vector directed azimuthally. If the system is also in hydrostatic
balance, and if we assume that fractional pressure perturbations are small com-
pared with fractional density perturbations, then we may approximate ∇p ≈ ρg,
where g is the modified gravitational acceleration, i.e. the gradient of the total
gravito–centrifugal potential. Equation (1.1) now becomes
∂
∂z (Ω
2r) ≈ 1ρ(g ×∇ρ) · eφ . (1.2)
So any axial gradient in the angular velocity must be accompanied by a horizon-
tal4 density gradient. In the absence of such gradients, the system reverts to a
“Taylor–Proudman state”, with no axial variation in Ω. (A simple example of
thermal-wind balance is the “Margules front” described in appendix D.)
In an axisymmetric fluid, neither pressure nor gravity exert any torque about the
axis of symmetry. So if we apply an axisymmetric torque to a fluid in a steady
balance of Coriolis, pressure and gravitational forces, we expect the fluid to find
4Throughout this dissertation, we will use “horizontal” to refer to the directions perpendic-
ular to g.
1.2 Rotating, stratified fluids 7
a new steady state in which the applied torque is balanced by a Coriolis torque.
Hence a retrograde applied torque, for example, will “gyroscopically pump” a flow
toward the rotation axis. Viewed in an inertial (i.e. non-rotating) frame of refer-
ence, the fluid spirals in toward the rotation axis; the orbital decay of satellites
is an analogous process. Since the fluid must conserve mass, the gyroscopically
pumped flow must form part of a meridional circulation (see appendix B). The
circulation transports angular momentum, and thus a retrograde torque applied
at one location in the fluid can gyroscopically pump a circulation that spins down
the entire fluid. Spin-down in a stirred cup of tea occurs in precisely this fash-
ion; in that case the meridional circulation is driven primarily by the retrograde
frictional torque at the bottom of the cup (e.g. Greenspan, 1968).
In the Earth’s stratosphere, breaking gravity and Rossby waves constitute, in
a time-averaged sense, an effective retrograde forcing. This forcing drives the
Brewer–Dobson and Murgatroyd–Singleton mean meridional circulations (MMCs)
(e.g. Andrews et al., 1987). However, in a stably stratified system such as the
stratosphere, buoyancy forces act to restore any vertically displaced fluid elements
to their original height. This produces oscillations known as internal gravity
waves, whose frequency is proportional to the “buoyancy frequency” N (see ap-
pendix C). For large N , buoyancy inhibits motions that deform the stratification
surfaces, producing a “horizontal stiffness”. The presence of stable stratification
therefore tends to inhibit the formation of meridional circulations in rotating sys-
tems. In the absence of a “thermal relaxation” mechanism (i.e. a mechanism that
diabatically returns the system to a reference-state temperature profile) there can
be no compromise between stable thermal stratification and meridional circula-
tions: either the circulation overturns the stratification surfaces or the stable
stratification shuts off the circulation.
When a thermal relaxation mechanism is present, however, meridional circula-
tions can co-exist with stable thermal stratification on lengthscales and timescales
for which the rate of relaxation matches the rate of advection of the stratifica-
tion. In that case, meridional flows driven persistently at a particular altitude
tend to “burrow” downward over time, increasing the vertical extent of the circu-
lation. This burrowing tendency gives rise to the principle of “downward control”
(Haynes et al., 1991), which states that the steady-state meridional mass flux at
8 1. Basic Solar Physics
any altitude is determined by the gyroscopic pumping at higher altitudes.
1.3 Magnetised fluids
Motions within an ionised fluid, or plasma, can induce a magnetic field. If the
plasma is highly collisional and non-relativistic then the evolution of the magnetic
field, B say, is described by the MHD induction equation, which states that the
magnetic field lines are advected by the fluid flow, and that the field diffuses
at a rate inversely proportional to σ, the electrical conductivity of the plasma.
Specifically,
∂B
∂t = ∇× (u×B− η∇×B) , (1.3)
where η = (4piσ)−1 is the magnetic diffusivity.5 The magnetic field influences the
dynamics of the plasma through its Lorentz force, FL = (4pi)−1(∇ × B) × B.
Since the magnetic field is solenoidal, i.e. ∇ · B = 0, the Lorentz force can be
written as the sum of gradient and curvature contributions,
FL = −
1
4pi∇(
1
2 |B|
2) + 14piB ·∇B , (1.4)
which are separately called the “magnetic pressure force” and “magnetic tension”.
In the limit of perfect electrical conductivity (i.e. zero magnetic diffusivity) the
field lines behave like elastic strings “frozen into” the fluid. Perturbations to the
field lines generate Alfve´n waves, which propagate along the field lines analogously
to tension waves on an elastic string (Alfve´n, 1942). The speed of propagation
is the Alfve´n speed, which is proportional to the strength of the magnetic field.
If the Alfve´n speed greatly exceeds the typical speeds of fluid motions within
the plasma then the field lines behave rigidly. Therefore a strong magnetic field
inhibits motions that perturb the field lines.
In the limit of strong magnetic field, only fluid motions that do not perturb the
field lines are permitted. This “field-line stiffness” is somewhat analogous to the
rotational stiffness described in §1.2. If the fluid is differentially rotating, then its
angular velocity must be constant along each magnetic field line; this is known
5We use Gaussian-cgs units throughout this dissertation.
1.4 Solar spin-down 9
as Ferraro’s law of isorotation, after Ferraro (1937). In the presence of weak
magnetic diffusion, the rigidity of the field lines becomes compromised, and fluid
may “leak” across the field lines on long timescales and small lengthscales.
In a fluid that features rapid rotation, stable stratification and strong magnetic
fields, motions of the fluid will be subject to all of the physical effects described
in this and the previous section. Some simple systems that demonstrate the
interplay between the various physical effects are described in appendices C and
D. A more complicated, but more immediately relevant, case study for many
of these effects is solar spin-down — i.e. the gradual loss of angular momentum
from the Sun over its main-sequence lifetime. This is described in detail in the
next section.
1.4 Solar spin-down
1.4.1 Magnetic braking
Solar spin-down owes its origin to the “solar wind” — that is, the material ejected
from the solar surface with sufficient energy to escape the Sun’s gravity — and
the interaction of the wind with the Sun’s external magnetic field (Schatzman,
1962). Close to the solar surface the Alfve´n speed exceeds the wind speed, and
so the magnetic field lines behave rigidly. Therefore the wind is constrained to
follow the field lines outward until it reaches the “Alfve´n radius”, at which the
wind speed overtakes the (outwardly decreasing) Alfve´n speed. Moreover, within
the Alfve´n radius each magnetic field line rotates with the same angular velocity
as its footpoint on the solar surface. So the matter ejected from the surface of
the Sun in the solar wind conserves its angular velocity (rather than its angular
momentum) until it reaches the Alfve´n radius. This implies that the magnetic
field exerts a prograde “Alfve´nic” torque on the solar wind, and a corresponding
retrograde Alfve´nic torque on the solar surface.
“Magnetic braking” from the Sun’s external magnetic field acts only on the out-
ermost layers of the Sun, since the tension in the field lines is overcome by the
turbulent motions deeper within the convection zone. However, these turbulent
10 1. Basic Solar Physics
convective motions themselves transport angular momentum, and thus the spin-
down of the solar surface is communicated throughout the convection zone on
a timescale comparable to the turnover time of the largest convective eddies,
which is about one month. We can crudely model this process by parametrising
the convective turbulence as a large “turbulent viscosity” within the convection
zone.
Within the stably stratified radiative envelope, the effective viscosity returns to
its smaller, microscopic value. The timescale for purely viscous transport of
angular momentum through this region is longer than the Sun’s lifetime. We
might therefore expect solar spin-down to be confined to the convection zone,
with the radiative envelope beneath rotating more rapidly. However, there are
several mechanisms that can exchange angular momentum between the convec-
tion zone and radiative envelope, even when viscosity is negligible. These are
discussed in the following sections.
1.4.2 Meridional circulations
If magnetic braking caused the solar convection zone to rotate more slowly than
the radiative envelope, then the time-averaged turbulent stresses at the base of
the convection zone would exert a drag on the top of the radiative envelope,
gyroscopically pumping mean meridional circulations (MMCs) by the process de-
scribed in §1.2. Spiegel (1972) described how, within a few rotation periods,
these “spin-down currents” would establish a Taylor–Proudman regime within
the outer part of the radiative envelope, wherein the rotation rate would match
that of the convection zone, say Ω. The radiative envelope’s stable stratification
would temporarily confine these MMCs to a layer of thickness ∼ (2Ω/N)R, where
R is the radius of the convective–radiative interface, and N is the buoyancy fre-
quency of the radiative envelope. Spiegel called this layer6 the “tachycline” (see
figure 1.2). On a longer timescale, thermal relaxation would allow the MMCs to
burrow deeper into the interior, as described in §1.2. Within the radiative env-
elope, thermal relaxation occurs through radiative diffusion, and so this burrow-
ing process is often called “radiative spreading”. Using κ to denote the radiative
6A layer of this kind is often called a Holton layer, after Holton (1965), although the vertical
lengthscale (2Ω/N)R is named the “Rossby height”, after Rossby (1938).
1.4 Solar spin-down 11
Figure 1.2: The convection zone’s slow rotation gyroscopically pumps spin-down currents
within a “tachycline” in the outer part of the radiative envelope. Taken from Spiegel (1972).
diffusivity, the timescale7 τ for the burrowing, or spreading, is
τ =
( N
2Ω
)2 R2
κ . (1.5)
If we take Ω to be the present rotation rate of the solar surface, then τ ≈ 1011
years — much longer than the age of the Sun. But in the early, faster rotating
Sun, the burrowing tendency would have been stronger, and this timescale would
have been 109 years or less. If no other angular momentum transport mechanism
were operating below the convection zone, then by now the surface spin-down
would have been communicated perhaps half way to the centre of the Sun (see
§2.4).
1.4.3 Internal gravity waves
As mentioned in §1.1.1, downward convective plumes close to the bottom of the
convection zone may well overshoot and penetrate some way into the stably strat-
ified radiative envelope. These overshooting plumes will excite internal gravity
7The physical mechanism that drives these circulations is closely analogous to that which
drives the Eddington–Vogt–Sweet circulation (e.g. Clark, 1975). Indeed, the timescale in (1.5)
is a (radially) local Eddington–Sweet timescale.
12 1. Basic Solar Physics
waves within the stably stratified interior. In fact, even in the absence of over-
shooting convection, fluctuating Reynolds stresses at the base of the convection
zone can excite internal gravity waves in the radiative envelope.
Internal gravity waves carry an angular momentum flux. In the absence of global
rotation, waves carrying westward and eastward angular momentum fluxes are
excited and dissipated equally. But the presence of significant global rotation
introduces asymmetry between these two types of waves, and therefore allows for
wave-induced angular momentum transport (Schatzman, 1993; Zahn et al., 1997).
In order for this transport to be effective in the solar interior, the waves must
be significantly dissipated within that region. Microscopic viscosity and radiative
diffusion are too small to produce the required dissipation, so it must be supposed
either that the interior is turbulent, with an effective turbulent diffusivity, or
that the waves steepen and break at some depth within the interior. For certain
parametrisations of the waves’ generation and dissipation, the resulting angular
momentum coupling between the convection zone and radiative envelope is found
to be significant on the timescale of solar evolution (see Charbonnel & Talon,
2007, for a review of the successes of this theory). The effect of internal gravity
waves on the rotation profile of the solar interior will be discussed further in
chapter 2.
1.4.4 The Ferraro constraint
If the solar interior contains a global-scale magnetic field, Bi say, then the pic-
ture of angular momentum transport changes dramatically (Mestel & Weiss, 1987;
Charbonneau & MacGregor, 1993). In that case angular momentum is communi-
cated along the field lines by Alfve´n waves, and the pattern of angular momentum
transport therefore depends on the topology of the field. Assuming that the field
is confined within the radiative interior, and that |Bi| & 10−2 gauss, Mestel &
Weiss (1987) argued that the field would wipe out any radial differential rotation
over the course of the Sun’s main-sequence lifetime. The spin-down of the con-
vection zone would thereby be efficiently communicated throughout the entire
Sun.
1.4 Solar spin-down 13
Charbonneau & MacGregor (1993) constructed a simple axisymmetric, time-
dependent computational model of spin-down in the solar interior in the presence
of a global-scale magnetic field. They assumed that the stable stratification of the
radiative interior would suppress any vertical motions within that region. Since
their system was axisymmetric, only rotational motions were therefore permit-
ted, and so Coriolis effects were absent. Furthermore, they neglected diffusion
of the poloidal components of the magnetic field, so that the poloidal field con-
figuration, which they took to be an axisymmetric dipole, could be statically
imposed throughout the simulation. Therefore the only equations solved were
the azimuthal components of the momentum and induction equations, and the
only forces considered were Lorentz and viscous torques.
Charbonneau & MacGregor (1993) parametrised the turbulence in the convec-
tion zone as a large turbulent viscosity, and introduced a distributed “angular
momentum sink” within the convection zone to represent magnetic braking by the
solar wind. Over the course of their simulations, turbulent viscosity maintained
near-uniform rotation within the convection zone, which was therefore spun down
at the same rate as the solar surface. The manner in which spin-down was
communicated to the interior was found to depend on the topology of the imposed
poloidal magnetic field Bi. The spin-down of the convection zone was rapidly
communicated along any field lines that connected the radiative envelope to the
convection zone, leading to a Ferraro state (recall §1.3). In cases for which
the field was located entirely below the base of the convection zone, spin-down
within the radiative interior was initially confined to a growing viscous boundary
layer at the interface with the convection zone. Once this boundary layer grew
sufficiently thick to encounter the deep magnetic field, a Ferraro state was again
established within the interior. In all the cases they considered, a quasi-steady
Ferraro rotation profile was established within approximately 107 years.
The effect of an interior magnetic field on angular momentum transport is dis-
cussed in further detail in chapter 3.
14
Chapter 2
The Sun’s Differential Rotation
2.1 Standard solar models
The internal structure of a star is, in most cases, uniquely determined by its
mass and composition (e.g. Ka¨hler, 1978). Hence the star’s future evolution is
also uniquely determined; this is known as the “Vogt–Russell theorem”. In one-
dimensional solar evolution models it is commonly assumed that the Sun was
compositionally homogeneous at the start of its main-sequence phase (see §5).
In principle, therefore, the internal structure of the Sun at any time on the main
sequence can be determined from its initial mass, and the initial mass fractions
of each of its chemical constituents. However, the effect of turbulence within the
Sun’s convection zone cannot be precisely quantified by any analytical formula,
and must instead be parametrised within one-dimensional solar evolution models.
The most commonly used parametrisation is the “mixing-length” treatment of
Bo¨hm-Vitense (1958).
Furthermore, most solar evolution models assume that the Sun’s mass is constant.
The necessary input parameters for these models are therefore the mixing-length
parameter α, the helium mass fraction Y , and the “metallicity” Z, which is
the mass fraction of all elements that are heavier than helium. The models are
then calibrated by adjusting these input parameters to match the Sun’s present
radius, luminosity and the ratio of Z/X at the solar surface, where X is the mass
fraction of hydrogen. The surface ratio Z/X is not directly observable, but can
15
16 2. The Sun’s Differential Rotation
be inferred from spectroscopic data, via a suitable model of the solar atmosphere
(e.g. Christensen-Dalsgaard, 2002, & refs).
Until fairly recently, surface measurements of the abundance ratio Z/X were
obtained using one-dimensional models of the solar atmosphere, and assumed
local thermal equilibrium (LTE) (e.g. Holweger & Mu¨ller, 1974). In such models,
the broadening and displacement of spectral lines by convective turbulence must
be parametrised, usually in terms of “microturbulence” and “macroturbulence”
(e.g. de Jager, 1972). More recent abundance measurements have been based on
three-dimensional models of convection at the solar surface (Stein & Nordlund,
1998), and incorporate many non-LTE effects. The consequences of such effects
for surface abundance measurements have been reviewed by Asplund (2005).
Arguably the main advantage of the new three-dimensional models over previous
models is that they contain no adjustable “turbulence” parameters (Asplund
et al., 2006), and therefore permit less ambiguity in the interpretation of the
results. The development of these models has led to significant revision of the
surface abundance measurements. The revised abundances of carbon, nitrogen
and oxygen are lower than previous estimates (e.g. Asplund et al., 2009), and
more closely reflect the abundances of the local interstellar medium (Grevesse
et al., 2007, & refs). The abundance ratio Z/X at the solar surface is now
estimated to be 0.0181, which is significantly lower than previous estimates (e.g.
Grevesse & Sauval, 1998, estimated 0.023).
The revised abundances have serious implications for solar evolution models, to
be discussed in the next section.
2.2 Helioseismology
In the 1970s, observations of the Sun’s surface oscillations demonstrated the
existence of acoustic (p) modes, excited by turbulence in the convection zone
(Appourchaux et al., 2010, & refs). These observations, taken together with
reasonable dynamical assumptions, allow certain properties of the solar interior
to be inferred. In particular, if the oscillations are assumed to be linear, adiabatic
perturbations to a spherically symmetric and hydrostatic background, then the
2.2 Helioseismology 17
radial profiles of pressure, density and sound speed can be inferred from the
oscillation spectrum (e.g. Christensen-Dalsgaard & Thompson, 2007, & refs).1
The sound speed is particularly well constrained (to within a few parts in 104)
since it is the principal factor determining the acoustic oscillation frequencies.
However, the vertical profiles of temperature and luminosity, for example, cannot
be inferred without invoking additional, thermodynamical assumptions.
Helioseismology provides a means of testing models of the Sun’s interior struc-
ture. In particular, it can be used to locate the base of the convection zone,
defined (in the context of helioseismology) to be the depth at which the mean
stratification changes from almost adiabatic to subadiabatic. This transition is
located at (0.713± 0.003)R (Christensen-Dalsgaard et al., 1991). Ideally, solar
evolution models should be able to reproduce this result with reasonable accu-
racy. Helioseismology can also be used to measure the helium mass fraction Y in
the convection zone, provided that the equation of state is known (Basu & Antia,
1995), which provides an additional test of solar evolution models.
Until recently, standard solar models were in close agreement with the results of
helioseismology (e.g. Christensen-Dalsgaard et al., 2009, & refs). However, solar
models that are calibrated to match the revised abundances described in §2.1
agree less well with helioseismology, primarily because the opacity of solar ma-
terial is sensitive to the abundance of heavy elements (e.g. Basu & Antia, 2004).
Recalibrating standard solar models to match the revised abundances leads to an
opacity decrease of around 25% near the base of the convection zone (Montalba´n
et al., 2006, & refs). As a result, the convection zone in these recalibrated models
is significantly thinner, by about 10Mm (Bahcall & Pinsonneault, 2004). These
recalibrated models also have a smaller helium mass fraction Y in the convec-
tion zone than is inferred from helioseismology, and a smaller sound speed in the
radiative envelope.
Many authors have sought a possible resolution to the discrepancies between the
recalibrated solar models and helioseismic results (see review in Montalba´n et al.,
2006). Since the most significant effect of the revised abundances is a reduction
1The inversion is achieved by iteratively refining a chosen reference model, and so the results
may depend, to some extent, on the reference model chosen. See Christensen-Dalsgaard &
Thompson (2007) for further details.
18 2. The Sun’s Differential Rotation
in opacity, the simplest resolution is simply to artificially increase the opacity
by an appropriate amount, over the required temperature range (e.g. Bahcall &
Serenelli, 2005; Christensen-Dalsgaard et al., 2009). However, there seems to be
no physical justification for such a large opacity increase (Badnell et al., 2005).
Since gravitational settling leads to a reduction in the surface abundance of heavy
elements (see chapter 5), an alternative resolution is to increase the settling ve-
locities used in solar models (e.g. Montalba´n et al., 2004). In this way, the present
surface abundances can be explained without significantly changing the primor-
dial abundances. However, this also reduces the abundance of helium in the
convection zone, again in contradiction with the helioseismic results.
The opacity of solar material is highly sensitive to the abundance of neon, whose
surface abundance cannot be measured directly by spectroscopy. Increasing the
neon abundance of standard solar models by a factor of three would largely
compensate for the reduction in other heavy-element abundances (Bahcall et al.,
2005). Whether such a high neon abundance is compatible with observations
remains controversial (Grevesse et al., 2007, & refs).
As mentioned in §2.1, the “convection zone” defined by helioseismology includes
any region in which entropy is well mixed, for example by convective overshoot.
Including an overshoot layer in the recalibrated solar evolution models, with a
thickness ≈ 10Mm, therefore yields a convection zone whose thickness and helium
abundance are in agreement with helioseismology. However, the sound speed
profile within the radiative envelope in these models still departs significantly
from the profile obtained seismologically (Montalba´n et al., 2006).
The present uncertainty regarding the Sun’s internal structure must be borne
in mind when constructing any model of its internal dynamics. Fortunately,
helioseismology is able to provide considerable information with relatively few
assumptions. In the next section we review the properties of the Sun’s interior
rotation that have been inferred from helioseismology.
2.3 Differential rotation 19
2.3 Differential rotation
Carrington’s sunspot observations of the 1860s revealed that the surface of the
Sun is differentially rotating, and that its angular velocity increases monotoni-
cally from pole to equator. Subsequent observations of sunspots and other mag-
netic tracers, as well as Doppler-shift measurements, allowed for more accurate
measurement of the surface rotation, and it was found that the pole-to-equator
variation is approximately 30% of the mean rotation rate.
With the advent of helioseismology, it became possible to infer the Sun’s interior
rotation rate from the “rotational splitting” of the acoustic frequency spectrum
(e.g. Christensen-Dalsgaard & Thompson, 2007, & refs). It was found that the
surface pattern of differential rotation persists qualitatively unchanged down to
the base of the convection zone, but that the radiative envelope beneath is in
approximately uniform rotation (see figure 2.1). The early helioseismic inferences
Figure 2.1: The angular velocity of the solar interior, adapted from Schou et al. (1998). The
radiative envelope (below the dashed line) rotates approximately uniformly, with an angular
velocity Ωi ≈ 2.7 × 10−6s−1. The convection zone (above the dashed line) exhibits significant
differential rotation.
were subsequently refined, and the transition from differential to uniform rotation
was found to occur across a thin shear-layer, which came to be known as the
“tachocline” (Spiegel & Zahn, 1992). Even with modern helioseismology, the
20 2. The Sun’s Differential Rotation
tachocline is too thin to be accurately resolved, but its thickness is thought to
be . 30Mm — less than 5% of the Sun’s radius R (e.g. Christensen-Dalsgaard
& Thompson, 2007; Howe, 2009, & refs). The tachocline appears to be centred
in the stably stratified radiative envelope, at 0.69R, but possibly straddles the
boundary with the convection zone, particularly in high latitudes (e.g. Basu &
Antia, 2003). In fact, Basu & Antia found that, whereas the bottom of the tacho-
cline is very nearly spherical, the top of the tachocline is prolate, although their
results are also consistent with a discontinuous change in the tachocline thickness
at the zero-shear latitudes ≈ ±30◦.
Superimposed on the convection zone’s differential rotation are the so-called “tor-
sional oscillations” first observed at the surface by Howard & LaBonte (1980).
These are bands of prograde and retrograde motion that propagate equatorward
in low latitudes and poleward in high latitudes, with period of approximately 11
years. They are believed to be closely connected to the Sun’s magnetic cycle (see
§2.5). Helioseismology reveals that these bands occupy at least the outer 60Mm
of the convection zone (e.g. Antia & Basu, 2000), and may extend all the way
to be the base of the convection zone in high latitudes (Vorontsov et al., 2002).
There is some evidence that the low latitude oscillation propagates upward, as
well as equatorward (Basu & Antia, 2003).
Howe et al. (2000) found evidence of similar torsional oscillations at the base of
the low-latitude convection zone, but with a period of 1.3 years. The angular
velocity variations were most pronounced either side of the tachocline, with one
band of oscillation at 0.63R, and another at 0.72R; the two bands oscillate
out of phase, such that the shear variations are largest in the tachocline. The
existence of tachocline oscillations has significant implications for theories of the
solar dynamo (see §2.5). However, subsequent studies have failed to confirm the
presence of these oscillations (e.g. Antia & Basu, 2000; Basu & Antia, 2003; Howe
et al., 2007). It is possible that the oscillations switched off in 2001, around the
time of maximum solar activity (Howe et al., 2007), and may reappear with the
new solar cycle.
The early measurements made using helioseismology suggested that the Sun’s
core, below about 0.2R, rotates significantly faster than the outer layers (e.g.
Claverie et al., 1981). However, subsequent studies have failed to yield consis-
2.3 Differential rotation 21
tent measurements for the rotation rate in this region. In addition, the rotation
rate close to the rotation axis is not well constrained by the observations (see
figure 2.1). All data gathered in the last three decades is consistent with uni-
form rotation throughout the radiative interior, with an interior angular velocity
Ωi ≈ 2.7× 10−6s−1 (see Howe, 2009, for a review).
Perhaps the most important observation concerning the Sun’s interior rotation
is that the average angular velocity over each spherical surface is roughly the
same, even within the differentially rotating convection zone. This demonstrates
that the spin-down of the solar surface by magnetic braking has been efficiently
communicated throughout the interior, and can help us to decide which (if any) of
the angular momentum transport mechanisms mentioned in §1.4 are predominant
in the solar interior.
However, helioseismology has also revealed some surprising features of the Sun’s
interior rotation. Prior to the advent of helioseismology, the existence of the ta-
chocline, i.e. the sharp transition from differential to uniform rotation, had not
been anticipated.2 In the following sections we review various attempts to answer
the following questions:
1. What is the origin of the convection zone’s differential rotation?
2. How is uniform rotation maintained in the radiative envelope?
3. Why is the tachocline so thin?
2.3.1 The influence of rotation
The differential rotation within the Sun’s convection zone must, in some fash-
ion, be driven by the turbulent convection within that region. Although a com-
plete description of this process is lacking, it is well known that weakly nonlin-
ear, rapidly rotating (quasi-geostrophic) convection in spherical geometry pro-
duces the kind of “equatorial acceleration” (i.e. latitudinal differential rotation)
2Although the tachocline bears a superficial resemblance to Spiegel’s tachycline (§1.4.2), the
tachycline was expected to thicken over time, and by the Sun’s present age would extend most
of the way to the Sun’s core — see §2.4.
22 2. The Sun’s Differential Rotation
observed at the solar surface. In that case, the phenomenon is attributed to
the “banana shape” of the convective rolls in low latitudes, which implies an
anisotropic Reynolds stress that transports angular momentum equatorward (e.g.
Busse & Hood, 1982). This behaviour survives into the nonlinear regime provided
that the convection is still “strongly influenced” by rotation (e.g. Gilman, 1977).
However, the rotation profile produced in such cases is usually subject to the
Taylor–Proudman constraint, i.e. the angular velocity contours are aligned with
the rotation axis. By contrast, the Sun’s angular velocity contours are more
“conical” (see figure 2.2).
Over the last 30 years many studies of rotating convection have sought to explain
the dynamics behind the convection zone’s differential rotation (for a more ex-
tensive review, see Miesch, 2005). Most of these fall into one of two categories:
direct numerical simulations and mean-field models. Here we will focus on direct
numerical simulations. As mentioned above, numerical simulations of rotating
Figure 2.2: The convection zone’s angular velocity contours. The left panel shows the nu-
merical model of Gilman & Miller (1981), in which the contours are aligned with the rotation
axis, particularly in low latitudes. The right panel shows the real Sun, as inferred from helio-
seismology, with conical contours inclined at about 30◦ to the rotation axis.
convection in spherical shells are able, in certain parameter regimes, to produce
equatorial acceleration similar to that observed for the Sun. However, the time-
averaged angular velocity contours in these simulations tend to be aligned with
2.3 Differential rotation 23
the rotation axis. The early computational models considered only Boussinesq
fluids, for reasons of computational ease, and therefore did not include the ef-
fects of compressibility. But such effects are expected to be significant in the
convection zone, particularly in its outer part, where the vertical density gradi-
ent is steepest. Compressible convection is characterised by narrow, concentrated
downflows and broad, weak upflows (Hurlburt et al., 1986; Cattaneo et al., 1991).
The downflows can take the form of isolated plumes or connected lanes, depend-
ing on the strength of the rotational influence (e.g. Brummell et al., 2002). This
suggests that compressible effects might significantly alter the angular momen-
tum transport properties of rotating convection (e.g. Glatzmaier et al., 2009, &
refs).
In order to more closely approximate the dynamics of the convection zone, most
modern computational models make use of the “anelastic” approximation, which
retains the energetics of low Mach number compressibility but filters out acoustic
waves (Gough, 1969). In these anelastic models, the quasi-geostrophic convection
rolls that characterise rapidly rotating, weakly nonlinear convection in Boussinesq
fluids are replaced by more complicated structures with distinct vertical asym-
metry. Nonetheless, rotating anelastic convection retains certain properties of
its Boussinesq counterpart. Provided that the rotational influence is sufficiently
strong, concentrated downflowing lanes with a large north–south extension are
found to form in low latitudes. Like the banana-shaped Boussinesq convection
rolls, these downflowing lanes are titled such as to transport angular momentum
equatorward (Miesch et al., 2000). Strongly rotating anelastic convection can
therefore produce the same equatorial acceleration as strongly rotating Boussi-
nesq convection. However, as with Boussinesq convection, the time-averaged
angular velocity contours are found to be roughly aligned with the rotation axis,
so the effects of compressibility alone do not overcome the Taylor–Proudman
constraint.
It is generally accepted that the Sun’s differential rotation is due to rotation-
ally induced anisotropies in the convective turbulence. However, it has not been
conclusively shown in which particular parameter ranges an equatorial accelera-
tion occurs, and whether the rotation profile in these parameter ranges is always
Taylor–Proudman. The strength of convection is usually measured in terms of
24 2. The Sun’s Differential Rotation
the Rayleigh number Ra. At sufficiently high Rayleigh number the effects of
baroclinicity overcome the Taylor–Proudman constraint (recall §1.2), and the
structure of convection is no longer geostrophic. But numerical simulations at
such high Rayleigh number typically exhibit an equatorial deceleration. Differ-
ential rotation of this kind is often called “anti-solar”, and may be the result of
global homogenisation of angular momentum by the convective turbulence (e.g.
Aurnou et al., 2007).
It is usually assumed, following Gilman (1977), that convection is “strongly in-
fluenced” by rotation when the “convective Rossby number” Roc is smaller than
unity. This is the Rossby number based on the timescale for the growth of con-
vective elements in the absence of dissipative effects, and can be expressed as
Roc =
(RaEk2
Pr
)1/2
, (2.1)
where Pr is the Prandtl number and Ek is the Ekman number. For fixed Prandtl
number, the condition Roc . 1 implies that the transitional Rayleigh number
Rat, below which convection is “strongly influenced” by rotation, is proportional
to Ek−2.
More recently, King et al. (2009) have argued that convectively unstable systems
will naturally adopt whichever form of convection allows the more efficient trans-
port of heat, which can be quantified by the Nusselt number, Nu. Laboratory
and numerical experiments suggest that, for fixed Prandtl number, the Nusselt
number in non-rotating (Ek = ∞) convection is proportional to Ra2/7, whereas
for quasi-geostrophic (Ek → 0) convection the Nusselt number is proportional
to Ra6/5Ek8/5. King et al. (2009) therefore argue that the transitional Rayleigh
number Rat must obey the relation
Ra2/7t ∝ Ra6/5t Ek8/5 (2.2)
⇒ Rat ∝ Ek−7/4. (2.3)
This implies that the transitional Nusselt number, Nut say, is proportional to
Ek−1/2. Furthermore they show that this dependence of Nut on Ek can be ex-
plained by a simple boundary layer argument. Whether the Prandtl number
2.3 Differential rotation 25
dependence of the transition can be similarly explained is unclear.
Any argument based only on parameter values, with no regard to either ge-
ometry or boundary conditions, is likely to be overly simplistic. But taken at
face-value, the scaling arguments above suggest that in the (narrow) parameter
range Ek−7/4  Ra  Ek−2 the structure of convection will be non-geostrophic
(i.e. non-Taylor–Proudman), and yet Coriolis effects will be significant on the
timescale of the convective motions. The existence of such a regime has not been
demonstrated, however.
2.3.2 Thermal-wind
If we assume that both hydrostatic balance and thermal-wind balance hold ro-
bustly within the solar interior, then we can use the rotation profile shown in
figure 2.1, together with equation (1.1), to calculate the density distribution,
and hence the temperature distribution, over every horizontal surface. From the
qualitative pattern of rotation shown in figure 2.1 we anticipate that, on each
horizontal surface, the temperature increases from equator to pole. In the Sun,
the pressure surfaces are very nearly spherical (more so than the density and
temperature surfaces), so we also expect to find “warm polar anomalies” on each
spherical surface. These anomalies will be most pronounced in the tachocline,
where axial variations in the angular velocity are greatest; Miesch (2005) esti-
mated that the base of the convection zone is roughly 5K warmer at the poles than
at the equator, although this estimate is sensitive to any errors in the rotation
profile obtained by helioseismology. By directly imposing a latitudinal entropy
gradient at the base of the convection zone in their simulations, corresponding to
a temperature increase of about 10K from equator to pole, Miesch et al. (2006)
were able to produce angular velocity contours in rough agreement with those
observed in the Sun. This suggests that the interaction between the convection
zone and the tachocline may play an essential role in determining the rotational
profile in the solar interior.
Miesch et al. (2006) suggested two different physical explanations for the pres-
ence of warm polar anomalies in the convection zone. The first, which they called
“thermal forcing”, supposes that the influence of rotation reduces the convective
26 2. The Sun’s Differential Rotation
heat flux at certain latitudes, and hence produces latitudinal variations in the
time-averaged temperature profile. Thermal forcing of this kind has been incor-
porated successfully into mean-field models (e.g. Kitchatinov & Ru¨diger, 1995) as
a means to overcome the Taylor–Proudman constraint. The other explanation,
which Miesch et al. called “mechanical forcing”, relies on the interaction between
the convection zone and the stably stratified tachocline. Whatever process drives
the convection zone’s differential rotation will also gyroscopically pump mean
meridional circulations (MMCs) that burrow into the tachocline (see §4.2.6 and
appendix B). These circulations are expected to be downwelling near the poles
(see §1.4.2 and further discussion in §2.4). Within the tachocline, these down-
welling flows will produce warm polar anomalies by adiabatic compression, i.e.
advection of the background entropy profile. This process has been observed in
the mean-field model of Rempel (2005).
Although Miesch et al. considered how thermal or mechanical forcing might drive
the convection zone and tachocline towards a state of thermal–wind balance,
they did not attempt to construct a complete picture of the balanced system.
However, if we assume that the anomalous polar heating required to explain the
Sun’s angular velocity profile does indeed originate in the tachocline, then it
is possible to make a strong statement about the MMCs in that region. Since
the stably stratified tachocline is a thermally relaxing system, departures from
local radiative equilibrium will tend to relax back toward zero. To maintain a
warm anomaly near the pole, there has to be persistent adiabatic compression
by downwelling. This point was recognised much earlier by Gough & McIntyre
(1998, see §4.2.1). In principle the required downwelling can be estimated from
the magnitude and structure of the warm anomaly, which in turn can be estimated
from the Sun’s rotation profile, as described above. In this fashion, Gough &
McIntyre estimated a downwelling velocity of 10−5cm s−1 in the polar tachocline.
To make a more precise estimate we would need a more accurate measurement
of the shear in the tachocline.
2.4 The tachocline 27
2.4 The tachocline
The most puzzling feature of the Sun’s rotation is the thinness of the tachocline.
As described in chapter 1, we expect any differences between the rotation rates
of the convection zone and the radiative envelope to gyroscopically pump mean
meridional circulations (MMCs) that burrow into the stably stratified interior.
These MMCs transport angular momentum, and so the burrowing would cause
the convection zone’s differential rotation to spread into the interior. For the
deep interior to remain in uniform rotation, the MMCs driven at the base of the
convection zone must somehow be prevented from burrowing into the radiative
envelope, and thereby confined to the tachocline.
The problem of confining the convection zone’s MMCs was first elucidated by
Spiegel & Zahn (1992), in the first paper on the tachocline. They described, using
a laminar, hydrodynamic model, how MMCs driven at the base of the convection
zone seek to establish a Taylor–Proudman state within the radiative envelope,
leading to a thickening of the tachocline on the local Eddington–Sweet timescale
(1.5). This is the same “radiative spreading” process described in chapter 1.
Spiegel & Zahn’s results were confirmed by the numerical model of Elliott (1997,
see figure 2.3); see also §4.2.6. If this spreading had been active throughout the
Sun’s main-sequence lifetime, including during the early part of the main sequence
when the Sun’s faster rotation made the burrowing tendency stronger, then the
MMCs and differential rotation would by now have burrowed approximately half-
way to the centre of the Sun. This result is not compatible with the helioseismic
data. We must therefore conclude that some additional mechanism, absent from
their laminar, hydrodynamic model, acts to stop the burrowing of the tachocline’s
MMCs, and thereby maintains the uniform rotation of the radiative envelope.
The tachocline’s MMCs, which are gyroscopically pumped by turbulent stresses
in the convection zone, can be prevented from burrowing only by compensating
gyroscopic pumping within the tachocline. That is, angular momentum must
be efficiently redistributed in the tachocline so as to precisely balance the anti-
frictional redistribution of angular momentum in the layers above. The (mathe-
matically) simplest mechanism that can redistribute angular momentum in this
way is a large horizontal viscosity, which was the mechanism invoked by Spiegel
28 2. The Sun’s Differential Rotation
Figure 2.3: Radiative spreading in the solar interior, from the model of Elliott (1997). The
left panel shows how deeply the convection zone’s prograde (solid contours) and retrograde
(dashed contours) differential rotation burrow after one solar age. The right panel shows the
MMCs that produce this burrowing, with dashed contours indicating anti-clockwise circulation.
& Zahn to halt radiative spreading in their model. They assumed that the com-
bination of latitudinal shear and stable stratification in the tachocline would pro-
duce predominantly horizontal hydrodynamic turbulence, i.e. turbulent motions
largely confined to horizontal surfaces. They further assumed that the turbulence
would redistribute angular momentum in the manner of a large horizontal eddy
viscosity (see also Zahn, 1992). However, several questions have subsequently
been raised regarding this turbulence model. Firstly, it is not clear whether the
shear in the tachocline is able to sustain turbulence. Many studies have examined
the possibility of hydrodynamic and magnetohydrodynamic (MHD) instability in
the tachocline (e.g. Gilman & Cally, 2007; Ogilvie, 2007; Parfrey & Menou, 2007;
Spruit, 1999, see appendix F). As yet, there seems to be no consensus as to which
instabilities are present, and whether they can lead to sustained turbulence in
the tachocline.
Secondly, the kind of horizontal, shear-driven, hydrodynamic turbulence proposed
by Spiegel & Zahn has no tendency to act like a viscosity. Instead, the horizontal
turbulent eddies act to locally homogenise potential vorticity (PV), which can
produce distinctly “anti-frictional” effects (e.g. McIntyre, 1994, 2003, & refs). In
the context of the tachocline, such turbulence would probably act to smooth any
large extrema in the PV distribution, bringing the system to a marginally stable
state without enforcing uniform rotation (Garaud, 2001).
2.4 The tachocline 29
On the other hand, horizontal MHD turbulence behaves rather differently from
its hydrodynamical counterpart, and might well be free from the anti-frictional
behaviour just described (McIntyre, 2003; Tobias et al., 2007). But even if hor-
izontal MHD turbulence can mix angular momentum locally, it would need to
act efficiently over all latitudes in the tachocline in order to prevent radiative
spreading. And even if horizontal turbulence is assumed to act at all latitudes
and depths within the radiative interior, an additional mechanism must then be
invoked to explain the lack of radial shear in that region.
Given the difficulties with using turbulence to explain the radiative envelope’s
uniform rotation, it seems worth considering whether the other mechanisms for
angular momentum transport mentioned in chapter 1 might provide a more plau-
sible explanation. For example, internal gravity waves can transport angular
momentum over large distances within the stably stratified interior (§1.4.3). If
the amplitude and spectrum of this wave generation were finely tuned over all
latitudes, then the resulting angular momentum transport could prevent the ta-
chocline’s MMCs from burrowing into the radiative envelope. However, this sit-
uation seems unlikely to be realised within the solar interior, particularly in the
early, faster rotating Sun when the burrowing tendency was stronger.
Nevertheless, several studies have advocated gravity waves as the explanation
for the Sun’s uniform interior rotation (e.g. Talon et al., 2002; Charbonnel &
Talon, 2007, & refs). These studies assume, following Zahn (1992), that horizon-
tal turbulence is able to enforce uniform rotation over each horizontal surface,
and therefore seek only to describe how gravity waves can act to reduce the inte-
rior’s radial shear. By assuming a particular spectrum for the wave generation,
which is based on the spectrum of turbulence in the lower convection zone, and
by introducing a large turbulent viscosity into any regions of radial shear, the
overall effect of angular momentum transport by gravity waves is found to be a
smoothing of the radial rotation profile. However, it is frequently argued (e.g.
Dintrans et al., 2005, & refs) that the main source of internal gravity waves in the
radiative envelope is convective overshoot, a source which is neglected in these
models. Overshooting convection typically generates a broad-band spectrum of
gravity waves (Hurlburt et al., 1986; Rogers & Glatzmaier, 2006), and so would
probably induce “anti-frictional” angular momentum transport (e.g. McIntyre,
30 2. The Sun’s Differential Rotation
1994; Gough & McIntyre, 1998; Rogers et al., 2008). Furthermore, these models
do not take into account Coriolis effects on either the generation or the propaga-
tion of the waves.
A model of gravity wave generation and dissipation in the radiative interior that
allows for latitudinal shear, as well as Coriolis effects, is currently being developed
(see Mathis, 2009), but results are not yet available. It is worth noting that
the presence of a global-scale magnetic field in the radiative interior would also
significantly alter both the generation and the propagation of the waves, as well
as their angular momentum transport properties.
Having recognised all of the difficulties associated with purely hydrodynamical
models of the tachocline, Gough & McIntyre (1998) concluded that the only
mechanism that can plausibly explain the uniform rotation of the radiative en-
velope is a global-scale interior magnetic field. The impact of such a magnetic
field on angular momentum transport was mentioned in §1.4.4, and is discussed
in detail in the next chapter. In chapter 4 we give a summary of the Gough &
McIntyre model.
2.5 Dynamo theory
The window into the Sun’s interior rotation provided by helioseismology has
significantly influenced theories of the solar dynamo and the solar cycle. In this
section we briefly review the historical development of dynamo theory as applied
to the Sun. More detailed discussions of dynamo theory can be found in the
books by Moffatt (1978) and Krause & Ra¨dler (1980).
As early as 1850 it was known that the number of sunspots on the Sun’s surface
follows an 11-year cycle, now known as the Schwabe cycle. Subsequently, Hale
(1908) and Hale et al. (1919) discovered that sunspots are the sites of strong
magnetic fields, and that the polarity of sunspot magnetic fields reverses at the
start of each 11-year cycle. The Sun therefore has a 22-year magnetic cycle, which
is now known as the Hale cycle. Also in 1919, Larmor proposed that the Sun’s
surface magnetic field is generated by a hydromagnetic dynamo mechanism.
2.5 Dynamo theory 31
The first significant development of dynamo theory was the proof by Cowling
(1933) that a dynamo magnetic field must necessarily be non-axisymmetric, which
rules out the possibility of simple axisymmetric dynamos. In an axisymmetric
system, the action of differential rotation can twist poloidal magnetic field lines
to generate a toroidal field, but there is no mechanism that can regenerate the
poloidal field from the toroidal field. A possible solution to the dynamo problem
was proposed by Parker (1955a). He noted that, in a rapidly rotating turbulent
system such as the Sun’s convection zone, Coriolis effects lead to cyclonic (i.e.
helical) motions within the fluid. He then showed that certain small-scale heli-
cal motions, in the presence of magnetic diffusion, can regenerate a large-scale
poloidal field from a large-scale toroidal field.
Parker’s idea, that small-scale turbulence can regenerate a large-scale poloidal
field, was developed into a formal mathematical theory by Steenbeck et al. (1966),
which came to be known as “mean-field electrodynamics”. Today, the generation
of poloidal field by small-scale turbulence, and the generation of toroidal field
by differential rotation, are known as the “alpha” (α) and “omega” (ω) effects
respectively, following the notation employed by Steenbeck & Krause (1969).
It is now generally accepted that the Sun’s 22-year magnetic cycle results from
a large-scale αω-dynamo process operating within the convection zone, and that
sunspots are the surface manifestation of the dynamo field. The Sun’s interior
differential rotation generates a strong toroidal magnetic field, which becomes
buoyantly unstable (Parker, 1955b) and rises to the surface in the form of “mag-
netic flux tubes”. Upon reaching the surface, the flux tubes extend into the
Sun’s atmosphere as coronal loops. Sunspots are the footpoints of these loops
on the solar surface. The strong magnetic field in sunspots inhibits convection,
and so sunspots are cooler, and hence darker, than the rest of the solar surface.
This model of sunspot formation also offers an explanation for many other ob-
served properties of sunspots (Parker, 1955b), including their bipolarity and near
east-west orientation.
However, it can be argued that any dynamo magnetic field within the bulk of
the convection zone will be brought to the surface, by a combination of magnetic
buoyancy and turbulent convection, on a timescale of months, rather than years
(Parker, 1975; Schmitt, 1993, & refs). This suggests that the toroidal component
32 2. The Sun’s Differential Rotation
of the dynamo field must be located in a boundary layer the base of the convection
zone (Spiegel & Weiss, 1980). Helioseismology provides additional support for
this view, since it shows that the Sun’s differential rotation is most pronounced
at the base of the convection zone, i.e. within the tachocline. Spiegel & Weiss
therefore proposed that the solar dynamo operates entirely within the stably
stratified overshoot layer beneath the convection zone, within which the alpha
effect arises from the helicity of overshooting convective plumes.
However, there are difficulties in reconciling sunspot observations with a model
of the solar dynamo based entirely at the base of the convection zone. Since very
few sunspots are observed at latitudes & 35◦, toroidal flux tubes must rise almost
vertically through the convection zone, rather than parallel to the Sun’s axis of
rotation (Parker, 1993, & refs). This requires the toroidal field strength to be
∼ 105 gauss at the base of the convection zone, in order for Lorentz forces to
dominate Coriolis forces during the buoyant rise of the flux tubes. However, the
magnetic energy associated with a field of this strength significantly exceeds the
kinetic energy associated with turbulent convective motions near the base of the
convection zone. The Lorentz force from the field is therefore expected to sup-
press any turbulence in its vicinity, and thus “quench” the alpha effect (Parker,
1993, & refs). Indeed, in a system with a large magnetic Reynolds number, such
as the convection zone, alpha-quenching is expected to occur well before the mag-
netic energy reaches equipartition with the turbulent kinetic energy (Cattaneo &
Hughes, 1996, & refs).
A possible resolution of the alpha-quenching problem was proposed by Parker
(1993), (see also Charbonneau & MacGregor, 1997). Parker proposed an “inter-
face dynamo” in which the alpha and omega effects occur in spatially distinct
regions, either side of the interface between the convection zone and radiative
envelope. In this model, the toroidal magnetic field is significantly weaker in
the region above the interface, as a result of turbulent magnetic diffusion within
the convection zone. Alpha quenching is therefore averted within this region.
Transport of poloidal and toroidal field across the interface is mediated by a
combination of laminar and turbulent diffusion. The solar cycle and torsional
oscillations can then be explained in terms of “dynamo waves” that propagate at
the interface between the alpha and omega regions.
2.5 Dynamo theory 33
An alternative resolution of the alpha-quenching problem can be found in the
model of Wang & Sheeley (1991), which is based on earlier work by Babcock
(1961) and Leighton (1969). In this “Babcock–Leighton” or “flux-transport”
dynamo model, generation of poloidal field3 is attributed both to the twisting of
rising flux tubes by Coriolis forces and to the redistribution of magnetic flux by
cyclonic motions and magnetic diffusion over the upper surface of the convection
zone. This is not strictly an alpha effect, because the poloidal field is generated
by large-scale motions rather than small-scale turbulence. Recent flux-transport
dynamo models have placed increasing emphasis on the role of the convection
zone’s mean meridional circulation in transporting the poloidal field from the
surface to the base of the convection zone (e.g. Dikpati & Charbonneau, 1999).
In these recent models, the turnover time of the convection zone’s meridional
circulation is the main factor determining the period of the solar cycle, and the
Sun’s torsional oscillations are attributed to meridional advection rather than a
dynamo-wave mechanism.
The Sun’s dynamo field, which reverses every 11 years, must remain disconnected
from the pseudo-steady interior magnetic field Bi needed to explain the interior’s
uniform rotation.4 This requires the existence of some mechanism that can con-
fine the interior magnetic field against diffusion, as will be discussed in detail in
chapter 4. If the interior field is sufficiently strong then it might still influence
the solar cycle, perhaps inducing a 22-year cyclicity in solar activity (Levy &
Boyer, 1982). There is some observational evidence for such cyclicity (Mursula
et al., 2001, & refs), although an interior magnetic field is not the only possible
explanation (e.g. Charbonneau et al., 2007). In chapter 3 we discuss the origin of
the Sun’s interior magnetic field and its angular momentum transport properties.
3The generation of toroidal field is still attributed to the omega effect, i.e. to the action of
differential rotation.
4Wale´n (1946) proposed that the solar cycle might arise from global Alfve´nic oscillations of
an interior magnetic field, but in that case the Sun’s differential rotation would also reverse
every 11 years — see §3.2.
34
Chapter 3
The Magnetic Interior
3.1 The Sun’s fossil field
In the neighbourhood of our solar system, the interstellar magnetic field has a
strength of roughly 10−6 gauss. The molecular cloud from which the Sun formed
in all likelihood contained a magnetic field of similar strength, which would have
been amplified during the Sun’s gravitational collapse. It is not known to what
extent this primordial field was modified by deep convection during the Sun’s
pre-main sequence Hayashi phase, but observations of T Tauri stars suggest that
surface field strengths of 103 gauss are typical (e.g. Bouvier et al., 2007).
Some component of the Sun’s primordial magnetic field was probably “fossilised”
in the deep interior following the retreat of the convection zone (e.g. Parker,
1981). Since the global-scale magnetic diffusion time for the Sun is roughly 1010
years, comparable to the Sun’s present age, it is likely that some global-scale
remnant of this field Bi still exists within the stably stratified interior, with a
field strength of perhaps a few hundred gauss (e.g. Mestel & Weiss, 1987, & refs).
However, only torque-free magnetic field configurations can persist on such long
timescales (e.g. Chandrasekhar & Prendergast, 1956), so the field configuration
in the Sun’s stably stratified interior cannot be deduced from magnetic diffusion
alone.
Recently, Braithwaite & Spruit (2004) numerically modelled the evolution of fos-
35
36 3. The Magnetic Interior
sil fields in Ap stars (see also Braithwaite & Nordlund, 2006). By allowing the
internal field to evolve from a random initial condition, they found that the star
established a roughly torque-free equilibrium on the Alfve´nic timescale — i.e. the
global travel time of an Alfve´n wave through the interior. The equilibrium config-
uration for the magnetic field was similar to that suggested by Prendergast (1956),
a roughly axisymmetric dipole with a toroidal component in a neighbourhood of
the neutral ring (see figure 3.1). The presence of the toroidal component is neces-
sary, since any purely poloidal field configuration would be unstable (Braithwaite
& Nordlund, 2006, & refs).
Figure 3.1: A possible configuration for the fossil magnetic field Bi inside an Ap star (Braith-
waite & Spruit, 2004). The field is a roughly axisymmetric dipole. The field in the insulating
region surrounding the star is purely poloidal. Below the surface of the star, the neutral ring
of the dipole is stabilised by a band of toroidal field, shown in blue. The right frame shows a
simplified schematic of the field configuration.
In the following sections we suppose that the Sun’s interior fossil field has a
similar configuration to that found by Braithwaite & Spruit, and describe the
effect of this field on angular momentum transport within the Sun.
3.2 The Ferraro constraint
As described in §1.3, in a system with a strong magnetic field the rotation must be
constant on each “Ferraro surface”, i.e. on the surface of revolution of each field
line. The presence of a magnetic field Bi within the Sun’s radiative interior can
3.2 The Ferraro constraint 37
therefore provide a natural explanation for the uniform rotation of that region.
Any initial deviation from a Ferraro state will deform the field lines, producing
Alfve´n waves. The simplest model of this process supposes that the field is
axisymmetric, and that the fluid motions are purely rotational, as in the model
of Charbonneau & MacGregor (recall §1.4.4). If we neglect viscosity and magnetic
diffusion, then each “Ferraro surface” oscillates independently of the others, and
with its own discrete set of harmonic frequencies, determined by the poloidal
Alfve´n speed within the surface (Plumpton & Ferraro, 1955). These global Alfve´n
modes, known as “torsional oscillations”,1 were first studied by Wale´n (1946) as
a possible explanation for the 22-year solar magnetic cycle (§2.5).
Since each Ferraro surface oscillates independently of its neighbours, the oscilla-
tions of neighbouring surfaces are generally not in phase, leading to large trans-
verse gradients in the azimuthal velocity and magnetic fields; this is known as
“phase mixing” (Charbonneau & MacGregor, 1993). Therefore the presence of
weak magnetic diffusion rapidly damps the oscillations, leading to a steady Fer-
raro state. In the absence of viscosity, there can be no angular momentum trans-
port between neighbouring Ferraro surfaces, so in general the steady state will be
differentially rotating. However, the system may be subject to various hydrody-
namic or MHD instabilities (see references in §2.4), leading to enhanced diffusion
of both magnetic field and angular momentum. Allowing for finite meridional
motions also introduces angular momentum coupling between the Ferraro sur-
faces, and may further enhance the damping of torsional oscillations (e.g. Mestel
& Weiss, 1987).
The picture just described applies if the magnetic field Bi is axisymmetric. For
a non-axisymmetric field, only uniform rotation is permitted in the steady state.
However, under some circumstances we would expect the field to be very nearly
axisymmetric. For example, if the field were initially quite weak, a few gauss or
less, and non-axisymmetric, then the action of differential rotation would wind up
a predominantly toroidal field with large meridional gradients in the toroidal field
strength. Magnetic diffusion would then rapidly eliminate the non-axisymmetric
component of the field; this process is known as “rotational smoothing” (Spruit,
1999, & refs).
1The torsional oscillations of the interior magnetic field should not be confused with the
torsional oscillations of the convection zone described in §2.3.
38 3. The Magnetic Interior
3.3 The magnetic confinement problem
As mentioned in §2.2, the early measurements made using helioseismology sug-
gested that the Sun’s core rotates at least twice as fast as the radiative envelope.
Mestel & Weiss (1987) argued, based on the picture put forward in §3.2, that the
presence of a magnetic field in the radiative interior is incompatible with such
rapid core-rotation, unless the magnetic field is very weak (10−2 gauss or less).
They suggested that the helioseismological measurements were in error, and that
in fact the entire region below the convection zone is locked into uniform rotation
by a global-scale fossil magnetic field Bi. They assumed that the magnetic field
would be confined to the radiative interior by “magnetic flux expulsion” from the
convection zone (Weiss, 1966). This mechanism is described in §4.1.
Ru¨diger & Kitchatinov (1997) constructed the first quantitative model of the in-
teraction between a magnetic field in the radiative envelope and the differential
rotation of the convection zone. Their model was similar to the spin-down model
of Charbonneau & MacGregor (1993), except that the convection zone’s differ-
ential rotation was directly imposed at the top of the radiative envelope. The
poloidal component of Bi, which they took to be an axisymmetric dipole, was im-
posed throughout the calculation, and meridional flows were forbidden. Ru¨diger
& Kitchatinov assumed that the turbulence in the convection zone would act as
a large eddy diffusivity, and thereby expel any magnetic field from that region
(see §4.1). They therefore only considered magnetic field configurations with no
direct magnetic coupling between the convection zone and radiative envelope. In
their simulations, the convection zone’s differential rotation spread into the ra-
diative envelope by viscous diffusion, and then proceeded to wind up the interior
field. A balance was achieved once the Lorentz torque from the magnetic field
became large enough to balance the viscous torque. In cases where the imposed
poloidal field was sufficiently strong (|Bi| & 10−4 gauss), the differential rotation
within the radiative envelope was confined to a thin shear-layer resembling the
tachocline. This shear-layer can be described by a simple boundary-layer theory
(Ru¨diger & Kitchatinov, 1997) in which we approximate the density ρ, viscosity
ν and magnetic diffusivity η as constant. Since the density is assumed constant
(ρ ≈ 0.21g cm−3, see appendix A), we may measure the magnetic field B in
units of Alfve´n speed, with 1 gauss corresponding to approximately 0.6cm s−1.
3.3 The magnetic confinement problem 39
The governing equations are the azimuthal components of the momentum and
induction equations, which can be expressed in cylindrical polar coordinates as
0 = 1rBi ·∇(rBφ) + ν
(
∇2 − 1r2
)
uφ , (3.1)
0 = rBi ·∇(uφ/r) + η
(
∇2 − 1r2
)
Bφ . (3.2)
Away from the poles, the dipolar field Bi is roughly horizontal below the base of
the convection zone. We therefore expect that
Bi ·∇ ∼ |Bi|/L , (3.3)
where L is the horizontal scale of the imposed differential rotation, which we take
to be the radius of the tachocline, 0.7R. Assuming that the thickness of the
steady-state shear-layer is ∆  L, we also expect that
∇2 ∼ 1/∆2 . (3.4)
Hence it is readily verified from (3.1) and (3.2) that
∆ ∼
(νηL2
|Bi|2
)1/4
, (3.5)
Bφ ∼ (ν/η)1/2uφ . (3.6)
So the thickness of the shear-layer decreases as the strength of the poloidal field
increases, and becomes . 0.05R when
|Bi| &
νη(0.7R)2
(0.05R)4
(3.7)
≈ 4× 10−7cm s−1. (3.8)
So a poloidal field with |Bi| & 10−6 gauss should be sufficient to confine the con-
vection zone’s differential rotation within a layer as thin as the tachocline. From
(3.6) we see that the magnitude of the toroidal field wound up in the shear-layer
is essentially independent of |Bi|. The typical magnitude of uφ is determined by
the imposed differential rotation at the top of the radiative envelope, which is
40 3. The Magnetic Interior
about 30% of the interior rotation rate Ωi, and so
uφ ∼ 0.3ΩiL . (3.9)
This leads to an estimate for Bφ of around 104 gauss — considerably larger than
the lower bound on |Bi| found above. However, a strongly toroidal field would
probably be unstable (see appendix F).
The model of Ru¨diger & Kitchatinov demonstrates that an interior magnetic field
can prevent the viscous spread of the convection zone’s differential rotation into
the interior. However, the radiative spreading process described in §2.4 is the
result of meridional flows seeking to establish a Taylor–Proudman state, and has
no connection with viscous diffusion. In the model of Ru¨diger & Kitchatinov
meridional flows are artificially suppressed, and so Coriolis effects are absent.
We might say that their model describes a non-rotating Sun. Nonetheless, their
model provides a simple illustration of how an interior magnetic field is able to
enforce uniform rotation, provided that the field is confined (i.e. predominantly
horizontal) beneath the convection zone.
Subsequently, MacGregor & Charbonneau (1999) considered a very similar model,
but where now the interior magnetic field Bi was unconfined, i.e. directly coupled
to the convection zone. They found that the convection zone’s differential rota-
tion propagated along the poloidal field lines, rapidly establishing a Ferraro state
of differential rotation within the interior (see figure 3.2). This demonstrates that
the geometry of the internal field is crucial to determining the pattern of rotation
in the interior. In order to explain the sharp transition between uniform and
differential rotation at the base of the convection zone, the magnetic field lines
must remain nearly horizontal within the outer part of the radiative envelope.
But even if the magnetic field is initially confined in this fashion, over the Sun’s
lifetime it can diffuse into the convection zone, whereupon the convection zone’s
differential rotation will be rapidly transferred to the interior. This can be called
the “magnetic confinement problem”, and has been illustrated by the recent
numerical simulations of Brun & Zahn (2006). They modelled the evolution
of a dipolar fossil magnetic field within the radiative envelope, using a three-
dimensional, anelastic numerical scheme. They modelled the radiative–convective
3.3 The magnetic confinement problem 41
Figure 3.2: The imposed poloidal magnetic field (left panel) and steady-state angular velocity
(right panel) in the model of MacGregor & Charbonneau (1999).
interface as a differentially rotating, impenetrable boundary; the differential ro-
tation imposed on this boundary was chosen to represent the differential rotation
of the lower convection zone. They found that, as the interior field diffused, it
connected to the differentially rotating outer boundary, and communicated this
differential rotation into the radiative envelope. This problem was particularly
conspicuous in high latitudes, where the convection zone’s slow rotation was
transmitted to an “apple-core” region surrounding the axis (see figure 3.3).
In order for the interior magnetic field to enforce uniform rotation within the
radiative envelope, a mechanism must be in place to confine the magnetic field,
holding the field lines horizontal below the base of the convection zone. Possible
mechanisms for achieving magnetic confinement are discussed in the next chapter.
42 3. The Magnetic Interior
Figure 3.3: If the interior magnetic field Bi (black lines) is allowed to diffuse into the convec-
tion zone, then the convection zone’s differential rotation is imprinted on the interior (Brun &
Zahn, 2006). The colours represent angular velocity (units are nHz).
Chapter 4
Magnetic Field Confinement
4.1 Magnetic pumping
The studies of Mestel & Weiss (1987) and Ru¨diger & Kitchatinov (1997) both
assumed that the Sun’s interior magnetic field would be confined to the radiative
interior by the turbulence in the convection zone. They argued that any magnetic
field in the convection zone would be characterised by the same small lengthscales
and timescales as the convective turbulence, and that therefore the time-averaged
mean field would reside entirely within the stably stratified interior.
A simple example of the interaction between convective turbulence and magnetic
fields was considered by Weiss (1966). He showed that turbulent eddies enhance
the rate of dissipation of magnetic fields, so that magnetic flux becomes concen-
trated at the boundaries between the eddies. This is often called “magnetic flux
expulsion”, i.e. expulsion of magnetic flux from the most turbulent regions of the
flow.
A related but separate concept is that of “topological pumping”, first proposed
by Drobyshevski & Yuferev (1974). Based on the structure of solar convection —
concentrated downwelling lanes separating weak upflowing regions — they argued
that horizontal flux tubes in the convection zone would be carried downward by
the convection (see also Tao et al., 1998, & refs). Recently, Tobias et al. (1998)
studied the interaction between turbulent compressible convection and a hori-
43
44 4. Magnetic Field Confinement
zontal layer of magnetic field (see also Dorch & Nordlund, 2001, & refs). They
performed cartesian-box simulations of a convectively unstable region overlying
a stably stratified region. Once the convection was fully developed, a horizontal
layer of magnetic field was inserted into the unstable region. They found that
much of this field was pumped into the stably stratified region by overshooting
convective plumes. They also found that this “magnetic flux pumping” operates
only when the convective turbulence is significantly compressible, and therefore
vertically asymmetric (recall §2.3.1). However, it is unclear from their simula-
tions whether the pumping results simply from the vertical anisotropy of the
turbulence, or whether the presence of the stably stratified layer is also required.
The concepts of flux expulsion and flux pumping are also present in mean-field
theory. In that context, mean-field transport is attributed to the antisymmet-
ric component of the “α-tensor” (e.g. Roberts & Soward, 1975), which can be
expressed in terms of an effective pumping velocity γ (Ra¨dler, 1968). Pumping
can arise from either inhomogeneous turbulence (so-called “diamagnetic pump-
ing”) or density gradients (e.g. Kitchatinov & Ru¨diger, 1992). In fact, inhomoge-
neous turbulence of itself cannot produce pumping; the turbulence must also be
anisotropic. Within the context of mean-field theory it is therefore straightfor-
ward to make a distinction between the flux expulsion described by Weiss (1966),
which occurs even for isotropic turbulence, and the flux pumping hypothesised
by Drobyshevski & Yuferev (1974). In numerical simulations this distinction may
be less clear.
Recently, Garaud & Rogers (2007) performed a numerical simulation of the inter-
action between a magnetic field Bi in the solar interior and the turbulence of the
convection zone. Their model was of a two-dimensional (non-rotating) “cylindri-
cal Sun” (see figure 4.1). They found that, as the initially confined magnetic field
diffused into the convection zone, it was quickly stretched and advected by the
turbulent convective eddies, so that no large-scale component survived within the
convection zone. However, the structure of the field within the radiative envelope
was similar to that found by Brun & Zahn (2006, see figure 3.3). A snapshot of
the magnetic field is shown in figure 4.1. After less than one global-scale magnetic
diffusion time, significant field-line connection is established with the convection
zone in the polar regions, although the field lines remain almost horizontal in
4.1 Magnetic pumping 45
middle and low latitudes. There is no indication of magnetic field being pumped
Figure 4.1: A cross-section through the “cylindrical Sun” of Garaud & Rogers (2007). Tur-
bulent eddies destroy any large-scale poloidal magnetic field within the convection zone.
out of the convection zone back into the radiative envelope. This may be because
the turbulence is not sufficiently compressible; the parameter regime described
by Tobias et al. (1998) cannot be matched in a global model.
Garaud & Rogers hypothesised that, although the interior magnetic field Bi is
able to connect to the convection zone, since the field lines follow an essentially
random path through that region they do not transmit any net torque to the
interior. However, since the effects of rotation could not be included in their two-
dimensional model, it was not possible to directly test this hypothesis. By analogy
with the results of Brun & Zahn (2006), it seems likely that the slow rotation
of the high-latitude convection zone would be communicated to an apple-core
region surrounding the axis in the radiative interior.
The question of whether magnetic flux pumping might act to confine the Sun’s
interior magnetic field has also been considered by Kitchatinov & Ru¨diger (2008),
in the context of an axisymmetric mean-field model. They parametrised the tur-
bulence in the convection zone using a combination of large turbulent magnetic
diffusivity and diamagnetic pumping. They found that the interior magnetic
field Bi could be confined, provided that the diamagnetic pumping was suffi-
46 4. Magnetic Field Confinement
ciently strong. The model of Kitchatinov & Ru¨diger is re-derived (and slightly
generalised) in appendix E. The main result is that strong downward pumping
at the base of the convection zone leads to an accumulation of horizontal mag-
netic field at the top of the radiative envelope. This result is at least superficially
similar to that of Tobias et al. (1998), although as mentioned above the same
behaviour was not observed in the model of Garaud & Rogers (2007).
Whether such strong diamagnetic pumping operates at all latitudes within the
convection zone, and whether it is directly related to the magnetic flux pumping
found by Tobias et al., remains unclear. Intuitively, we would expect pumping
to be most effective in middle and low latitudes, where the interior field is hor-
izontal, as in the model of Tobias et al. (1998). The pumping efficiency might
be further reduced in high latitudes by Coriolis effects (Tobias et al., 2001). We
might therefore tentatively conclude that diamagnetic pumping can confine the
horizontal field in middle and low latitudes, but that another mechanism may be
required in high latitudes. The importance of confining the magnetic field in the
high latitude regions is highlighted by the results of Brun & Zahn (figure 3.3).
4.2 Mean meridional circulations (MMCs)
4.2.1 The Gough & McIntyre (GM98) model
The problem of magnetic field confinement was first discussed by Gough & McIn-
tyre (1998, hereafter GM98). They proposed that the Sun’s fossil magnetic field
Bi is confined to the radiative interior by the tachocline’s MMCs, which hold the
field in the tachocline in advection–diffusion balance; that is, the upward diffusion
of the field is balanced by downward advection by the tachocline’s MMCs. As dis-
cussed in §2.3.2, they argued that the pattern of differential rotation observed in
the convection zone and tachocline imply the existence of persistent downwelling
in high latitudes. They estimated the magnitude of the downwelling velocity U
to be ∼ 10−5cm s−1 in the stably stratified polar tachocline. Since the magnetic
diffusivity η in the tachocline is ≈ 410cm2s−1 (see appendix A), the magnetic
advection–diffusion lengthscale δη with this velocity is η/U ≈ 0.4Mm, which is
much thinner than the tachocline. They reasoned that the high-latitude down-
4.2 Mean meridional circulations (MMCs) 47
welling therefore confines the interior magnetic field within a magnetic boundary
layer of thickness δη at the base of the high-latitude tachocline. Henceforth we
will refer to this boundary layer as the “magnetic confinement layer”.
The GM98 model comprises three distinct regions:
• Above the magnetic confinement layer, the tachocline is laminar,1 field-free
and in thermal-wind balance.
• Within the magnetic confinement layer, the poloidal field is held in advective–
diffusive balance. A toroidal field component, wound up by differential rota-
tion, exerts a prograde Lorentz torque on the fluid, turning the downwelling
flow equatorward, and thereby confining the meridional circulation.
• Below the confinement layer, the radiative envelope is stagnant and uni-
formly rotating, held in isorotation by the interior magnetic field Bi, which
GM98 assumed to be an axial dipole.
GM98 did not provide a complete model of magnetic confinement. In fact, they
modelled only a localised, high-latitude region, sufficiently far from the pole that
curvature terms could be neglected and the interior magnetic field assumed hori-
zontal. Moreover, they modelled only the deepest extent of the confinement layer,
which can be approximated as a linear perturbation to the stagnant, uniformly
rotating interior. Their model can therefore be described in terms of the MAC-
wave theory set out in appendix C. Following GM98 we neglect viscosity and
take the interior magnetic field Bi to be uniform. The steady-state temperature
perturbation T ′ is then described by the equation
[
(Bi ·∇)4 + η2(2Ωi ·∇)2∇2
]
∇2T ′ − ηN
2
κ (Bi ·∇)
2∇2HT ′ = 0 (4.1)
(cf. (C.20)–(C.21) in appendix C). Here, Ωi is the constant angular velocity
of the interior, κ is the thermal diffusivity, N is the buoyancy frequency and
∇2H is the horizontal Laplacian. Again following GM98, we assume that the
system is axisymmetric, and that there is a single latitudinal lengthscale, say L,
1Recently, McIntyre (2007) presented a significantly revised version of the GM98 model in
which the bulk of the tachocline is subject to Tayler instability — see appendix G.
48 4. Magnetic Field Confinement
which is determined by the global geometry of the tachocline (GM98 estimated
L ≈ 0.15R ≈ 100Mm). We therefore replace ∇2H → −1/L2 in (4.1). Since the
interior field Bi is assumed horizontal, we may also replace (Bi·∇)2 → −|Bi|2/L2.
For convenience we choose to measure the strength of the field in terms of the
“horizontal” Elsasser number
ΛH =
|Bi|2
2Ωiη
, (4.2)
which is a dimensionless number. If ΛH = 1 then |Bi| ≈ 8 × 10−2 gauss. Since
δη  L, we may make the boundary-layer approximation ∇2 ≈ ∂2/∂z2, where z
is the vertical coordinate, so that (4.1) becomes
(
ΛH +
L4
ΛH
∂4
∂z4
) ∂2
∂z2T
′ = N
2
2Ωiκ
T ′ . (4.3)
In GM98 it was assumed that thermal-wind balance is not modified by the Lorentz
force in the confinement layer. As a result of this assumption, rather than (4.3)
they obtained the equation
L4
ΛH
∂6
∂z6T
′ = N
2
2Ωiκ
T ′ . (4.4)
This equation has three independent solutions that decay with depth, and three
that grow with depth; the second three violate the assumptions made in the
linearisation, and are therefore rejected. The remaining three solutions all decay
exponentially with depth on the vertical lengthscale
( 2Ωiκ
ΛHN2L2
)1/6
L ≈ Λ−1/6H Mm . (4.5)
By matching this vertical lengthscale to the magnetic advection–diffusion length-
scale δη = η/U , GM98 were able to directly relate the downwelling velocity U to
the strength of the interior magnetic field, which here is quantified by ΛH . Using
their estimates U ≈ 10−5cm s−1 and L ≈ 100Mm they deduced that |Bi| ≈ 1
gauss.
4.2 Mean meridional circulations (MMCs) 49
By comparing (4.3) and (4.4), we find that the results of GM98 are valid only if
Λ2H 
N2L2
2Ωiκ
≈ 1012. (4.6)
This leads to an upper bound on |Bi| of approximately 50cm s−1, or 80 gauss.
The results of GM98 are therefore valid, although the margin of error is rather
small. In fact, the full equation (4.3) has the same essential property as (4.4): it
has three independent solutions that decay with depth, and three that grow with
depth. The GM98 picture is therefore qualitatively unaffected by the extra term
in (4.3).
4.2.2 The polar regions
In GM98 it was assumed that, in order to confine the tachocline’s MMCs, the
interior magnetic field Bi needed to be accurately horizontal. Near the poles,
where the assumed dipolar field is not horizontal, they speculated that the MMCs
would burrow into the radiative envelope, forming “polar pits” (see figure 4.2).
Such polar pits would have significant implications for lithium burning, as will
be described in chapter 6.
Figure 4.2: An interior dipole Bi (red lines) confined by the tachocline’s MMCs (black lines
with arrows). Near the poles, the MMCs burrow into the radiative envelope, forming polar pits.
Adapted from Gough & McIntyre (1998). The thickness of the tachocline (shown in green) has
been exaggerated.
50 4. Magnetic Field Confinement
To go beyond these speculations, we must formulate a complete and self-consistent
model of the magnetic confinement layer that includes the polar regions. This
is the ultimate goal of this dissertation. As a first step towards achieving this
goal, it is instructive to consider how the regime described by (4.3) changes as
we move closer to the pole. The first of GM98’s approximations to break down
will be that of horizontal magnetic field. Once the slope of the interior field lines
becomes comparable to the aspect ratio of the confinement layer, it is no longer
valid to neglect the vertical component of Bi. Indeed, within some neighbour-
hood of the pole Bi must be approximately vertical. Within this neighbourhood
we may approximate (Bi ·∇)2 ≈ B2iz∂2/∂z2, where Biz is the vertical component
of Bi. Equation (4.1) then becomes
∂2
∂z2
[(
ΛV
∂2
∂z2 +
1
ΛV
∇2
)
L2∇2T ′ + N
2
2Ωiκ
T ′
]
= 0 (4.7)
where we have replaced ∇2H → −1/L2 and defined a “vertical” Elsasser number
ΛV =
B2iz
2Ωiη
. (4.8)
Equation (4.7) supports two independent axisymmetric solutions with non-exponential
dependence on z. One corresponds to Ferraro differential rotation; the other cor-
responds to vertical flows along the field lines, the signature of GM98’s polar
pits. The flow perturbs the thermal stratification surfaces, thereby establishing
a thermal wind in the radiative envelope. In the laminar, inviscid system de-
scribed by equation (4.7) there is nothing that can limit the degree of Ferraro
differential rotation, nor the depth to which the flow along the field lines can
burrow. However, as discussed in §3.2, in the real Sun any differential rotation
in the polar regions would be limited by hydrodynamic and MHD instabilities,
which introduce angular momentum coupling between the Ferraro surfaces, and
allow the flow to leak across the field lines. We would therefore expect the polar
pits, and their associated thermal winds, to extend only a finite distance into the
radiative envelope,2 as indicated in figure 4.2.
2Since helioseismology cannot provide an accurate measurement of the Sun’s angular velocity
profile close to the rotation axis (see figure 2.1), differential rotation within this region cannot
be ruled out by present solar observations.
4.2 Mean meridional circulations (MMCs) 51
In addition to the two solutions just described, equation (4.7) supports four
boundary-layer (∇2 ≈ ∂2/∂z2  1/L2) solutions, which are described by the
equation
(
ΛV + Λ−1V
)
L2 ∂
4
∂z4T
′ + N
2
2Ωiκ
T ′ = 0 . (4.9)
As in GM98, solutions that grow exponentially with depth are not physically rel-
evant. This leaves two independent solutions, both of which decay exponentially
with depth on the vertical lengthscale
( 2Ωiκ
N2L2 (ΛV + Λ
−1
V )
)1/4
L . (4.10)
Equating this to the lengthscale δη, and again assuming that U ≈ 10−5cm s−1, we
deduce that (ΛV + Λ−1V ) ≈ 240. This implies that there are two possible values
for the vertical field strength |Biz|: either |Biz| ∼ 1 gauss or |Biz| ∼ 10−3 gauss.
4.2.3 Gyroscopic pumping
GM98 gave a qualitative description of the “gyroscopic pumping” process that
drives the downwelling in the high-latitude tachocline. Their argument can be
summarised as follows:
1. The turbulent stresses within the convection zone3 impart a systematic
retrograde forcing at high latitudes, causing that region to rotate slower
than the rest of the convection zone.
2. The confined interior magnetic field Bi holds the radiative envelope in uni-
form rotation, and thus constrains the high-latitude tachocline to rotate
faster than the layers above.
3. The retrograde forcing at the base of the high-latitude convection zone,
unable to spin down this region, instead gyroscopically drives a poleward
flow.
3And also any turbulent stresses within the tachocline — see McIntyre (2007) and ap-
pendix G.
52 4. Magnetic Field Confinement
4. The converging poleward flow establishes a meridional circulation within
the tachocline, which tries to burrow into the radiative envelope in the
manner described by Haynes et al. (1991) and Spiegel & Zahn (1992). The
resulting MMC is downwelling in the polar regions (see figure 2.3).
A similar process pumps the “spin-down currents” in Spiegel’s tachycline model
(§1.4.2). An analogy can also be made with the “mechanical forcing” described
by Miesch et al. (2006), but the gyroscopic pumping described by GM98 persists
long after the system has established thermal-wind balance. Some more recent
studies of gyroscopically pumped meridional circulations are discussed in §4.2.6.
The gyroscopic pumping described by GM98 relies on the interior magnetic field
being confined below the base of the convection zone. If, instead, the field lines
protrude into the convection zone in high latitudes, then the retrograde forcing
in that region will be transmitted directly into the radiative envelope along the
unconfined field lines, as in the model of Brun & Zahn (2006).
4.2.4 Global field confinement
GM98 modelled only linear perturbations to the interior magnetic field Bi, and
therefore did not explicitly demonstrate field confinement. Subsequent studies by
Garaud (2002) and Garaud & Garaud (2008) have sought to achieve field con-
finement by MMCs in a global, nonlinear, axisymmetric model of the radiative
interior. Garaud & Garaud (2008) also offered a two-fold explanation for the lack
of field-confinement in the model of Brun & Zahn (2006). Firstly, by imposing a
no-penetration boundary condition at the top of the radiative envelope, Brun &
Zahn inhibited the generation of MMCs in their model. Secondly, the viscosity
used by Brun & Zahn was too large in relation to the thermal diffusivity, signifi-
cantly weakening the burrowing tendency of the MMCs. In a field-free analogue
of the Brun & Zahn model, radiative spreading would be halted by viscosity be-
fore the MMCs could burrow more than ∼ 0.2R into the radiative envelope. As
a result, the MMCs in Brun & Zahn’s model were too weak to significantly offset
the outward diffusion of the magnetic field.
In Garaud & Garaud’s model, differential rotation and meridional flows are im-
4.2 Mean meridional circulations (MMCs) 53
posed at the top of the radiative envelope. They found, as expected, that the
MMCs in the radiative envelope distorted the magnetic field lines on lengthscales
comparable to the magnetic advection–diffusion lengthscale η/U , where U is a
typical value of the imposed downwelling at the top of the radiative envelope. In
particular, imposing downwelling in high latitudes was found to reduce the mag-
netic coupling between the interior field Bi and the differentially rotating outer
boundary. However, since they only considered cases where the downwelling van-
ished at the poles, the magnetic field in all cases was found to be unconfined in
the polar regions (see figure 4.3).
Figure 4.3: Left panel: An interior dipole field Bi partially confined by MMCs. Right panel:
The steady-state angular velocity; the colour axis runs from 278 nHz (black) to 463 nHz (red).
Taken from Garaud & Garaud (2008).
The model of Garaud & Garaud (2008) employs a Newton-Raphson-Kantorovich
(NRK) numerical algorithm to find solutions of the steady-state equations, rather
that evolving the time-dependent equations. With this method, convergence to
the steady state is usually more rapid than with a time-dependent algorithm.
However, this method provides no information regarding the dynamical conver-
gence to the steady state, and failure to find a steady state does not imply that no
steady state exists. Furthermore, any steady state found will be highly sensitive
to the choice of boundary conditions, and may also depend on the ansatz used
54 4. Magnetic Field Confinement
as a starting point for the algorithm. The discussion of gyroscopic pumping in
§4.2.3 suggests that two distinct steady states may exist in high latitudes: one
with a magnetic field confined by MMCs, and one with an unconfined field that
short-circuits the gyroscopic pumping of MMCs.
4.2.5 Burrowing
Recently, a number of authors have claimed that the tachocline’s MMCs are
unable to burrow significantly into the radiative envelope, even in the absence
of a magnetic field (e.g. Gilman & Miesch, 2004; Ru¨diger et al., 2005), and that
the tachocline can therefore be explained simply as an Ekman layer.4 Were these
claims true, it would undermine the majority of the work that has been reviewed
in the preceding chapters. In particular, it would invalidate the model of Spiegel
& Zahn (1992), which is the starting point for most models of the tachocline.
Detailed criticisms of these claims have been made elsewhere (e.g. McIntyre, 2007;
Garaud & Brummell, 2008) and will not be repeated here. Instead, two recent
studies of the gyroscopic pumping of MMCs by Garaud & Brummell (2008) and
Garaud & Acevedo Arreguin (2009) are summarised in §4.2.6.
However, based on the assumption that MMCs driven in the convection zone and
tachocline cannot penetrate deeper into the radiative envelope than the Ekman
length δEk = (ν/2Ωi)1/2, Kitchatinov & Ru¨diger (2006) have constructed a model
for the interaction between the tachocline’s MMCs and the fossil field Bi in
the radiative envelope. Using the tachocline’s microscopic viscosity value ν ≈
30cm2s−1 (see appendix A), the Ekman scale is δEk ≈ 3× 103 cm ∼ 10−8R. It
might be thought that such a small-scale flow would have negligible impact on
the diffusion of the interior field. In particular, since the microscopic magnetic
diffusivity η in the tachocline is larger than the viscosity ν by a factor η/ν ≈ 15,
the magnetic field diffuses across a layer of thickness δEk in much less than one
solar rotation period. However, if the horizontal shear flow within the Ekman
layer is sufficiently large, then it can “kink” the field lines that cross it. Indeed,
if the Ekman layer is assumed to be infinitely thin, so that its latitudinal flow
has the form of a delta function in radius, then the flow will produce a vertical
4Or Ekman–Hartmann layer if a magnetic field is present.
4.2 Mean meridional circulations (MMCs) 55
discontinuity in the latitudinal magnetic field. This is the situation considered
by Kitchatinov & Ru¨diger. In their model, the turbulence in the convection zone
is assumed to act in the manner of a large magnetic diffusivity, so that the field
outside the radiative envelope can be approximated as a vacuum field, for which
the latitudinal and vertical components are necessarily of similar magnitude. By
imposing an equatorward flow ∼ 103cm s−1 in a layer of thickness ∼ δEk at the
top of the radiative envelope, they found that the magnetic field lines crossing
the layer were sufficiently kinked that the field lines near the top of the radiative
envelope were nearly horizontal. In this way, the interior field could be confined
to the radiative envelope.
The shear associated with Kitchatinov & Ru¨diger’s equatorward flow is ∼ 1s−1,
which implies a gradient Richardson number less than 10−5. It is unlikely that
such an extreme shear flow would remain stable within the radiative envelope.
However, the model relies on the shear layer remaining laminar, or at most weakly
turbulent, so that only microscopic magnetic diffusion is present within this layer.
A larger magnetic diffusivity would not permit the required kinking of the mag-
netic field lines. Yet the magnetic diffusivity is assumed to be infinite in the
convection zone immediately above the shear layer, in order to prevent large
horizontal magnetic fields from accumulating in that region.
It should also be noted that the latitudinal flow assumed by Kitchatinov &
Ru¨diger is artificially imposed (via mass sources at the pole and correspond-
ing sinks at the equator) and is not derived within the model. Furthermore, the
Coriolis force arising from this flow is neglected, although it would be expected to
make a leading-order contribution to the total force balance. The Lorentz forces
arising from the kinking of the field lines in the Ekman layer are also not included
in the model.
4.2.6 Recent studies of gyroscopic pumping
The radiative spreading, or burrowing, described by Spiegel & Zahn (1992, see
§2.4) arises because of the departures from geostrophic balance implied by the
tachocline’s rotational shear. Such departures imply, via thermal-wind balance,
a horizontal temperature gradient, which must be maintained against thermal
56 4. Magnetic Field Confinement
relaxation by the advection of entropy by meridional circulations. In the absence
of viscous or magnetic torques, the Coriolis torque on the meridional circulation
cannot be balanced, and as a result the meridional circulations burrow into the
interior, spreading the differential rotation of the convection zone. Radiative
spreading is essentially a hyperdiffusive process, and so the tachocline in Spiegel
& Zahn’s model thickens as t1/4. More precisely, if the rotation rate in the tacho-
cline varies on a horizontal scale of order L, say, then the tachocline thickness ∆
grows like
∆ ∼ L(t/τ)1/4 (4.11)
where
τ =
( N
2Ωi
)2 L2
κ , (4.12)
cf. equation (1.5). Here, Ωi is the uniform rotation rate below the tachocline.
As radiative spreading brings the radiative interior progressively closer to a
geostrophically balanced Taylor–Proudman state, so the meridional flow speed
decreases. Eventually the flow would shut off completely, but only after a time
 τ , which is much longer than the age of the Sun.
In the presence of some additional physics, such as viscosity or magnetic fields,
that is capable of supporting a finite torque, the burrowing process can be halted
sooner, and finite meridional circulations can persist even in the steady state. In
the absence of radiative spreading, the presence of viscosity, ν, would itself cause
the tachocline thickness to grow like
∆ ∼ (νt)1/2. (4.13)
Equating (4.11) and (4.13) we deduce that viscous diffusion overtakes radiative
spreading, thereby shutting off the burrowing of the meridional circulations, when
∆ ∼ L
2
(ντ)1/2 (4.14)
∼
(κ
ν
)1/2
(2Ω
N
)
L (4.15)
(see Clark, 1973). In oceanography, the lengthscale (4.15) is usually called the
4.2 Mean meridional circulations (MMCs) 57
Lineykin length,5 after Lineykin (1955).
In a recent study, Garaud & Brummell (2008, hereafter GB08) have analysed the
steady-state structure of meridional flows in the radiative interior in the presence
of viscosity. Their model is similar to that of Spiegel & Zahn in that the flows,
which are assumed to be linear perturbations to a hydrostatic, stably stratified
background, are driven by boundary conditions imposed at the interface between
the radiative and convective zones. The particular sets of boundary conditions
they consider differ from those of Spiegel & Zahn, however. Specifically, they
impose the three components of velocity at the top of the radiative envelope,
together with one of two choices for the thermal boundary conditions. The first
choice imposes that the temperature perturbation has no derivative normal to
the radiative–convective interface,6 and the second choice imposes that the per-
turbations do not contribute to the total (advective plus diffusive) flux of heat
across the radiative–convective interface.
GB08 first consider a cartesian model, with the rotation axis aligned with gravity,
and with constant buoyancy frequency N and constant angular velocity Ωi. In
this situation, the importance of stratification in shaping the steady state is
quantified by the dimensionless parameter σ, say, where
σ =
(ν
κ
)1/2 N
2Ωi
(4.16)
(e.g. Barcilon & Pedlosky, 1967). Assuming that all fields have the same sinu-
soidal horizontal dependence, GB08 show that any steady-state solution can be
expressed as a superposition of eight linearly independent modes. In the param-
eter regime appropriate to the solar interior, four of these modes have the small
vertical lengthscale characteristic of Ekman layers, δEk = (ν/2Ωi)1/2. Of the
remaining four modes, two have vertical scales equal to the imposed horizontal
scale, L say, and two (the Lineykin modes) have the vertical scale given by (4.15),
i.e. σ−1L. With parameter values appropriate to the tachocline (see appendix A),
we have σ−1 ≈ 5. The Lineykin modes therefore penetrate considerably deeper
5Garaud & Acevedo Arreguin (2009) refer to this as the “thermo-viscous” lengthscale. We
will not use this term here, in order to avoid confusion with the thermo-viscous lengthscale
found by Tassoul & Tassoul (1986) in a very similar problem.
6Although GB08 refer to this as a “perfectly conducting” boundary condition, in fact it
serves to fix the diffusive heat flux across the boundary.
58 4. Magnetic Field Confinement
than the Ekman modes for any L that is larger than δEk.
If the domain is assumed to be infinitely deep then, for all the sets of boundary
conditions GB08 considered, virtually all of the meridional flow in the steady
state resides within an Ekman layer below the radiative–convective interface.
They found that the deep meridional flows were strongest in cases where vertical
flows were imposed at the interface, and when the second choice of temperature
boundary conditions was employed. Even in these cases, the vertical flows be-
neath the Ekman layer were found to be weaker than the vertical flows imposed
at the interface by a factor
σH
L ≈
H
5L (4.17)
where L is again the imposed horizontal scale and H = g/N2 is the scale height of
the background potential temperature. Furthermore, mass conservation implies
that the latitudinal flow within the thin Ekman layer is much faster than the
vertical flow, by a factor L/δEk, whereas the latitudinal flow within the deep
interior is slightly slower than the vertical flow, by a factor σ ≈ 0.2.
If we take the limit ν → 0, holding all other quantities fixed, then σ → 0. In this
limit all of the downwelling meridional flow imposed at the radiative–convective
interface is diverted into an Ekman layer of vanishing thickness, and the flow
within the bulk of interior vanishes. This is at first surprising, as it seems to
suggest that the presence of viscosity is essential for driving deep meridional
flows, contrary to the results of Spiegel & Zahn (1992). This apparent paradox is
resolved when we remember that the analysis of GB08 applies only to the steady
state, in which the burrowing of meridional circulations has either extended to
infinite depth and then shut off, or else has been halted by viscosity.
If the domain is now assumed to be of finite depth, and if the lower boundary is
taken to be impenetrable, then there is a limit on the depth to which meridional
circulations can burrow, and so we would expect stronger circulations in the
steady state. GB08 verify that this is indeed the case. In particular, if the depth
of the domain is taken to be the same as the imposed horizontal scale, L, then
the vertical component of the deep meridional circulation increases by a factor of
σ−1 ≈ 5, and the latitudinal component increases by a factor of σ−2 ≈ 25.
GB08 then consider a more realistic model of the radiative interior, in axisymmet-
4.2 Mean meridional circulations (MMCs) 59
ric spherical geometry, with background thermodynamic profiles from a standard
one-dimensional solar model. The linearised equations are solved using the same
NRK algorithm used by Garaud & Garaud (2008, see §4.2.4). They find that
the scaling laws derived from the cartesian model agree very well with the results
of the more realistic model. The agreement is particularly good in cases where
the deep meridional circulations are weak. Physically, this is because only the
deep circulations are affected by the geometrical curvature and the non-uniform
background profiles used in the spherical model.
The results of GB08 demonstrate that imposing strong flows at the radiative–
convective interface does not guarantee the existence of strong steady-state merid-
ional circulations within the radiative interior. Indeed, in the absence of viscosity
ν, any meridional circulation would have an unbalanced Coriolis torque, and so
the steady-state meridional circulation must vanish in the limit ν → 0. The
choice of temperature boundary conditions is also found to significantly affect
the steady state.
However, GB08 note that by considering only the linearised steady-state equa-
tions they have overlooked many possible types of behaviour that might otherwise
arise. For example, sufficiently strong meridional circulations within the stably
stratified interior can overturn the stratification surfaces. This kind of nonlin-
ear behaviour is beyond the scope of the linearised equations, and might be
missed completely by a steady-state analysis. (For further discussion, see McIn-
tyre (2007).) Similarly, angular momentum transport by meridional circulations
can lead to nonlinear departures from an initial state of uniform rotation. In-
deed, in several of the cases presented by GB08 the steady state has regions of
counterrotation! Clearly such solutions should not be taken literally.
In a subsequent study, Garaud & Acevedo Arreguin (2009, hereafter GAA09)
have considered meridional flows driven by Reynolds stresses within the convec-
tion zone, rather than by conditions imposed at the radiative–convective interface.
They model the convection zone as an unstratified layer, with a large (turbulent)
thermal diffusivity, overlying the stably stratified radiative interior. Within the
convection zone they parametrise the turbulent Reynolds stresses in terms of
a “generalised Rayleigh friction” that drives the system towards a state with
prescribed differential rotation and vanishing meridional flow within the con-
60 4. Magnetic Field Confinement
vection zone (cf. appendix B, in which we assume a purely azimuthal friction).
This friction can be regarded as a generalisation of the Darcy friction used by
Bretherton & Spiegel (1968) to model Reynolds stresses in the convection zone.
In the model of GAA09, the azimuthal component of the friction gyroscopically
pumps meridional circulations in the convection zone, which penetrate into the
interior; in this respect their model has more in common with those of Haynes
et al. (1991) and McIntyre (2007) than with that of Spiegel & Zahn (1992).
As in GB08, they first consider a cartesian model, and then a more realistic
spherical model. However, due to the complexity of modelling the interface be-
tween the convection zone and the radiative interior, even the cartesian model
can only be solved analytically in certain limiting cases. The cartesian model is
further simplified by assuming that the differential rotation towards which the
convection zone is driven is linear in altitude and sinusoidal in latitude.
With realistically small values for the diffusivities, GAA09 find that meridional
circulations within the convection zone are typically stronger than any meridional
circulations within the deep radiative interior by many orders of magnitude. In
“weakly stratified” cases, defined by σ  1, the strength of the deep merid-
ional circulations is found to depend on the choice of boundary conditions at the
bottom of the cartesian domain. In particular, cases with a no-slip lower bound-
ary condition typically have significantly stronger circulations than cases with
a stress-free lower boundary condition. This behaviour is most evident in cases
with no stratification, i.e. σ = 0 (see figure 4.4). But even in these cases, the
circulations within the radiative interior are still much weaker than those within
the convection zone (again, see figure 4.4).
At first sight, it is perhaps surprising that the interior meridional circulations are
so weak even in the absence of stratification. However, it should be remembered
that such circulations arise because of departures from geostrophic balance (i.e.
departures from a Taylor–Proudman state). In figure 4.4 it is evident that the
steady state below the convection zone is very nearly a Taylor–Proudman state,
and therefore it should be no surprise that the steady-state meridional circula-
tions are weak in that region. It is also evident that the presence of a no-slip
lower boundary breaks the Taylor–Proudman constraint, and therefore allows for
stronger steady-state meridional circulations. GAA09 note that “the dynamics
4.2 Mean meridional circulations (MMCs) 61
Figure 4.4: Steady-state meridional circulations and differential rotation in an unstratified,
cartesian model of the solar interior (Garaud & Acevedo Arreguin, 2009). Colours represent
angular velocity, and solid/dotted lines represent anti-clockwise/clockwise meridional stream-
lines. The left and right panels show the cases with stress-free and no-slip lower boundaries
respectively. The convection zone is located in the region 0.7 6 z 6 1.0.
of the lower boundary layer entirely control the mass flux through the system”.7
In “strongly stratified” cases, i.e. cases with σ  1, the deep meridional circula-
tions are confined to a thin Lineykin boundary layer below the convection zone,
with an aspect ratio ∼ σ. Nonetheless, the steady-state rotation profile is still
found to be roughly constant on cylinders, as a result of viscous diffusion over
the entire domain. Since the meridional flows do not extend to the bottom of
the domain, the choice of lower boundary conditions has little influence on the
structure of the steady state in these cases.
GAA09 note that the parameter σ depends not only on the stratification strength
N but also on the rotation rate Ωi. So at the time of the Sun’s arrival on the
main sequence, when the Sun rotated more rapidly, the burrowing tendency would
have been stronger. We expect this statement to apply to the transient meridional
flows as well as to the steady-state meridional flows (§2.4).
Finally, GAA09 note that viscosity is not the only mechanism capable of over-
coming the Taylor–Proudman constraint and sustaining a finite torque. An in-
7This statement on first reading seems to directly contradict the principle of “downward
control” described by Haynes et al. (1991, see §1.2). However, downward control operates
only in non-Boussinesq systems, i.e. systems with meridional circulations encompassing several
density scale heights.
62 4. Magnetic Field Confinement
terior magnetic field can also serve these two purposes, and thereby prevent the
burrowing of meridional circulations while also permitting finite meridional circu-
lations to persist in the steady state, as described by Gough & McIntyre (1998).
However, a complete model of these processes would also need to describe the
confinement of the magnetic field, and so there is no direct analogy between the
models of Gough & McIntyre and GAA09.
4.3 Summary
The observations made in the preceding sections highlight the need for a fully
self-consistent model of magnetic confinement and angular momentum balance
within the tachocline. From the various studies of the solar interior reviewed in
chapters 2–4 we may draw the following conclusions:
• In order to confine the convection zone’s differential rotation, and maintain
uniform rotation within the radiative envelope, angular momentum must be
transported from the low-latitude tachocline to the high-latitude tachocline.
• An interior magnetic field Bi can naturally produce this angular momentum
transport, provided that the field is confined, i.e. provided that the field
lines remain nearly horizontal within the tachocline.
• If the interior magnetic field is assumed to be an axisymmetric dipole, then
magnetic flux pumping by overshooting convection may be able to confine
the field in middle and low latitudes.
• Field confinement in high latitudes can probably only be achieved by the
tachocline’s MMCs, which, fortunately, are expected to be downwelling at
those latitudes.
In Part III we construct a model of the magnetic confinement layer responsible
for confining the interior magnetic field in high latitudes. First, in Part II, we
discuss certain observations regarding the Sun’s composition, and the physical
implications of these observations. These discussions prove to have a significant
bearing on our confinement-layer model.
Part II
Solar Composition
63
Chapter 5
Helium Settling
5.1 The µ-choke
The molecular cloud from which the Sun formed comprised mostly hydrogen and
helium, plus traces of heavier elements1 (collectively known as “metals”). When
the Sun first formed, its composition by mass was roughly 72% hydrogen, 27%
helium and 1% metal (e.g. Asplund et al., 2009).
During its pre-main-sequence Hayashi phase, the entire Sun was in a state of tur-
bulent convection. The turbulence maintained an approximately uniform compo-
sition throughout the interior. At the end of the Hayashi phase the centre of the
Sun became convectively stable, and the base of the convection zone began to re-
treat towards its current location (e.g. Iben, 1965). The retreat of the convection
zone took approximately 3× 107 years, and concluded at about the same time as
the Sun arrived on the main sequence. This timescale is sufficiently short that
the Sun was still “well mixed”, i.e. compositionally homogeneous, at the start of
the main sequence, with a uniform helium mass fraction Y ≈ 0.27.
In the absence of interactions between the various chemical species inside the Sun
(including free electrons) each species would attain separate hydrostatic equilib-
rium on a dynamical timescale, with the hydrogen nuclei located further from the
gravitational centre of the Sun, on average, than the helium nuclei, for example.
1Including lithium and beryllium, which are discussed in the next chapter.
65
66 5. Helium Settling
But collisions between the various species reduce the rate at which this com-
positional segregation occurs, and under certain conditions can even reverse the
direction of the compositional flux. The slow segregation that occurs in the pres-
ence of collisions is usually called either “gravitational settling” or “diffusion”,
and has a timescale much longer than the lifetime of the Sun (see §5.3 below).
The compositional changes produced by gravitational settling may therefore be
regarded as linear perturbations from the perspective of stellar evolution.
The terms “settling” and “diffusion” are used interchangeably in the literature.
In this dissertation, we will use “settling” to refer to a compositional flux that
is proportional to the concentration, and “diffusion” to refer to a compositional
flux that is proportional to the concentration gradient, although this distinction
is somewhat artificial.
Within the Sun’s core, the helium mass fraction Y increases steadily during
the main sequence as a result of hydrogen fusion (see figure 5.1). Within the
convection zone the composition remains roughly uniform throughout the main-
sequence phase, because of turbulent mixing, but as heavy elements settle out of
the convection zone a vertical compositional gradient develops in the outer part
of the radiative envelope. The resulting gradient in the mean molecular weight µ
contributes significantly, in today’s Sun, to the total stable stratification within
that region. Since the formation of this layer is principally due to gravitational
settling of helium, it is usually known as the “helium settling layer”. The reality
of the helium settling layer is strongly indicated by both standard solar-evolution
models and helioseismology, notwithstanding the compositional uncertainties de-
scribed in §2.2 (e.g. Christensen-Dalsgaard et al., 1993; Ciacio et al., 1997; Elliott
& Gough, 1999; Christensen-Dalsgaard & Thompson, 2007, again see figure 5.1).
Because the rate of diffusion of heavy elements through hydrogen is typically very
small compared to thermal diffusion, compositional stratification is much more
effective at inhibiting vertical motions than thermal stratification (recall §1.2).
This near-impermeability of compositionally stratified regions strongly inhibits
MMCs, and has been called the “µ-choke” (Mestel & Moss, 1986).
The effects of compositional stratification were neglected in the tachocline model
5.1 The µ-choke 67
Figure 5.1: The helium mass fraction Y in a one-dimensional model of the solar interior,
plotted as a function of radius R in units of R (J. Christensen-Dalsgaard, personal communi-
cation). Nuclear fusion produces helium in the region R . 0.2R. The helium settling layer is
visible below the convection zone in the region 0.6R . R . 0.7R.
of GM98. Subsequently, McIntyre (2007) presented a significantly revised version
of the GM98 model that included the effects of compositional stratification. In
the revised model, the helium settling layer acts as a barrier against the MMCs
responsible for confining the interior field.2 Furthermore, the stable compositional
stratification within the helium settling layer holds the µ isopleths horizontal to
within “a very tiny fraction of a megametre”. Therefore the tachocline’s MMCs
are unable to dig the polar pits conjectured in §4.2.2. Moreover, the MMCs do not
decay exponentially with depth beneath the magnetic confinement layer in the
manner described in §4.2.1, but terminate abruptly at a “tachopause” marking
the top of the helium settling layer.
These conclusions, which have significant implications for lithium burning (chap-
ter 6), rely on the existence of the helium settling layer and the compositional
stratification therein. Our aim in this chapter is to delineate the physical pro-
cesses involved in the formation of the helium settling layer, and to estimate its
contribution to the total buoyancy frequency. In this way, we can test the robust-
2Spiegel & Zahn (1992) also considered the effect of compositional stratification, as a possible
means to halt radiative spreading in their model of the tachocline. They concluded that radia-
tive spreading is not hindered by the presence of compositional stratification. This conclusion
is in error, as explained in §8.5 of McIntyre (2007).
68 5. Helium Settling
ness of the arguments made in McIntyre (2007) against the assumptions made
by solar evolution models.
5.2 The settling equation
We consider a simple one-dimensional model of solar evolution. For simplicity, we
suppose that the Sun contains only hydrogen and helium, i.e. we neglect heavier
elements. The settling and diffusion of helium within the radiative envelope,
viewed in an inertial frame, can be described by the equation
∂
∂t(ρR
2Y ) = − ∂∂R
(
ρR2(uR − v)Y − ρR2χ
∂Y
∂R
)
(5.1)
(e.g. Proffitt & Michaud, 1991). Here, ρ(R, t) is the density, uR(R, t) is the
radial velocity, χ(R, t) is the diffusivity of helium through hydrogen and v(R, t)
is the “settling velocity”. Both χ and v are calculated from inter-species collision
integrals, and are related by the formula
v = χ
( kp
Hp
+ kTHT
)
(5.2)
where Hp and HT are the pressure and temperature scale heights, and kp and kT
are the corresponding Soret factors, which depend only on the local composition.
The scale heights and Soret factors can be calculated to a high degree of accuracy,
but the formulae for χ and v depend on the Coulomb logarithm, and consequently
are less tightly constrained (e.g. Michaud & Proffitt, 1993). Therefore we will
make only order-of-magnitude statements in what follows. Approximate values
for all these quantities at 0.7R are given in appendix A. We note that v is always
positive, and so the settling flux is always directed downward, as we would expect,
whereas the diffusive flux may be directed either up or down.
Within the nuclear core, a source term must be introduced into equation (5.1) to
represent the production of helium by nuclear fusion. The changes in the core’s
composition drive the main-sequence evolution of the Sun. However, the com-
positional changes within the outer radiation zone are small, as can be seen in
figure 5.1, and play no role in the Sun’s main-sequence evolution to a first ap-
5.2 The settling equation 69
proximation. In describing the formation of the helium settling layer, the density
ρ, the pressure p and the temperature T may therefore be assumed independent
of Y in (5.1), and regarded as known functions of R and t. This assumption also
allows us to approximate the Soret factors in (5.2) as constant, rather than as
functions of Y . Similarly we may neglect the Y -dependence of χ and v, so that
the settling and diffusion coefficients are also known functions of R and t. In
particular, since the Coulomb logarithm is approximately constant throughout
the solar interior, we have
χ ∝ T
5/2
ρ (5.3)
(Michaud & Proffitt, 1993). Under these assumptions (5.1) becomes a linear PDE
for the helium mass fraction Y within the radiative envelope.
Equation (5.1) assumes purely radial motions, and therefore neglects any compo-
sitional mixing arising from turbulence or MMCs. In stellar evolution modelling,
these mixing processes are usually parametrised in terms of a turbulent diffusivity
(see discussion in §5.5). In this study, for simplicity, we will assume that tur-
bulent mixing is confined to the convection zone, where it maintains a perfectly
homogeneous composition throughout the Sun’s main-sequence life. This treat-
ment of convective mixing is equivalent to assuming an infinitely large turbulent
diffusivity within the convection zone.
The location of the convection zone does not remain fixed in space during the
main sequence, because the Sun’s outer layers expand over time. However, the
convection zone’s lower boundary is approximately fixed in the Lagrangian de-
scription. It is therefore convenient to rewrite equation (5.1) in terms of new
coordinates M and t, where M(R, t) represents the mass contained within a
sphere of radius R at time t. Then
(∂M
∂R
)
t
= 4piρR2 (5.4)
and
(∂M
∂t
)
R
= −4piρR2uR . (5.5)
70 5. Helium Settling
It follows that the “mass coordinate” M is a Lagrangian coordinate, since
(∂R
∂t
)
M
= −
(∂M
∂t
)
R
/(∂M
∂R
)
t
(5.6)
= uR . (5.7)
Written in terms of the new coordinates M and t, equation (5.1) becomes
∂Y
∂t =
∂
∂M
(
4piρR2vY + (4piρR2)2χ ∂Y∂M
)
. (5.8)
So, in the Lagrangian description, the mass fraction Y is a conserved quantity.
Integrating equation (5.8) between fixed mass coordinates, say M = M1 and
M = M2, we obtain
d
dt
∫ M2
M1
Y dM =
[
4piρR2vY + (4piρR2)2χ ∂Y∂M
]M2
M1
(5.9)
which quantifies how the total mass of helium between concentric material sur-
faces changes due to settling and diffusive fluxes over those surfaces.
Though the evolution equation (5.8) does not hold within the convection zone,
the principle of helium conservation can still be applied there. Neglecting any
loss of helium due to outflows from the solar surface we may therefore express
the conservation of helium in the convection zone as
d
dt
∫ M
MRZ
Y dM = −
(
4piρR2vY + (4piρR2)2χ ∂Y∂M
)∣
∣
∣
∣
M=MRZ
(5.10)
where MRZ ≈ 0.975M is the mass coordinate at the top of the radiation zone,
and M is the total solar mass, both of which we approximate as constant. The
right-hand side of equation (5.10) is to be evaluated immediately below the base
of the convection zone, where χ takes a microscopic (non-turbulent) value. Our
assumption that turbulent mixing homogenises the chemical composition within
the convection zone implies that Y can be taken outside the integral on the left-
hand side of (5.10). We thus obtain an equation for the evolution of Y within
5.3 The formation of the helium settling layer 71
the convection zone:
MCZ
dY
dt = −
(
4piρR2vY + (4piρR2)2χ ∂Y∂M
)
∣
∣
∣
∣
M=MRZ
(5.11)
where MCZ = M−MRZ is the mass of the convection zone. This equation serves
as an outer boundary condition when solving (5.8) within the radiative envelope.
5.3 The formation of the helium settling layer
Since the Sun’s chemical composition is assumed to be uniform at the start of the
main sequence, there is no diffusive flux anywhere in the interior at that time.
Therefore the initial evolution within the bulk of the radiative envelope (i.e. away
from the core and convection zone) arises entirely from spatial variations in the
settling rate:
∂Y
∂t = Y
∂
∂M (4piρR
2v) (5.12)
∼ −Y vHρ
, (5.13)
where Hρ is the density scale height. This process has a timescale of order
Hρ/v ∼ 60Gyr — considerably longer than the Sun’s main sequence lifetime —
and produces a slight compositional gradient within the bulk of the radiation
zone (see figure 5.1, between 0.4R and 0.6R). The formation of the settling
layer is a separate process that is driven by the boundary condition (5.11). For
the rest of this section we will neglect the effects of spatial or temporal variations
in ρR2, v and χ on the formation of the settling layer. The importance of such
effects will be considered in §5.3.3.
Formation of the settling layer, at least in the early stages of the main sequence,
is driven by helium settling out of the convection zone. From (5.11) we see that
the helium mass fraction Y at the top of the radiative envelope initially decreases
exponentially on the timescale tCZ given by
tCZ =
MCZ
v
∂R
∂M
∣
∣
∣
∣
MRZ
(5.14)
72 5. Helium Settling
which is of the order of 30Gyr. This depletion of helium from the convection zone
produces a steep gradient in Y at the top of the radiative envelope. Diffusion of
helium back into the convection zone reduces this gradient, producing a diffusive
tail in Y in the outer part of the radiative envelope whose thickness increases as
(χt)1/2. This diffusive tail is the nascent helium settling layer. During these early
stages, the compositional gradient in the settling layer is roughly ∂ lnY/∂r ∼
−(t/tCZ)/(χt)1/2 ∝ −t1/2. The subsequent evolution of the settling layer depends
crucially on the relative magnitudes of tCZ and the settling–diffusion time, tSD =
χ/v2.
5.3.1 The solar case
In the outer part of the Sun’s radiative envelope, the settling–diffusion time tSD is
approximately 7Gyr, which is significantly shorter than the timescale for helium
depletion in the convection zone, tCZ ≈ 30Gyr. As the helium settling layer
steepens, the upward diffusive flux of helium into the convection zone increases,
while the downward settling flux remains roughly constant. After a time of order
tSD, the upward diffusive flux is still smaller than the downward settling flux,
by a factor (χ/v)∂ lnY/∂R ∼ (χ/v)(tSD/tCZ)/(χtSD)1/2 ∼ (tSD/tCZ). To a first
approximation we may therefore neglect the diffusive flux of helium into the con-
vection zone throughout this period. The upper boundary condition (5.11) can
therefore be approximated as
dY
dt = −
Y
tCZ
at M = MRZ . (5.15)
For times t . tSD, the helium settling layer thickens like (χt)1/2. When t ∼ tSD
the thickness of the settling layer becomes comparable to the settling–diffusion
length, χ/v ≈ 0.02R. At this time the divergence of the settling flux in (5.8)
becomes of the same order as the divergence of the diffusive flux. For times
t & tSD the depth of the settling layer therefore increases roughly linearly in
time, like vt, and the compositional gradient therefore converges to a constant
5.3 The formation of the helium settling layer 73
value,
∂ lnY
∂R
∣
∣
∣
∣
MRZ
≈ −(t/tCZ)/(vt) = −1/(vtCZ) (5.16)
⇔ ∂Y∂M
∣
∣
∣
∣
MRZ
≈ − YMCZ
. (5.17)
In fact, this implies that the settling flux out of the convection zone exceeds
the diffusive flux into the convection zone even for times > tSD. Throughout
the evolution the abundance of helium in the convection zone therefore decreases
approximately exponentially, but on a long timescale tCZ ≈ 30Gyr, and so helium
depletion in the convection zone occurs at a roughly constant rate during the main
sequence.
5.3.2 The super-solar case
Stars that are slightly more massive (and therefore hotter and less opaque) than
the Sun have smaller convection zones (e.g. Pinsonneault et al., 2001). For such
stars the available “reservoir” of helium in the convection zone is smaller, and the
timescale tCZ is therefore shorter. A star of mass 1.1M, for example, would have
tCZ  tSD. In that case the upward diffusive flux of helium into the convection
zone quickly comes into balance with the downward settling flux, after a time of
roughly t2CZ/tSD  tCZ. From this time onwards the upper boundary condition
(5.11) can be approximated as
0 ≈
(
vY + 4piρR2χ ∂Y∂M
)
∣
∣
∣
∣
M=MRZ
(5.18)
and the formation of the helium settling layer proceeds as though the convection
zone were infinitely thin. From (5.18) we deduce that
∂ lnY
∂R
∣
∣
∣
∣
MRZ
≈ −v/χ . (5.19)
74 5. Helium Settling
5.3.3 The effect of spatio-temporal variations
The formation of the settling layer is affected quantitatively by the variations of
ρR2, v and χ in both space and time, which have been neglected in the analysis
above. The effect of these variations can be estimated by a simple experiment in
which we evolve equation (5.8) with static profiles of ρR2, v and χ. To this end
we have written a numerical code to timestep (5.8) using profiles for ρR2, v and
χ calculated from a one-dimensional model of the solar interior (J. Christensen-
Dalsgaard, personal communication). The numerical scheme employed is identi-
cal to that described in Proffitt & Michaud (1991), except that the settling flux
is calculated using a first-order upwind difference (as opposed to their centred
difference). This gives greater numerical stability, with no loss in overall accu-
racy. In figure 5.2 we plot the resulting evolution of Y , from an initially uniform
profile of Y = 0.27, over 5Gyr on the main sequence. The formation of the
0.2
0.3
0.4
0.5
0.6
0.7
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
1
5
Figure 5.2: Plots of Y against R/R after 1, 2, 3, 4 and 5Gyr on the main sequence.
helium settling layer closely follows the analytical predictions, suggesting that
spatial variations in the profiles of ρR2, v and χ do not significantly influence the
process. We note in particular that the abundance of helium in the convection
zone decreases approximately linearly with time, indicating that the downward
settling flux out of the convection zone dominates the upward diffusive flux into
5.4 The compositional gradient 75
the convection zone throughout the period simulated. The plots in figure 5.2
are also in rough agreement with those of Ciacio et al. (1997), where the full
nonlinear, time-dependent equations were employed (see figure 5.3). The most
conspicuous difference — aside from the lack of nuclear fusion in our model — is
due to the gradual shift in the location of the convection zone mentioned in §5.2.
Nonetheless, both models produce comparable estimates for the compositional
gradient within the helium settling layer, so the effect of temporal variations in
the profiles of ρR2, v and χ also has little impact on the formation of this layer.
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.2
0.3
0.4
0.5
0.6
0.1 Gyr
1
2
3
4.57
Figure 5.3: Plots of Y against R/R at various stages during main sequence. Taken from
Ciacio et al. (1997).
5.4 The compositional gradient
We can use the scaling arguments in §5.3.1 to estimate the contribution from
compositional stratification, Nµ say, to the buoyancy frequency N at the top of
the radiative envelope. In this region both hydrogen and helium are fully ionised,
and so the mean molecular weight µ, neglecting contributions from heavier ele-
ments, depends only on Y :
µ = 12− 54Y
(5.20)
76 5. Helium Settling
(e.g. Stro¨mgren, 1938). Hence
d lnµ = 54µdY (5.21)
and so, assuming an ideal gas equation of state,
N2µ = −g
∂ lnµ
∂R = −
5
4gµ
∂Y
∂R . (5.22)
Using the upper bound for ∂ lnY/∂R given by (5.16) we obtain an estimate for Nµ
of roughly 10−3s−1. This is in good agreement with the value Nµ ≈ 0.5×10−3s−1
computed by McIntyre (2007, §8.5) from information in Christensen-Dalsgaard
& Thompson (2007).
5.5 Summary
The discussions above indicate that the helium settling layer is a robust feature of
the solar interior, i.e. that its existence is not conditional on the assumptions or
parametrisations made in solar-evolution models. In particular, the presence of
additional turbulent mixing below the base of the convection zone cannot prevent
the formation of the settling layer. If the turbulent mixing is very efficient, then
it can be regarded as an extension of the mixing in the convection zone, and will
therefore only cause the settling layer to form deeper inside the radiative envelope.
Less efficient turbulent mixing, on the other hand, can make the settling layer
thicker, and thereby reduce the compositional gradient below the convection zone.
For example, if we suppose that turbulence below the convection zone produces an
fifty-fold increase in the effective helium diffusivity χ, then the effective settling–
diffusion time tSD also increases by a factor of fifty. This puts the Sun into the
regime discussed in §5.3.2 (i.e. tSD  tCZ). The compositional gradient at the
top of the radiative envelope is then reduced by a factor tCZ/tSD ≈ 0.1 from the
value estimated in §5.4, meaning that N2µ is about ten times smaller. A larger
turbulent diffusivity would produce a greater reduction in Nµ, but could not be
reconciled with helioseismology (Christensen-Dalsgaard et al., 1993).
As described by McIntyre (2007), the presence of the helium settling layer has
5.5 Summary 77
significant implications for the convection zone’s lithium abundance. These are
discussed in detail in the next chapter.
We note that the helium settling layer, though an effective barrier against MMCs,
does not nullify the arguments in favour of an interior magnetic field put forward
in chapter 2. In the absence of such a magnetic field, the burrowing MMCs
described by Spiegel & Zahn (1992) would ventilate the outer part of the radiative
envelope, thickening the tachocline and causing the helium settling layer to form
much deeper than is observed. Furthermore, the region below the helium settling
layer, inaccessible to MMCs driven from the convection zone, would probably
not be spun down as efficiently, leading to a rapidly rotating core. This point is
discussed further in the next chapter.
78
Chapter 6
Depletion of Lithium and
Beryllium
Observations of solar-type stars indicate that the surface lithium abundance de-
creases with age on the main sequence (e.g. Soderblom et al., 1990), although
not all stars fit the general trend (e.g. Gaidos, 1998). The Sun seems to be a
particularly extreme example of lithium depletion; the lithium to hydrogen ratio
in the solar convection zone is now less than 1% of its primordial value (Vauclair
et al., 1978; Asplund et al., 2009). The surface beryllium abundance is less well
determined by observations (e.g. Pinsonneault & Delahaye, 2009) but is believed
to be at least half of the primordial value.
Many mechanisms for lithium depletion in stars have been proposed, including
gravitational settling (Michaud, 1986), convection-zone dilution caused by mass
loss (Weymann & Sears, 1965), and various sources of internal mixing (e.g. Garc´ıa
Lo´pez & Spruit, 1991; Charbonnel et al., 1992; Fritts et al., 1998; Israelian et al.,
2009, & refs). Nearly all1 of the proposed mechanisms rely on the fact that 7Li
is destroyed (or “burned”) by proton capture at around 2.5 × 106 K (Michaud
& Charbonneau, 1991). This is the temperature of the solar interior at a depth
of about 0.05R below the bottom of the convection zone (e.g. Christensen-
1The gravitational settling mechanism proposed by Michaud (1986) is an exception. This
mechanism is not thought to be particularly relevant to the Sun, since it would tend to deplete
lithium and beryllium equally, but it may contribute to the “lithium dip” observed in F stars
(Boesgaard & Tripicco, 1986).
79
80 6. Depletion of Lithium and Beryllium
Dalsgaard et al., 1992). If material from this depth is efficiently mixed into the
convection zone then the surface depletion of lithium can be explained. If the
mixed region extends significantly beyond 0.1R below the bottom of the convec-
tion zone, where the temperature exceeds 3 × 106 K, then the convection zone’s
beryllium will also be burned.2 The surface abundance observations therefore
impose tight constraints on the degree of mixing within the solar interior.
Arguably the simplest process that can produce the required mixing is main-
sequence mass loss (Weymann & Sears, 1965). As matter is lost from the solar
surface, the Sun’s internal structure adjusts, and layers below the convection zone
that were previously hot enough to burn lithium become cooler. If the rate of mass
loss is sufficiently fast, then the cool, lithium-poor layers may become convectively
unstable, leading to “lithium dilution” in the convection zone. However, the mass-
loss rate would need to be quite finely tuned in order to dilute the convection
zone’s lithium by a factor of 100 without greatly diluting the beryllium as well
(Swenson & Faulkner, 1992).
The mean meridional circulations (MMCs) described in chapter 2 can transport
lithium, as well as angular momentum, within the radiative envelope. This is
another potential source of mixing, which may be further enhanced by shear in-
stabilities and turbulence triggered by the angular momentum transport. This
combination of processes is often called “rotational mixing” (e.g. Charbonnel
et al., 1992). The early models of rotational mixing sought to explain both the
spin-down of the solar interior and the depletion of lithium in the convection
zone simultaneously. However, the limits on interior mixing implied by the sur-
face beryllium abundance make it difficult to explain solar spin-down by this
mechanism. Furthermore, the production of helium in the Sun’s core establishes
strong, stable compositional stratification deep within the radiative interior (fig-
ure 5.1) which practically forbids the propagation of MMCs into that region
(Mestel, 1953). Therefore models based entirely on rotational mixing tend to
have a rapidly rotating core, which contradicts the current measurements from
helioseismology.
Another potential source of turbulence and mixing in the solar interior arises
2Although rapid beryllium burning occurs only at temperatures exceeding 3.4× 106 K (e.g.
Christensen-Dalsgaard et al., 1991).
81
from internal gravity waves. As described in §2.4, gravity waves generated at
the bottom of the convection zone can transport angular momentum within the
interior. This transport can induce MMCs and turbulence, and hence lead to
mixing (e.g. Fritts et al., 1998).
Using a combination of rotational mixing and gravity wave mixing, Charbonnel
& Talon (2005) constructed a model that can simultaneously explain spin-down
and lithium depletion within the solar interior. In their model, differential rota-
tion within the interior is assumed to drive turbulence that maintains constant
angular velocity on each spherical surface (see §2.4). This turbulence also dissi-
pates internal gravity waves generated by the convection zone’s turbulent stresses,
leading to spin-down. At the same time, rotational mixing below the convection
zone leads to lithium depletion. Spin-down by gravity waves is an essential part
of this picture — without a mechanism for reducing the vertical shear in the
interior, their model would produce too much mixing, and therefore too much
lithium depletion.
However, the angular momentum transport produced by internal gravity waves is
sensitive to the manner in which the waves are generated and dissipated (§2.4).
Furthermore, the presence of a fossil magnetic field in the radiative envelope
would strongly inhibit the MMCs required for rotational mixing, as described
in chapter 3. In the GM98 model, it was argued that MMCs would be able to
burrow into the radiative envelope only at the polar “weak spots” of the dipo-
lar magnetic field (figure 4.2). To produce lithium depletion in the convection
zone, these “polar pits” would need to be well ventilated and deep enough to
burn lithium. However, in this model the effects of compositional stratification
were neglected. In the revised model presented by McIntyre (2007) the existence
of polar pits was ruled out by the µ-choke from the the helium settling layer.
Instead, it was conjectured that the outer 0.1R of the radiative envelope is ven-
tilated by a combination of MMCs and MHD turbulence (see appendix G), and
that the helium settling layer forms below this ventilated region. At the time,
the uncertainties in the helioseismic evidence described in §2.2 permitted such
speculation (Christensen-Dalsgaard & di Mauro, 2007). But more recent work
seems to confirm that the top of the helium settling layer is in fact located sig-
nificantly higher than 0.6R (Christensen-Dalsgaard & Gough, in preparation).
82 6. Depletion of Lithium and Beryllium
This suggests that the Sun’s lithium must have been depleted in the early stages
of the main sequence, before the formation of the helium settling layer. (This
hypothesis has some support from observations, e.g. Chen et al. (2001).)
In Part III we will develop a self-consistent model of magnetic field confinement
in the high-latitude tachocline that allows us to test some of the speculation in
GM98 and McIntyre (2007). In Part IV we will discuss the implications of our
model for solar lithium depletion.
As a postscript to this review of lithium depletion, we note that the theories
outlined above all need to be rather finely tuned in order to produce the observed
99% depletion of lithium in the Sun’s convection zone. Such fine tuning would
not be necessary if it could be shown that lithium is actively produced in the
outer layers of the Sun, for example by spallation in solar flares (e.g. Christian
& Mathioudakis, 2000). In that case the observations could be explained even
if all the convection zone’s lithium were destroyed early in the main sequence
(Guzik et al., 1987). However, most “autogenesis” theories predict a higher ratio
of 6Li/7Li than is observed (Feast, 1966; Soderblom et al., 1990).
Part III
The Magnetic Confinement Layer
83
Chapter 7
The High-Latitude Tachocline
7.1 The structure of the magnetic field
Based on the arguments put forward in Parts I and II, we now attempt to con-
struct a global model of angular momentum balance and magnetic field con-
finement in the solar tachocline. Our focus will be on today’s Sun, as inferred
from helioseismology. However, in Part IV we will consider the implications of
our model for the early main-sequence Sun, with particular regard to lithium
depletion.
Figure 7.1 is a sketch suggesting the simplest possible structure for the time-
averaged magnetic field in the tachocline. The inner, grey sphere represents the
top of the helium settling layer. In high latitudes, the field lines emerging from
the helium settling layer are swept aside by the tachocline’s downwelling MMC.
The latitudinal shear in the tachocline winds up these field lines, transmitting a
prograde “Alfve´nic torque” to high latitudes, as required to confine the down-
welling MMC. We presume that the time-averaged effect of convective overshoot
and other turbulent processes in the tachocline is to hold the field lines horizon-
tal in middle and low latitudes, as described in §4.1. We presume also that the
turbulent processes act in the manner of a magnetic diffusivity, preventing the
field lines from being wound up arbitrarily tightly by the latitudinal shear.
Assuming that the high-latitude confinement layers at either pole behave in a
85
86 7. The High-Latitude Tachocline
Figure 7.1: Sketch of the time-averaged magnetic field in the tachocline. Poloidal magnetic
field lines (red) emerge from the helium settling layer (inner sphere) in high latitudes, and are
wound up by differential rotation in the tachocline, acting against eddy diffusion. A prograde
Alfve´nic torque is transmitted from low to high latitudes along these field lines. The cutaway
outer sphere, which represents the bottom of the convection zone, is shaded according to angular
velocity. The dashed lines indicate the latitudes at which the rotation of the convection zone
matches that of the interior, Ωi = 2.7× 10−6s−1.
similar fashion, the field lines must emerge from the helium settling layer with
no toroidal component Bφ. Otherwise an Alfve´nic torque would be exerted on
the top of the helium settling layer at either pole, producing torsional oscilla-
tions within the apple-core region between the two confinement layers. A similar
topological constraint arises in the simulations of Braithwaite & Spruit (2004),
although in their case it results from the assumption of perfectly insulating bound-
aries (see figure 3.1 in §3.1). There must still be a nonvanishing Bφ component
deeper within the interior, in order to stabilise the neutral ring of the dipolar
field (again see figure 3.1).
7.2 The structure of the tachocline
It is apparent from the discussions in chapter 4 that the high-latitude regions are
the most critical to the magnetic confinement problem. Our aim in the following
7.3 Angular momentum balance 87
chapters is to construct a fully self-consistent model of the high-latitude magnetic
confinement layer, based on the global picture presented in figure 7.1. We seek
to describe the complete vertical structure of the confinement layer, including
how it matches on to the mean-field-free tachocline above, and on to the Ferraro-
constrained radiative envelope below. However, we do not attempt to describe
in detail the turbulent processes within the convection zone and tachocline that
pump the downwelling MMC responsible for confining the field. To do so would
require an accurate description of the (magneto)hydrodynamic turbulence within
the convection zone and tachocline, which is beyond the scope of this study. The
existence of downwelling in the high-latitude tachocline is strongly supported
by the arguments in §2.3.2 and §4.2.1 (see also appendix G and discussions in
Part IV).
Since our focus is on today’s Sun, our model must describe how the structure of
the confinement layer is affected by the presence of the compositionally stratified
helium settling layer in the outer part of the radiative envelope. As described in
chapter 5, the µ-choke from the settling layer will severely limit the magnitude
and scale of any MMCs within that region. The downwelling MMC that confines
the magnetic field in high-latitudes must therefore feed into an equatorward flow
immediately above the helium settling layer. In the absence of helium diffusivity,
we would expect a discontinuity in the vertical compositional gradient at the
interface between the magnetic confinement layer and helium settling layer, as in
the model of McIntyre (2007). Here, we shall allow for finite helium diffusivity,
and so the change in the compositional gradient will be resolved over a thin
“helium sublayer”.
Figure 7.2 shows a vertical cross-section through the confinement layer in a region
surrounding the north pole.
7.3 Angular momentum balance
The most crucial role of the interior magnetic field is to supply angular momentum
to the high-latitude confinement layer, and thereby to prevent the tachocline’s
MMCs from burrowing into the radiative envelope. Consider the overall angular
88 7. The High-Latitude Tachocline
−5 −4 −3 −2 −1 0 1 2 3 4 5−1
0
1
2
3
4
5 u
B
Ωi
tachocline
confinement
layer
helium
sublayer
helium
settling
layer
Figure 7.2: The magnetic confinement layer near the north pole. The field strength falls off
exponentially with altitude. The streamlines with arrows show the downwelling responsible for
the confinement. Compositional stratification is indicated by shading. With typical parameters
the confinement layer is a fraction of a megametre thick. The plot is from a numerical solution
of the confinement-layer equations (§7.4). The horizontal and vertical axes are colatitude and
altitude in units of the magnetic advection–diffusion length δη = η/U .
momentum balance in the cylindrical volume illustrated in figure 7.2, which we
denote as V . Adopting cylindrical polar coordinates (r, φ, z), angular momentum
balance in the non-rotating frame can be expressed as
∫
∂V
rρuφu · dS =
1
4pi
∫
∂V
rBφB · dS+
∫
∂V
r2ρν∇
(uφ
r
)
· dS , (7.1)
i.e. the net flux of angular momentum into or out of the volume V must be bal-
anced by Maxwell and viscous stresses exerted over the boundary, ∂V . Provided
that departures from uniform rotation remain small, the system will be mag-
netostrophic, i.e. Coriolis and Lorentz forces will dominate viscous and relative
inertial forces. We may then approximate uφ ≈ Ωir, so that (7.1) becomes
∫
∂V
r2ρΩiu · dS ≈
1
4pi
∫
∂V
rBφB · dS . (7.2)
For any flow qualitatively like that in figure 7.2 the left-hand side of (7.2) is
positive — i.e. the flow produces a net export of angular momentum. To main-
tain balance, angular momentum must be supplied by the magnetic field, via an
Alfve´nic torque over the boundary ∂V . In principle this Alfve´nic torque can be
supplied either over the bottom of the domain, or over its periphery, or a com-
7.4 The equations and parameters 89
bination of the two. However, as discussed in §7.1, any Alfve´nic torque exerted
on the top of the helium settling layer in high-latitudes will lead to torsional
oscillations in the apple-core region between the two polar confinement layers.
Therefore, in the steady state, the Alfve´nic torque over the bottom of the do-
main V must vanish. This is equivalent to the condition that Bφ = 0 over the
bottom of the domain. This boundary condition ensures that the Alfve´nic torque
on the confinement layer comes exclusively from the confinement layer’s sideways
connection to low latitudes,1 via field lines like those sketched in figure 7.1.
7.4 The equations and parameters
Within the high-latitude cylindrical domain shown in figure 7.2 we may take
gravity to be uniform and anti-parallel to the rotation axis. We impose uniform
downwelling of magnitude U through the top of the domain, and an axial dipolar
magnetic field underneath to represent the interior magnetic field Bi. We an-
ticipate that the characteristic vertical scale of the magnetic confinement layer
will be the magnetic advection–diffusion length δη = η/U , where the magnetic
diffusivity2 η ≈ 4.1 × 102cm2s−1. If the imposed downwelling U has magnitude
∼ 10−5cm s−1 then this lengthscale is ∼ 0.4Mm, which is far smaller than the lo-
cal pressure scale height, Hp ≈ 60Mm. We may therefore employ the Boussinesq
approximation in modelling the confinement layer.
We work in a frame rotating with the same angular velocity as the interior,
Ωi = 2.7× 10−6s−1. We nondimensionalise the Boussinesq equations using δη as
the lengthscale, U as the scale for the velocity field u, and δη/U as the timescale,
∼ 105yr if U ∼ 10−5cm s−1. We measure the magnetic field B in units of Alfve´n
speed, 0.6cm s−1 per gauss at a tachocline density of 0.2g cm−3, and nondimen-
sionalise B with respect to (2Ωiη)1/2 ≈ 0.05cm s−1. The significance of this scaling
will be explained in chapter 8. We suppose that the thermal and compositional
stratifications are approximately uniform within the helium settling layer, shaded
1If the interior field were assumed to be a quadrupole, rather than a dipole, then part of the
Alfve´nic torque could be transmitted from low to high latitudes via the deep radiative interior.
However, it is not known whether a quadrupolar interior field can be dynamically stable, so
here we consider only the simpler, dipolar case illustrated in figure 7.1.
2Values for all the relevant dimensional parameters in the tachocline are listed in appendix A.
90 7. The High-Latitude Tachocline
in figure 7.2. Writing Tˆ and µˆ for the fractional perturbations of temperature
and mean molecular weight, we define the buoyancy frequencies NT and Nµ to
be exactly constant at the bottom of the domain,
∂Tˆ
∂z
∣
∣
∣
∣
∣
bottom
= N
2
T δη
g = const., (7.3)
∂µˆ
∂z
∣
∣
∣
∣
bottom
= −
N2µδη
g = const., (7.4)
where z is now the dimensionless axial or altitude coordinate. At the top of
the helium settling layer we have NT ≈ 0.8 × 10−3s−1, and Nµ ≈ 0.5 × 10−3s−1
(appendix A). In place of Tˆ and µˆ we define dimensionless rescaled quantities T
and µ by
N2T δη
g T = Tˆ , (7.5)
N2µδη
g µ = µˆ . (7.6)
The thermal and compositional diffusivities are denoted respectively by κ ≈
1.4 × 107 cm2s−1 and χ ≈ 0.9 × 101 cm2s−1. We then arrive at the following
dimensionless Boussinesq equations,
RoDuDt + ez × u = −∇p + αT Tez − αµµez
+ (∇×B)×B+ Ek∇2u (7.7)
0 = ∇ · u (7.8)
∂B
∂t = ∇× (u×B) +∇
2B (7.9)
0 = ∇ ·B (7.10)
DT
Dt =
κ
η∇
2T (7.11)
Dµ
Dt =
χ
η∇
2µ (7.12)
where D/Dt = ∂/∂t+u ·∇, the Lagrangian or material derivative, and where ez
is the unit vector parallel to the rotation axis. We have defined four dimensionless
7.4 The equations and parameters 91
coefficients Ro, Ek, and
αT =
N2T δ2η
2Ωiη
, αµ =
N2µδ2η
2Ωiη
; (7.13)
Ro and Ek are Rossby and Ekman numbers defined by
Ro = U/δη2Ωi
, Ek = ν2Ωiδ2η
= νηRo (7.14)
where ν is the kinematic viscosity, the sum of molecular and radiative contribu-
tions, ≈ 2.7×101cm2s−1 (appendix A). In (7.12) we have neglected gravitational
settling, an excellent approximation in virtue of the short timescale of the con-
finement-layer dynamics relative to the Sun’s lifetime.
For U ∼ 10−5cm s−1 the Rossby number is tiny, Ro ∼ 0.5×10−7, and the Ekman
number is smaller still, because ν/η ≈ 0.7 × 10−1. To excellent approximation,
therefore, the flows under consideration will be magnetostrophic. In (7.7) the
Coriolis force will be balanced against the combined pressure-gradient, buoyancy
and Lorentz forces:
ez × u = −∇p + αT Tez − αµµez + (∇×B)×B . (7.15)
For simplicity we restrict attention to axisymmetric steady states. Then the
azimuthal components of (7.15) and its curl are respectively
ur =
1
rB ·∇(rBφ) , (7.16)
∂uφ
∂z = αT
∂T
∂r − αµ
∂µ
∂r +
1
r
∂
∂z (B
2
φ)− rB ·∇
( [∇×B]φ
r
)
(7.17)
where r is the dimensionless perpendicular distance from the rotation axis, and
where suffixes (r, φ, z) denote vector components in cylindrical polar coordinates.
Equation (7.16) represents the local torque balance about the rotation axis, after
multiplication by r (cf. equation (7.2)). It describes how the retrograde Coriolis
torque from the equatorward flow is balanced by the prograde Lorentz torque
from the confined magnetic field. Equation (7.17) represents low Rossby num-
ber thermal-wind balance generalised to include compositional gradients and the
92 7. The High-Latitude Tachocline
Lorentz force-curl. In the upper part of figure 7.2, where the magnetic field and
compositional stratification are both negligible, this balance becomes the stan-
dard thermal-wind balance
∂uφ
∂z = αT
∂T
∂r (7.18)
(cf. equation (1.2)).
7.5 The limit of infinite stratification
Equation (7.17) describes how the Coriolis and Lorentz force-curls act to tilt the
stratification surfaces. Now if the stable stratification is sufficiently strong, i.e. if
αT and αµ are both sufficiently large, then the tilting will be only slight. With
U ∼ 10−5cm s−1 we have
αT ∼ 5× 1011 (7.19)
and αµ ∼ 2× 1011. (7.20)
So we anticipate that both the thermal and compositional stratification surfaces
will be nearly “flat”, meaning gravitationally horizontal, within the confinement
layer. This flatness will be more precisely quantified in chapters 8 and 9 using
scale analysis. In fact, the analysis in chapter 9 shows that, to a first approx-
imation, the helium sublayer can be modelled as an impermeable, frictionless,
horizontal lower boundary.
If we take the asymptotic limit in which αT , αµ →∞, with u and B both finite,
then we may seek solutions in which the stratification surfaces are perfectly flat,
i.e. T = T (z) and µ = µ(z). It then proves possible to solve the steady-state con-
finement-layer equations analytically, as a system of coupled ordinary differential
equations, as we now show. We will refer to these analytical solutions as “perfectly
flat”.
In the limit of perfect flatness, (7.11) and (7.12) imply for steady flow that
uz =
κ
η
d
dz ln
dT
dz =
χ
η
d
dz ln
dµ
dz . (7.21)
7.5 The limit of infinite stratification 93
Therefore uz is a function of z alone, uz = uz(z). From (7.8) it then follows
that ur is r times a function of z alone, on the assumption of regularity at the
pole r = 0. We say that the poloidal velocity field is “horizontally self-similar”,
since it is invariant under a uniform stretching of the horizontal coordinates.3
The induction equation (7.9) then permits a steady poloidal magnetic field that
is horizontally self-similar in the same sense. Equations (7.7)–(7.10) then reduce
to
rur = B ·∇(rBφ) (7.22)
0 = 1r
∂(rur)
∂r +
duz
dz (7.23)
ru ·∇Bφr = rB ·∇
uφ
r +
(
∇2 − 1r2
)
Bφ (7.24)
uz
dBz
dz = Bz
duz
dz +
d2Bz
dz2 (7.25)
0 = 1r
∂(rBr)
∂r +
dBz
dz . (7.26)
This is a complete set because perfect flatness implies that the meridional compo-
nents of (7.7) are not needed; more precisely, equation (7.17) becomes degenerate
in the limit where αT , αµ →∞ and ∂T/∂z, ∂µ/∂z → 0. The azimuthal compo-
nent of (7.7) is replaced by (7.16) and multiplied by r to give (7.22), expressing
the balance between Coriolis and Lorentz torques. Equations (7.23)–(7.26) cor-
respond to (7.8)–(7.10); equations (7.11) and (7.12) have no further role, beyond
their connection to the downwelling expressed by (7.21).
Equations (7.22)–(7.26) represent a slight generalisation of the confinement-layer
equations presented by Wood & McIntyre (2007, hereafter WM07). In particular,
equation (7.24) permits more general toroidal magnetic and differential-rotation
fields than were considered by WM07. Perhaps surprisingly, these equations
admit a family of solutions in which the function uz(z) is arbitrary except for
certain restrictions on its asymptotic behaviour as z → ±∞. In particular,
we require uz(z) → −1 as z → +∞ and uz(z) → 0 as z → −∞. The necessary
restrictions will shortly be made more precise. Solutions can be most conveniently
3This horizontal self-similarity is an indication that the confinement layer’s intrinsic hori-
zontal lengthscale becomes infinite in the asymptotic limit αT , αµ →∞. This is confirmed by
the scale analysis in chapters 8 and 9.
94 7. The High-Latitude Tachocline
constructed by specifying a suitable uz(z) at the start even though, in virtue of
(7.21), we could instead choose to specify either T (z) or µ(z).
With uz(z) specified, we can find Bz(z) immediately by solving the vertical com-
ponent (7.25) of the induction equation as a linear ordinary differential equation,
assuming that Bz vanishes far above the confinement layer (z → +∞) and that
it matches on to the interior dipolar magnetic field structure beneath (z → −∞).
The interior dipole has the same horizontally self-similar structure as the con-
finement layer, with components satisfying Biφ = 0, Bir/r = constant, and Biz a
linear function of z consistent with (7.26). Even though the balance in (7.25) is
not simple advective–diffusive, we find that Bz decays upward like exp(−z).
The radial components of u and B can be found directly from their vertical
components, by using (7.23) and (7.26) and assuming regularity at the pole r = 0:
ur = −
r
2
duz
dz , (7.27)
Br = −
r
2
dBz
dz . (7.28)
So once we have Bz(z) we can calculate Br from (7.28), and then the toroidal field
Bφ from (7.22) by using (7.27) and taking advantage of the hyperbolic character
of the operator B ·∇. By calculating Bφ in this way, we ensure that the Lorentz
torque balances the Coriolis torque along each magnetic field line. Requiring
that Bφ(r, z) → 0 for all r as z → −∞ leads to the following, unique solution of
(7.22):
Bφ = Bz
∫ z
−∞
ur
B2z
dz . (7.29)
For any ur profile that decays exponentially as z → −∞, this solution for Bφ,
and with it the Maxwell stress and Alfve´nic torque, will also decay exponentially
as z → −∞. The expression (7.29) then shows that Bφ has the same self-similar
functional form as ur and Br, namely r times a function of z alone.
To ensure that Bφ decays aloft, as z → +∞, it is sufficient to assume that
ur = O(exp(−γz)) as z → +∞ (7.30)
7.5 The limit of infinite stratification 95
where γ > 1, implying that uz(z) ∼ −1 + O(exp(−γz)). The three cases γ > 2,
γ = 2, and 2 > γ > 1 need separate consideration. When γ > 2, the only
case considered in WM07, the integral in (7.29) converges to a constant plus
O(exp(−(γ−2)z)) as z → +∞. That in turn means that Bφ decays upward like
exp(−z). When γ = 2, the integral in (7.29) asymptotes to a linear function of
z, and Bφ decays upward like z exp(−z). When 2 > γ > 1, the integral in (7.29)
increases upward like exp((2−γ)z), but Bφ still decays upward, like exp(−(γ−1)z).
In all these cases it is clear from (7.29) that Bφ does not have exactly the same
z-dependence aloft as does Bz, as might have been expected from advective–
diffusive balance with advection by constant downwelling uz = −1. Having the
same z-dependences would make the right-hand side of (7.22) vanish. Hence it is
the more or less subtle departures from advective–diffusive balance, including the
contribution to Bφ from the twisting of field lines by the differential rotation uφ,
that enable the Lorentz torque on the right of (7.22) to support the flow ur at all
altitudes.
The uφ field that does the twisting can be calculated next, from (7.24) and the
condition that the interior rotates uniformly, uφ → 0 as z → −∞. Again this
calculation depends on the hyperbolic character of the operator B · ∇. When
Bφ is given by (7.29) we have, uniquely,
uφ =
∫ z
−∞
(
uz
∂Bφ
∂z −
∂2Bφ
∂z2
) dz
Bz
(7.31)
showing that uφ is also r times a function of z alone. That is, the differential
rotation is shellular solid rotation. In the three cases γ > 2, γ = 2, and 2 > γ > 1,
the behaviours of uφ as z → +∞ are respectively uφ ∼ constant, uφ ∼ ±z, and
uφ ∼ ± exp((2−γ)z).
Cases with negative shear aloft, ∂uφ/∂z < 0 — especially the last case, with
exponentially-increasing negative shear aloft — are suggestive of a possible way
to match upwards to the observed, much stronger negative shear in the bulk of
the tachocline. By using (7.29) to eliminate Bφ from ∂/∂z of (7.31), then using
(7.25) to eliminate d2Bz/dz2 and (7.27) to eliminate duz/dz, we find
∂uφ
∂z =
uzur
B2z
− 1B2z
∂ur
∂z −
2ur
r
∫ z
−∞
ur
B2z
dz ≈ (γ−1) urB2z
as z →∞ . (7.32)
96 7. The High-Latitude Tachocline
The asymptotic behaviour on the right comes from the first two terms in the
exact expression. The third term involving the integral is smaller by a factor
O(exp(−γz)). Asymptotically, therefore, the sign of the shear ∂uφ/∂z aloft is
the same as the sign of ur aloft. We can therefore find solutions that match on
to the strong negative tachocline shear provided that there is an exponentially
weak poleward mass flux above the confinement layer. However, a more precise
description of such matching may well need to take account of other sensitivities
aloft (see discussions in §§10.3, 10.5), including a sensitivity to slight departures
from magnetostrophic balance due to small-scale MHD turbulence (e.g. Parfrey
& Menou, 2007, & refs.), as well as slight departures from perfect flatness.
Although the structure aloft is sensitive to the precise conditions as z → +∞, the
solutions (7.29) and (7.31) show that, by contrast, the structure of the rest of the
confinement layer is relatively insensitive. An illustrative solution is plotted in
figure 7.3. In this case the interior field Bi has Br/r = 1. The downwelling profile
uz(z) was adapted from the numerical solution shown in figure 7.2, as described
in §10.5, and has γ = 2.24. For other γ values the bottom half of each plot
looks qualitatively the same, whereas the top half may in some cases be wildly
different; again see §10.5.
Some three-dimensional streamlines and magnetic field lines corresponding to
the solution in figure 7.3 are plotted in figure 7.4, visualising how the prograde
Lorentz torque on the right of (7.22), associated with field-line curvature, balances
the retrograde Coriolis torque on the left of (7.22) and satisfies the overall torque
balance (7.2).
7.6 Finite stratification
The perfectly flat confinement-layer equations (7.22)–(7.26) are formally valid
only in the limit Ro, Ek → 0 and αT , αµ → ∞. Allowing for finite αT or αµ
introduces non-flat corrections to the solutions. In chapters 8 and 9 we consider
the scalings and dynamics within the confinement layer and helium sublayer
respectively, in order to estimate the magnitude of these corrections.
The full equations (7.7)–(7.12), with finite Ro, Ek, αT and αµ, are solved numeri-
7.6 Finite stratification 97
−1 0 1 2
−1
0
1
2
3
4
5
6
uz
ur
r
uφ
r
z
Velocity Field
0 1 2 3
−1
0
1
2
3
4
5
6
BzBr
r
Bφ
r
z
Magnetic Field
Figure 7.3: Vertical profiles from a perfectly flat solution of the confinement-layer equations
in the limit αT , αµ → ∞. The downwelling profile uz(z), solid curve on left, was chosen to
match the downwelling profile from the numerical solution shown in figure 7.2. For numerical
reasons, small adjustments were made to this profile in the “slippery” region z > 1.5δη —
see §10.5. In (7.30) the decay constant γ = 2.24, and the uφ profile therefore approaches a
constant like exp(−0.24z).
Figure 7.4: Streamlines (left) and magnetic field lines (right) in a perfectly flat solution of the
confinement-layer equations, whose vertical profiles are displayed in figure 7.3. The peripheral
shading indicates the helium settling layer and helium sublayer.
98 7. The High-Latitude Tachocline
cally in chapter 10. The numerical solutions verify the results of chapters 8 and 9,
and provide additional insights. In particular, they confirm that the thermal and
compositional stratification surfaces are very nearly flat within the confinement
layer and sublayer, provided that the stable stratification is realistically strong.
The vertical cross-section shown in figure 7.2 is taken from a numerical solution;
the slight concavity of the helium surfaces, for example, is barely visible.
The numerical solutions of chapter 10 also verify that the flows within the confine-
ment layer and sublayer are magnetostrophic to excellent approximation. How-
ever, magnetostrophic balance becomes more delicate in the upper part of the
flow, wherein the Lorentz and Coriolis forces become vanishingly small. This
lends a certain “slipperiness” to the solutions, one symptom of which is the
abovementioned sensitivity of the perfectly flat solutions to the exact value of γ.
7.7 Summary
The results of section 7.5 may be summarised as follows:
• In the strongly stratified, magnetostrophic limit, the confinement-layer equa-
tions reduce to (7.22)–(7.26), which can be solved for any chosen profile of
uz(z).
• In order for the magnetic field |B| to be confined, i.e. to decay exponentially
as z → +∞, we require only that uz(z) ∼ −1 +O(exp(−γz)) as z → +∞,
with γ > 1.
• If uz(z) → 0 as z → −∞ then there is a unique solution with uniform
rotation and vanishing Alfve´nic torque below the confinement layer.
It is therefore apparent that plausible confinement-layer solutions can be ob-
tained with only minor restrictions on the downwelling profile uz(z). The physi-
cal interpretation of the considerable freedom this affords becomes clearer when
considering the boundary conditions for the numerical solutions (§10.3).
Chapter 8
The Magnetic Confinement Layer
Consider the scaling regime within the bulk of the confinement layer, well above
the helium sublayer. Because the photon mean free path makes the thermal
diffusivity relatively large, κ/η ∼ 3×104, the confinement-layer flow only weakly
perturbs the background thermal stratification. Consistently with (7.3) and (7.5)
we define the temperature perturbation T ′ so that T = z + const. + T ′. Then,
sufficiently close to the pole, the leading-order balance in the steady-state thermal
equation (7.11) involves only the vertical component uz of u,
uz =
κ
η∇
2T ′ . (8.1)
From this we estimate the magnitude of the thermal anomaly to be T ′ ∼ η/κ  1.
By contrast, the diffusivity of helium is relatively small, χ/η ≈ 2× 10−2, so the
diffusive flux of helium up into the bulk of the confinement layer is negligible. The
confinement layer is therefore well ventilated, i.e. the mean molecular weight µ
within the bulk of the confinement layer is close to being constant. Compositional
stratification is important only in the helium sublayer and the helium settling lay-
er, as indicated by the shading in figure 7.2 and further discussed in chapter 9.
The aim in this chapter is to estimate, using scale analysis, the characteristic
horizontal lengthscale of the confinement-layer flow, and hence to quantify the
“flatness” of the confinement layer. For the moment we do not assume axisym-
metry. We therefore begin with the full set of equations (7.7)–(7.12). We neglect
99
100 8. The Magnetic Confinement Layer
the αµ term in (7.7), since we expect the effects of compositional stratification
to be negligible within the confinement layer, and therefore equation (7.12) plays
no further role in the dynamics. We also use (8.1) in place of (7.11). Next,
we assume that the confinement-layer flow is perfectly magnetostrophic, i.e. we
neglect the O(Ro) and O(Ek) terms in (7.7). After taking the curl of (7.7) to
eliminate the pressure force, the steady-state confinement-layer equations are
−∂u∂z = −αT ez ×∇T
′ + B ·∇(∇×B)− (∇×B) ·∇B (8.2)
0 = ∇ · u (8.3)
0 = ∇× (u×B) +∇2B (8.4)
0 = ∇ ·B (8.5)
uz =
κ
η∇
2T ′ . (8.6)
We now assume that the horizontal scale of the confinement-layer flow is much
larger than the vertical lengthscale δη, by a factor F  1, say. We call F the
“flatness number”, since it measures the characteristic aspect ratio of the confine-
ment-layer flow. From (8.3) and (8.5) we deduce that the horizontal components
of the velocity and magnetic fields must be larger than the respective vertical
components by a factor of F . We therefore split all vector fields and operators
into their vertical and horizontal components:
u = uzez + uH , (8.7)
B = Bzez +BH , (8.8)
∇ = ∂∂zez + ∇H , (8.9)
∇2 = ∂
2
∂z2 +∇
2
H , (8.10)
and assume that
uz ∼ 1, |uH | ∼ F, (8.11)
Bz ∼ 1, |BH| ∼ F, (8.12)
∂
∂z ∼ 1, |∇H | ∼ 1/F. (8.13)
By substituting definitions (8.7)–(8.10) into (8.2)–(8.6), and neglecting contribu-
101
tions O(F−2), we obtain a set of boundary-layer equations for the confinement-
layer flow. These boundary-layer equations have ∂2/∂z2 in place of ∇2 in (8.4)
and (8.6), and the boundary-layer equivalent of (8.2) is
− ∂u∂z = −αT ez ×∇T
′ + B ·∇(∇×BH)− (∇×BH) ·∇B . (8.14)
We now define a set of “horizontally stretched” coordinates, with respect to
which both the horizontal and vertical lengthscales of the confinement-layer flow
are order unity. This is equivalent to rescaling all horizontal lengthscales and
horizontal vector-field components by a factor of F . In the stretched coordinates,
the vorticity equation (8.14) becomes
− ∂u∂z = −
αT
F 2ez ×∇T
′ + B ·∇(∇×BH)− (∇×BH) ·∇B , (8.15)
but, crucially, the other equations are invariant. By comparing (8.14) and (8.15)
we see that increasing the thermal stratification in the confinement layer, and
thus increasing the parameter αT , has the same effect as horizontally stretching
the confinement-layer flow (boundary conditions permitting). Indeed, since each
term in (8.15) is O(1), we deduce that the characteristic flatness of the confine-
ment layer is F 2 ∼ αTT ′. From (8.6) we have T ′ ∼ η/κ, and so
F 2 ∼ αT (η/κ) . (8.16)
However, (8.11)–(8.13) implicitly assume that the two horizontal components of
the vector fields u and B are of similar magnitude, and that Bz has a dimension-
less magnitude of order unity. In Wood & McIntyre (2010) a more careful scale
analysis was performed in a neighbourhood of the pole, assuming axisymmetry.
In that case, the azimuthal and colatitudinal components of u and B are gener-
ally not of similar magnitude. If the dimensional magnitude of Bz at the bottom
of the confinement layer is denoted as B, then
Bz ∼ Λ1/2 (8.17)
and uruφ
∼ BrBφ
∼ Λ (8.18)
where Λ is the Elsasser number based on the vertical component of the magnetic
102 8. The Magnetic Confinement Layer
field:
Λ = B
2
2Ωiη
. (8.19)
Because of axisymmetry, the thermal stratification surfaces are necessarily hori-
zontal at the pole. Wood & McIntyre (2010) showed that tilting of the thermal
stratification surfaces first becomes dynamically significant at a dimensionless
distance ∼ rT from the pole, where
r2T = αT (η/κ)min(Λ,Λ−1) . (8.20)
Wood & McIntyre’s rT is therefore analogous to the flatness number F derived
above. With NT ≈ 0.8 × 10−3s−1, U ∼ 10−5cm s−1, and Λ ∼ 1, we have rT =
F ∼ 4 × 103. With these numbers, which imply δη ∼ 0.4Mm, the dimensional
colatitudinal distance at which the tilting of the thermal stratification surfaces
first becomes significant is rT δη ∼ 1500Mm, roughly 2R away from the pole.1
Of course our cylindrical model with its assumption of uniform downwelling will
itself cease to apply, almost certainly, before we get as far as rT δη from the pole.
At some colatitude the downwelling must give way to upwelling. The confine-
ment-layer regime cannot then apply even qualitatively. Instead, the interior
magnetic field lines are free to advect and diffuse upward until they encounter
the magnetic flux pumping associated with the convective overshoot layer, as
assumed in figure 7.1.
With Λ ∼ 1 and the other numbers just given, the interior field Bi has dimen-
sional orders of magnitude typified by Bir ∼ 30cm s−1, equivalently 50 gauss, at
colatitude ∼ 30◦, or a third of those values at colatitude ∼ 10◦. Such |Bi| values
are far above the threshold, more like 10−2 gauss, for the field strength required
to enforce the Ferraro constraint in the interior (chapter 3). However, a sepa-
rate question concerns the possible range of field strengths, and MHD turbulent
intensities in the bulk of the tachocline, required to support the time-averaged
picture suggested in figure 7.1. There, field strengths greater than 10−2 gauss
might well be needed (e.g. appendix G).
1Even when the tilting is significant, the slopes of the thermal stratification surfaces are
far smaller than the geometrical aspect ratio r−1T , partly because of the parabolic shape of the
stratification surfaces and partly because T ′ is only a small perturbation to the background
stratification.
103
The range of interior field strengths accommodated by the confinement-layer
regime is determined by (8.20). The condition for the regime to apply at least
qualitatively within, say, 10◦ colatitude of the poles is rT δη & 90Mm. We can use
(8.20) together with realistic NT and diffusivity values to write this condition as
max(Λ,Λ−1) . 3× 102
( U
10−5cm s−1
)−4
. (8.21)
So for U ∼ 10−5cm s−1 the regime applies over a range of about five decimal
orders of magnitude in Λ. If the downwelling is weaker then a much wider range
of Λ values may become possible.
104
Chapter 9
The Helium Sublayer
9.1 Sublayer scalings
The helium sublayer marks the transition between the compositionally well-ven-
tilated confinement layer and the nearly impermeable, compositionally stratified
helium settling layer. The simplest description of the sublayer arises from neglect-
ing both viscosity and helium diffusivity, as in McIntyre (2007). The sublayer is
then infinitely thin, frictionless and impermeable, and the downwelling flow in the
confinement layer feeds into a “slip-flow” immediately above the sublayer. The
magnetic field lines crossing the sublayer are kinked by this slip-flow, to an extent
that is limited by magnetic diffusion. In fact, the sublayer in this description is
an example of the “magnetic Margules front” described in appendix D, except
that the slope of the sublayer is controlled by compositional rather than thermal
stratification. Hence there is a discontinuity in the mean molecular weight µ
across the sublayer.
To estimate the magnitude of the discontinuity in µ, and hence the slope of the
sublayer, we must go beyond this simple description and allow for finite helium
diffusivity χ. In order to simplify the analysis, we continue to neglect viscosity,
and assume that
χ  η . (9.1)
We also assume for the moment that the magnetic field Bi below the sublayer is
105
106 9. The Helium Sublayer
uniform and vertical. We may then equate the dimensional value of the vertical
field strength to B in the definition (8.19) of the Elsasser number Λ. We then
have
Bi = Λ
1
2ez. (9.2)
We also make the boundary-layer approximation, ∇2 ≈ ∂2/∂z2.
Allowing for finite χ values does not significantly alter the dynamics within the
confinement layer, provided that the sublayer remains nearly flat and that χ  η.
Therefore the timescale of the flow within the sublayer is inherited from the con-
finement layer. In particular, the flow in the sublayer has the same strain rate as
the confinement layer, U/δη. The thickness of the sublayer, δχ say, is determined
by a balance between advection and diffusion of helium in (7.12). The strain rate
U/δη must therefore match the helium diffusion rate χ/δ2χ. Since U/δη = η/δ2η,
it follows that
δχ/δη ∼ (χ/η)1/2  1 . (9.3)
With realistic solar values for the tachocline diffusivities, we have χ/η ≈ 0.02 (see
appendix A). The numerical solutions in chapter 10 confirm that this provides
sufficient scale separation for the sublayer to be regarded as distinct from the
confinement layer (see figure 7.2).
The magnetic diffusion rate η/δ2χ in the sublayer greatly exceeds the helium dif-
fusion rate χ/δ2χ. The flow within the sublayer can therefore induce only a small
departure B−Bi = B′, say, from the interior field Bi. In figures 7.2 and 7.4 the
field lines are hardly deflected as they cross the sublayer. We therefore analyse the
sublayer as a perturbation to the state with u = 0 and B = Bi, and assume that
the perturbation to Bi is linear. Under these assumptions, equations (7.7)–(7.12)
9.1 Sublayer scalings 107
become
RoDuDt + ez × u = −∇p + αT Tez − αµµez + (∇×B
′)×Bi (9.4)
0 = ∇ · u (9.5)
∂B′
∂t = ∇× (u×Bi) +
∂2B′
∂z2 (9.6)
0 = ∇ ·B′ (9.7)
DT
Dt =
κ
η
∂2T
∂z2 (9.8)
Dµ
Dt =
χ
η
∂2µ
∂z2 . (9.9)
Since κ  χ, the sublayer flow produces a negligible perturbation to the helium
settling layer’s uniform thermal stratification NT . We may therefore assume that
the thermal buoyancy force αT Tez within the sublayer arises entirely from the
uniform background thermal stratification. Since this force is curl-free, it can
be incorporated into the pressure field, and will therefore be neglected for the
remainder of this section.
On the timescale δη/U of the strain flow we may approximate magnetic diffusion
across the sublayer as instantaneous, and therefore neglect the left-hand side of
the induction equation (9.6). Substituting Bi = Λ
1
2ez into (9.6), we then find
that
0 = Λ
1
2
∂u
∂z +
∂2B′
∂z2 . (9.10)
Since both B′ and u vanish beneath the sublayer, we can integrate (9.10) to give
0 = Λ
1
2 u + ∂B
′
∂z . (9.11)
Using (9.11), we can now write the Lorentz force in the momentum equation (9.4)
as
(∇×B′)×Bi = (∇×B′)× (Λ
1
2ez) (9.12)
= −Λ
1
2 ∇B′z + Λ
1
2
∂B′
∂z (9.13)
= −Λ
1
2 ∇B′z − Λu . (9.14)
The term −Λ12 ∇B′z is curl-free, and can therefore be incorporated into the pres-
108 9. The Helium Sublayer
sure force. The term −Λu has the form of a Darcy or Rayleigh drag, showing
that the sublayer behaves like a porous medium on the timescale set by the strain
flow. The µ-choke together with the sublayer’s flatness and thinness act to keep
the flow nearly horizontal, compelling it to push past, and slightly deflect, the
magnetic field lines spanning the sublayer. Similar behaviour occurs in Hartmann
layers (e.g. Debnath, 1973), although in that case viscosity also contributes to
the balance of forces.
Neglecting terms of order Ro, we obtain the following set of magnetostrophic
equations for the sublayer flow:
ez × u = −∇p− αµµez − Λu (9.15)
0 = ∇ · u (9.16)
Dµ
Dt =
χ
η
∂2µ
∂z2 . (9.17)
A careful scaling analysis, presented in Wood & McIntyre (2010), confirms that
equations (9.15)–(9.17) accurately describe the dynamics of the sublayer in a
neighbourhood of the pole, even if the interior magnetic field is an axisymmetric
dipole, rather than the uniform vertical field assumed here. The same method
used to determine the “flatness” of the confinement layer in chapter 8 can now
be applied to the sublayer. Wood & McIntyre (2010) showed that the tilting
of the compositional stratification surfaces becomes dynamically significant at a
dimensionless distance from the pole ∼ rµ, where
r2µ = αµ(χ/η)min(Λ,Λ−1) . (9.18)
Comparing this with (8.20) we find
rµ
rT
∼ NµNT
(κχ)1/2
η
≈ 20 , (9.19)
and so the sublayer is “even flatter” than the confinement layer, as was asserted
in §7.5.
9.2 The subtail 109
9.2 The subtail
Within the lower extremity of the sublayer, or “subtail”, the MMCs decay expo-
nentially with depth as the µ-choke takes hold. Within this subtail, the helium
settling layer suffers only small perturbations to its otherwise uniform composi-
tional stratification, ∂µ/∂z = −1, and we may therefore perform a linear analysis
similar to those in §4.2.1 and §4.2.2. We denote the perturbation to µ by µ′. In
the steady state, equation (9.17) may then be approximated as
− uz =
χ
η
∂2µ′
∂z2 . (9.20)
We now take the curl and double-curl of (9.15), and apply (9.16), in order to
obtain
− ∂u∂z = αµez ×∇µ
′ − Λω, (9.21)
−∂ω∂z = αµ∇
2µ′ez − αµ∇
∂µ′
∂z + Λ∇
2u , (9.22)
where ω = ∇× u is the vorticity. Taking the vertical components of (9.21) and
(9.22), and making the boundary-layer approximation ∇2 ≈ ∂2/∂z2, we obtain
− ∂uz∂z = −Λωz (9.23)
−∂ωz∂z = αµ∇
2
Hµ′ + Λ
∂2uz
∂z2 (9.24)
where ∇2H = ∇2 − ∂2/∂z2 is the horizontal Laplacian. We can combine (9.23)
and (9.24) with (9.20) to yield a single equation for the perturbation µ′,
αµ(η/χ)∇2Hµ′ =
(
Λ + Λ−1
) ∂4µ′
∂z4 . (9.25)
In dimensional variables, this equation is
N2µ
2Ωiχ
∇2Hµ′ =
(
Λ+ Λ−1
) ∂4µ′
∂z4 , (9.26)
which is closely analogous to equation (4.9) in §4.2.2, once we identify ∇2H =
−1/L2. In §4.2.2, the horizontal scale L was determined by the global geometry
110 9. The Helium Sublayer
of the tachocline. Here, we expect the subtail to inherit the horizontal scale of the
sublayer. Since (9.19) implies that the sublayer is even flatter than the confine-
ment layer, the sublayer in turn inherits the horizontal scale of the confinement
layer, rT δη. If we replace ∇2H → −1/(rT δη)2 in the dimensional subtail equation
(9.26), and apply definitions (9.3) and (9.18), we obtain
δ4χ
∂4µ′
∂z4 ≈ −
r2µ
r2T
µ′ . (9.27)
Therefore the vertical scale of the subtail, δ` say, is smaller than the vertical scale
of the bulk of the sublayer δχ, by a factor
δ`
δχ
=
(rT
rµ
)1/2
≈ 0.2 . (9.28)
In the sublayer analysis of §9.1 it was assumed that viscous forces are negligible.
However, this approximation may break down within the subtail, if the length-
scale δ` is sufficiently small. If we allow for finite viscosity, then equation (9.26)
becomes
N2µ
2Ωiχ
(
Λ− δ2Ek
∂2
∂z2
)
∇2Hµ′ =
[
1 +
(
Λ− δ2Ek
∂2
∂z2
)2
]
∂4µ′
∂z4 (9.29)
where δEk = (ν/2Ωi)1/2 is the Ekman lengthscale. So viscous effects are negligible
within the subtail only if the Darcy friction from the magnetic field dominates
the fluid viscosity, i.e. if
Λ  δ2Ek/δ2` (9.30)
= Ek (rµ/rT )(η/χ) (9.31)
∼ 10−6. (9.32)
This implies a lower bound on the vertical field strength B of approximately
10−4 gauss.
Chapter 10
The Numerical Model
10.1 Numerical complexity
Unfortunately it is not possible to solve the complete set of equations (7.7)–
(7.12) analytically. To go beyond the perfectly flat analytical solutions described
in §7.5 we must compute solutions numerically. One way to gauge the complexity
of this task is to count the number of spatial and temporal derivatives in these
equations. The order of the complete system can be determined by considering
small-scale perturbations to an idealised background. In the absence of compo-
sitional stratification, such perturbations are governed by equations (C.20) and
(C.21) of appendix C. In (C.21), the term with the highest number of derivatives
is
( ∂
∂t − κ∇
2
)( ∂
∂t − ν∇
2
)2 ( ∂
∂t − η∇
2
)2
∇2 (10.1)
from which we deduce that the system is 12th order in space and 5th order in
time. Allowing for compositional stratification introduces an additional factor of
( ∂
∂t − χ∇
2
)
(10.2)
and so the system becomes 14th order in space and 6th order in time.
The large number of spatial derivatives relative to temporal derivatives arises
from the presence of diffusion in each of the prognostic equations, which render
111
112 10. The Numerical Model
the equations “almost parabolic” on small spatial scales. However, the system
remains partially elliptic because of the incompressibility constraint (7.8), which
leads to the ∇2 factor in (10.1). The system also supports hyperbolic behaviour,
for example through Alfve´n wave propagation.
From the discussions in chapters 8 and 9 it is clear that the confinement-layer
flow incorporates a wide range of spatial and temporal scales. Resolving all of
these scales in a numerical model is computationally demanding. Since we are
primarily interested in steady states, it would be desirable to solve the steady-
state equations directly, as in the model of Garaud & Garaud (2008). However,
the numerical scheme employed by Garaud & Garaud is not guaranteed to con-
verge to a solution when applied to a nonlinear system such as (7.7)–(7.12), and
does not provide any information regarding the dynamical convergence towards
the steady state (see §4.2.4).
An alternative approach is to neglect any terms in equations (7.7)–(7.12) that
we expect to be small. In particular, we might choose to neglect the O(Ro)
and O(Ek) terms in (7.7), thereby imposing magnetostrophic balance directly.
However, this approach leads to pathological behaviour at small scales, for reasons
discussed in appendix H.
Instead, we have chosen to solve the full set of time-dependent equations numeri-
cally. Although this is computationally demanding, we can be confident that the
system will eventually converge to a steady state, provided that the confinement
layer is dynamically stable. Timestepping the equations in a physically mean-
ingful way also provides some insight into the dynamics of the confinement layer
(§I.1).
Because of the computational complexity of the task, we have not attempted to
solve the global tachocline problem. We consider only the high-latitude confine-
ment layer illustrated in figure 7.2. Fitting the high-latitude confinement layer
into a global tachocline model remains a challenge for the future, and will require
the quantification of turbulent processes in the tachocline and convection zone.
In particular, a global model would need to include realistic descriptions of mag-
netic flux pumping and turbulent magnetic diffusion in the bulk of the tachocline.
Also crucial to the global picture is the turbulent gyroscopic pumping of polar
10.2 The numerical scheme 113
downwelling, and compensating upwelling in lower latitudes, together with the
resulting global-scale pattern of heat flow and the feedback on NT within the
tachocline.
The numerical scheme employed is described briefly in the next section, with
further technical details in appendix H. Several consistency tests have been per-
formed on the numerical solutions, the results of which are set out in appendix I.
10.2 The numerical scheme
We have written a numerical code to solve the axisymmetric version of (7.7)–
(7.12) in a cylinder of radius rd, say. The computational domain is taken to be a
rectangular grid in r and z, with regular grid intervals ∆r and ∆z, as illustrated
in figure 10.1.
r
z
∆r
∆z
r = 0 r = rd
Figure 10.1: The computational domain used for numerical solution.
Since we are solving a local problem, we are required to specify boundary con-
ditions (listed in §10.3). This inevitably involves artificial choices. In order to
minimise the number of boundary conditions required by the numerical scheme,
we have chosen to evaluate all spatial derivatives using low order finite differences.
In particular, all first-order spatial derivatives are evaluated using two-point, one-
114 10. The Numerical Model
sided (first-order) finite differences, and all second-order spatial derivatives are
evaluated using three-point, centred (first-order) finite differences. The direc-
tions of the one-sided finite differences are chosen to ensure numerical stability
at the grid scale (further detail in appendix H). The prognostic equations are
timestepped explicitly using a (first-order) Euler scheme. Since our primary in-
terest is in the steady state, there would be no benefit from using a more accurate,
higher-order timestepping scheme.
As explained in appendix H, the vertical grid interval ∆z must be small enough to
resolve the helium sublayer accurately, and the timestep ∆t must be small enough
to resolve thermal diffusion at this scale. In terms of the dimensionless units1
defined in chapter 7, this imposes the following upper bound on the dimensionless
timestep ∆t:
∆t < χ/κ . (10.3)
We anticipate that the system will converge to a steady state on the timescale of
magnetic diffusion across the confinement layer, which is order unity in the nondi-
mensionalised system. The disparity between the diffusion rates of helium and
temperature is therefore the principal factor determining the number of timesteps
required for convergence to the steady state.
10.3 Boundary conditions
We impose a dipolar magnetic field structure underneath the confinement layer,
with Br ∝ r, and a uniform downwelling of dimensional magnitude U from the
field-free region aloft. In the example shown in figure 7.2 the numerical domain
was defined by 0 6 r 6 5, i.e. rd = 5, and −1 6 z 6 6, one dimensionless unit
taller than shown in the figure. We imposed uz = −1 at z = 6 and Br/r = 1 at
z = −1. The imposed value of Br/r determines the characteristic magnitude of
the poloidal magnetic field components in the confinement layer, as can be seen
from the analytical solutions of §7.5 (e.g. dashed curve in figure 7.3). This in
turn determines the magnitude of the Elsasser number Λ defined by (8.19). In
1We use dimensionless units for the remainder of this chapter, except where explicitly stated
otherwise.
10.3 Boundary conditions 115
fact,
Λ ∼ B2z (10.4)
∼ (Br/r)2 . (10.5)
From the analytical solutions of §7.5, we anticipate that the poloidal components
of the magnetic field will decay with altitude like exp(−z) above the confine-
ment layer. We therefore impose that ∂Bz/∂z = −Bz at z = 6, to simulate the
presence of such an exponential tail in the region above the numerical domain.
Since Bz is exponentially small at the top of the domain, we do not expect the
precise choice of boundary condition there to significantly influence the steady
state. This expectation is confirmed in appendix I. We also impose that Bφ = 0
at z = 6, to ensure that the Maxwell stress over the top of the domain vanishes.
As shown in §7.5, the bulk of the confinement layer is relatively insensitive to
conditions within the field-free region aloft, and in particular to the vertical shear
∂uφ/∂z aloft. In the numerical solutions, the vertical shear aloft is tied to the
temperature distribution via the thermal-wind equation (7.18), and is therefore
sensitive to the upper boundary condition for the temperature. This behaviour
is not seen in the perfectly flat solutions derived in §7.5, because equation (7.18)
becomes degenerate in the perfectly flat limit where αT →∞ and ∂T/∂r → 0.
The “true” upper boundary condition for temperature can only be determined
by solving the global-scale heat flow problem, which is beyond the scope of the
present study. For simplicity, we impose T = const. and ∂uφ/∂z = 0 at z = 6,
which is consistent with equation (7.18) and the imposed uniform downwelling,
and also ensures that no viscous torque is exerted on the top of the domain.
Matching the confinement layer on to the observed negative vertical shear in the
bulk of the tachocline is beyond the scope of this study, for reasons mentioned in
§7.5 (but see discussions in Part IV).
At the periphery of the domain, the artificial cylindrical surface r = rd = 5,
the numerical algorithm requires us to specify the thermal and compositional
stratification profiles T (z) and µ(z), in a manner consistent with scalings in the
confinement layer and helium sublayer (chapters 8 and 9). In this way we ar-
tificially fix the altitude of the helium sublayer, which in a global model would
116 10. The Numerical Model
be determined rather by the competition between the burrowing of MMCs and
the diffusion of the interior field (see Part IV). The z-origin in our local, Boussi-
nesq model is arbitrary, and has been chosen to coincide approximately with the
altitude of the helium sublayer.
The peripheral T (z) and µ(z) profiles strongly influence the velocity field, which
is tightly linked to the two stratifications by equations (7.11) and (7.12). In
the limit of infinitely strong stratification described in §7.5, equations (7.11)
and (7.12) reduce to (7.21), and so T (z) and µ(z) can no longer be specified
independently. In §7.5, T (z) and µ(z) were both determined by (7.21), up to
boundary conditions, as soon as uz(z) was specified. The T (z) and µ(z) profiles
obtained in this manner can be used as peripheral boundary conditions for the
numerical solutions.
In addition to T (z) and µ(z), the numerical solutions require a third vertical
profile to be specified at the periphery. An obvious choice is the vertical profile
of Maxwell stress, BrBφ, which represents the field lines’ connection to lower
latitudes and the Alfve´nic torque exerted therefrom. As discussed in §7.3, we
expect the Maxwell stress over the periphery to be prograde within the confine-
ment layer, and to vanish both above and below the confinement layer. However,
in practice, we found it more convenient to specify Bφ(z) at the periphery, since
Bφ decays less rapidly with height above the confinement layer than does the
Maxwell stress, and is therefore less susceptible to numerical truncation errors.
The freedom to specify Bφ(z) is also present in the analytical framework of §7.5,
but was given up in order to ensure the vanishing of the Alfve´nic torque as
z → −∞. This condition uniquely determined Bφ(z) via the expression (7.29).
Also allowed by the analytical framework is the freedom to specify uφ(z), which we
similarly gave up in order to ensure solid rotation as z → −∞, thus determining
uφ(z) via (7.31).
Within the analytical framework both uφ and Bφ can be determined from their
values at any single location on each magnetic field line. This is because of
the Alfve´nic coupling along the field lines, expressed by the B · ∇ operator in
equations (7.22) and (7.24). There is slightly less freedom within the numerical
framework. The time-dependence, in the equations solved numerically, replaces
10.3 Boundary conditions 117
static Alfve´nic coupling by Alfve´nic wave propagation, requiring a single bound-
ary condition at either end of each field line. This is analogous to the need for
a single boundary condition at either end of an elastic string in motion. Having
already chosen to specify Bφ(z) at the periphery r = 5 we must therefore spec-
ify uφ(r) at the bottom of the domain, i.e. at z = −1. Alternatively, we could
choose to specify Bφ(r) at z = −1 and uφ(z) at r = 5. Although both choices
are mathematically valid, we found it more convenient to specify Bφ(z) at the
periphery, for reasons discussed in §10.4.
At z = −1 we impose uφ = 0, i.e. that the interior is in solid rotation, as expected
from the global picture. However, care must be taken to ensure that an Ekman
layer does not form at the bottom of the domain. Uniform rotation must therefore
be imposed indirectly, via the azimuthal component of the induction equation
(7.9), rather than via the azimuthal component of the momentum equation (7.7).
As shown in appendix D, any discontinuity in ∂Bφ/∂z over a horizontal surface
must be balanced by a discontinuity in uφ. More precisely,
−
[∂Bφ
∂z
]
= Bz [uφ] (10.6)
in dimensionless units (cf. equation (D.8)) where square brackets indicate a dis-
continuity. By imposing a discontinuity in ∂Bφ/∂z across the lower boundary of
the numerical domain, with amplitude equal to the value of Bzuφ immediately
below the lower boundary, we thereby ensure that uφ vanishes immediately above
the lower boundary. Using +/− to denote conditions above/below the boundary
z = −1, this boundary condition can be expressed as
[∂Bφ
∂z
]+
−
= Bz uφ|− . (10.7)
As the system converges to a steady state, we expect uφ to converge to zero below
z = −1, so that the discontinuity in ∂Bφ/∂z vanishes. Indeed, in appendix D we
show that such a discontinuity can be maintained only by a discontinuity in T
(or µ), which is ruled out by the presence of temperature and helium diffusion.
At z = −1, we further impose the conditions (7.3) and (7.4), equivalently ∂T/∂z =
1 and ∂µ/∂z = −1. We also impose uz = 0, i.e. that there is no meridional flow
118 10. The Numerical Model
into the interior.
The remaining boundary conditions are used to promote smoothness in the so-
lutions, and thereby to minimise spurious boundary effects in the steady state.
The complete set of boundary conditions is shown in figure 10.2. The boundary
conditions shown between braces are only needed to evaluate the terms of order
Ro and Ek in the momentum equation (7.7), and so do not significantly affect
the steady state. Since the system is 14th order in both z and r we require 7
boundary conditions on each boundary. We also need to impose regularity at the
coordinate singularity (see appendix H for details).
r
z
uz = −1,
{∂ωφ
∂z = 0
}
,
{∂uφ
∂z = 0
}
,
∂Bz
∂z = −Bz , Bφ = 0, T = 0, µ = 0
uz = 0,
{∂ωφ
∂z = 0
}
,
{
(∇2 − r−2)uφ = 0
}
,
Br/r = 1,
[∂Bφ
∂z
]+
−
= Bzuφ|− ,
∂T
∂z = 1,
∂µ
∂z = −1
∂
∂r
(ur
r
)
= 0,
{ ∂
∂r
(ωφ
r
)
= 0
}
,
{ ∂
∂r
(uφ
r
)
= 0
}
,
∂Bz
∂r = 0,
Bφ = Bφ(z),
T = T (z),
µ = µ(z)
Figure 10.2: The boundary conditions used for numerical solution.
10.4 The numerical solutions
Computing limitations preclude a perfect match to the real Sun’s parameter
values. They also require a slight modification to (7.7)–(7.12), explained in ap-
10.4 The numerical solutions 119
pendix H, in which artificial horizontal diffusivities νH , χH are introduced in
order to maintain numerical stability. However, from the scale analyses outlined
in chapters 8 and 9 we may identify the conditions most essential for achieving a
qualitatively similar parameter regime — that is, qualitatively similar to a regime
with a perfect parameter match to the real Sun. These essential conditions are:
1. The Rossby number Ro should be small in comparison with unity, so that
the steady state is close to magnetostrophic.
2. The thermal diffusivity κ should be large in comparison with the magnetic
diffusivity η, so that the confinement-layer flow only weakly perturbs the
background thermal stratification.
3. The confinement layer and helium sublayer should both be reasonably flat,
at least within the numerical domain r 6 rd. With (9.19) in mind, we also
take rµ > rT , where rT > rd > 1.
4. The helium diffusivity χ should be small in comparison with the magnetic
diffusivity η, so that the helium sublayer is thinner than, and therefore
distinct from, the magnetic confinement layer.
5. The viscosity ν should be small enough that viscous effects are negligible
throughout the confinement layer and sublayer. This condition implies a
lower bound on Λ, given by (9.31), or equivalently an upper bound on Ek.
Leaving νH and χH aside for the moment (see Appendix H) we can characterise
the system by seven dimensionless parameters, including the Elsasser number Λ,
which enters through the boundary condition for Br at z = −1, as described in
§10.3. Table 10.1 presents the other six dimensionless parameters, with nominal
solar values alongside the values used in the numerical solution presented in
figure 7.2 (and also in figures 10.3, 10.4 and 10.5).2 The last column echoes
aspects of the qualitative parameter conditions just stated. The nominal solar
values are based on the parameter values quoted in appendix A, and also assume
U = 10−5cm s−1.
2In figures 7.2 and 10.3 the uppermost part of the numerical domain, 5 6 z 6 6, is not
shown.
120 10. The Numerical Model
Table 10.1: Parameter values and conditions; see text.
Dim’less Nominal Value for Condition
parameter solar value num. sol’n
Ro 5× 10−8 10−2  1
κ/η 3× 104 102  1
αT (η/κ) 2× 107 50 > max(Λ,Λ−1)a
(rµ/rT )2 3× 102 2 > 1
χ/η 2× 10−2 2× 10−2  1
Ek 3× 10−9 2× 10−4  Λ(rT/rµ)(χ/η)
a This condition follows from rT > 1, with rT defined by (8.20).
Figure 10.3 shows plots of the steady-state streamlines and magnetic field lines
from the numerical solution whose parameter values are listed in table 10.1, and
whose meridional cross-section was presented in figure 7.2. Figure 10.4 shows the
vertical profiles of uz, ur/r, uφ/r, Bz, Br/r and Bφ/r on the rotation axis, from
the same numerical solution.
Figure 10.3: A numerical solution of the confinement-layer equations with r2T ≈ 14 and
Elsasser number Λ ≈ 3.5. Other parameter values are given in the second-last column of
Table 10.1. From the same solution as shown in figure 7.2.
From figure 10.4 we see that Λ = B2iz|z=0 ≈ 3.5. The nominal range of Λ for
which solutions are nearly flat within dimensionless radius r is determined by
(8.20) with rT > r. With the parameter values in the second-last column of
Table 10.1, rT > r is equivalent to
max(Λ,Λ−1) < 50/r2 . (10.8)
10.4 The numerical solutions 121
−1 0 1 2
−1
0
1
2
3
4
5
6
uz
ur
r
uφ
r
z
Velocity Field
0 1 2 3
−1
0
1
2
3
4
5
6
BzBr
r
Bφ
r
z
Magnetic Field
Figure 10.4: Vertical profiles from a numerical solution of the confinement-layer equations
with r2T ≈ 14 and Elsasser number Λ ≈ 3.5. Other parameter values are given in the second-last
column of Table 10.1. From the same solution as shown in figure 7.2.
Since Λ ≈ 3.5, (10.8) is violated near the edge of the domain, r = rd = 5.
Nevertheless, as seen for instance in figure 7.2, the solution remains fairly close
to being flat. This shows that the confinement-layer structure is robust even at
colatitudes for which the analysis of §7.5 breaks down.
The numerical solutions confirm that the confinement layer and helium sublayer
are magnetostrophically balanced (see plots in §I.2). However, magnetostrophic
balance holds less robustly in the upper part of the domain. In this upper region
the Coriolis and Lorentz forces are exponentially small, and the effects of small
but finite Ro and Ek are therefore relatively more significant. This leads to
difficulties when comparing the numerical solutions directly with the analytical
solutions of §7.5, as will be described in §10.5.
For the numerical solution shown in figures 7.2, 10.3 and 10.4, the vertical profiles
of T , µ and Bφ at r = 5 were initially taken from an analytical solution of the kind
described in §7.5. The resulting steady-state meridional flow, and the Coriolis
torque it exerts on each field line, cannot be precisely known in advance. So
the steady state found, with this choice of Bφ, will generally include a non-zero
profile of Alfve´nic torque exerted at the bottom of the computational domain,
122 10. The Numerical Model
z = −1. This is inconsistent with the picture presented in figure 7.1, in which
Bi is dipolar and no Alfve´nic torque is exerted on the top of the helium settling
layer. We therefore make iterative adjustments to the Bφ(z) profile at r = 5,
aiming to make the Alfve´nic torque at z = −1 vanish in the steady state.
It was found that making small changes in the peripheral Bφ(z) profile typically
leads to a new steady state in which the meridional velocity within the confine-
ment layer is essentially unchanged. Therefore the steady-state poloidal magnetic
field is also essentially unchanged. Since the numerical solutions are very nearly
magnetostrophic, this observation suggests a simple algorithm for determining
the required adjustment to Bφ(z). Equation (7.16) implies that any change in
rBφ at r = 5 is transmitted to the bottom of the domain along the poloidal
field lines. This is the same “static Alfve´nic coupling” described in §10.3. So
the required adjustment to Bφ(z) at r = 5 can be determined directly from the
steady-state Bφ(r) profile at z = −1.
In practice, slight departures from magnetostrophic balance mean that this ad-
justment procedure never works perfectly at the first attempt, but nonetheless
the Alfve´nic torque at the bottom of the domain generally decreases with each
successive adjustment to Bφ(z), and after several such adjustments is close to
zero. In figure 10.5 we display the steady-state profile of rBφ in a vertical cross-
section through the same solution shown in figure 7.2. Figure 10.5 also displays
some of the poloidal field lines in the confinement layer. It is apparent that the
Alfve´nic torque on the bottom of the domain is close to zero, and that rBφ is
roughly constant along each magnetic field line in the region below the helium
sublayer, where the left-hand side of (7.16) vanishes.
10.5 Upper-domain “slipperiness”
On the rotation axis, where the profiles in figure 10.4 were taken, the stratification
surfaces are flat for any finite αT and αµ. If the solution were in perfect mag-
netostrophic balance then we could use its uz(z) profile to calculate the other
field components on the axis by the procedure for constructing perfectly flat
solutions described in §7.5. But, as described above, the numerical solutions are
10.5 Upper-domain “slipperiness” 123
 
 
−5 0 5
−1
0
1
2
3
4
5
0.5
1
1.5
2
Figure 10.5: The steady-state profile of √|rBφ| in a vertical cross-section through the same
numerical solution shown in figure 7.2. Poloidal magnetic field lines are shown in black.
generally not in perfect balance, especially toward the upper part of the domain.
In particular, the numerical uz(z) will generally not conform to the decay law
(7.30) as z increases. So the perfectly flat solution obtained by this process
cannot precisely match the numerical solution, even on the axis. Indeed, such
a perfectly flat solution will often exhibit wild deviations from the numerical
solution toward the upper part of the domain. There, the delicate balance of
terms gives the dynamics a certain “slipperiness”, as already evidenced by the
upper-domain sensitivity to values of the decay constant γ in (7.30).
To enable a meaningful comparison between the numerical and analytical solu-
tions, we are therefore compelled to make small adjustments to uz(z) in the upper
domain, to make it conform to (7.30), before using it to compute a perfectly flat
solution. In the case shown here the required adjustment to uz(z) is very small
indeed. The solid uz curves on the left of figures 7.3 and 10.4 are practically
indistinguishable.
124
Part IV
Conclusions
125
Chapter 11
Discussion
The confinement-layer solutions described in Part III represent the first self-
consistent model of magnetic confinement and angular momentum balance in the
high-latitude tachocline. However, in order to formulate a complete tachocline
theory, we must describe how the high-latitude confinement layer fits into the
global picture suggested by figure 7.1.
Figure 11.1 sketches the way in which the confinement layer might fit into its
surroundings near the bottom of the polar tachocline (cf. figure 7.2). At the
periphery of the polar downwelling region, the field lines (solid) emerge from the
confinement layer on their way to lower latitudes. They will tend to splay out
and slant upward as they exit the downwelling region. The extra-polar tacho-
cline is therefore characterised by stronger magnetic fields with greater vertical
lengthscales. The magnetic Reynolds number, which is of order unity in the con-
finement layer, increases with colatitude, and the equatorward MMC at the base
of the tachocline begins to slant upward, flowing approximately along the field
lines.
At even greater colatitudes, the field lines continue to rise through the tachocline
until they encounter the convection zone’s overshoot layer, where they are held
horizontal by turbulent magnetic flux pumping. On the way we must expect
turbulent eddy fluxes to become increasingly important, decoupling the MMC’s
upwelling streamlines from the time-averaged field lines, and also limiting the
winding-up of the field lines by the tachocline’s latitudinal shear (see figure 7.1).
127
128 11. Discussion
Figure 11.1: Sketch of the magnetic confinement layer and its immediate surroundings at
the bottom of the high-latitude tachocline. Close to the pole the interior magnetic field (solid
lines) is confined by the downwelling MMC (dashed streamlines). The vertical scale has been
greatly exaggerated.
At all latitudes the stratification surfaces (shown dotted in figure 11.1) remain
nearly flat. However, the region ventilated by the MMC (unshaded in figure 11.1)
is deepest at the poles, in a manner reminiscent of the “polar pits” mentioned
in chapter 4. In today’s Sun, the ventilated region extends down only as far as
the top of the helium settling layer, indicated by dark shading. But in the early
Sun, before the formation of the settling layer, the tachocline’s MMCs might
well have burrowed much deeper into the radiative envelope (Wood & McIntyre,
2010). During the early main sequence, when the Sun rotated more rapidly, the
burrowing tendency would have been stronger (§2.4),1 and the polar pits might
perhaps have been many tens of megametres deep. Ventilation to greater depths
during the first gigayear or so of the Sun’s main-sequence evolution might prove
to be the long-sought explanation for the Sun’s observed lithium and beryllium
abundances (chapter 6).
The depth to which the MMCs burrow into the radiative envelope is determined
by the structure of the internal field Bi as well as by the strength of the gyroscopic
pumping in the layers above. In our local model of the high-latitude confine-
ment layer these factors are fixed by the boundary conditions, and so we cannot
quantitatively predict the depth of ventilation below the convection zone. But if
this ventilation is the dominant mechanism for solar lithium depletion, then the
polar pits would need to extend approximately 0.1R into the radiative envelope
1The interior field Bi also would have been stronger, but probably by only a modest fraction,
since the global-scale magnetic diffusion time is longer that the Sun’s present age.
129
(see chapter 6 and further discussion in chapter 12).
The global tachocline model that would be needed to test, and to begin to quan-
tify, the foregoing speculations would have to describe
1. the precise way in which turbulent stresses in the convection zone and tacho-
cline gyroscopically pump the polar downwelling responsible for confining
Bi;
2. the global-scale distribution of temperature and heat flow that fits in with
these MMCs;
3. the turbulent magnetic flux pumping by convective overshoot that we as-
sume confines Bi in extra-polar latitudes;
4. the extent to which the winding-up of the time-averaged toroidal field in
extra-polar latitudes (figure 7.1) is limited by turbulent eddy fluxes;
5. the reaction of the overlying turbulent layers to all of the above, espe-
cially the deficit in the convection zone’s differential rotation governing the
torques exerted from above, whether via gyroscopically-pumped MMCs or
via turbulent stresses in the bulk of the tachocline, or via both.
Progress on these formidable problems will of course depend on finding suitable
ways to model the turbulent processes. In chapter 12 we suggest ways in which
the high-latitude model presented in Part III could be extended to incorporate
some of the processes described above.
130
Chapter 12
Future Work
12.1 Gyroscopic pumping
In the confinement-layer model presented in Part III, the downwelling required
to confine Bi is imposed as an upper boundary condition. As described in §4.2.3,
downwelling in the high-latitude tachocline is driven by gyroscopic pumping in
the layers above, and is therefore closely connected to the vertical shear in those
layers. The case of gyroscopic pumping by a Tayler–Spruit dynamo (Spruit,
2002) was described by McIntyre (2007), whose arguments are summarised in
appendix G. If the angular momentum transport produced by the dynamo is
parametrised as a turbulent viscosity νT, then the vertical component of the
MMC driven in high-latitudes is
uz = νT
d lnΩ
dz (12.1)
where Ω is the absolute angular velocity and z is the vertical coordinate.
By introducing a turbulent viscosity into the upper layers of our confinement-
layer model, it should be possible to incorporate gyroscopic pumping directly,
rather than as a boundary condition. In order to drive downwelling, the angular
velocity Ω must decrease with altitude in the upper part of the domain. In the
numerical solution presented in chapter 10, constant temperature was imposed
at the top of the domain, and so the angular velocity has no vertical variation
131
132 12. Future Work
aloft (see §10.3). With more general temperature boundary conditions it should
be possible to obtain solutions with negative shear aloft, as is observed in the
high-latitude tachocline. The analysis in §7.5 suggests that the structure of the
confinement layer will not be strongly affected by the presence of vertical shear
and turbulent viscosity in the layers above.
12.2 Lithium burning
To verify the picture suggested in figure 11.1 we need to extend our confine-
ment-layer model horizontally, beyond the polar downwelling region. However,
to determine whether the ventilated polar pits can extend deep enough to burn
lithium, we will need to construct a global model of magnetic field confinement,
perhaps similar to that of Garaud & Garaud (2008, see §4.2.4). The analysis
performed in chapter 8 might allow us to develop their model to achieve field
confinement in the polar regions.
If the polar pits can be shown to extend deep enough to burn lithium, we must
then address how long they can persist, before being choked off by the formation of
the helium settling layer. Incorporating compositional stratification into global
tachocline models is problematic because of the tiny diffusivities involved, but
some insight might be gained from our local high-latitude model. By varying
the degree of compositional stratification in our numerical solutions, we hope
to describe in detail the formation of the helium sublayer, and the onset of the
µ-choke.
Part V
Appendices
133
Appendix A
Parameter Values
135
136 A. Parameter Values
Table A.1: Tachocline parameters at 0.7R, and global solar param-
eters. All values are taken from Gough (2007) except Nµ, which is
taken from McIntyre (2007, see chapter 5).
density ρ 0.21 g cm−3
pressure p 6.7× 1013 g cm−1s−2
temperature T 2.3× 106 K
sound speed c 2.3× 107 cm s−1
gravitational acceleration g 5.4× 104 cm s−2
density scale height Hρ 0.12R
pressure scale height Hp 0.08R
pressure Soret factor kp 2.9
temperature Soret factor kT 2.6
adiabatic index γ 1.665
thermal buoyancy frequency NT 8× 10−4 s−1
compositional buoyancy frequency Nµ 5× 10−4 s−1
total buoyancy frequency N 9× 10−4 s−1
magnetic diffusivity η 4.1× 102 cm2s−1
kinematic viscosity ν 2.7× 101 cm2s−1
thermal diffusivity κ 1.4× 107 cm2s−1
helium diffusivity χ 8.7 cm2s−1
total solar mass M 1.99× 1033 g
total solar radius R 6.96× 1010 cm
interior angular velocity Ωi 2.7× 10−6 s−1
Appendix B
Meridional Circulation
Mean meridional circulations (MMCs) arise naturally in rotating fluids that are
driven away from uniform rotation. As a simple example, we consider here an
incompressible fluid in approximately uniform rotation with angular velocity Ω.
In the frame rotating with the fluid, the momentum equation is then
Du
Dt + 2Ω× u = −
1
ρ∇p + g + ν∇
2u , (B.1)
where D/Dt is the Lagrangian or material derivative. We define the operator
〈.〉 as an average over time and longitude, and use it to decompose the veloc-
ity and pressure fields u and p into steady, axisymmetric and fluctuating, non-
axisymmetric parts: u = 〈u〉+ u′ and p = 〈p〉+ p′. By taking the average of the
momentum equation (B.1) we obtain the following equation for 〈u〉:
〈u〉 ·∇ 〈u〉+ 2Ω× 〈u〉 = −1ρ∇ 〈p〉+ g + ν
〈
∇2u
〉
− 〈u′ ·∇u′〉 . (B.2)
The final two terms on the right-hand side incorporate all the angular momen-
tum transport processes, including hydrodynamic turbulence and wave breaking,
that can act to drive the system away from uniform rotation. To properly de-
scribe these processes would require solving the complete set of nonlinear, three-
dimensional fluid equations. Here, we instead adopt a crude parametrisation,
replacing these two terms with a “generalised Rayleigh friction” that forces the
fluid towards a state of differential rotation. Using cylindrical polar coordinates
137
138 B. Meridional Circulation
(r, z, φ) aligned with the rotation axis, we write this “anti-frictional” force1 as
α(u˜φ − 〈uφ〉) eφ, where 〈uφ〉 is the averaged azimuthal velocity in the rotating
frame, and u˜φ is the steady, axisymmetric differential rotation that the system
is forced towards. The parameter α is a relaxation rate, which quantifies how
strongly the system is driven away from uniform rotation, and might be inter-
preted physically in terms of the turnover rate of the local convective eddies.
Provided that the forcing produces only small departures from uniform rotation
(i.e. u˜φ/r  |Ω|) and is not too strong (α . |Ω|) we may neglect the nonlinear
term on the left-hand side of B.2. We then have
2Ω× 〈u〉 = −1ρ ∇〈p〉+ g + α(u˜φ − 〈uφ〉) eφ . (B.3)
The azimuthal components of (B.3) and its curl are
2|Ω| 〈ur〉 = α(u˜φ − 〈uφ〉), (B.4)
2|Ω| ∂∂z 〈uφ〉 = 0. (B.5)
Equation (B.5) is an explicit demonstration of the Taylor–Proudman theorem,
and equation (B.4), taken together with boundary conditions, describes the gy-
roscopic pumping of the meridional circulation. For example, if the domain is
bounded by axisymmetric impermeable surfaces z = z1(r) and z = z2(r) then
mass conservation implies that
∫ z2
z1 dz 〈ur〉 = 0 for all r. Hence, by integrating
(B.4) between z = z1 and z = z2, we find that
∫ z2
z1
dz (u˜φ − 〈uφ〉) = 0 (B.6)
⇒ 〈uφ〉 =
∫ z2
z1 dz u˜φ
z2 − z1
(B.7)
i.e. 〈uφ〉 is the axial average of u˜φ.
So the anti-frictional forcing produces only as much differential rotation as is
permitted by the Taylor–Proudman theorem. The residual forcing serves to drive
the meridional circulation. The only effect of increasing the strength of the forcing
(i.e. increasing α) is to drive a faster circulation.
1We adopt here a similar parametrisation to that of Garaud & Acevedo Arreguin (2009),
except that here the forcing acts only in the azimuthal direction.
139
As an illustrative example, we now apply these results to the case of a spherical
shell with the same aspect ratio as the Sun’s convection zone. We choose u˜φ
to mimic the observed pattern of differential rotation. Since the results do not
depend sensitively on the precise form of u˜φ we adopt a simple “hyperbolic”
approximation for the solar angular velocity contours:
u˜φ/r = (ar)2 − (bz)2 (B.8)
where a and b are constants. In Figure B.1 we show u˜φ/r for a particular choice
of a and b, together with the resulting profile of 〈uφ/r〉. Streamlines of the merid-
ional flow are shown in Figure B.2; the shape of the streamlines is independent
of both a and b.
Figure B.1: Left panel: the differential rotation u˜φ/r, chosen to approximate the rotation of
the convection zone. Right panel: the differential rotation 〈uφ/r〉 in our model.
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Figure B.2: Streamlines of the meridional flow in our model.
140
Appendix C
MAC Waves
We consider perturbations to a Boussinesq, magnetised fluid, with uniform grav-
itational acceleration g. We suppose that the unperturbed flow is uniform rota-
tion with angular velocity Ω, and adopt a frame of reference in which the fluid
is at rest. We further assume uniform thermal stratification and an unperturbed
background magnetic field of the form
B = B0 + 12J0 × x , (C.1)
where the vectors B0 and J0 are both uniform. Using primes to denote the
perturbations, the governing equations can be written in the form
( DDt − ν∇
2)ω′ + 12J0 × J
′ = (2Ω+ ω′) ·∇u′ + (B+B′) ·∇J′
− (J0 + J′) ·∇B′ + g ×∇T ′ (C.2)
∇ · u′ = 0 (C.3)
( DDt − η∇
2)B′ + 12J0 × u
′ = (B+B′) ·∇u′ (C.4)
∇ ·B′ = 0 (C.5)
( DDt − κ∇
2)T ′ = N
2
|g|2g · u
′ (C.6)
where ω′ = ∇×u′ is the vorticity, J′ = ∇×B′ is the electric current in suitable
units and T ′ represents fractional perturbations to the background temperature
141
142 C. MAC Waves
profile. The terms involving J0 on the left-hand sides of (C.2) and (C.4) arise
from the curvature of the background magnetic field B. The operator
D
Dt =
∂
∂t + u
′ ·∇ (C.7)
is the material derivative.
If the perturbations are sufficiently small then they can be described by the
linearised versions of (C.2), (C.4) and (C.6), which are
( ∂∂t − ν∇
2)ω′ + 12J0 × J
′ = (2Ω ·∇)u′ + (B ·∇)J′
− (J0 ·∇)B′ + g ×∇T ′ (C.8)
( ∂∂t − η∇
2)B′ + 12J0 × u
′ = (B ·∇)u′ (C.9)
( ∂∂t − κ∇
2)T ′ = N
2
|g|2g · u
′ . (C.10)
It proves useful to also consider the curls of (C.8) and (C.9), which are
−( ∂∂t − ν∇
2)∇2u′ − 12J0 ×∇
2B′ = (2Ω ·∇)ω′ − (B ·∇)∇2B′ − (J0 ·∇)J′
+ g∇2T ′ − g ·∇∇T ′ (C.11)
( ∂∂t − η∇
2)J′ + 12J0 × ω
′ = (B ·∇)ω′ . (C.12)
In the particular case where J0 and Ω are both aligned with g, the solutions
of the linearised equations are remarkably simple. Defining u′g, ω′g, B′g and J ′g
as the “downward” components of u′, ω′, B′ and J′ (i.e. the components in the
direction of g) we then need consider only the downward components of our linear
143
equations, which are
( ∂∂t − ν∇
2)ω′g = (2Ω ·∇)u′g + (B ·∇)J ′g − (J0 ·∇)B′g (C.13)
( ∂∂t − η∇
2)B′g = (B ·∇)u′g (C.14)
−∇2( ∂∂t − ν∇
2)u′g = (2Ω ·∇)ω′g −∇2(B ·∇)B′g
− (J0 ·∇)J ′g + |g|∇2HT ′ (C.15)
( ∂∂t − η∇
2)J ′g = (B ·∇)ω′g (C.16)
( ∂∂t − κ∇
2)T ′ = N
2
|g|u
′
g (C.17)
where ∇2H = ∇2 − |g|−2(g ·∇)2 is the horizontal Laplacian. These five equations
form a closed set and, moreover, all the differential operators commute. For
example,
[B ·∇, 2Ω ·∇] = (Ω× J0) ·∇ (C.18)
= 0. (C.19)
We can therefore combine all five equations into a single equation
M[u] = 0 (C.20)
where u represents any of the dependent variables {u′g, ω′g, B′g, J ′g, T ′}, and M is
the differential operator
( ∂∂t − κ∇
2)
[
( ∂∂t − ν∇
2)( ∂∂t − η∇
2)− (B ·∇)2
]2
∇2
+( ∂∂t − κ∇
2)
[
(2Ω ·∇)( ∂∂t − η∇
2)− (J0 ·∇)(B ·∇)
]2
+( ∂∂t − η∇
2)
[
( ∂∂t − ν∇
2)( ∂∂t − η∇
2)− (B ·∇)2
]
N2∇2H . (C.21)
The unstratified case has been studied by several authors (e.g. Malkus, 1967;
Craik, 1988; Dritschel, 1991). Malkus (1967) considered the case with B0 = 0 in
a spherical geometry, and hence this particular field configuration is often called
the “Malkus field”. Craik (1988) and Dritschel (1991) considered more general
144 C. MAC Waves
field configurations, and showed that certain solutions of (C.20) are also exact
solutions of the full, nonlinear equations.
In the special case where J0 = 0 (i.e. the unperturbed magnetic field is uniform)
equation C.20 applies even if Ω is not aligned with g. Furthermore, the system
supports plane wave solutions, known as “magneto–Archimedes–Coriolis” (MAC)
waves, with spatio–temporal dependence of the form exp(ik · x − iωt), provided
that the wave vector k and frequency ω satisfy the dispersion relation
0 = −ω(κ)
[
ω(ν)ω(η) − (B · k)2
]2 |k|2
+ω(κ)ω(η)ω(η)(2Ω · k)2
+ω(η)
[
ω(ν)ω(η) − (B · k)2
]
N2|kH |2 . (C.22)
In (C.22) we have defined
kH = k− |g|−2(g · k)2, (C.23)
ω(κ) = ω + iκ|k|2, (C.24)
ω(η) = ω + iη|k|2, (C.25)
and ω(ν) = ω + iν|k|2. (C.26)
The special cases of Alfve´n waves, internal gravity waves and Coriolis waves can
all be deduced from (C.22).
Appendix D
The Magnetic Margules Front
We consider a rotating, magnetised fluid within which there is a temperature dis-
continuity across a surface of negligible curvature; we call this surface a “front”.
We suppose that the baroclinicity associated with the temperature discontinuity
is supported by discontinuities in the components of the velocity u and electric
current J that lie within the front. We may describe this situation by considering
linear perturbations to a uniformly rotating system containing a uniform mag-
netic field. Adopting the same notation as in appendix C, the relevant equations
are (C.8) and (C.9), together with
J′ = ∇×B′ . (D.1)
Neglecting viscosity, and setting J0 = 0, the steady-state versions of (C.8) and
(C.9) are
0 = (2Ω ·∇)u′ + (B ·∇)J′ + g×∇T ′ (D.2)
and −η∇2B′ = (B ·∇)u′. (D.3)
Defining the normal to the front as the unit vector n, we decompose the vector
fields u′ and J′ into their components parallel and perpendicular to n:
u′ = u′‖n+ u′⊥ (D.4)
and J′ = J ′‖n+ J′⊥. (D.5)
145
146 D. The Magnetic Margules Front
After making a similar decomposition for the components of g, Ω and B, we can
use (D.1)–(D.3) to derive the following “jump conditions”:
[J′⊥] = n× [(n ·∇)B′⊥] (D.6)
0 = 2Ω‖[u′⊥] +B‖[J′⊥] + g⊥ × n[T ′] (D.7)
−η [(n ·∇)B′⊥] = B‖[u′⊥]. (D.8)
where the notation [.] refers to the discontinuity in a quantity across the front.
The jump conditions can be combined into a single equation relating the discon-
tinuity in T ′ to the discontinuities in [u′⊥]:
0 = 2Ω‖[u′⊥]− 2Ω‖Λ‖n× [u′⊥] + g⊥ × n[T ′] (D.9)
where for notational convenience we have defined an Elsasser number
Λ‖ =
B2‖
2Ω‖η
. (D.10)
The vector product of (D.9) with n is
0 = 2Ω‖n× [u′⊥] + 2Ω‖Λ‖[u′⊥] + g⊥[T ′]. (D.11)
Eliminating n× [u′⊥] between (D.9) and (D.11) we find
0 = 2Ω‖(1 + Λ2‖)[u′⊥] + g⊥ × n[T ′] + Λ‖g⊥[T ′]. (D.12)
In an axisymmetric system, the front must be a surface of revolution about the
rotation axis. Then g and n both lie within the meridional plane, and the az-
imuthal component of (D.12) is
0 = 2Ω‖(1 + Λ2‖)[u′φ] + (g⊥ × n) · eφ[T ′]. (D.13)
If Ω is anti-parallel with g, we deduce that the front has an angle of inclination,
ε say, given by the formula1
tan ε = −(1 + Λ2‖)
2Ω
g
[u′φ]
[T ′] . (D.14)
1This is an implicit formula, since in general Λ‖ depends on the angle ε.
147
This situation is illustrated in figure D.1. For a purely hydrodynamic system
(Λ‖ = 0), equation (D.14) reduces to the well-known Witte–Margules equation
for the slope of a Margules front.
Ω
g
n
ε
Figure D.1: Meridional cross-section through an axisymmetric front with inclination angle ε
and unit normal n.
148
Appendix E
Magnetic Pumping
In the presence of weakly non-linear, non-rotating, unstratified, “quasi-isotropic”
small-scale turbulence, and in the limit of large turbulent viscosity and magnetic
diffusivity, a weak, large-scale magnetic field experiences diamagnetic pump-
ing with a characteristic pumping velocity γ, say. If the turbulence is three-
dimensional, then
γ = −12∇ηe (E.1)
(e.g. Kitchatinov & Ru¨diger, 1992), where ηe(x) is the effective magnetic diffu-
sivity, the sum of microscopic and turbulent contributions. If the turbulence is
two-dimensional, then
γ = −∇ηe (E.2)
(Zel’dovich, 1957). The evolution of the large-scale field B, in the absence of
large-scale fluid motions, is then described by the mean-field induction equation
∂B
∂t = ∇× (γ ×B− ηe∇×B). (E.3)
For a pumping velocity of the form
γ = −c∇ηe, (E.4)
where c is constant, (E.3) can be written more compactly as
∂B
∂t = −∇× (η
1−c
e ∇× (ηceB)). (E.5)
149
150 E. Magnetic Pumping
The strength of the pumping can be measured in terms of the coefficient c, with
c = 0 corresponding to no pumping.
Following Kitchatinov & Ru¨diger (2008), we consider a three-layer solar model.
We assume that the effective magnetic diffusivity ηe takes a small value, η0 say,
within the radiative interior, but a much larger value, η1 say, within the convection
zone. Both η0 and η1 are assumed to be constant, implying a discontinuity in ηe at
the radiative–convective interface. The region outside the Sun is approximated as
a vacuum, wherein ∇×B = 0. We assume that B is poloidal and axisymmetric,
so that it can be defined in terms of a scalar potential A:
B = ∇× (Aeφ) (E.6)
where eφ is the unit vector directed azimuthally. The induction equation can
then be “uncurled” and written in the form
∂A
∂t = rη
1−c
e ∇ ·
(ηce
r2 ∇(rA)
)
(E.7)
where r is the perpendicular distance from the axis of symmetry. The discontinu-
ity in ηe at the base of the convection zone implies a delta function in the pumping
velocity γ, and leads to the following matching conditions at the interface with
the radiative envelope:
[A] = 0, (E.8)
[ηcen ·∇(rA)] = 0, (E.9)
where n is the normal to the interface. These matching conditions imply that
the vertical field n ·B is continuous at the interface, but that the horizontal field
n×B is discontinuous, and is smaller on the convection-zone side of the interface
by a factor (η0/η1)c.
We are interested in the long-time behaviour of the system, i.e. the behaviour
on timescales comparable to the global magnetic diffusion time in the radiative
interior. On such timescales we may assume that magnetic diffusion within the
151
convection zone is instantaneous, since η1  η0. Hence (E.7) becomes
∇ ·
( 1
r2∇(rA)
)
= 0 (E.10)
within the convection zone, which is equivalent to ∇ × B = 0. The location
of the solar surface is now immaterial, since we have a vacuum field everywhere
outside the radiative interior.
After a sufficiently long time, the magnetic field in the radiative interior will
become dipolar, since the dipolar component has the longest decay time. If the
radiative–convective interface is approximated as spherical, then the vacuum field
outside the radiative interior must be
A = r/|x|3 (E.11)
after this time. Hence n×B and n ·B are of similar magnitude at the base of the
convection zone, regardless of the pumping strength c. At the top of the radiative
envelope the matching conditions (E.8) and (E.9) then imply that the horizontal
field n×B is typically larger than the vertical field n ·B by a factor of (η1/η0)c.
Hence if (η1/η0)c  1 then the magnetic field is confined, i.e. predominantly
horizontal, below the base of the convection zone.
However, we note that solutions of equation (E.3) are generally not force-free.
Indeed, the Lorentz force will be particularly large at the base of the convection
zone, where the vertical discontinuity in the horizontal field implies the existence
of a current sheet. A complete and self-consistent model of magnetic confinement
by turbulent pumping would need to take into account the dynamical effects of
this force, including its effect on the turbulence and hence the pumping velocity
γ.
152
Appendix F
The Tayler Instability
In this section, we re-derive a simple case of Tayler instability described by Spruit
(1999). (Spruit considered more general toroidal field configurations than we do
here.) We consider perturbations to a uniformly rotating, stably stratified fluid
containing a magnetic field of the form
B = 12J0 × x, (F.1)
where J0 is uniform. We suppose that gravity g and the fluid’s angular velocity Ω
are both uniform and parallel to J0. As shown in appendix C, linear perturbations
in such a fluid can be described by the equation
M[u] = 0 (F.2)
where M is the operator
( ∂∂t − κ∇
2)
[
( ∂∂t − ν∇
2)( ∂∂t − η∇
2)− (B ·∇)2
]2
∇2
+( ∂∂t − κ∇
2)
[
(2Ω ·∇)( ∂∂t − η∇
2)− (J0 ·∇)(B ·∇)
]2
+( ∂∂t − η∇
2)
[
( ∂∂t − ν∇
2)( ∂∂t − η∇
2)− (B ·∇)2
]
N2∇2H . (F.3)
We suppose that the magnetic field is weak, in the sense that |J0|  |Ω|, and
that the system is “heavily stratified”, i.e. |Ω|  N . In the absence of thermal
153
154 F. The Tayler Instability
and magnetic diffusion, the system would be stable. In the presence of thermal
or magnetic diffusion, there are unstable modes that grow on a timescale τ ∼
|Ω|/|J0|2 (Pitts & Tayler, 1985). In order for diffusion to destabilise the motions,
they must have a small vertical scale, say δ. On the other hand, the motions
must have a large horizontal scale, in order to avoid significantly perturbing the
stratification surfaces. For such motions we may approximate ∇2 ≈ ∂2/∂z2,
where z is the vertical coordinate.
In the radiative envelope, where κ  η  ν, Tayler instability is possible only if
δ lies in the range
κ/δ2  τ−1 & η/δ2  ν/δ2 (F.4)
(Spruit, 1999). In that case we may neglect viscosity, and assume that thermal
diffusion is instantaneous, so that ( ∂∂t − κ∇2) → (−κ∇2) in (F.3). Since the
timescale for the instability τ is much longer than the rotation period, we also
make the magnetostrophic approximation, and therefore neglect all terms in (F.3)
containing the factor ( ∂∂t − ν∇2).
Defining cylindrical coordinates (r, φ, z) centred on the rotation axis, we seek
solutions to (F.2) of the form
u = uˆ(r) exp(inz + imφ− iωt) , (F.5)
where n is real and m is integer. Applying the approximations listed above, (F.2)
then becomes
(∇2H + L−2)u = 0 , (F.6)
with
L−2 = iκn
4|J0|2
N2(ω + iηn2)
[
(
1 + 4Ω(ω + iηn
2)
|J0|2m
)2
− m
2
4
]
. (F.7)
Necessary conditions for instability are that L2 > 0 and Imω > 0. After writing
ω as the sum of (dimensionless) real and imaginary parts
ω = |J0|
2
4Ω (ω˜R + iω˜I) (F.8)
we take the real and imaginary parts of (F.7), treating L2 as real, which leads to
the following equations for ω˜I and L:
(
ω˜I +
4Ωηn2
|J0|2
)2
= m2
[ 1
4m2 − (1 + ω˜R/m)2
1 + 2m/ω˜R
]
(F.9)
N2L−2
8Ωκn4 = −
( 1
m +
1
ω˜R
)[ 1
4m2 − (1 + ω˜R/m)2
1 + 2m/ω˜R
]1/2
. (F.10)
For a mode of the form (F.5) to be unstable, the right-hand sides of (F.9) and
(F.10) must be positive, which requires that m = 1 and −12 < ω˜R < 0.
If there are no further constraints on n and L then there are always unstable
modes. If, however, the domain is of finite horizontal extent, then there will be
an upper bound on L, say L0. This in turn implies a lower bound on the vertical
wavenumber n, which can be deduced from (F.10). But n cannot be too large,
or else the instability will be damped by magnetic diffusion, as can be seen from
(F.9). In fact, for a given value of ω˜R, instability is possible if and only if
(4Ωηn2
|J0|2
)2
<
1
4 − (1 + ω˜R)2
1 + 2/ω˜R
, (F.11)
with n given by
N2L−20
8Ωκn4 = −(1 + 1/ω˜R)
[ 1
4 − (1 + ω˜R)2
1 + 2/ω˜R
]1/2
. (F.12)
By considering all ω˜R in the range (−12 , 0), we deduce that the necessary and
sufficient condition for instability is
2Ωη2N2L−20
κ|J0|4
< maxω˜ ∈ (−12 , 0)
{
−(1 + 1/ω˜)
[ 1
4 − (1 + ω˜)2
1 + 2/ω˜
]3/2}
(F.13)
≈ 0.05 (F.14)
⇒ |J0|4 &
40Ωη2N2
κL20
. (F.15)
155
156 F. The Tayler Instability
Appendix G
The Tayler–Spruit Dynamo
We consider a horizontal layer of stably stratified rotating fluid, threaded by
a weak, approximately homogeneous magnetic field. If we impose vertical dif-
ferential rotation within this layer, then the weak field will be wound up into a
roughly axisymmetric toroidal field, in accordance with the principle of rotational
smoothing described in §3.2. Assuming that the vertical shear is not too large,
the toroidal field will continue to grow in strength, but remain stable, until the
Tayler-instability criterion (F.15) is exceeded. The instability generates a verti-
cal field component. If the imposed shear winds up this vertical field component
into new toroidal field sufficiently quickly, the instability can operate continually,
resulting in a dynamo (Spruit, 2002). Spruit argued that the dynamo would op-
pose the imposed shear through its Maxwell stress, yielding an effective turbulent
viscosity νT. By assuming that the vertical and azimuthal field components in
the saturated dynamo are correlated in the manner of the most unstable Tayler
mode, Spruit estimated this viscosity to be
νT ∼ Ωδ2 (G.1)
in a neighbourhood of the rotation axis, where δ is the vertical scale of the
instability. He also estimated the minimum vertical shear required to sustain
dynamo action as
∣
∣
∣
∣
dΩ
dz
∣
∣
∣
∣
∼ η/δ3 (G.2)
157
158 G. The Tayler–Spruit Dynamo
where z is the vertical coordinate. In the idealised scenario described in ap-
pendix F, the vertical scale of the Tayler instability can be deduced from (F.12):
δ−4 ∼ n4 ∼ N
2
ΩκL20
. (G.3)
McIntyre (2007) suggested that a Tayler–Spruit dynamo operating within the
high-latitude tachocline might gyroscopically pump the downwelling required to
confine the Sun’s interior magnetic field. The downwelling that can be pumped by
this mechanism can be estimated using Spruit’s estimate (G.1) for the turbulent
viscosity. We define cylindrical polar coordinates (r, φ, z) centred on the rotation
axis. Within a neighbourhood of the north pole, the angular velocity Ω and
density ρ within the tachocline may be assumed to depend only on height z.
We assume that the effective viscosity νe is also a function of z. Within the
bulk of the tachocline, we have νe = νT, with νT given by (G.1). Within the
laminar magnetic confinement layer at the base of the tachocline, the effective
viscosity drops to its microscopic value, which we shall take to be zero here.
The latitudinal flow pumped by the tachocline’s Tayler–Spruit dynamo can be
estimated by balancing the azimuthal Coriolis force against the divergence of the
turbulent stress:
2Ωρur =
d
dz
(
ρrνe
dΩ
dz
)
. (G.4)
Since νe(z) increases steeply with height within the tachocline, against a back-
ground of negative shear dΩ/dz, we anticipate that the gyroscopically pumped
flow will be poleward, with ur ∝ −r. Now applying mass conservation, ∇·(ρu) =
0, we deduce that
d(ρuz)
dz = −2ρur/r (G.5)
= − ddz
(
ρνe
d lnΩ
dz
)
. (G.6)
From the principle of downward control (Haynes et al., 1991) we expect the
converging poleward flow to turn downward. Assuming that there is no gyroscopic
pumping in the layers above the tachocline, the vertical flow uz within the bulk
of the tachocline must vanish. We can then integrate (G.6) vertically to calculate
the downwelling pumped into the top of the confinement layer. Neglecting the
159
small changes in density within the tachocline, we find
uz = νT
d lnΩ
dz . (G.7)
Following McIntyre (2007) we suppose that the turbulent stresses within the ta-
chocline prevent the vertical shear from growing significantly beyond the thresh-
old value given by (G.2). We can then apply the estimates (G.1) and (G.2) to
the downwelling formula (G.7) to obtain
uz = η/δ . (G.8)
The estimate (G.1) for the turbulent viscosity assumes that the magnetic field
generated by the Tayler–Spruit dynamo has the large-scale horizontal structure
associated with Tayler instability. However, Denissenkov & Pinsonneault (2007)
have argued that such a large-scale field would not be sustained in the nonlinear
regime, and that (G.1) therefore overestimates the angular momentum transport
that can be produced by a Tayler–Spruit dynamo. In that case, we should regard
(G.8) as an upper bound for the downwelling that can be pumped by a Tayler–
Spruit dynamo operating at threshold within the tachocline.
160
Appendix H
The Numerical Scheme
H.1 Equations
We wish to solve a suitable version of equations (7.7)–(7.12) in axisymmet-
ric cylindrical polar coordinates. We introduce streamfunctions Ψ and A, i.e.
azimuthal vector-potential components, for the poloidal velocity and magnetic
fields, such that
uz =
1
r
∂(rΨ)
∂r and ur = −
∂Ψ
∂z , (H.1)
Bz =
1
r
∂(rA)
∂r and Br = −
∂A
∂z , (H.2)
guaranteeing that the fields are divergence-free. In the fluid-dynamical literature
rΨ is sometimes called the Stokes streamfunction. The azimuthal vorticity ωφ
and electric current Jφ are related to Ψ and A by
ωφ = −
(
∇2 − r−2
)
Ψ (H.3)
and Jφ = −
(
∇2 − r−2
)
A . (H.4)
For numerical reasons described below, we introduce anisotropic viscosity and
helium diffusivity with dimensionless horizontal components νH and χH . So
161
162 H. The Numerical Scheme
equations (7.7)–(7.12) are replaced by
Ro 1r
D(ruφ)
Dt =
∂Ψ
∂z +
1
rB ·∇(rBφ)
+ Ek
[ ∂2
∂z2 + νH
(
∇2H − r−2
)
]
uφ (H.5)
Ro
[
rD(ωφ/r)Dt +
∂(ruφ, uφ/r)
∂(z, r)
]
= ∂uφ∂z − αT
∂T
∂r + αµ
∂µ
∂r
+ rB ·∇(Jφ/r) +
∂(rBφ, Bφ/r)
∂(z, r)
+ Ek
[ ∂2
∂z2 + νH
(
∇2H − r−2
)
]
ωφ (H.6)
rD(Bφ/r)Dt = rB ·∇(uφ/r) +
(
∇2 − r−2
)
Bφ (H.7)
1
r
D(rA)
Dt =
(
∇2 − r−2
)
A (H.8)
DT
Dt =
κ
η∇
2T (H.9)
Dµ
Dt =
χ
η
[ ∂2
∂z2 + χH∇
2
H
]
µ . (H.10)
Because of the axisymmetric cylindrical geometry, the lack of spatial periodic-
ity, and the wide range of spatial scales inherent in the polar confinement-layer
problem, a spectral or pseudospectral code would be unsuited to the task of solv-
ing these equations. Instead, a simple finite-difference code has been written in
cylindrical polar coordinates, with an Eulerian grid regularly spaced in r and z
at intervals ∆r and ∆z. The outer boundary of the computational domain is
at r = rd. The inner boundary is at r = 2∆r, i.e. two grid intervals from the
coordinate singularity at the rotation axis. For reasons of symmetry and good
behaviour near the coordinate singularity, the finite differencing is done by locally
approximating the fields Ψ/r, uφ/r, ωφ/r, A/r, Bφ/r, T , and µ as functions that
are linear in z and in r2 over a single grid interval. This ensures that the error is
O(∆r) even for small r. Field values for r < 2∆r are obtained by extrapolation
from r = 3∆r and r = 2∆r, again assuming linear functional dependence on r2.
With the parameter values given in table 10.1, the dimensionless helium-sublayer
and Ekman-layer thicknesses are δχ/δη = (χ/η)1/2 ≈ 0.14 and δEk/δη = Ek1/2 ≈
0.01 respectively. We have chosen a vertical grid interval ∆z = 0.01, dimen-
H.1 Equations 163
sionally 0.01δη, which is small enough to resolve the helium sublayer accurately.
This ∆z is too large to resolve any Ekman layers. However, Ekman layers are
prevented from becoming significant by careful choice of the code representing
the boundary conditions (see §10.3).
An explicit Eulerian timestepping scheme is used to evolve the system. The
timestep ∆t must be small enough to resolve thermal diffusion at the grid scale,
which is the fastest process at this scale and therefore determines the Courant–
Friedrichs–Lewy condition. With the parameter values given in table 10.1, this
constraint becomes ∆t . (η/κ)(∆z)2 = 10−2 × (0.01)2 = 10−6, dimensionally
10−6δη/U or 10−4(2Ωi)−1. The system typically takes several domain-scale mag-
netic diffusion times to reach a steady state, and multiple iterations of the pe-
ripheral Bφ(z) profile are required to achieve a steady state with vanishing Bφ(r)
at the bottom. To make the computation feasible, in a domain wide enough to
accommodate noticeable tilting effects, we have used rd = 5, dimensionally 5δη,
and a horizontal grid interval ∆r = 0.1, dimensionally 0.1δη, larger than the
vertical grid interval ∆z by a factor of 10. To ensure that the diffusive terms in
(H.5), (H.6) and (H.10) are dominant at the grid scale, we increased the horizon-
tal viscosity and helium diffusivity each by a factor of 10, i.e. νH = χH = 10. We
have verified, in smaller computational domains, that the coarser horizontal res-
olution and the anisotropic diffusion do not qualitatively affect the steady state
of the system (see §§I.4 and I.5). The increased horizontal viscosity introduces
stronger angular momentum coupling between neighbouring Ferraro surfaces, and
typically brings the system’s steady state closer to uniform rotation below the
confinement layer. In this sense the artificial viscosity acts rather like the hy-
pothesized shear-induced turbulence mentioned in §3.2.
At each timestep, the azimuthal vorticity ωφ is updated and the streamfunction Ψ
then computed from (H.3) by inverting the operator ∇2 − r−2, approximated us-
ing centred differences. The inversion is performed iteratively, using a successive-
overrelaxation method described in Press et al. (1986). During the early evo-
lution, when the dynamics is dominated by timescales not much longer than
the timestep ∆t (see §I.1), many such iterations are required, at each timestep,
to achieve convergence. At later times the same degree of convergence can be
achieved with far fewer iterations. Since we are interested only in the ultimate
164 H. The Numerical Scheme
steady state, we can tolerate a larger error in the inversion during the early
evolution. Further details of the numerical code are given in §H.2.
As anticipated from scaling arguments, the steady state is found to be close to
magnetostrophic balance (see §I.2). It might be thought that imposing magne-
tostrophic balance throughout the evolution, as first suggested by Taylor (1963),
would filter out all the fast oscillations, including inertial or epicyclic oscillations,
and thereby allow larger timesteps to be used. However, the imposition of magne-
tostrophic balance leads to pathological behaviour at small scales (Walker et al.,
1998). Far from eliminating or slowing the fast oscillations, the imposition of
balance exacerbates the problem, for reasons explained in §H.3.
H.2 Finite-differences
We consider grid-scale perturbations in a region of the numerical domain away
from the axis, so that curvature terms may be neglected. We linearise equa-
tions (H.5)–(H.9) about a background state with a uniform magnetic field B and
a uniform vertical temperature gradient dT/dz = 1. For simplicity we ignore
gradients in µ. (The terms in our equations involving µ will be computed in the
same fashion as the terms involving T .) We may also assume that the background
state is at rest, provided that the grid intervals ∆r and ∆z are small enough that
diffusion dominates advection at the grid scale.
Using primes to denote the perturbations, the linearised equations are
Ro ∂∂tu
′
φ −
∂
∂z
(1)
Ψ′ =B ·∇(1)B′φ + Ek
[ ∂2
∂z2 + νH∇
2
H
]
u′φ (H.11)
Ro ∂∂tω
′
φ −
∂
∂z
(2)
u′φ =− αT
∂
∂r
(1)
T ′ +B ·∇(2)J ′φ + Ek
[ ∂2
∂z2 + νH∇
2
H
]
ω′φ (H.12)
∂
∂tB
′
φ =B ·∇(3)u′φ +∇2B′φ (H.13)
∂
∂tA
′ =B ·∇(4)Ψ′ +∇2A′ (H.14)
∂
∂tT
′ + ∂∂r
(2)
Ψ′ =κη∇
2T ′ (H.15)
ω′φ =−∇2Ψ′ (H.16)
J ′φ =−∇2A′ (H.17)
H.2 Finite-differences 165
where the operators ∂/∂z, B · ∇, etc. represent finite difference operators. The
Laplacian-type operators are all evaluated using three-point, centred (first-order)
finite differences. Operators with superscripts in brackets are evaluated using
two-point, one-sided (first-order) finite differences. We would like to choose the
directions of these one-sided differences such that small scales in the system
behave as MAC waves (see appendix C).
To simplify the notation, we define new operators
Ro ∂∂t
(ν)
= Ro ∂∂t − Ek
[ ∂2
∂z2 + νH∇
2
H
]
, (H.18)
∂
∂t
(η)
= ∂∂t −∇
2 , (H.19)
and ∂∂t
(κ)
= ∂∂t −
κ
η∇
2 . (H.20)
The linear equations (H.11)–(H.17) can be combined into a single equation
M[u] = 0, (H.21)
where u represents any of the fields {Ψ′, u′φ, A′, B′φ, T ′}, and M is the operator
∂
∂t
(κ)
[
Ro ∂∂t
(ν) ∂
∂t
(η) −B ·∇(1)B ·∇(3)
] [
Ro ∂∂t
(ν) ∂
∂t
(η) −B ·∇(2)B ·∇(4)
]
∇2
+ ∂∂t
(κ) ∂
∂z
(1) ∂
∂z
(2) ∂
∂t
(η) ∂
∂t
(η)
+ ∂∂t
(η)
[
Ro ∂∂t
(ν) ∂
∂t
(η) −B ·∇(1)B ·∇(3)
]
αT ∂∂r
(1) ∂
∂r
(2) (H.22)
(cf. equation (C.21)). In the numerical scheme, the one-sided difference operators
∂
∂z
(1) and ∂∂z
(2) are evaluated in opposite directions (i.e. one is “up” and one is
“down”) so that
∂
∂z
(1) ∂
∂z
(2)
= ∂
2
∂z2 (H.23)
where ∂2∂z2 is a first-order-accurate centred difference. Similarly, ∂∂r
(1) is evaluated
in the opposite direction to ∂∂r
(2), B · ∇(1) is evaluated in the opposite direction
to B ·∇(3), and B ·∇(2) is evaluated in the opposite direction to B ·∇(4). These
choices ensure that the system’s small-scale behaviour mimics that of true MAC
waves as closely as possible. The remaining degrees of freedom in choosing the
166 H. The Numerical Scheme
Table H.1: The directions chosen for the one-sided finite differences. “Up-
ward”, for example, means that the derivative at each grid point is evaluated
using the difference between the value at that grid point and the value at
the grid point immediately above.
Finite Direction Finite Direction
difference chosen difference chosen
∂
∂z
(1) upward ∂∂z
(2) downward
∂
∂r
(1) outward ∂∂r
(2) inward
B ·∇(1) upward and inward B ·∇(3) downward and outward
B ·∇(2) downward and outward B ·∇(4) upward and inward
directions for the differences have been used to allow the most straightforward
implementation of the boundary conditions (§10.3). For example, since we specify
a peripheral boundary condition for T , it is convenient to define the difference
operator ∂∂r
(1) in (H.12) as an “outward” directed difference. The directionalities
chosen for all one-sided finite differences are listed in table H.1.
The code is written in C and currently runs on a single processor only. The source
code is available online1 or by request.2
H.3 Magnetostrophic balance
Because the confinement-layer problem involves only steady or nearly-steady flow,
it is natural to try to save computational resources by filtering out the fast os-
cillations (e.g. Taylor, 1963). These include Alfve´n waves, gravity waves, iner-
tia/Coriolis/epicyclic waves and the various hybrid types (see appendix C). Such
filtering is familiar, and often effective, in many other problems involving stiff
differential equations. A well known example is that of fluid flow in non-MHD
fluid systems with strong rotation and stable stratification. The standard “quasi-
geostrophic equations” result from small-Ro approximations that filter out inertia
and gravity waves, as well as sound waves, allowing relatively large time steps.
1www.damtp.cam.ac.uk/user/tsw25/cl-code/
2t.wood.02 "at" cantab.net
H.3 Magnetostrophic balance 167
Such filtering turns out, however, to be ineffective in the confinement-layer prob-
lem. Indeed — at first sight paradoxically — the imposition of magnetostrophic
balance leads to pathological behaviour in the following sense. Far from elimi-
nating fast oscillations, it introduces spurious modes of oscillation that are even
faster, as shown by Walker et al. (1998) in the context of the terrestrial dynamo
problem. Following Walker et al., we show how the pathology can be understood
through an idealised analysis of the fast oscillations, first in the unfiltered and
then in the filtered equations.
The reason for the pathology is the interplay between the Coriolis and Lorentz
forces. Stratification N2 is relatively unimportant, as will be shown shortly. We
therefore start with the linear theory of MC (magneto–Coriolis) waves, i.e. small
plane-wave disturbances to an unstratified, incompressible fluid with uniform
rotation Ω and a uniform magnetic field B. Neglecting viscosity and magnetic
diffusivity, we find the well-known dispersion relation
ω2 − (B · k)2 = ±2Ω · kω/|k| , (H.24)
where ω is the frequency and k is the wavevector, both dimensional here. If we
take the limit of rapid rotation, |Ω| → ∞, then for most choices of k the four
roots of this dispersion relation are asymptotically
ω ∼ ±2Ω · k|k| , (H.25)
and ω ∼ ±(B · k)
2
2Ω · k |k| . (H.26)
The modes corresponding to (H.25) are inertial waves — in this context some-
times called “fast MC waves” — and those corresponding to (H.26) are “slow
MC waves”. By imposing magnetostrophic balance we neglect relative fluid ac-
celerations, which corresponds to dropping the ω2 term from the left-hand side
of (H.24). The dispersion relation then becomes
ω = ±(B · k)
2
2Ω · k |k| , (H.27)
so imposing magnetostrophic balance eliminates the two “fast” branches (H.25)
of the full dispersion relation (H.24).
168 H. The Numerical Scheme
However, not all modes of the full dispersion relation (H.24) have the asymptotic
behaviour described by (H.25) and (H.26). Even in the presence of rapid rotation,
there are always some modes whose k values satisfy
|B · k|  |2Ω · k|/|k| as |Ω| → ∞ , (H.28)
and such modes behave like Alfve´n waves, i.e. ω ∼ ±B · k. Imposing mag-
netostrophic balance removes the mechanism for Alfve´n wave propagation, and
must therefore alter the behaviour of these modes. From (H.27) we find that
these modes now have a very high frequency. In fact if we fix 2Ω · k and allow
|k| → ∞ then ω grows as |k|3. The conclusion is that even in a rapidly rotating
system some modes of the full dispersion relation always feel the Lorentz force
more strongly than the Coriolis force, and these modes become ill-behaved under
the assumption of magnetostrophic balance. A numerical scheme that imposes
magnetostrophic balance will therefore be ill-conditioned.
If we introduce stratification N2 then (H.27) becomes
ω = ± B · k2Ω · k
[
(B · k)2|k|2 +N2(|k|2 − k2V )
]1/2 , (H.29)
where kV is the vertical component of k. Therefore the presence of stratification
serves only to increase the frequency of the ill-behaved modes, and thereby to
exacerbate the problem.
Appendix I
Numerical Test Cases
I.1 Transient adjustment
The numerical scheme timesteps equations (H.5)–(H.10) forwards in time from a
given initial condition. We typically take the initial condition to be a previously
computed numerical or analytical confinement-layer solution. The behaviour at
early times is characterised by “sloshing” of the flow through the confinement
layer, leading to undulations of the helium sublayer on timescales of the order
of the inverse buoyancy frequency. Figure I.1 shows the meridional streamlines
in a vertical cross-section through the computational domain during this tran-
sient adjustment phase. An animated version of this figure is available online at
www.damtp.cam.ac.uk/user/tsw25/adjust.avi. The animation shows the evo-
lution of the poloidal streamlines over a dimensionless time interval of 5× 10−2,
dimensionally 50/N , where N is the total buoyancy frequency at the bottom of
the domain. The meridional flow typically settles into a quasi-steady state after a
dimensionless time much less that unity. The magnetic field, on the other hand,
typically takes a dimensionless time of order unity to achieve a steady state.
169
170 I. Numerical Test Cases
Figure I.1: Meridional streamlines in a vertical cross-section through the confinement layer
during its transient adjustment towards a steady state.
I.2 Magnetostrophic balance
The evolution of angular momentum in the numerical scheme is described by
equation (H.5). Since both Ro and Ek are much smaller than unity, we anticipate
that the solutions will be close to magnetostrophic balance, with
0 ≈ 1r
∂Ψ
∂z +
1
r2B ·∇(rBφ) . (I.1)
In figure I.2 we plot the first and second terms on the right-hand side of (I.1) in a
vertical cross-section through the same numerical solution1 presented in figure 7.2.
In figure I.3 we plot the sum of these two terms, which would vanish if the solution
were in perfect magnetostrophic balance. We see that both terms in (I.1) are
largest within the confinement layer, and both are close to zero above and below
the confinement layer. As expected, their sum is close to zero throughout the
domain, indicating that the solution is very nearly magnetostrophically balanced.
The largest absolute departures from magnetostrophic balance are localised next
to the vertical outer boundary, where we impose profiles of T and µ taken from
a perfectly flat, analytical confinement-layer solution (§7.5). The largest relative
departures from magnetostrophic balance occur above and below the confinement
layer, where both terms in (I.1) are small. Below the confinement layer, the most
significant departures from magnetostrophic balance are associated with viscous
1Figure I.2 shows the entire numerical domain, whereas the cross-section in figure 7.2 was
truncated at z = 5.
I.3 Numerical boundaries 171
angular-momentum exchange between the poloidal magnetic flux surfaces.
 
 
−4 −2 0 2 4−1
0
1
2
3
4
5
−1.6
−1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
0
 
 
−4 −2 0 2 4−1
0
1
2
3
4
5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Figure I.2: Contours of the first and second terms in (I.1), in a vertical cross-section through
a numerical solution. Both terms are largest in magnitude within the confinement layer.
 
 
−4 −2 0 2 4−1
0
1
2
3
4
5
−0.05
−0.04
−0.03
−0.02
−0.01
0
0.01
Figure I.3: Contours of the total right-hand side of (I.1), representing departures from mag-
netostrophic balance.
I.3 Numerical boundaries
I.3.1 The upper boundary
In the numerical solution presented in figure 7.2, the boundary condition
∂Bz/∂z = −Bz (I.2)
172 I. Numerical Test Cases
was imposed at the top of the domain (see §10.3) in order to mimic the presence of
an exponential tail in the layers above. Contours of Bz in a vertical cross-section
through that solution are presented in figure I.4. To verify that the solution is not
sensitive to this particular choice of boundary condition, we have also computed
a solution with the condition Bz = 0, instead of (I.2), imposed at the top of the
domain, but with all other boundary conditions unchanged. Figure I.4 presents
contours of the difference in Bz between the two solutions. As expected, the effect
of the change in the upper boundary condition is localised to the region close to
the upper boundary, and is small in magnitude.
 
 
−5 0 5−1
0
1
2
3
4
5
0.5
1
1.5
2
2.5
3
3.5
 
 
−5 0 5−1
0
1
2
3
4
5
−4
−3.5
−3
−2.5
−2
−1.5
−1
−0.5
0x 10
−3
Figure I.4: The effect of changing the upper boundary condition on Bz . The left panel shows
contours of Bz obtained when matching on to an exponentially decaying profile above (see
§10.3). The right panel shows the change in Bz when, instead, a vanishing field is imposed
above.
I.3.2 Boundary locations
Since the boundaries of the numerical domain are artificial, we should check that
their exact locations do not significantly affect the solutions. In figure I.5 we
show the result of shifting the domain vertically. The left-hand plot shows the
streamlines, magnetic field lines and compositional stratification in a vertical
cross-section through the same solution presented in figure 7.2. The right-hand
plot shows an equivalent cross-section through a solution in which the numerical
domain was shifted downward by a dimensionless distance of 0.5, i.e. 0.5δη, via a
suitable translation of the peripheral boundary conditions. The solution is found
to be entirely robust to the vertical shift of the numerical domain.
I.4 Horizontal resolution 173
−5 0 5
−1
0
1
2
3
4
5
−5 0 5
−1
0
1
2
3
4
5
Figure I.5: Shifting the numerical domain downward does not significantly affect the solution.
The velocity streamlines and magnetic field lines are shown in blue and red respectively. The
shading indicates compositional stratification, ∂µ/∂z.
I.4 Horizontal resolution
As mentioned in appendix H, the constraints on the vertical grid interval ∆z
are stronger than the constraints on the horizontal grid interval ∆r, due to the
thinness of the helium sublayer. The numerical solution presented in figure 7.2
has ∆z = 0.01 and ∆r = 0.1. To check that the relatively coarse horizontal reso-
lution does not qualitatively affect the numerical results, we have also computed
a solution with ∆r = 0.05, in a domain with cylindrical radius rd = 2.5. The pe-
ripheral boundary profiles T (z), µ(z) and Bφ(z) in this solution were taken from
the original, wider solution at r = 2.5. In figure I.6 we plot contours of Bφ/r
from both solutions. Although the two solutions are not identical, the differences
between the two are small.
 
 
−5 0 5−1
0
1
2
3
4
5
0
0.05
0.1
0.15
0.2
 
 
−5 0 5−1
0
1
2
3
4
5
0
0.05
0.1
0.15
0.2
Figure I.6: The influence of horizontal resolution on the numerical solutions. The left panel
shows contours of Bφ/r in a solution with ∆r = 0.1 and rd = 5. The right panel shows contours
of Bφ/r in a solution with ∆r = 0.05 and rd = 2.5.
174 I. Numerical Test Cases
I.5 Horizontal diffusion
The increased horizontal diffusivities νH and χH , introduced for numerical sta-
bility, undoubtedly have an effect on the numerical solutions. We would like to
quantify this effect, in order to ensure that the solutions do not depend signifi-
cantly on the values of νH and χH chosen. In figure I.7 we present the rotation
profile uφ/r in a vertical cross-section through the numerical solution described
in §I.4, which has ∆r = 0.05, rd = 2.5 and νH = χH = 10. We also present the
rotation profile in a solution that has νH = χH = 5 but is otherwise identical.
The most noticeable differences between the two solutions are in the upper part of
the domain, where magnetostrophic balance becomes delicate (§I.2) and diffusion
of angular momentum therefore becomes increasingly significant. Nonetheless,
changing the horizontal diffusivities does not qualitatively affect the structure of
the confinement layer.
 
 
−5 0 5−1
0
1
2
3
4
5
0
0.5
1
1.5
2
 
 
−5 0 5−1
0
1
2
3
4
5
0
0.5
1
1.5
2
Figure I.7: The rotation profile uφ/r in two numerical solutions. The left panel has νH =
χH = 10; the right panel has νH = χH = 5. The two solutions have identical boundary
conditions.
Bibliography
Alfve´n, H. 1942 Existence of electromagnetic-hydrodynamic waves. Nature
150, 405–406.
Andrews, D. G., Holton, J. R. & Leovy, C. B. 1987 Middle atmosphere
dynamics . New York: Academic Press.
Antia, H. M. & Basu, S. 2000 Temporal variations of the rotation rate in the
solar interior. Astrophys. J. 541, 442–448.
Appourchaux, T., Belkacem, K., Broomhall, A.-M., Chaplin, W. J.,
Gough, D. O., Houdek, G., Provost, J., Baudin, F., Boumier, P.,
Elsworth, Y., Garc´ıa, R. A., Andersen, B. N., Finsterle, W.,
Fro¨hlich, C., Gabriel, A., Grec, G., Jime´nez, A., Kosovichev, A.,
Sekii, T., Toutain, T. & Turck-Chie`ze, S. 2010 The quest for the solar
g modes. A&A Rev. 18, 197–277.
Asplund, M. 2005 New light on stellar abundance analyses: departures from
LTE and homogeneity. Annu. Rev. Astron. Astrophys. 43, 481–530.
Asplund, M., Grevesse, N. & Sauval, A. J. 2006 The new solar abundances
— Part I: the observations. Communications in Asteroseismology 147, 76–79.
Asplund, M., Grevesse, N., Sauval, A. J. & Scott, P. 2009 The chemical
composition of the Sun. Annu. Rev. Astron. Astrophys. 47, 481–522.
Aurnou, J., Heimpel, M. & Wicht, J. 2007 The effects of vigorous mixing
in a convective model of zonal flow on the ice giants. Icarus 190, 110–126.
Babcock, H. W. 1961 The topology of the Sun’s magnetic field and the 22-year
cycle. Astrophys. J. 133, 572–587.
175
176 Bibliography
Badnell, N. R., Bautista, M. A., Butler, K., Delahaye, F., Mendoza,
C., Palmeri, P., Zeippen, C. J. & Seaton, M. J. 2005 Updated opacities
from the Opacity Project. MNRAS 360, 458–464.
Bahcall, J. N., Basu, S. & Serenelli, A. M. 2005 What is the neon
abundance of the Sun? Astrophys. J. 631, 1281–1285.
Bahcall, J. N. & Pinsonneault, M. H. 2004 What do we (not) know the-
oretically about solar neutrino fluxes? Phys. Rev. Lett. 92 (12), 121301.
Bahcall, J. N. & Serenelli, A. M. 2005 How do uncertainties in the surface
chemical composition of the Sun affect the predicted solar neutrino fluxes?
Astrophys. J. 626, 530–542.
Barcilon, V. & Pedlosky, J. 1967 A unified linear theory of homogeneous
and stratified rotating fluids. J. Fluid Mech. 29, 609–621.
Basu, S. & Antia, H. M. 1995 Helium abundance in the solar envelope. MN-
RAS 276, 1402–1408.
Basu, S. & Antia, H. M. 2003 Changes in solar dynamics from 1995 to 2002.
Astrophys. J. 585, 553–565.
Basu, S. & Antia, H. M. 2004 Constraining solar abundances using helioseis-
mology. Astrophys. J. 606, L85–L88.
Boesgaard, A. M. & Tripicco, M. J. 1986 Lithium in the Hyades Cluster.
Astrophys. J. 302, L49–L53.
Bo¨hm-Vitense, E. 1958 U¨ber die Wasserstoffkonvektionszone in Sternen ver-
schiedener Effektivtemperaturen und Leuchtkra¨fte. Mit 5 Textabbildungen.
Zeitschrift fur Astrophysik 46, 108–143.
Bouvier, J., Alencar, S. H. P., Harries, T. J., Johns-Krull, C. M.
& Romanova, M. M. 2007 Magnetospheric accretion in classical T Tauri
stars. In Protostars and Planets V (ed. B. Reipurth, D. Jewitt, & K. Keil),
pp. 479–494. Tuscon: Univ. Arizona Press.
Braithwaite, J. & Nordlund, A˚. 2006 Stable magnetic fields in stellar in-
teriors. A&A 450, 1077–1095.
BIBLIOGRAPHY 177
Braithwaite, J. & Spruit, H. C. 2004 A fossil origin for the magnetic field
in A stars and white dwarfs. Nature 431, 819–821.
Bretherton, F. P. & Spiegel, A. E. 1968 The effect of the convection zone
on solar spin-down. Astrophys. J. 153, L77–L80.
Brummell, N. H., Clune, T. L. & Toomre, J. 2002 Penetration and over-
shooting in turbulent compressible convection. Astrophys. J. 570, 825–854.
Brun, A. S. & Zahn, J.-P. 2006 Magnetic confinement of the solar tachocline.
A&A 457, 665–674.
Busse, F. H. & Hood, L. L. 1982 Differential rotation driven by convection
in a rapidly rotating annulus. Geophys. Astrophys. Fluid Dyn. 21, 59–74.
Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A. & Hurl-
burt, N. E. 1991 Turbulent compressible convection. Astrophys. J. 370, 282–
294.
Cattaneo, F. & Hughes, D. W. 1996 Nonlinear saturation of the turbulent
α effect. Phys. Rev. E 54, 4532–4535.
Chandrasekhar, S. & Prendergast, K. H. 1956 The equilibrium of mag-
netic stars. Proceedings of the National Academy of Science 42, 5–9.
Charbonneau, P., Beaubien, G. & St-Jean, C. 2007 Fluctuations in
Babcock–Leighton dynamos. II. Revisiting the Gnevyshev–Ohl rule. Astrophys.
J. 658, 657–662.
Charbonneau, P. & MacGregor, K. B. 1993 Angular momentum transport
in magnetized stellar radiative zones. II. The solar spin-down. Astrophys. J.
417, 762.
Charbonneau, P. & MacGregor, K. B. 1997 Solar interface dynamos. II.
Linear, kinematic models in spherical geometry. Astrophys. J. 486, 502–520.
Charbonnel, C. & Talon, S. 2005 Influence of gravity waves on the internal
rotation and Li abundance of solar-type stars. Science 309, 2189–2191.
178 Bibliography
Charbonnel, C. & Talon, S. 2007 On internal gravity waves in low-mass
stars. In Unsolved Problems in Stellar Physics: A Conference in Honor of
Douglas Gough (ed. R. J. Stancliffe, J. Dewi, G. Houdek, R. G. Martin, &
C. A. Tout), American Institute of Physics Conference Series, vol. 948, pp.
15–26.
Charbonnel, C., Vauclair, S. & Zahn, J.-P. 1992 Rotation-induced mixing
and lithium depletion in galactic clusters. A&A 255, 191–199.
Chen, Y. Q., Nissen, P. E., Benoni, T. & Zhao, G. 2001 Lithium abun-
dances for 185 main-sequence stars: Galactic evolution and stellar depletion of
lithium. A&A 371, 943–951.
Christensen-Dalsgaard, J. 2002 Helioseismology. Reviews of Modern
Physics 74, 1073–1129.
Christensen-Dalsgaard, J. & di Mauro, M. P. 2007 Diffusion and he-
lioseismology. In Stellar Evolution and Seismic Tools for Asteroseismology —
Diffusive Processes in Stars and Seismic Analysis (ed. C. W. Straka, Y. Le-
breton, & M. J. P. F. G. Monteiro), EAS Publications Series , vol. 26, pp.
3–16.
Christensen-Dalsgaard, J., di Mauro, M. P., Houdek, G. & Pijpers,
F. 2009 On the opacity change required to compensate for the revised solar
composition. A&A 494, 205–208.
Christensen-Dalsgaard, J., Gough, D. O. & Thompson, M. J. 1991
The depth of the solar convection zone. Astrophys. J. 378, 413–437.
Christensen-Dalsgaard, J., Gough, D. O. & Thompson, M. J. 1992 On
the rate of destruction of lithium in late-type main-sequence stars. A&A 264,
518–528.
Christensen-Dalsgaard, J., Proffitt, C. R. & Thompson, M. J. 1993
Effects of diffusion on solar models and their oscillation frequencies. Astrophys.
J. 403, L75–L78.
BIBLIOGRAPHY 179
Christensen-Dalsgaard, J. & Thompson, M. J. 2007 Observational re-
sults and issues concerning the tachocline. In The Solar Tachocline (ed.
D. W. Hughes, R. Rosner, & N. O. Weiss), p. 53.
Christian, D. J. & Mathioudakis, M. 2000 Is Lithium produced in active
stellar atmospheres? In NOAO Proposal ID #2000B-0165 , p. 165.
Ciacio, F., degl’Innocenti, S. & Ricci, B. 1997 Updating standard solar
models. A&AS 123, 449–454.
Clark, A. 1973 The linear spin-up of a strongly stratified fluid of small Prandtl
number. J. Fluid Mech. 60, 561–580.
Clark, A. 1975 Time scales in stellar spin-down. Me´moires Socie´te´ Royale des
Sciences de Lie`ge 8, 43–46.
Claverie, A., Isaak, G. R., McLeod, C. P., van der Raay, H. B. &
Roca Cortes, T. 1981 Rapid rotation of the solar interior. Nature 293,
443–445.
Cowling, T. G. 1933 The magnetic field of sunspots. MNRAS 94, 39–48.
Craik, A. D. D. 1988 A class of exact solutions in viscous incompressible
magnetohydrodynamics. Proc. R. Soc. Lond. A 417, 235–244.
de Jager, C. 1972 Micro- and macroturbulent motions and the velocity spec-
trum of the solar photosphere. Sol. Phys. 25, 71–80.
Debnath, L. 1973 On Ekman and Hartmann boundary layers in a rotating fluid.
Acta Mechanica 18, 333–341.
Denissenkov, P. A. & Pinsonneault, M. 2007 A revised prescription for
the Tayler–Spruit dynamo: magnetic angular momentum transport in stars.
Astrophys. J. 655, 1157–1165.
Dikpati, M. & Charbonneau, P. 1999 A Babcock–Leighton flux transport
dynamo with solar-like differential rotation. Astrophys. J. 518, 508–520.
Dintrans, B., Brandenburg, A., Nordlund, A˚. & Stein, R. F. 2005
Spectrum and amplitudes of internal gravity waves excited by penetrative con-
vection in solar-type stars. A&A 438, 365–376.
180 Bibliography
Dorch, S. B. F. & Nordlund, A˚. 2001 On the transport of magnetic fields
by solar-like stratified convection. A&A 365, 562–570.
Dritschel, D. G. 1991 Generalized helical Beltrami flows in hydrodynamics
and magnetohydrodynamics. J. Fluid Mech. 222, 525–541.
Drobyshevski, E. M. & Yuferev, V. S. 1974 Topological pumping of mag-
netic flux by three-dimensional convection. J. Fluid Mech. 65, 33–44.
Elliott, J. R. 1997 Aspects of the solar tachocline. A&A 327, 1222–1229.
Elliott, J. R. & Gough, D. O. 1999 Calibration of the thickness of the solar
tachocline. Astrophys. J. 516, 475–481.
Feast, M. W. 1966 Lithium in β Hydri and other stars. MNRAS 134, 321–326.
Ferraro, V. C. A. 1937 The non-uniform rotation of the Sun and its magnetic
field. MNRAS 97, 458–472.
Fritts, D. C., Vadas, S. L. & Andreassen, O. 1998 Gravity wave excita-
tion and momentum transport in the solar interior: implications for a residual
circulation and lithium depletion. A&A 333, 343–361.
Gaidos, E. J. 1998 Nearby young solar analogs. I. Catalog and stellar charac-
teristics. PASP 110, 1259–1276.
Garaud, P. 2001 Latitudinal shear instability in the solar tachocline. MNRAS
324, 68–76.
Garaud, P. 2002 Dynamics of the solar tachocline — I. An incompressible
study. MNRAS 329, 1–17.
Garaud, P. & Acevedo Arreguin, L. 2009 On the penetration of meridional
circulation below the solar convection zone. II. Models with convection zone,
the Taylor–Proudman constraint, and applications to other stars. Astrophys.
J. 704, 1–16.
Garaud, P. & Brummell, N. H. 2008 On the penetration of meridional
circulation below the solar convection zone. Astrophys. J. 674, 498–510.
BIBLIOGRAPHY 181
Garaud, P. & Garaud, J.-D. 2008 Dynamics of the solar tachocline — II.
The stratified case. MNRAS 391, 1239–1258.
Garaud, P. & Rogers, T. 2007 Solar rotation. In Unsolved Problems in Stel-
lar Physics: A Conference in Honor of Douglas Gough (ed. R. J. Stancliffe,
J. Dewi, G. Houdek, R. G. Martin, & C. A. Tout), American Institute of
Physics Conference Series, vol. 948, pp. 237–248.
Garc´ıa Lo´pez, R. J. & Spruit, H. C. 1991 Li depletion in F stars by internal
gravity waves. Astrophys. J. 377, 268–277.
Gilman, P. A. 1977 Nonlinear dynamics of Boussinesq convection in a deep
rotating spherical shell. I. Geophys. Astrophys. Fluid Dyn. 8, 93–135.
Gilman, P. A. & Cally, P. S. 2007 Global MHD instabilities of the tachocline.
In The Solar Tachocline (ed. D. W. Hughes, R. Rosner, & N. O. Weiss), p.
243.
Gilman, P. A. & Miesch, M. S. 2004 Limits to penetration of meridional
circulation below the solar convection zone. Astrophys. J. 611, 568–574.
Gilman, P. A. & Miller, J. 1981 Dynamically consistent nonlinear dynamos
driven by convection in a rotating spherical shell. ApJS 46, 211–238.
Glatzmaier, G., Evonuk, M. & Rogers, T. 2009 Differential rotation in
giant planets maintained by density-stratified turbulent convection. Geophys.
Astrophys. Fluid Dyn. 103, 31–51.
Gough, D. 2007 An introduction to the solar tachocline. In The Solar Tachocline
(ed. D. W. Hughes, R. Rosner, & N. O. Weiss), p. 3.
Gough, D. O. 1969 The anelastic approximation for thermal convection. J. At-
mos. Sci. 26, 448–456.
Gough, D. O. & McIntyre, M. E. 1998 Inevitability of a magnetic field in
the Sun’s radiative interior. Nature 394, 755–757.
Greenspan, H. P. 1968 The theory of rotating fluids. Cambridge: University
Press.
182 Bibliography
Grevesse, N., Asplund, M. & Sauval, A. J. 2007 The solar chemical com-
position. Space Sci. Rev. 130, 105–114.
Grevesse, N. & Sauval, A. J. 1998 Standard solar composition.
Space Sci. Rev. 85, 161–174.
Guzik, J. A., Willson, L. A. & Brunish, W. M. 1987 A comparison between
mass-losing and standard solar models. Astrophys. J. 319, 957–965.
Hale, G. E. 1908 On the probable existence of a magnetic field in sun-spots.
Astrophys. J. 28, 315–343.
Hale, G. E., Ellerman, F., Nicholson, S. B. & Joy, A. H. 1919 The
magnetic polarity of sun-spots. Astrophys. J. 49, 153–178.
Haynes, P. H., McIntyre, M. E., Shepherd, T. G., Marks, C. J. &
Shine, K. P. 1991 On the ‘downward control’ of extratropical diabatic circu-
lations by eddy-induced mean zonal forces. J. Atmos. Sci. 48, 651–680.
Holton, J. R. 1965 The influence of viscous boundary layers on transient mo-
tions in a stratified rotating fluid. Part I. J. Atmos. Sci. 22, 402–411.
Holweger, H. & Mu¨ller, E. A. 1974 The photospheric barium spectrum —
Solar abundance and collision broadening of BA II lines by hydrogen. Sol. Phys.
39, 19–30.
Howard, R. & LaBonte, B. J. 1980 The Sun is observed to be a torsional
oscillator with a period of 11 years. Astrophys. J. 239, L33–L36.
Howe, R. 2009 Solar interior rotation and its variation. Living Reviews in Solar
Physics 6, 1.
Howe, R., Christensen-Dalsgaard, J., Hill, F., Komm, R., Schou, J.,
Thompson, M. J. & Toomre, J. 2007 Temporal variations in solar rotation
at the bottom of the convection zone: The current status. Advances in Space
Research 40, 915–918.
Howe, R., Christensen-Dalsgaard, J., Hill, F., Komm, R. W.,
Larsen, R. M., Schou, J., Thompson, M. J. & Toomre, J. 2000 Dy-
BIBLIOGRAPHY 183
namic variations at the base of the solar convection zone. Science 287, 2456–
2460.
Hurlburt, N. E., Toomre, J. & Massaguer, J. M. 1986 Nonlinear com-
pressible convection penetrating into stable layers and producing internal grav-
ity waves. Astrophys. J. 311, 563–577.
Iben, Jr., I. 1965 Stellar evolution. I. The approach to the main sequence.
Astrophys. J. 141, 993–1018.
Israelian, G., Delgado Mena, E., Santos, N. C., Sousa, S. G., Mayor,
M., Udry, S., Dom´ınguez Cerden˜a, C., Rebolo, R. & Randich, S.
2009 Enhanced lithium depletion in Sun-like stars with orbiting planets. Nature
462, 189–191.
Ka¨hler, H. 1978 The Vogt-Russell theorem, and new results on an old problem.
In The HR Diagram — The 100th Anniversary of Henry Norris Russell (ed.
A. G. D. Philip & D. S. Hayes), IAU Symposium, vol. 80, pp. 303–311.
King, E. M., Stellmach, S., Noir, J., Hansen, U. & Aurnou, J. M. 2009
Boundary layer control of rotating convection systems. Nature 457, 301–304.
Kitchatinov, L. L. & Ru¨diger, G. 1992 Magnetic-field advection in inho-
mogeneous turbulence. A&A 260, 494–498.
Kitchatinov, L. L. & Ru¨diger, G. 1995 Differential rotation in solar-type
stars: revisiting the Taylor-number puzzle. A&A 299, 446–452.
Kitchatinov, L. L. & Ru¨diger, G. 2006 Magnetic field confinement by merid-
ional flow and the solar tachocline. A&A 453, 329–333.
Kitchatinov, L. L. & Ru¨diger, G. 2008 Diamagnetic pumping near the base
of a stellar convection zone. Astron. Nachr. 329, 372–375.
Krause, F. & Ra¨dler, K.-H. 1980 Mean-field magnetohydrodynamics and
dynamo theory . Oxford: Pergamon Press.
Larmor, J. 1919 How could a rotating body such as the Sun become a magnet?
Rep. Brit. Assoc. Adv. Sci. p. 159.
184 Bibliography
Leighton, R. B. 1969 A magneto-kinematic model of the solar cycle. Astrophys.
J. 156, 1–26.
Levy, E. H. & Boyer, D. 1982 Oscillating dynamo in the presence of a fossil
magnetic field — The solar cycle. Astrophys. J. 254, L19–L22.
Lineykin, P. C. 1955 On the determination of the thickness of the baroclinic
layer in the ocean. Dokl. Akad. Nauk. SSSR 101, 461–464.
MacGregor, K. B. & Charbonneau, P. 1999 Angular momentum trans-
port in magnetized stellar radiative zones. IV. Ferraro’s theorem and the solar
tachocline. Astrophys. J. 519, 911–917.
Malkus, W. V. R. 1967 Hydromagnetic planetary waves. J. Fluid Mech. 28,
793–802.
Mathis, S. 2009 Transport by gravito-inertial waves in differentially rotating
stellar radiation zones. I — Theoretical formulation. A&A 506, 811–828.
McIntyre, M. 1994 The quasi-biennial oscillation (QBO) : some points about
the terrestrial QBO and the possibility of related phenomena in the solar in-
terior. In The Solar Engine and its Influence on Terrestrial Atmosphere and
Climate (ed. E. Nesme-Ribes), p. 293.
McIntyre, M. E. 2003 Solar tachocline dynamics: eddy viscosity, anti-friction,
or something in between? In Stellar astrophysical fluid dynamics (ed. Thomp-
son, M. J. & Christensen-Dalsgaard, J.), pp. 111–130. Cambridge: University
Press.
McIntyre, M. E. 2007 Magnetic confinement and the sharp tachopause. In
The Solar Tachocline (ed. D. W. Hughes, R. Rosner, & N. O. Weiss), p. 183.
Mestel, L. 1953 Rotation and stellar evolution. MNRAS 113, 716–745.
Mestel, L. & Moss, D. L. 1986 On mixing by the Eddington–Sweet circula-
tion. MNRAS 221, 25–51.
Mestel, L. & Weiss, N. O. 1987 Magnetic fields and non-uniform rotation in
stellar radiative zones. MNRAS 226, 123–135.
BIBLIOGRAPHY 185
Michaud, G. 1986 The lithium abundance gap in the Hyades F stars — The
signature of diffusion. Astrophys. J. 302, 650–655.
Michaud, G. & Charbonneau, P. 1991 The lithium abundance in stars.
Space Sci. Rev. 57, 1–58.
Michaud, G. & Proffitt, C. R. 1993 Particle transport processes. In IAU
Colloq. 137: Inside the Stars (ed. W. W. Weiss & A. Baglin), Astronomical
Society of the Pacific Conference Series, vol. 40, pp. 246–259.
Miesch, M. S. 2005 Large-scale dynamics of the convection zone and tachocline.
Living Reviews in Solar Physics 2, 1.
Miesch, M. S., Brun, A. S. & Toomre, J. 2006 Solar differential rotation
influenced by latitudinal entropy variations in the tachocline. Astrophys. J.
641, 618–625.
Miesch, M. S., Elliott, J. R., Toomre, J., Clune, T. L., Glatzmaier,
G. A. & Gilman, P. A. 2000 Three-dimensional spherical simulations of
solar convection. I. Differential rotation and pattern evolution achieved with
laminar and turbulent states. Astrophys. J. 532, 593–615.
Moffatt, H. K. 1978 Magnetic field generation in electrically conducting fluids.
Cambridge: University Press.
Montalba´n, J., Miglio, A., Noels, A., Grevesse, N. & di Mauro,
M. P. 2004 Solar model with CNO revised abundances. In SOHO 14 Helio-
and Asteroseismology: Towards a Golden Future (ed. D. Danesy), ESA Special
Publication, vol. 559, pp. 574–576.
Montalba´n, J., Miglio, A., Theado, S., Noels, A. & Grevesse, N.
2006 The new solar abundances — Part II: the crisis and possible solutions.
Communications in Asteroseismology 147, 80–84.
Mursula, K., Usoskin, I. G. & Kovaltsov, G. A. 2001 Persistent 22-year
cycle in sunspot activity: Evidence for a relic solar magnetic field. Sol. Phys.
198, 51–56.
186 Bibliography
Ogilvie, G. I. 2007 Instabilities, angular momentum transport and magnetohy-
drodynamic turbulence. In The Solar Tachocline (ed. D. W. Hughes, R. Rosner,
& N. O. Weiss), p. 299.
Parfrey, K. P. & Menou, K. 2007 The origin of solar activity in the
tachocline. Astrophys. J. 667, L207–L210.
Parker, E. N. 1955a Hydromagnetic dynamo models. Astrophys. J. 122, 293–
314.
Parker, E. N. 1955b The formation of sunspots from the solar toroidal field.
Astrophys. J. 121, 491–507.
Parker, E. N. 1975 The generation of magnetic fields in astrophysical bodies.
X — Magnetic buoyancy and the solar dynamo. Astrophys. J. 198, 205–209.
Parker, E. N. 1981 Residual fields from extinct dynamos. Geophys. Astrophys.
Fluid Dyn. 18, 175–196.
Parker, E. N. 1993 A solar dynamo surface wave at the interface between
convection and nonuniform rotation. Astrophys. J. 408, 707–719.
Pedlosky, J. 1979 Geophysical fluid dynamics. New York: Springer-Verlag.
Pinsonneault, M. H. & Delahaye, F. 2009 The solar heavy element abun-
dances. II. Constraints from stellar atmospheres. Astrophys. J. 704, 1174–1188.
Pinsonneault, M. H., DePoy, D. L. & Coffee, M. 2001 The mass of
the convective zone in FGK main-sequence stars and the effect of accreted
planetary material on apparent metallicity determinations. Astrophys. J. 556,
L59–L62.
Pitts, E. & Tayler, R. J. 1985 The adiabatic stability of stars containing
magnetic fields. IV — The influence of rotation. MNRAS 216, 139–154.
Plumpton, C. & Ferraro, V. C. A. 1955 On toroidal magnetic fields in the
Sun and stars. Astrophys. J. 121, 168–174.
Prendergast, K. H. 1956 The equilibrium of a self-gravitating incompressible
fluid sphere with a magnetic field. I. Astrophys. J. 123, 498–507.
BIBLIOGRAPHY 187
Press, W. H., Flannery, B. P., Teukolsky, S. A. & Vetterling, W. T.
1986 Numerical recipes. The art of scientific computing . Cambridge: University
Press.
Proffitt, C. R. & Michaud, G. 1991 Diffusion and mixing of lithium and
helium in population II dwarfs. Astrophys. J. 371, 584–601.
Ra¨dler, K. H. 1968 On the electrodynamics of conducting fluids in turbu-
lent motion. II. Turbulent conductivity and turbulent permeability. Zeitschrift
Naturforschung Teil A 23, 1851–1860.
Rempel, M. 2005 Solar differential rotation and meridional flow: The role of a
subadiabatic tachocline for the Taylor–Proudman balance. Astrophys. J. 622,
1320–1332.
Roberts, P. H. & Soward, A. M. 1975 A unified approach to mean field
electrodynamics. Astron. Nachr. 296, 49–64.
Rogers, T. M. & Glatzmaier, G. A. 2006 Angular momentum transport by
gravity waves in the solar interior. Astrophys. J. 653, 756–764.
Rogers, T. M., MacGregor, K. B. & Glatzmaier, G. A. 2008 Non-linear
dynamics of gravity wave driven flows in the solar radiative interior. MNRAS
387, 616–630.
Rossby, C. G. 1938 On the mutual adjustment of pressure and velocity distri-
butions in certain simple current systems, II. J. Mar. Res. 1, 239–263.
Ru¨diger, G. & Kitchatinov, L. L. 1997 The slender solar tachocline: a
magnetic model. Astron. Nachr. 318, 273–279.
Ru¨diger, G., Kitchatinov, L. L. & Arlt, R. 2005 The penetration of
meridional flow into the tachocline and its meaning for the solar dynamo. A&A
444, L53–L56.
Schatzman, E. 1962 A theory of the role of magnetic activity during star for-
mation. Annales d’Astrophysique 25, 18–29.
Schatzman, E. 1993 Transport of angular momentum and diffusion by the
action of internal waves. A&A 279, 431–446.
188 Bibliography
Schmitt, D. 1993 The solar dynamo. In The Cosmic Dynamo (ed. F. Krause,
K. H. Ra¨dler, & G. Ru¨diger), IAU Symposium, vol. 157, pp. 1–12.
Schou, J., Antia, H. M., Basu, S., Bogart, R. S., Bush, R. I., Chitre,
S. M., Christensen-Dalsgaard, J., di Mauro, M. P., Dziembowski,
W. A., Eff-Darwich, A., Gough, D. O., Haber, D. A., Hoeksema,
J. T., Howe, R., Korzennik, S. G., Kosovichev, A. G., Larsen,
R. M., Pijpers, F. P., Scherrer, P. H., Sekii, T., Tarbell, T. D., Ti-
tle, A. M., Thompson, M. J. & Toomre, J. 1998 Helioseismic studies of
differential rotation in the solar envelope by the Solar Oscillations Investigation
using the Michelson Doppler Imager. Astrophys. J. 505, 390–417.
Soderblom, D. R., Oey, M. S., Johnson, D. R. H. & Stone, R. P. S.
1990 The evolution of the lithium abundances of solar-type stars. I — The
Hyades and Coma Berenices clusters. AJ 99, 595–607.
Spiegel, E. A. 1972 A history of solar rotation. NASA Special Publication 300,
61.
Spiegel, E. A. & Weiss, N. O. 1980 Magnetic activity and variations in solar
luminosity. Nature 287, 616–617.
Spiegel, E. A. & Zahn, J.-P. 1992 The solar tachocline. A&A 265, 106–114.
Spruit, H. C. 1999 Differential rotation and magnetic fields in stellar interiors.
A&A 349, 189–202.
Spruit, H. C. 2002 Dynamo action by differential rotation in a stably stratified
stellar interior. A&A 381, 923–932.
Steenbeck, M. & Krause, F. 1969 Zur Dynamotheorie stellarer und plan-
etarer Magnetfelder I. Berechnung sonnena¨hnlicher Wechselfeldgeneratoren.
Astron. Nachr. 291, 49–84.
Steenbeck, M., Krause, F. & Ra¨dler, K.-H. 1966 A calculation of the
mean electromotive force in an electrically conducting fluid in turbulent motion,
under the influence of Coriolis forces. Zeitschrift Naturforschung Teil A 21, 369.
Stein, R. F. & Nordlund, A˚. 1998 Simulations of solar granulation. I. General
properties. Astrophys. J. 499, 914–933.
BIBLIOGRAPHY 189
Stro¨mgren, B. 1938 On the helium and hydrogen content of the interior of the
stars. Astrophys. J. 87, 520–534.
Swenson, F. J. & Faulkner, J. 1992 Lithium dilution through main-sequence
mass loss. Astrophys. J. 395, 654–674.
Talon, S., Kumar, P. & Zahn, J.-P. 2002 Angular momentum extraction
by gravity waves in the Sun. Astrophys. J. 574, L175–L178.
Tao, L., Proctor, M. R. E. & Weiss, N. O. 1998 Flux expulsion by inho-
mogeneous turbulence. MNRAS 300, 907–914.
Tassoul, J.-L. & Tassoul, M. 1986 On the nature of the core-envelope bound-
ary layer in a slowly rotating star. Geophys. Astrophys. Fluid Dyn. 36, 303–315.
Taylor, J. B. 1963 The magneto-hydrodynamics of a rotating fluid and the
Earth’s dynamo problem. Proc. R. Soc. Lond. A 274, 274–283.
Tobias, S. M., Brummell, N. H., Clune, T. L. & Toomre, J. 1998
Pumping of magnetic fields by turbulent penetrative convection. Astrophys. J.
502, L177–L180.
Tobias, S. M., Brummell, N. H., Clune, T. L. & Toomre, J. 2001
Transport and storage of magnetic field by overshooting turbulent compressible
convection. Astrophys. J. 549, 1183–1203.
Tobias, S. M., Diamond, P. H. & Hughes, D. W. 2007 β-Plane magnetohy-
drodynamic turbulence in the solar tachocline. Astrophys. J. 667, L113–L116.
Vauclair, S., Vauclair, G., Schatzman, E. & Michaud, G. 1978 Hydro-
dynamical instabilities in the envelopes of main-sequence stars — Constraints
implied by the lithium, beryllium, and boron observations. Astrophys. J. 223,
567–582.
Vorontsov, S. V., Christensen-Dalsgaard, J., Schou, J., Strakhov,
V. N. & Thompson, M. J. 2002 Helioseismic measurement of solar torsional
oscillations. Science 296, 101–103.
Wale´n, C. 1946 Ark. Mat. Astr. Fys. 30A, no. 15.
190 Bibliography
Walker, M. R., Barenghi, C. F. & Jones, C. A. 1998 A note on dynamo
action at asymptotically small Ekman number. Geophys. Astrophys. Fluid Dyn.
88, 261–275.
Wang, Y.-M. & Sheeley, Jr., N. R. 1991 Magnetic flux transport and the
Sun’s dipole moment — New twists to the Babcock–Leighton model. Astrophys.
J. 375, 761–770.
Weiss, N. O. 1966 The expulsion of magnetic flux by eddies. Proc. R. Soc.
Lond. A 293, 310–328.
Weymann, R. & Sears, R. L. 1965 The depth of the convective envelope on
the lower main sequence and the depletion of lithium. Astrophys. J. 142, 174.
Wood, T. S. & McIntyre, M. E. 2007 Confinement of the Sun’s interior
magnetic field: Some exact boundary-layer solutions. In Unsolved Problems in
Stellar Physics: A Conference in Honor of Douglas Gough (ed. R. J. Stan-
cliffe, J. Dewi, G. Houdek, R. G. Martin, & C. A. Tout), American Insti-
tute of Physics Conference Series, vol. 948, pp. 303–308. Corrigendum at
www.atm.damtp.cam.ac.uk/people/mem/papers/SQBO/solarfigure.html#chirality.
Wood, T. S. & McIntyre, M. E. 2010 Confinement of the Sun’s interior
magnetic field by laminar magnetostrophic flow. Submitted to the Journal of
Fluid Mechanics .
Zahn, J.-P. 1992 Circulation and turbulence in rotating stars. A&A 265, 115–
132.
Zahn, J.-P., Talon, S. & Matias, J. 1997 Angular momentum transport by
internal waves in the solar interior. A&A 322, 320–328.
Zel’dovich, Ya. B. 1957 The magnetic field in the two-dimensional motion of
a conducting turbulent fluid. Soviet Journal of Experimental and Theoretical
Physics 4, 460.