Meridional Circulation in Solar and Stellar Convection Zones

8.3 MB·English

by  Nicholas A. Featherstone

Draft version January 27, 2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
MERIDIONAL CIRCULATION IN SOLAR AND STELLAR CONVECTION ZONES
Nicholas A. Featherstone
JILA,UniversityofColorado,Boulder,CO80309-0440
Mark S. Miesch
HighAltitudeObservatory,NationalCenterforAtmosphericResearch,Boulder,CO,80307-3000,USA:[email protected]
Draft version January 27, 2015
5
1
ABSTRACT
0
2 We present a series of 3-D nonlinear simulations of solar-like convection, carried out using the
Anelastic Spherical Harmonic (ASH) code, that are designed to isolate those processes that drive and
n
shapemeridionalcirculationswithinstellarconvectionzones. Thesesimulationshavebeenconstructed
a
soastospanthetransitionbetweensolar-likedifferentialrotation(fastequator/slowpoles)and“anti-
J
solar’ differential rotation (slow equator/fast poles). Solar-like states of differential rotation, arising
6
whenconvectionisrotationallyconstrained, arecharacterizedbyaverydifferentconvectiveReynolds
2
stress than anti-solar regimes, wherein convection only weakly senses the Coriolis force. We find that

the angular momentum transport by convective Reynolds stress plays a central role in establishing
]
R the meridional flow profiles in these simulations. We find that the transition from single-celled to
multi-celled meridional circulation profiles in strong and weak regimes of rotational constraint is
S
linked to a change in the convective Reynolds stress, a clear demonstration of gyroscopic pumping.
.
h Latitudinalthermalvariationsdifferbetweenthesedifferentregimes,withthoseinthesolar-likeregime
p conspiringtosuppressasinglecellofmeridionalcirculation,whereasthecoolpolesandwarmequator
- established in the anti-solar states tend to promote single-celled circulations. Though the convective
o
angular momentum transport becomes radially inward at mid-latitudes in anti-solar regimes, it is
r
t the meridional circulation that is primarily responsible for establishing a rapidly-rotating pole. We
s
conclude with a discussion of how these results relate to the Sun, and suggest that the Sun may lie
a
[ near the transition between rapidly-rotating and slowly-rotating regimes.

1
v 1. INTRODUCTION to dynamo modeling rests on the availability of suffi-
1 cientobservationaldatatoconstrainthemeanflowsthat
Mean (longitudinally-averaged) flows play an essential
0 are employed. However, this observational data is lim-
roleinallrecentmodelsofthesolaractivitycycle. Differ-
5 ited. Global helioseismic inversions provide a reliable
ential rotation (the mean zonal flow) generates toroidal
6 measure of the time-averaged angular velocity profile
field from poloidal field and amplifies it, tapping the ki-
0 Ω(r,θ)throughoutmostofthesolarconvectionzone(CZ,
neticenergyofshear. Thisisthewell-knownΩ-effectand
.
1 it is generally held responsible for the observed dispar- Thompson et al. 2003; Howe 2009), but as discussed in
0 itybetweenthetoroidalfluxinferredfrombipolaractive §2.2 relatively little is known about the MC. Even less
5 is known about mean flows in other stars (see §7). Yet,
regions and the much smaller poloidal flux inferred from
1 stars provide an invaluable test of dynamo models; any
high-latitude magnetograms.
: viablemodelofthesolaractivitycyclemustalsoaccount
v The role of the mean meridional circulation (MC) in
for cyclic magnetic activity in solar-like stars.
i the operation of the solar dynamo is more controversial
X Without some theoretical guidance as to what merid-
but potentially no less profound. As discussed in §2.1,
ional and zonal flows one should expect to find in solar
r mean-field dynamo models have shown that the ampli-
a and stellar convection zones, kinematic dynamo models
tude and structure of the MC can have a substantial
cannotcompletelyaccountforobservedpatternsinsolar
influence on cycle properties, and in particular, it may
and stellar magnetic activity. Furthermore, theoretical
largely regulate the cycle period.
guidance is needed to help motivate and interpret ongo-
Yet, despite the central role that mean flows play in
ing observational investigations of the solar meridional
solar and stellar dynamo models, their existence is often
circulation(§2.2)whichseektoprobemoregenerallythe
takenforgranted. Mostmean-fielddynamomodelspub-
internal dynamics of the Sun and the origins of solar
lished in the literature are kinematic in the sense that
magnetic activity.
the velocity field is specified so the origin of mean flows
In a previous paper, (Miesch et al. 2012, hereafter Pa-
liesoutsidethescopeofthemodel. Non-kinematicmean-
per 1), we sought such theoretical guidance based only
fieldmodelsdoexistbutmanytakeaheuristicapproach,
on the observed structure of solar mean flows and ro-
introducingparameterizationsforconvectivemomentum
bust dynamical balances obtained from the fundamen-
and energy transport that produce mean flows with pre-
tal magnetohydrodynamic (MHD) equations. From this
conceived amplitudes and profiles (e.g. Rempel 2006).
alone we were able to identify gyroscopic pumping as
Theutilityofthesekinematicandheuristicapproaches
the principle physical mechanism that maintains the so-
lar MC. This work is summarized in §2. However, fur-
[email protected]
2 Featherstone & Miesch
ther progress requires an understanding of turbulent so- et al. 2012), it is clear that the meridional flow plays an
lar convection and in particular the convective angular important role in redistributing the magnetic flux that
momentum transport. To this end, in this paper, we threads through the solar surface.
have initiated a series of global convection simulations These observations and models of surface flux trans-
designed to elucidate the nature of convective transport porthelptojustifytheprominentrolethattheMCplays
inrotating,stratified,sphericalshellsandthemeanflows in many recent dynamo models of the solar cycle. This
it produces. We focus here on solar-like stars with con- role is most dramatic in a class of models known as flux-
vectiveenvelopesbutsimilardynamicsarelikelytooccur transport (FT) dynamo models, in which the equator-
in other types of stars with different convection zone ge- ward advection of toroidal flux by the meridional flow
ometries. in the CZ and tachocline accounts for the observed mi-
In§2wesurveyhowthisseriesofsimulationsfitswithin gration of sunspots toward the equator over the course
the context of previous theoretical and observational of a cycle, known as the solar butterfly diagram (Wang
work and in §3 we discuss the motivation and method- & Sheeley 1991; Dikpati & Charbonneau 1999; Ku¨ker
ology behind how the numerical experiments were set et al. 2001; Dikpati & Gilman 2006; Yeates et al. 2008;
up. We then present an overview of our results in §4, Munoz-Jaramillo et al. 2009; Jouve & Brun 2007). By
identifying two distinct dynamical regimes that exhibit regulating the solar butterfly diagram, the MC also reg-
qualitatively different mean flows. As in previous stud- ulates the cycle period.
ies (see §2.3), the dynamical regimes are distinguished In a subset of flux-transport dynamo models that op-
by the degree of rotational influence on the convection, erate in the so-called advection-dominated regime, the
as quantified by the Rossby number Ro = U/(2Ω L), MC plays another important role. In practice, most ex-
0
whereU andLarecharacteristicconvectivevelocityand isting advection-dominated flux-transport (FT) dynamo
lengthscalesandΩ ismeanrotationrate. Fastrotators models are also Babcock-Leighton (BL) dynamo models
0
produce solar-like Ω profiles (fast equator, slow poles) inwhichthemainsourceofpoloidalfluxliesinthebuoy-
with multi-celled MC profiles while slow rotators pro- ant destabilization, emergence and subsequent dispersal
duce anti-solar Ω profiles (slow equator, fast poles) with of toroidal flux (Dikpati & Gilman 2009; Charbonneau
single-celledMCprofiles(polewardflowintheupperCZ 2010). Mean-field parameterizations of the BL mecha-
andequatorwardinthelowerCZ).Wedemonstratethat nism typically take the form of a poloidal source term
these slow and fast regimes can be achieved not only by confined to the surface layers (Dikpati & Charbonneau
changing the rotation rate but also by changing the me- 1999; Rempel 2006; Munoz-Jaramillo et al. 2010, e.g.).
chanical and thermal dissipation, which in turn impacts Thus, in order for the dynamo to operate, there must
the buoyancy driving. In §5 we interpret the resulting be a coupling mechanism that transports magnetic flux
MC profiles as a consequence of the changing nature of from the poloidal source region near the surface to the
the convective Reynolds stress, within the context of gy- regionwherethetoroidalfieldisgenerated,whichgener-
roscopic pumping. We address the establishment and allyliesinthelowerCZ.Inadvection-dominatedBL/FT
significance of thermal gradients in §6. dynamo models, this coupling mechanism is the MC.
Which regime is the Sun in? Since it has a solar-like Changes in the amplitude and structure of the merid-
Ω profile by definition, one would expect it to be classi- ional flow in FT dynamo models can have a dramatic
fied as a rapid rotator. However, the situation is more effect on the length and amplitude of magnetic cycles.
subtle than this. In §7 we demonstrate that the tran- Faster meridional flows or multiple cells in latitude can
sition between solar and anti-solar differential rotation shorten the cycle period and influence the polar field
is intimately linked to the MC and may be sensitive to strength while slower flows can lengthen cycles and sup-
the complex dynamics occurring in the boundary layers, press activity, potentially inducing a grand minimum
namely the tachocline and the near-surface shear layer. (Dikpati & Charbonneau 1999; Jouve & Brun 2007;
The Sun lies very close to this transition, which may Schrijver & Liu 2008; Jiang et al. 2009; Karak 2010).
imply that it has somewhat atypical mean flow profiles. Multiplecellsinradiuscanhaveanevengreaterimpact,
Thisinturnhasimportantimplicationsforsolardynamo disrupting the operation of the dynamo and potentially
modelsandtheirapplicationtootherstars. Theseissues reversing the sense of the butterfly diagram (Jouve &
tooareaddressedin§7andourprincipleresultsandcon- Brun 2007).
clusions are summarized in §8. The profound role of the MC in contemporary solar
dynamo models has motivated intense, ongoing observa-
tionaleffortstoprobethemeridionalflowstructureinthe
2. THEORETICALANDOBSERVATIONALBACKGROUND
deepCZ.However,thisisaformidablechallenge,largely
2.1. Meridional Circulation and the Solar Dynamo
because the flow is expected to be so weak; typical am-
The potential significance of the MC to the opera- plitudes used in FT dynamo models are only 2-3 m s−1
tion of the solar dynamo is apparent from the observed near the base of the CZ. Thus, observational constraints
evolution of the radial magnetic flux in the solar pho- are currently limited but may improve with more data,
tosphere. In particular, the cyclic polarity reversals of as we now address.
polar fields observed at the solar surface are often at-
2.2. Observational Context
tributed to the poleward transport of the trailing mag-
netic flux in emerging active regions due to the com- One of the most notable triumphs of helioseismology
bined action of the meridional flow and turbulent dif- hasbeenthemappingoftheinternalsolarrotationprofile
fusion (Sheeley 2005; Baumann et al. 2006; Schrijver & Ω as a function of radius r and colatitude θ (Thompson
Liu 2008). Though the timing of the most recent north et al. 2003; Howe 2009). These helioseismic rotation in-
polar reversal in 2012 challenges this paradigm (Shiota versions span most of the solar convection zone, though
Meridional Circulation in Stars 3
they become less reliable near the poles and in the deep ies suggested that this poleward flow may extend signif-
interior due both to the coarser resolution of the inver- icantly deeper, spanning at least the upper half of the
sion kernels and the smaller amplitude of the intrinsic CZ (Braun & Fan 1998; Chou & Dai 2001). However,
frequency splitting near the rotation axis. They reveal a more recent observations and inversions based on longer
monotonic decrease of Ω by about 30% from equator to time series, new analysis techniques, and new data from
pole that persists throughout most of the CZ, with lit- HMI/SDO have begun to cast doubt on this conclusion.
tle radial variation apart from two boundary layers near Hathaway (2012b) argued based on surface feature
the bottom and top, known as the tachocline and the tracking that the return equatorward flow required for
near-surface shear layer respectively. mass conservation may occur as shallow as 0.93R. In-
The most well-established result for the solar MC, ferring subsurface flows with feature tracking relies on
meanwhile,istheexistenceofasystematicpolewardflow the hypothesis that larger horizontal scales in the pho-
of10-20ms−1 atthesurfaceatlatitudesbelowabout± tosphere trace convective structures that extend deeper
50-60◦. However,differenttechniquesoftengivedifferent into the CZ and may thus be advected by deeper mean
results as to the subsurface structure, the high-latitude flows (Hathaway 2012a). Hathaway (2012a,b) tests this
structure, and even the low-latitude amplitude. For ex- hypothesis by comparing inferred zonal flows to helio-
ample, measurements of the meridional flow based on seismic rotational inversions and then applies it to infer
the correlation tracking of magnetic features generally themeridionalflow. Resultssuggestaflowreversalfrom
give somewhat lower flow speeds than those based on poleward to equatorward at a depth of roughly 50 Mm
Doppler shifts and local helioseismic inversions (Ulrich (r ∼ 0.93, where R is the solar radius). This is notably
2010). Dikpati et al. (2010b) attribute this difference to shallow compared to earlier estimates obtained from lo-
theeffectsofturbulentdiffusion. Muchofthevariationin cal helioseismic inversions.
amplitude estimates arises from the intrinsic variability More recent helioseismic investigations, fueled by the
oftheflowitself,whichisknowntochangewiththesolar high-qualitydatafromHMI/SDOandbytheavailability
cycle (Ulrich 2010; Basu & Antia 2010; Hathaway 2010; of long time series from GONG and SOHO/MDI1, sug-
Hathaway & Rightmire 2011). However, these solar cy- gest that this picture may be changing. Time-distance
cle variations can be understood as perturbations about helioseismic inversions by Zhao et al. (2013), taking into
a persistent background circulation (poleward near the account systematic center-to-limb effects on the inver-
surfaceatlow-midlatitudes)thatareinducedbysurface sions, appear to support this conclusion, showing evi-
magnetism and its impact on radiative cooling. In this dence for multiple reversals in the latitudinal flow com-
paperweareinterestedonlyinthepersistentbackground ponent, the shallowest of which being 0.91 R (see also
circulation so we disregard the solar cycle variations. Mitra-Kraev & Thompson 2007). However, such a shal-
Athigherlatitudestheresultsarelessclear. Therecent lowreversalisnotfoundinglobalinversionsoflow-order
Doppler analysis by Ulrich (2010) suggests the possible oscillation modes where the meridional flow enters as
presenceofahigh-latitudecounter-cellwithequatorward higher-order effect (Schad et al. 2012).
surface flow at latitudes above about 60◦. This feature Given the high stakes from the perspective of solar
is most apparent near solar minimum and may be in- and stellar dynamo theory, the current uncertainty with
termittent. Indeed, Dikpati et al. (2010a) argue that regard to the meridional flow structure in the deep con-
the absence of such a high-latitude counter cell through vection zone, and the vibrancy and promise of ongoing
much of solar cycle 23 may account for the relatively observational efforts, further theoretical insight is sorely
deepminimumof2008-2009andthedelayedonsetofso- needed. Here we seek such insight from convection sim-
lar cycle 24. However, evidence for such a high-latitude ulations, after a brief overview of the fundamental pro-
counter-cell based on feature tracking and helioseismic cesses involved (§2).
inversions has so far been lacking, with the former sug-
gestingthatpolewardflowmaypersistallthewaytothe 2.3. Dynamical Models
poles (Rightmire-Upton et al. 2012) and the latter being
As noted above, kinematic mean-field dynamo models
thus far inconclusive (Zhao et al. 2012). Further insight
cannot provide any insight into the origin and structure
into thisquestionshould be possibleinthe comingyears
oftheMC.Forthis,dynamicalmodelsarerequired. Such
withtheavailabilityoflongertimeseriesfromtheHelio-
modelsmustaddressthetransportofenergyandmomen-
seismic Magnetic Imager (HMI) on the Solar Dynamics
tum by convection and the subsequent nonlinear feed-
Observatory (SDO), which has the higher spatial resolu-
backs involving the MC, differential rotation, and ther-
tion needed to minimize limb effects, and also with the
mal stratification (Miesch & Toomre 2009). This can be
advent of future missions such as Solar Orbiter, which
doneeitherexplicitlythroughglobal3-Dconvectionsim-
will have a high-latitude vantage point.
ulations or implicitly through non-kinematic mean-field
Of these three principle techniques for measuring the
models that include parameterizations for the turbulent
meridional flow (Doppler measurements, feature track-
transport.
ing, and local helioseismic inversions), only helioseismic
Global convection simulations have revealed two dis-
inversionsandfeaturetrackingcanprobethemeridional
tinct dynamical regimes characterized by qualitatively
flow below the surface. While the observational consen-
differentdifferentialrotation(DR)profiles(Gilman1976,
sus is that there is a systematic poleward flow of about
1977; Glatzmaier & Gilman 1982; Aurnou et al. 2007;
10-20 m s−1 at latitudes below about 60◦ in the near-
Brun & Palacios 2009; Matt et al. 2011; Gastine et al.
surfacelayersoftheSun,thesubsurfacestructureismore
controversial. The first studies based on SOHO/MDI 1 the Global Oscillations Network Group (GONG) and the
data by Giles et al. (1997) suggested that the poleward Michelson Doppler Imager (MDI) on the Solar and Heliospheric
flow extended down to at least 0.96R. Subsequent stud- Observatory(SOHO)
4 Featherstone & Miesch
2013; Guerrero et al. 2013a; Gastine et al. 2014; K¨apyl¨a port. In such models the MC profile is largely deter-
et al. 2014). At strong rotational influence Ro (cid:46) 1 mined by departures from thermal wind balance due to
the differential rotation is “solar-like” in the sense that the meridional Reynolds stress, which is typically mod-
the equator spins faster than the poles. In other words eled as a turbulent diffusion (Kitchatinov 2012; Dikpati
∆Ω≡Ω −Ω >0 where Ω is the equatorial rotation 2014). This can lead to nonlocal driving of the MC from
eq p eq
rateattheouterboundaryofthedomain,andΩ issome the upper and lower boundary layers (Kitchatinov 2012)
p
measure of the rotation rate at high latitudes, poleward ortohigh-latitudecounter-cellsforsteepdensitystratifi-
of 70◦. For the weakly rotating regime, Ro(cid:38)1, the dif- cations(Dikpati2014). However,inthiscase,theformof
ferential rotation is referred to as “anti-solar”, meaning the MC profile is sensitive to the nature of the turbulent
∆Ω<0. Theanti-solarregimehasoftenbeenattributed viscosity profile that is imposed. This alone breaks the
to the tendency for the convection to conserve angular degeneracy of the thermal wind balance equation with
momentumlocally, spinningupasmotionsapproachthe respect to the meridional flow (see §5.1 below).
rotation axis (Gilman 1977; Aurnou et al. 2007). Whenthefulltime-dependentmean-fieldequationsare
Though these regimes can be achieved by varying solved, the driving of the meridional flow is fundamen-
the rotation rate Ω (or its non-dimensional equivalent), tallydifferent. HeretheMCprofileisdeterminedmainly
several studies have demonstrated that the transition by the non-diffusive component of the convective angu-
from solar to anti-solar differential rotation can also be lar momentum transport, which is typically modeled by
achieved by keeping Ω (or its non-dimensional equiva- the Λ-effect in the zonal momentum equation (Rempel
lent) fixed and changing the convective driving through 2005, 2007). This is the process of gyroscopic pumping
the Rayleigh number (Gastine et al. 2014; K¨apyl¨a et al. that will be described in detail in §5.1 and that is re-
2014). This tends to increase Ro by decreasing the con- sponsible for the establishment of MC in our convection
vective turnover time L/U. In other words, convective simulations, as demonstrated in §5.2.
velocityscalesincreasewhilelengthscalesdecrease. Near
the solar/anti-solar transition, Gastine et al. (2014) and 3. DESCRIPTIONOFNUMERICALCONVECTION
K¨apyl¨aetal.(2014)alsofindthattheDRexhibitsbista- EXPERIMENTS
bility,withbothsolarandanti-solarregimespossiblede- We present a series of numerical experiments designed
pending on initial conditions and hysteresis. to examine solar-like convection and its resulting mean
Many of these studies reported a qualitative change in flows under conditions that span a range of rotational
theMCprofileinthetworegimes(Mattetal.2011;Guer- constraints and convective efficiency. We model the so-
rero et al. 2013a; Gastine et al. 2013, 2014; K¨apyl¨a et al. lar convection zone using the anelastic approximation
2014). The anti-solar regime (Ro < 1) often exhibits (Gough 1969; Gilman & Glatzmaier 1980). Such an
a single circulation cell per hemisphere, with poleward approach is appropriate in deep stellar interiors where
flow in the upper CZ and equatorward flow in the lower plasma motions are subsonic and perturbations to ther-
CZ. Meanwhile, the solar regime (Ro>1) often exhibits modynamic variables are small compared to their mean,
multiple cells in latitude and radius, aligned with the horizontally averaged values. As thermodynamic per-
rotation axis outside the tangent cylinder2. turbations are small, thermodynamic variables are lin-
Recenthigh-resolutionmodelshavealsobeguntocap- earized about a spherically symmetric reference state
ture aspects of the solar near-surface shear layer (NSSL) with density ρ¯, pressure P¯, temperature T¯. and specific
wherein small-scale convection in the surface layers es- entropy S¯. Fluctuations about this state are denoted as
tablishes an inward Ω gradient (Guerrero et al. 2013b; ρ, P,T,andS.Intheuniformly-rotatingreferenceframe
Hotta et al. 2014). As suggested by Miesch & Hindman of the star, the anelastic equations under these assump-
(2011), the NSSL in these simulations is established by tions may be expressed as
aradiallyinwardangularmomentumtransport,ismain-
tained by meridional Reynolds stress, and is intimately ∇·(ρ¯v)=0, (1)
linked to a poleward meridional flow. This will be dis-
cussed further in §7.2 below.
The solar and anti-solar regimes in global convection (cid:20) (cid:21)
Dv
simulations stand in stark contrast to mean-field models ρ¯ +2Ω ×v =−∇P +ρg−∇·D, (2)
whichtypicallyproducesolar-likeDRprofiles(e.g.Ku¨ker Dt 0
et al. 2011; Kitchatinov 2012; Kitchatinov & Olemskoy
and
2012). Anti-solar profiles generally only occur if there
is a strong, baroclinic, single-celled meridional flow in-
duced by additional factors such as polar starspots or DS
ρ¯T¯ =∇·[κρ¯T¯∇S+κ ρ¯c ∇(T¯+T)]
tidal forcing from a binary companion (Kitchatinov & Dt r p
Ru¨diger 2004; Ru¨diger & Kitchatinov 2007). (cid:20) (cid:21) (3)
1
Mean-field models of solar-like stars generally exhibit +2ρ¯ν e e − (∇·v)2 .
ij ij 3
single-celled MC profiles but their origin and structure
vary depending on the modeling approach. One ap-
The velocity v expressed in spherical coordinates is
proach is to seek steady-state solutions of the mean-
v = (v ,v ,v ) relative to a frame rotating at constant
field equations with heuristic or model-based prescrip- r θ φ
angular velocity Ω , g is the gravitational acceleration,
tions for the convective Reynolds stress and heat trans- o
and D is the viscous stress tensor given by
2 The tangent cylinder is a cylindrical surface parallel to the (cid:20) 1 (cid:21)
rotationaxisandtangenttothebaseoftheCZ D =−2ρ¯ν e − (∇·v)δ , (4)
ij ij 3 ij
Meridional Circulation in Stars 5
wheree isthestrainratetensor. Thekinematicviscos- adopted for each simulation such that
ij
ityisdenotedbyν andthethermaldiffusivitybyκ. The
radiative diffusivity is indicated by κr, cp is the specific v = ∂(vθ/r) = ∂(vφ/r) =0| . (8)
heat at constant pressure, and D/Dt is the advective r ∂r ∂r r=rbot,rtop
derivative. κ and c remain constant in time.
r p The radial entropy gradient is forced to vanish at the
This set of equations is closed by assuming the ther-
lower boundary of the convection zone, and the entropy
modynamic fluctuations satisfy the linear relations
perturbationsareforcedtovanishattheupperboundary,
ρ P T P S with
ρ¯ = P¯ − T¯ = γP¯ − cp, (5) ∂∂Sr =0|r=rbot, S =0|r=rtop. (9)
assuming the ideal gas law
Thus, there is no diffusive entropy flux across the lower
P¯ =Rρ¯T¯, (6) boundary. Heat is transported across this boundary by
the radiative heat flux instead, ∝κ ∇T, which drops to
where R is the gas constant. r
nearzeroattheupperboundary(seeFig.10below). The
We evolve the anelastic equations using the anelastic
vanishing of the entropy gradient at the lower boundary
spherical harmonic (ASH) code, details of which may be
implies that this is, by definition, the base of the CZ,
foundinCluneetal. (1999)andBrunetal. (2004). ASH
withnoovershootregion. Theentropygradientsteepens
solves the 3-D MHD equations in a rotating, spherical
towardthetopofthedomain,wheremostoftheconvec-
shell using a pseudo-spectral approach. Spherical har-
tive driving occurs. This will be discussed further in §6.
monics are employed in the horizontal dimension, and
Withthevalueofentropypinnedattheupperboundary,
4th-order finite-differences in the radial direction. ASH
the entropy gradient there is free to vary with time. As
employs a poloidal/toroidal decomposition for the mass
the simulation equilibrates, this gradient attains a sta-
fluxvectortosatisfythedivergence-freeconstraintofeq.
tisticallysteadyvaluethatissufficienttoexpelonesolar
1, with
luminosity through the upper boundary via thermal dif-
fusion.
ρ¯v =∇×∇×(Wer)+∇×(Zer). (7) A summary of our input parameter space for those
cases rotating at the solar rate of Ω =2.6×10−6 s−1 is
Here W and Z are the poloidal and toroidal stream- (cid:12)
provided in Table 3. We vary the viscosities by a fac-
functions respectively, and are the quantities that are
tor of four and the thermal diffusivities by a factor of
explicitly evolved in time (as opposed to the individual
two for these simulations. The naming convention is
componentsofv). Theunitvectorintheradialdirection
such that the letter denotes the level of thermal diffu-
is indicated by e .
r sion, with A corresponding to the highest thermal diffu-
The suite of convection zone models that we exam-
sivity of 3.2×1013 cm2 s−1 and C corresponding to the
ine are identical in all respects save for the diffusivities
lowest level (1.6×1013 cm s−1). The numbers 0–3 indi-
and the rotation rate. The reference state temperature,
cate the level of viscous diffusion, with 0 corresponding
pressure, and functional form of κ are derived from a
r to the lowest kinematic viscosity (2×1012 cm s−1), and
one-dimensionalsolarstructuremodel(Brunetal.2002).
3 corresponding to the highest value of 8×1012 cm s−1.
While the specific heat c of the one-dimensional struc-
p CaseB2,whichliesinthecenterofourparameterspace,
ture model varies weakly throughout the solar convec-
serves as the point around which rotation rate is varied
tion zone, we have adopted a constant value of c , taken
p for the results presented in §4. We have also provided
from the base of that model. This requires us to itera-
a summary of the resulting energy balances, as well as
tivelysolveforthereferencestatedensityρ¯, usingthe1-
their relative ratios, for these cases in Table 2. There,
Dmodeldensityasaninitialguess,untiladensityprofile
we report on the globally-averaged (over the full spher-
that satisfies both hydrostatic balance and our equation
ical shell) convective kinetic energy density (CKE), the
ofstateisreached. Oursimulatedconvectionzonesspan
energyassociatedwiththedifferentialrotation(DRKE),
from 0.72 R to 0.965 R , extending to 25 Mm below
(cid:12) (cid:12) andtheenergyassociatedwiththemeridionalcirculation
the solar surface, and encompassing roughly 3.5 density
(MCKE). We define these as
scale heights. Convective perturbations about this refer-
ence state remain small throughout the simulation, and
thuswechoosenottoupdatethereferencestate,holding CKE = 1ρ(cid:2)(v −(cid:104)v (cid:105))2+(v −(cid:104)v (cid:105))2+(v −(cid:104)v (cid:105))2(cid:3),
2 r r θ θ φ φ
it fixed throughout each calculation.
(10)
As the microscopic viscous and thermal diffusivities
appropriate for plasma motions within the convection MCKE = 1ρ(cid:2)(cid:104)v (cid:105)2+(cid:104)v (cid:105)2(cid:3), (11)
zonearebeyondthecapabilitiesofcurrentcomputational 2 r θ
models, our models should be viewed as large eddy sim-
and
ulations. Weadoptaradialprofileforν andκwitheach
quantitydiminishingindepthasν ∼ρ−1/2. Simulations DRKE = 12ρ(cid:2)(cid:104)vφ(cid:105)2(cid:3), (12)
are thus more diffusive near the surface where we ex-
pect the unresolved scales of convective motion to play where the brackets denote averages over longitude.
a prominent role in the transport of heat and momen- Each model discussed in this paper has been evolved
tum. This formulation does not, however, capture any for several thousand days, with elapsed simulated time
non-dissipative behavior of the unresolved scales. ranging from 7,000 days to 14,000 days depending on
Impenetrable and stress free boundary conditions are theindividualcase. Thiscorrespondstoseveralhundred
6 Featherstone & Miesch
4.1. Identification of Mean Flow Regimes
TABLE 1
Variable Input Parameters Wefindthatthestructureofmeridionalcirculationre-
Case ν κ Pr Ro Ra/105 R∗ (cid:96) alized in our convection zone models is intimately linked
top top max tothe differentialrotationestablished. This relationship
A1 4 32 1/8 0.15 1.79 2.61 340
ismostreadilyillustratedthroughaseriesofexperiments
A2 6 32 3/16 0.13 1.19 2.61 170
wherein the stellar rotation rate is varied. Using Case
A3 8 32 1/4 0.11 0.89 2.61 85
B2asourstartingpoint,wehaveexploredrotationrates
B0 2 24 1/12 0.16 6.35 3.48 680
ranging from 2Ω to 0.75Ω .
(cid:12) (cid:12)
B0.5 3 24 1/8 0.16 4.24 3.48 340 We plot the resulting differential rotation (upper row)
B1 4 24 1/6 0.15 3.18 3.48 340 andmeridionalcirculation(lowerrow)profilesforthisse-
B2 6 24 1/4 0.13 2.12 3.48 170 riesofsimulationsinFig. 1Atthehigherrotationrates,
B3 8 24 1/3 0.12 1.59 3.48 170 most notably 2Ω and Ω , the resulting differential ro-
(cid:12) (cid:12)
C1 4 16 1/4 0.24 7.14 5.22 340 tation is decidedly solar-like, possessing a fast equator
C2 6 16 3/8 0.20 4.76 5.22 340 and retrograde poles (see §2.3). We note that contours
C3 8 16 1/2 0.12 3.57 5.22 170 in these simulations are noticeably cylindrical with re-
Note. — Input parameters for all cases. Viscous and spect to the Sun. This may be attributed to the lack of
thermal diffusivities (ν and κ respectively) are quoted in an underlying stable zone which is thought to promote
units of 1012 cm2 s−1. The Prandtl number, Pr, is con- conical Ω profiles by inducing baroclinic torques (Rem-
stant with depth. The Rossby number Ro= v˜/(2ΩL) is pel 2005; Balbus & Schaan 2012). At rotation rates less
quoted at mid-convection zone. The Rayleigh number Ra= than Ω , the solar-like differential rotation transitions
(cid:12)
(∂ρ/∂S)(∂S/∂r)gL4/(ρνκ) is quoted at the top of the sim- to a state characterized by a weakly retrograde equator
ulation. The length scale L is taken to be the depth of the and rapidly prograde poles. These prograde polar states
convection zone (1.72×1010 cm) in both cases. The modi- tend to possess a much stronger equator-to-pole rota-
fied Rayleigh Number R∗ = g|dSc/dr|/(cpΩ2) is quoted at tional contrast than their solar-like counterparts. Simi-
mid-convection zone. lar results were also reported by other authors, includ-
ing Gastine et al. (2013, 2014), Guerrero et al. (2013a),
K¨apyl¨a et al. (2014).
TABLE 2
As the lower row of Fig. 1 is traversed from left to
Energy Balance right,circulatorypatternsassumetwodistinctmorpholo-
Case CKE DRKE DRKE/CKE MCKE MCKE/CKE gies and exhibit a clear correlation with the structure of
A1 2.09 2.94 1.4091 0.019 0.0091 the differential rotation. In the fast-equator cases, the
A2 1.46 2.38 1.6241 0.012 0.0079 low latitudes are characterized by a counter cell of cir-
A3 1.22 1.77 1.4447 0.010 0.0080 culation occurring near mid-convection zone. High lati-
B0 4.39 0.59 0.1345 0.032 0.0073 tudes in these cases are largely characterized by a single
cell of circulation (poleward at the surface), with some
B0.5 3.52 1.87 0.5305 0.026 0.0075
hints of multi-cellular behavior at high latitudes for the
B1 2.99 1.73 0.5777 0.024 0.0079
casewiththefastestrotation. Astherotationrateofthe
B2 2.23 1.79 0.8024 0.017 0.0078
system decreases, the high latitude cell extends to lower
B3 1.96 1.31 0.6721 0.016 0.0084
latitudes, ultimately supplanting the counter cell in the
C1 4.66 9.92 2.1285 0.054 0.0115
equatorial regions for rotation rates lower than Ω . The
C2 3.73 6.04 1.6202 0.044 0.0117 (cid:12)
shift from multi-cellular MC profiles in the solar regime
C3 2.90 0.18 0.0625 0.025 0.0085
to single-celled profiles in the anti-solar regime has also
Note. — Globallyaveragedenergybalanceasrealizedin beennotedbyotherauthors(Guerreroetal.2013a;Gas-
each run. Energy densities have been averaged over 10 rota- tine et al. 2014; K¨apyl¨a et al. 2014).
tion periods at the end of each simulation, and are reported
in units of 106 erg cm−3
4.2. Transition Between Mean Flow Regimes at Fixed
rotational and convective overturning times. All simu-
Ω
lations have been well equilibrated viscously and ther- 0
mally. As Prandtl numbers are always less than unity, The series of simulations illustrated in Fig. 1 possess
the viscous diffusion timescale is the most limiting of acommonmid-corermsvelocityamplitudeofabout100
these two timescales. With the exception of case B0, m s−1. The transition in mean-flow regimes spanned by
which has been run for one viscous diffusion time, all those simulations is thus most naturally characterized
simulations have been run for at least two viscous dif- by a transition in timescale ratios, namely the ratio of
fusion times. We find that case B0 is transitioning to the convective overturning time to the rotation period.
an anti-solar differential rotation rate after one viscous This is quantified by the Rossby number, introduced in
diffusion time of evolution. As we are able to explore §1. We can similarly span this transition by holding Ω
anti-solar regimes adequately with the more computa- fixed and varying the amplitude of the convective veloc-
tionally tractable cases C1 and C2, we have chosen not ities. For systems with a fixed luminosity, such as those
to pursue the evolution of B0 further, but include it in under consideration here, this is most naturally accom-
Tables 1 and 2, as well as Figure 3, to help delineate the plishedbyvaryingthetransportcoefficientsν andκ,and
pointoftransitionbetweensolarandanti-solarbehavior. through them the Rayleigh number Ra. In so doing, we
can also explore how our transport parameters ν and
4. DISTINCTMEANFLOWREGIMESFORFASTAND κ contribute to the balances realized in the slowly- and
SLOWROTATORS rapidly-rotatingregimesbysamplingarangeofresulting
Meridional Circulation in Stars 7
Fig. 1.—(a-e)Differentialrotationprofiles(Ω−Ωframe)forcasesrotatingatdifferentfractionsofthesolarrotationrateΩ0 (indicated)
achievedbyspinningupandspinningdowncaseB2,whichrotatesatthesolarrate. Rotationratedecreasesfromlefttoright,andunitsare
indicatedinnHz. (f-j)Meridionalcirculationprofilescorrespondingtoeachofthedifferentialrotationprofilesintheupperrow. Redtones
denotecounter-clockwisemotion,andbluetones,clockwise. Profileshavebeenaveragedintimeoverroughly200daysineachinstance. As
theSunisspundown,thepolarregionsdevelopprograderotation,andtheequatorialregionsretrograde. Similarly,multi-cellularprofiles
ofmeridionalcirculationtransitionintosingle-celledprofilesatlowrotationrates.
Reynolds numbers. contrastinmodelssuchasours,analogoustothetemper-
Othershaveshownthatthetransitionbetweentheso- ature contrast in RBC systems, is not held fixed, but is
larandanti-solarrotationregimescanbeachievedbyfix- anoutputquantity, muchastheReynoldsnumber. This
ingtherotationrateandchangingtheconvectivedriving entropy contrast is built entirely in the upper boundary
by varying the Rayleigh number (Gilman 1977; Aurnou layer, resulting from the entropy gradient realized there
et al. 2007; Guerrero et al. 2013a; Gastine et al. 2014; and the width of that boundary layer. Convective driv-
K¨apyl¨a et al. 2014). The key parameter in determin- ingisthusaccomplishedprimarilyintheupperboundary
ing the transition is the Rossby number. However, since layersofthesesystems,asituationmoreanalogoustothe
this is an output of the simulation rather than an input, Sun than RBC flow where driving occurs at the upper
the transition is often described in terms of the modified and lower boundaries.
Rayleigh number R∗ (see Table 3), which is a measure Whileentropyisfixedatthetop, theentropygradient
of the relative strengths the buoyancy and Coriolis force is implicitly determined by the solar luminosity, being
(Gastine et al. 2013). The simulations presented in this forced to satisfy, in a time-averaged sense, the relation
sectionwillbeusedinsections5toexplorethestructure
∂s L
and origin of the meridional circulation, a topic that has ρ¯T¯κ | = (cid:12) . (13)
received less attention in previous studies. While our ∂r r=rtop 4πrt2op
demonstration of the solar/anti-solar transition is itself
WethuschoosetodefinetheourRayleighnumbersusing
not new, we find it useful to revisit some aspects of this
∂s/∂r at the upper boundary because this quantity is a
transition as they also bear on the resulting meridional
well-defined simulation input. The consequence is that
circulation. our Rayleigh number now varies as κ−2 instead of the
It is worth mentioning that the Rayleigh number in more traditional κ−1.
systemssuchasthese, whereentropyisheldfixedatone
The variation of differential rotation for each case is
boundary, and the entropy gradient fixed at the other,
illustrated in Fig. 2. There, profiles of angular velocity,
possesses an implicit dependence on κ not found in clas-
averaged over several overturning times, are laid out in
sical Rayleigh-B´enard convection (RBC). The entropy
a grid, with κ increasing along the vertical axis and ν
8 Featherstone & Miesch
Fig. 2.— Differentialrotationprofiles(Ω−Ω(cid:12))forthea,b,andc seriesofsimulations. Casenamesareindicatedandarearrangedsuch
thatviscosityν increasesfromlefttoright,andthermaldiffusivityκincreasesfrombottomtotop. Themostlaminarcase(A3)appears
intheupperright,andthemostturbulentcase(C1)appearsinthelowerleft. Profileshavebeenaveragedintimeoverroughly200days
in each instance. A clear transition between solar-like and anti-solar rotation occurs between the middle and bottom rows (i.e. between
theB andC seriescases.)
increasing along the horizontal axis. The upper row is Rayleigh number noted above.
comprised by simulations from Series A, and the lower As noted above, the key parameter in determining the
rowbySeriesCsimulations. ConvectionmodelsfromSe- transitionbetweenmeanflowregimesistheRossbynum-
ries A and B each exhibit solar-like differential rotation, ber Ro. Fig. 3 shows the relative differential rotation,
whereasseriesCtendstoexhibitanti-solarrotation,with ∆Ω/Ω =(Ω −Ω )/Ω ,inallsimulationsversusRo,
0 eq pole 0
case C3 (lower right) possibly representing a transition exhibiting a clear transition at R ∼ 0.17. This transi-
o
point between the fast equator and fast pole regimes. tional value is lower than the R ∼ 1 value quoted in
o
Variation of the differential rotation with ν is present, previous papers (e.g. Gastine et al. 2013, 2014). We at-
but the effect is greatly reduced with respect to that of tributethisdiscrepancytotheambiguityindefiningRo.
κ. This can be attributed to the ν−1κ−2 scaling of the Here we use the depth of the CZ as the length scale L
Meridional Circulation in Stars 9
0.5 One might ask if the transition in mean-flow regimes
is purely a result of rotational constraint, or if the level
of turbulence might also play some role. In particular,
by admitting much smaller spatial scales into the sim-
0.0
ulation, are the correlations of v and v necessary to
r φ
establishafastequatorbeingdisrupted? Inotherwords,
might the Reynolds number R = UL/ν play a role in
e
ΩΩ/ -0.5 athdedsietiqounetsotiotnhsewReahyalevieghexnteunmdbeedr?theInsimorudleartiotonseoxfpFloirge.
∆
2toincludethecaseB0.5. Thiscasepossessesaκidenti-
caltotheotherB-seriescases, butaν thatislowerthan
-1.0 any of of the C-series. Case B0.5 was initiated by lower-
ingthediffusivityinthematurecaseB1,whichhaditself
been evolved for a total of two viscous diffusion times.
An overview of the flows in Case B0.5 is shown in Fig.
-1.5
6. Theconvectivepatternsaresmall-scaleinnature,like
0.10 0.15 0.20 0.25 C3, but do not exhibit the retrograde tilting character-
Ro istic of the patterns in this case. Case B0.5 possesses
a decidedly solar-like differential rotation, but, surpris-
ingly, seems to exhibit a single circulation cell within
Fig. 3.— Pole-to-equatordifferentialrotationcontrast∆Ω/Ωas
each hemisphere. Moreover, the flow amplitudes indi-
afunctionofRossbynumberRofortheA−Cseriesofsimulations.
Transitional case B0.5 is indicated with a hollow diamond. The cated near the surface are similar to those of case C3,
slightlymoreturbulentcaseB0,whichistransitioningtoananti- as are the mid-convection zone velocities (not shown).
solar state, is indicated by the red triangle. Ro was computed at CaseB0.5doespossessasomewhatlowerRayleighnum-
mid-convection zone in all cases. Equatorial and polar rotation
ber, and a different Prandtl number than case C3, how-
rateswerecalculatedbyaveragingzonalvelocityinbothlongitude
and depth, and then over a latitude range of 0◦N-20◦N latitude ever, and so the differences between these two cases are
and 70◦N-90◦N respectively. Anti-solar cases tend to arise when likely to be subtle ones. It would seem that case B0.5 is
Roisgreaterthanabout0.163. just straddling the transition between the solar-like and
anti-solar regimes of differential rotation.
and the rms value of the non-axisymmetric velocity in
Gastine et al. (2014) and Ka¨pyl¨a et al. (2014) have
the mid-CZ v(cid:48) as the velocity scale U. By contrast,
rms recently shown that the behavior of the differential ro-
(Gastine et al. 2014) use the volume-averaged value of
tationnearthesolar/anti-solartransitioncandependon
v(cid:48) for U and a value of L based on the peak of the
rms the history of the simulation, representing a type of hys-
kinetic energy spectrum. This tends to shift U higher
teresis. For a given transitional Ro, the simulation may
and L lower, resulting in a larger value of Ro. We refer
exhibiteithersolaroranti-solarDR.Thoughwehaveob-
thereaderto(Gastineetal.2014)foradiscussionofhow
served unusual mean flows near the transition with case
thetransitionalvalueofRodependsonpreciselyhowRo
B0.5, we have not searched explicitly for this bistability
is defined.
andhysteresis. ItispossiblethatcaseB0.5itselfisslowly
The meridional circulation profiles corresponding to
undergoing a transition from solar-like to anti-solar dif-
thesimulationsshowninFig.2areshowninFig.4. They
ferential rotation. However, this case has been run for
exhibitacorrelationwiththedifferentialrotationsimilar
18,000 days, about a factor of four longer than the time
tothatdiscussedin§4.1. Single-celledprofilesofcircula-
scaleforviscousdiffusionacrosstheconvectionzone,and
tionappearinalloftheC-seriescases(i.e. thosesystems
we see no signs yet that it is undergoing such a transi-
with fast poles), while those of the A-B series are multi-
tion. The globally-averaged energy densities associated
cellular in the equatorial regions.
withthedifferentialrotation(DRKE)andmeridionalcir-
Convectioninthesolarandanti-solarregimesexhibits
culation(MCKE)forcaseB0.5areshowninFig. 7. The
noticeably different morphologies as well. Fig. 5 il-
mean flows have varied by as much as 30% over this in-
lustrates the radial velocities in the upper and mid-
terval, but with no apparent long-term trends. We will
convectionzonesofcaseC1(a,c)andcaseA3(b,d). Case
discuss the significance of case B0.5 further in §7.
C1, which possess a anti-solar-like differential rotation,
also exhibits velocity amplitudes that are roughly twice
that of A3. Elongated convective cells in the equatorial 5. MAINTENANCEOFMEANFLOWS
regions tend to tilt in the negative φ-direction for case In this section we argue that the distinct mean flow
C1. In addition to lower velocity amplitudes, Case A3 regimesdescribedin§4arisefromachangeinthenature
possesses much more prominent columnar structuring of oftheconvectiveangularmomentumtransport. Further-
its convection, with a tilting evident in the positive φ- more, we argue that it is the convective angular momen-
direction. tum transport that largely determines the mean merid-
An additional effect is evident when comparing the ional flow profiles in simulations as well as stars.
flowsofcaseA3tothoseofcaseC1,namelythatthecon-
vection exhibits a much broader range of spatial scales. 5.1. Dynamical Balances and Gyroscopic Pumping
Case C3 is thus both more turbulent and more strongly
In an anelastic (low Mach number) system, the mass
driven.
flux ρv is divergenceless and the time evolution of the
mean meridional flow ρ(cid:104)v (cid:105) is governed by the zonal
4.3. A Transitional Regime m
vorticity equation, which may be written as (Miesch &
10 Featherstone & Miesch
Fig. 4.— Contoursofmeridionalcirculationstreamlines, withstreamfunctionunderlay, forcasesintheA-C seriesasindicatedwithin
each panel. The layout and time averaging are as in Fig. 2. Red tones denote counter-clockwise flow, and blue tones clockwise. As with
caseswheretherotationrateisvaried,multi-cellularcirculationsarisewheneversolar-likedifferentialrotationprofilesarepresent.
Hindman 2011, hereafter MH11) dinal θˆ components of a vector. We have focused here
onthetwodominantcontributorstothemeridionalforc-
∂ ∂Ω2 g ∂(cid:104)S(cid:105) ing, namely the inertia of the differential rotation and
(cid:104)ω (cid:105)=λ − +... , (14)
∂t φ ∂z rC ∂θ the baroclinicity of the mean stratification, which make
P
up the first and second terms on the right-hand-side of
where ω = (∇×v )·φˆ is the zonal component of the (14) respectively. These are sufficient to illustrate the
φ m
mean vorticity, λ = rsinθ and z = rcosθ are the radial essential dynamics. We have neglected the meridional
and axial coordinates in a (rotating) cylindrical coordi- components of the Reynolds stress ((cid:104)ρ¯v(cid:48) v(cid:48) (cid:105); primes
m m
nate system, and Ω = Ω +λ−1(cid:104)v (cid:105) is the total angu- denote fluctuations about the longitudinal mean), the
0 φ
lar velocity, including uniform and differential rotation baroclinicity of thermal fluctuations, the viscous stress,
components. We have adopted the convention that a and quadratic terms in vm. All of these are thought to
subscript m on a vector denotes the radial rˆand latitu- be negligible in simulations and stars. We have also ne-