Table Of ContentACCEPTEDFORPUBLICATIONINTHEASTROPHYSICALJOURNALON2012/01/30
PreprinttypesetusingLATEXstyleemulateapjv.08/22/09
HOWFARCANMINIMALMODELSEXPLAINTHESOLARCYCLE?
R.D.SIMITEV
SchoolofMathematicsandStatistics,UniversityofGlasgow,Glasgow,G128QW,UK
AND
F.H.BUSSE
InstituteofPhysics,UniversityofBayreuth,95440Bayreuth,Germany
AcceptedforpublicationintheAstrophysicalJournalon2012/01/30
2
1 ABSTRACT
0 Aphysicallyconsistentmodelofmagneticfieldgenerationbyconvectioninarotatingsphericalshellwith
2 a minimumofparametersis appliedto the Sun. Despite its unrealisticfeaturesthe modelexhibitsa number
n ofpropertiesresemblingthoseobservedontheSun. Themodelsuggeststhatthelargescalesolardynamois
a dominatedbyanon-axisymmetricm=1componentofthemagneticfield.
J Subjectheadings:Convection–Magnetohydrodynamics–Sun:dynamo–sunspots
0
3
1. INTRODUCTION fermuchfromunity in the highlyturbulentsolar convection
] zone, it must be realized that numerical simulations of con-
Duringthepasttwodecadestheunderstandingofthe pro-
R
vection driven dynamos depend sensitively on the parame-
cessofmagneticfieldgenerationbymotionsofanelectrically
S ters P and P near their values of unity. In particular it has
conductingfluidhasprogressedrapidlyandincreasinglyreal- m
. been found that a transition between two different types of
h isticmodelsofdynamosoperatinginplanetaryinteriorshave
dynamos occurs which is hysteretic as a function of P and
p beenachieved.Modelsofthegeodynamo,oftheancientMar-
P (Simitev&Busse 2009). It thus appears to be of inter-
- tiandynamoorofMercury’sdynamoareprimarilyrestricted m
o esttostudysimpleconvectiondrivendynamosindependence
bythelackofknowledgeofcertainpropertiesoftherespec-
r on these two parametersin addition to the Rayleigh number
t tiveplanetarycore.ThedynamoprocessoperatingintheSun,
s ontheotherhand,isstillsubjecttocontroversieswithseveral andtherotationparameter.Morerealisticsolardynamomod-
a els involve many more parameters and because of the high
competingproposals.Thedensityvariationsofseveralorders
[ demandfor numericalresolutionthe opportunitiesfor com-
ofmagnitudeinthesolarconvectionzoneandthecomplica-
prehensivestudiesofparameterdependencesarerestricted.
1 tionsintroducedbythecompressibilityofthefluidhavepre-
v ventedsofarafullyconvincingmodelofthesolarcycle. For the present study a numerical model for convection
7 Inthissituationitcanbehelpfultoconsidertheproperties drivendynamosthatpreviouslyhasbeenemployedforinves-
2 of a most simple physically consistent model of convection tigations motivated by problemsof the geodynamo(see, for
2 driven dynamos in a rotating spherical shell in order to see example,Simitev&Busse2005)hasbeenmodifiedforsolar
6 applications. A muchthinnersphericalshellis assumedand
whichpropertiesmightresemblethoseofthesolarcycle.As-
. no-slipconditionsattheinnerboundaryandstress-freecondi-
1 sumingagivenradiusratioofthesphericalfluidshellthemin-
tionsattheouterboundaryareused.TheBoussinesqapproxi-
0 imum number of dimensionless parameters is four, namely
mationhasbeenretainedeventhoughitisnotappropriatefor
2 theRayleighnumberR,therotationparameterτ, thePrandtl
any realistic solar dynamo model and causes an unrealistic
1 numberPdenotingtheratioofviscousandthermaldiffusiv-
dependenceofthesolardifferentialrotationondepth.Never-
: ity,andthemagneticPrandtlnumberP denotingtheratioof
v m theless,aswehopetodemonstrate,therearepropertiesofthe
viscousandmagneticdiffusivity.
i dynamosolutionsthatresemblethoseobservedontheSun.
X The success of numericalsimulationsof the generationof
magnetic fields by turbulent fluid motions relies of the con-
r 2. MATHEMATICALFORMULATION
a ceptofeddydiffusivitieswhichrepresentthediffusiveactions
of small scale turbulence that remains unresolvedin the nu- We consider a spherical fluid shell rotating about a fixed
merical schemes. The ratios of those eddy diffusivities are verticalaxis.Thebasicequationspermitastaticsolutionwith
expressed by the values of P and P used in the numerical thetemperaturedistribution
m
computations. FrequentlytheassumptionP=Pm=1ismade T =T +∆Tηr- 1(1- η)- 2,
basedontheargumentthattheturbulentdiffusionaffectsthe S 0
transportof heat, momentumand magnetic flux in the same wherer denotesthedistancefromthecenterofthespherical
way.Thisassumptionisnotcorrect,however,sincethetrans- shellmeasuredinmultiplesoftheshellthicknessd. Theratio
ports of scalar and vector quantities are affected in different ofinnertoouterradiusoftheshellisdenotedbyη.∆T isthus
ways (Kaysetal. 2005; Yousefetal. 2003). Moreover, the thetemperaturedifferencebetweenthe twoboundaries. The
diffusivitiesintheabsenceofanymotionoftendifferbymany gravityfield is givenby g=- dγr.In additionto d, the time
ordersofmagnitude,-inparticularwhentheradiativetrans- d2/ν,thetemperatureν2/γαd4,andthemagneticfluxdensity
port is taken into account in the heat diffusivity -, such that
ν(µ̺)1/2/d are used asscales for the dimensionlessdescrip-
even highly turbulentfluids do notfully erase differencesin
tion of the problemwhereν denotesthe kinematic viscosity
theeffectivediffusivities.
of the fluid, κ its thermal diffusivity, ̺ its density, and µ its
Even if it is assumed that ratios of the effective diffu-
magnetic permeability. The equations of motion for the ve-
sivities for heat, momentum and magnetic flux do not dif- locityvectoru,theheatequationforthedeviationΘfromthe
2 SimitevandBusse
R
105
10 7 1.2 1.41.61.8
5 4 3 0.8 1
2 0.6
1.5
Pm P
Figure1. (Coloronline)Convection-drivendynamosasafunctionofR,PandPmforτ=2000.Decayingdynamosareindicatedbyblackdots,MDdynamosare
indicatedbybluediamonds,FDdynamosareindicatedbyredstars,quadrupolardynamosareindicatedbygreentriangles,andcoexistingMDandFDdynamos
areindicatedbypinksquares.
static temperaturedistribution,andtheequationofinduction boundaries,
forthemagneticfluxdensityBarethusgivenby
v=∂ v=w=Θ=0 atr=r ≡η/(1- η) (9)
r i
∂tu+u·∇u+τk×u=- ∇π+Θr+∇2u+B·∇B, (1) v=∂r2rv=∂r(w/r)=Θ=0 atr=ro≡1/(1- η). (10)
∇·u=0, (2)
For the magnetic field, electrically insulating boundariesare
P(∂Θ+u·∇Θ)=(Rηr- 3(1- η)- 2)r·u+∇2Θ, (3) assumedsuchthatthepoloidalfunctionhmustbematchedto
t
thefunctionh(e), whichdescribesthe potentialfieldsoutside
∇·B=0, (4)
thefluidshell
∇2B=P (∂B+u·∇B- B·∇u), (5)
m t g=h- h(e)=∂ (h- h(e))=0 atr=r andat r=r . (11)
r i o
where∂ denotesthe partialderivativewith respectto timet
t Throughoutthispaperweshalluseη=0.65.
andwherealltermsintheequationofmotionthatcanbewrit-
Thenumericalintegrationproceedswithapseudo-spectral
ten as gradients have been combined into ∇π. The Boussi-
methodasdescribedbyTilgner(1999)whichisbasedonan
nesq approximation is assumed in that the density ̺ is re-
expansion of all dependentvariables in spherical harmonics
gardedasconstantexceptinthegravitytermwhereitstemper-
fortheθ,ϕ-dependences,i.e.
aturedependence,givenbyα≡- (∂̺/∂T)/̺=const.,istaken
intoaccount. The RayleighnumbersR, theCoriolisnumber v=XVlm(r,t)Plm(cosθ)exp(imϕ) (12)
τ,thePrandtlnumberPandthemagneticPrandtlnumberP
m l,m
aredefinedby
andanalogousexpressionsforw,Θ,handg. HerePmdenotes
l
αγ∆Td4 2Ωd2 ν ν theassociatedLegendrefunctions. Forther-dependenceex-
R= , τ = , P= , P = , (6)
m pansionsinChebychevpolynomialsareused.Forthecompu-
νκ ν κ λ
tationstobereportedinthefollowingaminimumof41col-
whereλisthemagneticdiffusivity.Becausethevelocityfield locationpointsintheradialdirectionandsphericalharmonics
u as well as the magnetic flux density B are solenoidalvec- uptotheorder128havebeenused.
torfields, thegeneralrepresentationintermsofpoloidaland Thenumericalmodeloutlinedaboveisessentiallythesame
toroidalcomponentscanbeused as has been used by Gilman&Miller (1981) and Gilman
(1983). At that time, 30 years ago, computer resources al-
u=∇×(∇v×r)+∇w×r , (7)
lowedtimeintegrationsonlyforafewselectedparameterval-
B=∇×(∇h×r)+∇g×r . (8) ues and the phenomenon of bistability (see below) was not
yet known. It thus seems appropriate to return to the most
Equations for v and w are obtained by multiplication of the
simplephysicallyconsistentformulationofconvectiondriven
curl2andofthecurlofequation(1)byr. Analogouslyequa- dynamosinrotatingsphericalfluidshellsasformulatedinthis
tions for h and g are obtained through the multiplication of section.
equation(5)andofitscurlbyr.
No-slip boundary conditions will be used at the inner 3. DYNAMOSOLUTIONSINDEPENDENCEONPANDPm
boundary and stress-free conditions are applied at the outer In choosing parameter values in our model we follow
boundary,whilethetemperatureisassumedtobefixedatboth Gilman (1983)inselectingavaluearound2000forτ which
Howfarcanminimalmodelsexplainthesolarcycle? 3
MD FD
1000
200 (a) (c)
3
ˇ
150 M¯t
100 Mdip100 MMˇpp
2
Mx 10 Nui 50 M¯t ¯MMˇˇtp
1 0 Mt,p
(b) (d)
600
¯
1 Et
0 400
0.2 0.4 0.6 0.8 111 1.2 1.4 1.6 1.8 222
P E
ˇ
Figure2. (Coloronline)Coexistingnonlineardynamosolutionsatidentical 200 Eˇt Eˇt
pTahreamfluecteturavtainlugespoinlotihdealcamseagηne=ti0c.6e5neτrg=y2c0o0m0p,oRn=en1t5M0ˇ0p0(0soalniddPci/rPclmes=)0a.n2d. 0 EE¯tp Eˇp
themeanpoloidalmagneticenergycomponentMp(emptycircles)arescaled 0.0 0.5 1.0 1t.5 2.0 2.5 0.0 0.5 1.0 1t.5 2.0 2.5
ontheleftordinate, whiletheNusseltnumberatr=ri (soliddiamonds)is
scaledontherightordinate. CasesstartedfrominitialconditionsontheFD Figure3. (Coloronline)Coexistingdistinctnonlineardynamosolutionsat
branch areshown inredandconnected bydashedlines, andthose started identical parameter values – a MD (left column (a,b)) and a FD dynamo
frominitialconditionsontheMDbranchareshowninblueandconnected (rightcolumn(c,d))bothatR=1.5×106,τ =2×104,P=1andPm=5.
bydottedlines. Panels(a,c)and(b,d)showtimeseriesofmagneticdipolarenergydensities
andkineticenergydensities, respectively. ThecomponentXp isshownby
ˇ ˇ
thicksolidblackline,whileXt,Xp,andXt areshownbyred,greenandblue
lines,respectively,andlabeled.XstandsforeitherMorE.
in the case of the Sun corresponds to an eddy viscosity ν =
1.6×108m2/s.ValuesforRhaveusuallybeenselectedwithin showninfigure3. Herethedefinitions
arangeofupthreetimesthevalueneededfordynamoaction. M¯ =h|∇×(∇¯h×r)|2i/2, M¯ =h|∇g¯×r|2i/2 (13)
Otherwise the dynamosbecometoo chaotic andregularfea- p t
turesdisappear.Ourstudyhasthusfocussedonthevariations Mˇ =h|∇×(∇ˇh×r)|2i/2, Mˇ =h|∇gˇ×r|2i/2 (14)
p t
ofdynamoswiththeparametersPandP .
m
havebeenusedwheretheangularbracketsindicatetheaver-
Infigure1theexistenceofsolutionsforconvectiondriven ¯
ageoverthefluidshellandwherehreferstotheazimuthally
dynamos in the (P,P ,R)–space is outlined. In order to
achievedynamoactionmtheRayleighnumberRmustbesuffi- averagedcomponentof h and ˇh is given by ˇh=h- ¯h. Anal-
cientlyhighsuchthatthemagneticReynoldsnumberR ex- ogous relationships hold for the kinetic energy densities for
m
ceeds a critical value. The latter is inverselyproportionalto whichM,h,garereplacedbyE,v,w.
the magnetic Prandtl number P . Since our study focusses Near the end points of the bistability regions transitions
m
on reasonably regular dynamoswe have found it difficultto from MD- to FD-dynamos and vice versa may, of course,
obtain dynamos for small values of P since dynamos will occurowingtoexceptionallylargefluctuations. Suchatran-
m
become highly chaotic in this case or even decay. As indi- sition was observed in the long-time dynamo integration of
catedinfigure1thereexistalargevarietyofdynamosranging Goudard&Dormy (2008). As their convection driven dy-
from dipolar to quadrupolarones. At higher Rayleigh num- namo changed from the MD-state to the FD-state these au-
bersthedynamosareoftenofmixedpolarity. Whilemostof thorsnotedasolarlikebehaviorinbutterfly-likediagramsin
thedynamosolutionshavebeenobtainedfortheradiusratio whichtheyhadplottedtheaxisymmetricϕ-componentofthe
η=0.65,somecasesobtainedforη=0.6havebeenaddedin magneticfield,Bϕ,asafunctionoflatitudeandtime.Asmust
the figure since it was found that the dynamos depend only beexpectedfromthedecreasewithdepthoftheirdifferential
weaklyonthisparameter. rotation(Yoshimura1975),thewavesofthebutterflydiagram
Of particular interest is the distinction between MD- propagatedtohigherlatitudesinsteadoflowerones. Itturns
dynamoswhicharedominatedbyastrongandnearlysteady out, however, that non-axisymmetric components of the az-
axisymmetricdipolarmagneticfieldandtheFD-dynamosfor imuthalmagneticfieldarecomparableinamplitudetotheax-
whichmagneticenergiesoftheaxisymmetriccomponentsare isymmetriccomponentB . Itthusappearsto beappropriate
ϕ
far less than those of the non-axisymmetriccomponentsand tostudythesedynamosinmoredetail.
which exhibit a cyclic behavior. Over a considerable extent In figure 4 a typical example of a nearly periodic FD-
of the parameter range these two types of dynamosexist si- dynamo is shown. In the four columns of plots the devel-
multaneously at the same values of the external parameters opmentofvariouscomponentsofthemagneticfieldisshown
asdemonstratedbytheexampleoffigure2. Itisremarkable in five equidistantsteps in time coveringabouthalf a period
thatthetwotypesofdynamosexhibitnearlythesameNusselt ofthe magneticcycle. Asis evidentfromthe radialcompo-
number. Thispropertyindicatesthatbothtypesaresimilarly nent B displayed in the third column, the dynamo appears
r
efficient in the use of the available buoyancy. For more de- to operatemainly in one of two meridionalhemispheres. In
tails on the phenomenonof bistability we refer to the paper the second column maxima of the horizontal magnetic field
by Simitev&Busse (2009). A comparison of the magnetic at the depth of 0.1 from the surface have been indicated as
andkineticenergydensitiesforthe two typesof dynamosis a measure of the likelihood of the appearance of sunspots.
4 SimitevandBusse
Figure4. (Coloronline)Approximatelyhalfaperiodofadynamooscillationinthecaseτ =2000,R=120000,P=1.2,Pm=4.5. Thefirstcolumnshows
meridionallinesofconstantBϕontheleftandpoloidalfieldlines,rsinθ∂h/∂θontheright. Thesecondcolumnshowslinesofconstanthorizontalmagnetic
field(B2ϕ+B2θ)1/2atr=0.9correspondingto- 0.9,- 0.8,- 0.7,0.7,0.8,0.9ofitsmaximumabsolutevaluewiththesigndeterminedbythesignofBϕ,andthe
thirdcolumnshowslinesofconstantBratr=ro. ThelastcolumnshowsRe(∂gm=1/∂θ)totheleftandIm(∂gm=1/∂θ)totheright. Thefiverowsareseparated
equidistantlyintimeby∆t=0.0224.
Howfarcanminimalmodelsexplainthesolarcycle? 5
80 thelowerplotoffigure5. Aphaseshiftbetweenthepatterns
oftheupperandthelowerplotcanalsobenoticedinthatthe
60 signchangeofB occursafterthatofB . Thisfeatureagrees
r ϕ
withthewellknownsolarphenomenonthatthechangeinthe
40
signofB occursnearthemaximumofthesunspotcycle.
r
20
de 4. CONCLUSION
u 0
atit Mostmodelsofthesolarcyclearebasedonessentiallyax-
l
-20 isymmetricdynamos,butthereareexceptionssuchasthedy-
namoproposedby Ruzmaikinetal. (1988) based on a mean
-40
field approach. It is not quite clear how strong the non-
-60 axisymmetriccomponentsofthesolartoroidalmagneticfield
with wavenumbersm=1 or m=2 are in relation to the ax-
-80
isymmetric component. Estimates range from 0.1 based on
74.15 74.2 74.25 74.3 74.35 sunspotnumbers(Usoskinetal.2007)to0.43basedonpow-
time erfulX-rayflares (Zhangetal. 2007). A main conclusionto
80 be derived from the dynamo model outlined in this paper is
that the large scale non-axisymmetriccomponentsare likely
60
to be an importantingredient of solar magnetism. The well
40 documentedexistenceofactivelongitudescertainlysupports
thisconclusion.
20 There are, of course, also major differences between so-
de larphenomenaandcertainfeaturesofthepresentBoussinesq
u 0
atit model. Someofthesearewellunderstoodasforexamplethe
l
-20 differenceinthepropagationofdynamowaves. Thisismost
probably caused by the fact that the model always exhibits
-40
adifferentialrotationdecreasingwithdepthforFD-typedy-
-60 namos,whiletheSunpossessesadifferentialrotationthatin-
creases with depth in the uppermost 15% of the convection
-80 zone.Forapreliminaryattempttotakethiseffectintoaccount
74.15 74.2 74.25 74.3 74.35 we refer to a recent companion paper by Simitev&Busse
time (2012).
Figure5. "Butterfly"diagrams: Bm=0+|Bm=1|sgn(Bm=0)(top)andBm=0+
ϕ ϕ ϕ r
|Bmr=1|sgn(Bmr=0)(bottom)areplottedasfunctionsoflatitudeandtimeinthe Theresearchof R.S. hasbeenpartlysupportedbythe UK
caseP=1.2,τ=2000,R=120000,Pm=4.5. RoyalSocietyunderResearchGrant2010R2. Theresearch
Since the underlying hemispherical structure hardly moves ofF.B.hasbeensupportedbyNASAGrantNNX09AJ85G.
the pattern of the second column seems to reflect the phe-
nomenonof the active longitudeof sunspots(see, for exam-
REFERENCES
ple, Usoskinetal. (2007) and referencestherein) or of solar
X-rayflares(Zhangetal.2007). Sincesometimesthem=2-
component of the magnetic field dominates instead of the Gilman,P.A.,1983,Astrophys.J.Suppl.,53,243
Gilman,P.A.,AndMiller,J.,1981,Astrophys.J.Suppl.,45,351.
m=1-component the phenomenon of two active longitudes
Goudard,L.,&Dormy,E.2008,Europ.Phys.Letts.,83,59001
180 degrees apart can also be observed in the dynamo sim-
Kays,W.,Crawford,M.,andWeigand,B.,ConvectiveHeatandMass
ulations. Persistent active longitudes have been detected in Transfer,2005(McGraw-Hill,4thedn)
varioustypesofactivestarsaswell(Usoskinetal.2007). Ruzmaikin,A.A.,Sokoloff,D.D.,&Starchenko,S.V.,1988,SolarPhysics,
Infigure5two"butterfly"diagramshavebeenplotted. Be- 115,5-15
SimitevR.,&BusseF.H.,2005,J.FluidMech.,532,365
causethestrengthoftheaxisymmetriccomponentofB alone
ϕ SimitevR.,&BusseF.H.,2010,Europ.Phys.Letts.,85,19001
isnotrepresentativeofthehorizontalmagneticfieldthatmay SimitevR.,&BusseF.H.,2012,PhysicaScripta,(inpress)
giverisetosunspotswehaveaddedintheupperplotthemain TilgnerA.,1999,Int.J.Num.Meth.Fluids,30,713
non-axisymmetriccomponentcharacterizedbytheazimuthal Usoskin,I.G.,Berdyugina,S.V.,Moss,D.,&Sokoloff,D.D.,2007,Adv.
wavenumberm=1.Whilethereissomewave-likemovement SpaceRes.40,951-958
Yousef,T.A.,Brandenburg,A.,andRüdiger,G.,2003,Astron.&
towardslowerlatitudes,themainpatternmovestowardshigh
Astrophys.,411,321–327
latitudes.Thismustbeexpectedsincethedifferentialrotation Yoshimura,H.,1975,Astrophys.J.,201,740
decreases with depth Yoshimura (1975). The movementto- Zhang,L.Y.,Wang,H.N.,Du,Z.L.,Cui,J.M.,andHe,H.,2007,Astron.&
wardhigherlatitudesisevenmorepronouncedinthepattern Astrophys.,411,321–327
oftheradialcomponentB ofthemagneticfieldasshownin
r
