Numerical simulations of the magnetorotational instability in protoneutron stars: I. Influence of buoyancy

2.2 MB·English

by  Jerome Guilet

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed27April2015 (MNLATEXstylefilev2.2)
Numerical simulations of the magnetorotational instability
in protoneutron stars : I. Influence of buoyancy
J´erˆome Guilet1,2 & Ewald Mu¨ller1
1 Max-Planck-Institut fu¨r Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany
2 Max Planck/Princeton Center for Plasma Physics
5
1
27April2015
0
2

r ABSTRACT
p
The magneto-rotational instability (MRI) is considered to be a promising mecha-
A
nism to amplify the magnetic field in fast rotating protoneutron stars. In contrast to

4 accretion discs, radial buoyancy driven by entropy and lepton fraction gradients is
2 expected to have a dynamical role as important as rotation and shear. We investigate
the poorly known impact of buoyancy on the non-linear phase of the MRI, by means

] of three dimensional numerical simulations of a local model in the equatorial plane
E
of a protoneutron star. The use of the Boussinesq approximation allows us to utilise
H
a shearing box model with clean shearing periodic boundary conditions, while taking
. into account the buoyancy driven by radial entropy and composition gradients. We
h
find significantly stronger turbulence and magnetic fields in buoyantly unstable flows.
p
- On the other hand, buoyancy has only a limited impact on the strength of turbulence
o andmagneticfieldamplificationforbuoyantlystableflowsinthepresenceofarealistic
r thermaldiffusion.Thepropertiesoftheturbulenceare,however,significantlyaffected
t
s in the latter case. In particular, the toroidal components of the magnetic field and
a
of the velocity become even more dominant with respect to the poloidal ones. Fur-
[
thermore, we observed in the regime of stable buoyancy the formation of long lived

2 coherent structures such as channel flows and zonal flows. Overall, our results sup-
v port the ability of the MRI to amplify the magnetic field significantly even in stably
6 stratified regions of protoneutron stars.
3
6 Key words: stars: magnetars – stars: neutron – supernovae: general – MHD – in-
7 stabilities – magnetic fields
0
.
1
0
5
1 INTRODUCTION et al. 2012; Mo¨sta et al. 2014). The magnetic field amplifi-
1
cationbytheMRItoastrengthof1015Gcouldfurthermore
:
v be a scenario to explain the formation of the most strongly
The impact of rotation and magnetic field on the dynamics
i magnetized neutron stars known as magnetars (Woods &
X of core collapse supernovae is still uncertain and represents
Thompson 2006, and references therein). The birth of such
r a promising avenue of research. In the case of fast rotat-
a ing progenitors giving birth to a protoneutron star (PNS) magnetars with rotation periods in the millisecond range is
a potential central engine for long gamma-ray bursts (e.g.
rotating with a period of a few milliseconds, the rotational
Duncan & Thompson 1992; Metzger et al. 2011) associated
energyrepresentsanimportantenergyreservoirthatcanbe
with supernovae that have extreme kinetic energies (also
efficiently tapped if strong magnetic fields are present. In
called ”hypernovae” or type Ic BL, e.g. Drout et al. 2011).
these conditions, the magnetorotational instability (MRI)
Thedelayedenergyinjectionduetothespindownofafast
hasbeenproposedasamechanismtoefficientlyamplifythe
rotating, highly magnetized neutron star has also been in-
magneticfield(Akiyamaetal.2003).Itsdevelopmentcould
vokedasanexplanationforsomesuperluminoussupernovae
impact the explosion by converting shear energy into ther-
like SN 2008 bi (Kasen & Bildsten 2010; Woosley 2010;
malenergy(Thompsonetal.2005)orbygeneratingastrong
Dessart et al. 2012; Nicholl et al. 2013; Inserra et al. 2013).
large-scale magnetic field which can extract the rotational
energy and launch a magnetorotational explosion (LeBlanc SincetheseminalworkofBalbus&Hawley(1991),the
& Wilson 1970; Bisnovatyi-Kogan et al. 1976; Mu¨ller & MRI has been thoroughly investigated in the context of ac-
Hillebrandt 1979; Symbalisty 1984; Moiseenko et al. 2006; cretiondiscs.ThephysicalconditionsinthePNSdifferfrom
Shibataetal.2006;Burrowsetal.2007;Dessartetal.2008; those prevailing in accretion discs in several aspects. The
Takiwaki et al. 2009; Takiwaki & Kotake 2011; Winteler mainsupportagainstgravityisprovidedbypressurerather
(cid:13)c 0000RAS
2 Guilet & Mu¨ller
than centrifugal forces, and as a consequence the rotation Choosingaboxsizesmallerthanthedensityscale-heighten-
profile can be non-Keplerian. A strong differential rotation sures that this local approach is approximately justified.
appropriateforMRIgrowthnonethelessresultsfromthecol- The paper is organized as follows. In Section 2, we de-
lapseoftheprogenitorcoretoaPNS(Akiyamaetal.2003; scribethephysicalandnumericalsetupconsidered.There-
Ottetal.2006).ThedependenceofMRIturbulenceonthe sults are then presented in Section 3 for the linear expo-
steepnessoftherotationprofilewasstudiedbyMasadaetal. nential growth phase of the MRI, and in Section 4 for the
(2012).Anotherconsequenceofthehighpressureisthatra- saturated non-linear phase that follows. Finally we discuss
dial buoyancy forces driven by entropy and lepton fraction the validity of our assumptions in Section 5 and draw con-
gradients have an important dynamical role. This stands in clusions in Section 6. A linear analysis is also presented in
contrasttothecaseofthinaccretiondiscswhereonlyverti- Appendix A and B, the results of which have been used in
cal buoyancy is dynamically relevant. The impact of buoy- Sections 3 and 4. Finally, a convergence test of the simula-
ancy on the linear growth of the MRI has been studied by tions is described in Appendix C.
Balbus&Hawley(1994);Menouetal.(2004);Masadaetal.
(2006,2007).Theyshowedthatthediffusionofheatandlep-
tons driven by neutrinos alleviates the otherwise stabilising
2 NUMERICAL SETUP
effect of a stable stratification, such that the MRI can still
growsufficientlyfast.Neutrinoradiationalsoinducesahigh 2.1 Governing equations
viscosity,whichcanlimitthegrowthoftheMRIiftheinitial
magnetic field is too low (Guilet et al. 2015). The simulations performed in this article are designed to
represent a small part of a fast rotating PNS located in the
Howthenon-linearphaseoftheMRIisaffectedbythe equatorial plane at a radius of around 20km. We describe
radialbuoyancyandneutrinoradiationismuchlessknown. thelocaldynamicsintheframeworkofaCartesianshearing
TheimpactofbuoyancyontheMRIhasbeenstudiedinnu- box (e.g. Goldreich & Lynden-Bell 1965). The coordinates
mericalsimulationsof”semi-global”models(Obergaulinger x, y, and z represent the radial, azimuthal and vertical di-
et al. 2009; Masada et al. 2015), but these simulations did rection,respectively,andthecorrespondingunitvectorsare
nottakeintoaccounttheimpactofneutrinosasasourceof e , e , and e . The angular frequency vector points in z
x y z
thermal diffusion and viscosity (although it has an impor- direction,i.e.Ω=Ωe ,whilegravityisin−xdirection,i.e.
z
tantimpact,seeGuiletetal.2015).Thepioneeringsimula- g = −ge . At the location of the box considered here, the
x
tions of Obergaulinger et al. (2009) considered a small part neutrinosareinthediffusiveregime(Guiletetal.2015),and
of a PNS but were not purely local in that they included theireffectsonthedynamicscanappropriatelybedescribed
global gradients, in particular of entropy and density. They by a viscosity ν as well as thermal and lepton number dif-
investigated the impact of buoyancy driven by a radial en- fusivities.Forthesakeofsimplicity,weassumeinthiswork
tropy gradient, but artefacts induced by the radial bound- that the thermal and lepton number diffusivities are equal
aryconditionspreventedanadequatestudyoftheimpactof anddenotedbyχ,whichallowsustodescribethebuoyancy
buoyancyonthenon-linearphaseoftheMRI.Indeed,soon associated with both entropy and lepton number gradients
after the MRI reached non-linear amplitudes they observed with the use of a single buoyancy variable (see below). For
a flattening of the entropy gradient with a strong gradient thesakeofconciseness,werefertoχasthethermaldiffusiv-
onlypresentinathinlayerneartheradialboundaries.The ity. We also considered the effects of a magnetic diffusivity
radially global (but vertically local) model of Masada et al. η in our study.
(2015) may not suffer from such boundary artefacts (and Asdiscussedintheintroduction,buoyancyisdescribed
couldinprincipledescribemeaningfullytheevolutionofthe in the framework of the Boussinesq approximation, i.e. the
global radial profiles), but the radial and azimuthal reso- buoyancy force is the only compressibility effect taken into
lution that they could afford was in our opinion much too account. The Boussinesq approximation holds if both the
coarse for an accurate study of the MRI. flowandAlfv´envelocityaremuchlessthanthesoundspeed
(incompressible flow), and if the density perturbations as-
Inordertoavoidtheseshortcomingsandtoreducethe
sociated with the buoyancy are small. The approximation
computationalcost(thusallowinganexplorationofparam-
further assumes a uniform background density ρ , i.e. one
eter space with 3D simulations), we used a more simplified 0
neglects the radial density gradient. The validity of these
setup.WeperformedshearingboxsimulationsintheBoussi-
assumptions will be assessed in Section 5.1.
nesqapproximation.Thisapproximationtakesintoaccount
In the shearing box approximation, the diffusive MHD
the buoyancy effects due to entropy/lepton fraction gradi-
Boussinesqequationsreadinaframerotatingwiththecen-
ents, but neglects other compressibility effects. As will be
tre of the box with an angular frequency Ω
argued a posteriori, this assumption is very well justified
since magnetic and kinetic energies remain small compared 1 1
∂ v+v·∇v = − ∇Π+ (∇×B)×B, (1)
to the thermal energy. Our approach retains the effects of t ρ 4πρ
0 0
global gradients of angular frequency, entropy and lepton −2Ω×v+2qΩ2xe −N2θe +ν∆v
x x
fraction through shear and buoyancy, but is purely local
∂ B = ∇×(v×B)+η∆B, (2)
in the sense that the volume averages of these quantities t
are kept fixed and that the density gradient is neglected. It ∂tθ+v·∇θ = vx+χ∆θ, (3)
therefore retains less global effects than the ”semi-global” ∇·v = 0, (4)
simulationsbyObergaulingeretal.(2009).Thisseemsnec-
∇·B = 0, (5)
essary,however,tohavecleanshear-periodicboundarycon-
ditionswhichavoidnumericalartefactsatthegridboundary. wherev andB aretheflowvelocityandthemagneticfield,
(cid:13)c 0000RAS,MNRAS000,000–000
MRI with buoyancy 3
respectively. The buoyancy variable to a dimensionless Prandtl number (characterizing the rela-
tive importance of viscosity and thermal diffusion)
g δρ
θ≡− (6)
N2 ρ0 Pr≡ν/χ=0.1 and 1, (13)
describesthedensityperturbationδρduetothejointeffect
respectively. Hence, the dimensionless Peclet number (char-
ofentropyandleptonnumberperturbations,andhasthedi-
acterizingtherelativeimportanceofadvectionandthermal
mensionofalength.Thegradientofthepressureperturba-
diffusion) has a value of
tion∇Πisobtainedfromtheconstraintofadivergence-free
flow field (equation 4). The quantity L2Ω
P ≡ z =100 and 1000 (14)
e χ
q≡−dlogΩ/dlogr (7)
for the high and low diffusion case, respectively.
measures the local shear (we assume q = 1.25 in all sim-
Theresistivityη isexpectedtobemanyordersofmag-
ulations). The square of the Brunt-Va¨isa¨la¨ frequency N is
nitude smaller than the viscosity ν, but numerical con-
defined by
straints do not allow one to perform numerical simulations
g (cid:20)∂ρ(cid:12) dS ∂ρ (cid:12) dY (cid:21) for such values of the resistivity. In our study, we used a
N2 ≡−ρ ∂S(cid:12)(cid:12)P,Ye dr + ∂Ye(cid:12)(cid:12)P,S dre , (8) resistivity η = 2.5×109cm2s−1, which corresponds to a
magnetic Prandtl numbercharacterizingtherelativeimpor-
where ρ, S, Ye, and P are the density, entropy, electron tance of viscosity and resistivity of
fraction,andpressure,respectively.Astationarysolutionof
these equations is: Pm ≡ν/η=4. (15)
v = −qΩxe , (9) This much too small value is a compromise between the
0 y
computing time constraints and the wish to simulate flows
B = B e , (10)
0 0 z
in the high magnetic Prandtl regime appropriate to a PNS
θ0 = 0. (11) (P ∼ 1013 ; Thompson & Duncan 1993; Masada et al.
m
Weinitialisedthenumericalsimulationswiththisstationary 2007). Note that this regime is unusual since the magnetic
solution, to which we imposed random velocity and mag- Prandtlnumberisordinarilyverysmallforstars,mostparts
neticfieldperturbations(asdescribedinmoredetailinnext of accretion discs and liquid laboratory metals. The high
subsection). value of the magnetic Prandtl number in a PNS is due to
theneutrinosinducingaverylargeviscositywhiletheresis-
tivityremainsverysmall.Thedependenceoftheresultson
2.2 Physical parameters and dimensionless themagneticPrandtlnumberwillbeinvestigatedinaforth-
numbers coming publication. The corresponding magnetic Reynolds
number,whichcharacterizestherelativeimportanceofmag-
We chose the physical parameters of the simulations to
netic advection to magnetic diffusion (resistivity), is
represent the fast rotating PNS model studied by Guilet
et al. (2015) at a radius of (cid:39) 20km with a density ρ = L2Ω
0 R ≡ z =4×103 (16)
2×1013gcm−3, a viscosity ν = 1010cm2s−1, and a rota- m η
tion angular frequency Ω = 103s−1. The box size was set
in our study.
to(L ,L ,L )=(4,4,1)km.Theradialwidthoftheboxis
x y z We considered a single value of the magnetic field
somewhatsmallerthanthedensityscaleheight∼8km,i.e.
strength,whose valueB =1.6×1013G lies within thevis-
the neglect of the background density gradient is roughly 0
cousregimedescribedbyGuiletetal.(2015)(seeleftpanel
justified (see Sections 5.1 and 5.3 for a discussion). The as-
oftheirFig.10).Adimensionlessnumbercharacterizingthe
pect ratio of the box, elongated in the horizontal directions
strength of the magnetic field is defined as
compared to the vertical direction, was chosen to ensure
that the box fits the most unstable parasitic modes grow- L2Ω2
β ≡ z =104, (17)
ing on MRI channel modes (Goodman & Xu 1994; Pessah v2
A
& Goodman 2009; Latter et al. 2009; Obergaulinger et al. √
2009; Pessah 2010). where vA = B/ 4πρ0 is the Alfv´en speed associated with
Given these values for the box size and for the viscos- the background vertical magnetic field. The impact of vis-
ity, the dimensionless Reynolds number (characterizing the cosityonthelineargrowthoftheMRIisdeterminedbythe
importanceofviscosityoverthelargestverticalscaleinthe viscous Elsasser number
box) is v2 R
E ≡ A = e =0.1. (18)
L2Ω ν νΩ β
R ≡ z =103. (12)
e ν SinceE <1,oursimulationsbelongtotheviscousregime,
ν
where the growth rate of the MRI is significantly reduced
According to Masada et al. (2007), the thermal diffu-
by neutrino viscosity compared to ideal MRI (e.g. Guilet
sivity χ is larger than the neutrino viscosity ν by a factor
etal.2015,andreferencestherein).Similarly,theimpactof
10−100. To investigate its impact on the results we chose
resistivity on the linear growth of the MRI is quantified by
two different values for the thermal diffusivity, namely a
”high” one χ = 1011cm2s−1, and a ten times lower one the resistive Elsasser number
χ = 1010cm2s−1 referred to as ”low diffusion” in the fol- v2 R
E ≡ A = m =0.4. (19)
lowing. These values of the thermal diffusivity correspond η ηΩ β
(cid:13)c 0000RAS,MNRAS000,000–000
4 Guilet & Mu¨ller
Contrary to what would be expected under realistic condi- method using the 2/3 dealiasing rule. The time integration
tions,theimpactofresistivityonthelinearMRIistherefore is performed using an implicit procedure for the diffusive
not negligible in our setup (since E < 1). Finally, the im- terms, while other terms use an explicit 3rd order Runge
η
pactofthermaldiffusiononthelineargrowthoftheMRIis Kutta scheme. snoopy is parallelized using both MPI and
quantified by the thermal Elsasser number OpenMP techniques. snoopy has been used in the past to
study the MRI (in the absence of buoyancy) by Lesur &
v2 P
E ≡ A = e =0.01 and 0.1, (20) Longaretti(2007),Longaretti&Lesur(2010),Rempeletal.
χ χΩ β
(2010),andLesur&Longaretti(2011).TheBoussinesqim-
forthehighandlowdiffusioncase,respectively.Thisshows plementationhasbeenusedtostudyverticalconvectionand
thatthermaldiffusionisexpectedtoplayanimportantrole the subcritical baroclinic instability in an accretion disc by
in the MRI linear growth, as will be detailed in Section 3. Lesur & Papaloizou (2010) and Lesur & Ogilvie (2010).
Furthermore, we varied the Brunt-Va¨isa¨la¨ frequency Most simulations presented in this paper were per-
(see equation 8), which characterizes the impact of buoy- formed using a standard grid resolution of (n ,n ,n ) =
x y z
ancy, performing a set of simulations with N2/Ω2 ∈ (256,128,64) zones. A convergence study, presented in Ap-
{−1.5,−1,−0.5,0,1,2,3,4,6,8,10} both for the low and pendix C, shows that this resolution is sufficient to resolve
high thermal diffusion case. thedissipativescalesandthattheresultsareconsistentwith
In order to trigger the growth of the MRI, we added those obtained at a twice higher resolution. All our simu-
velocity and magnetic field perturbations to the stationary lations with N2/Ω2 (cid:62) −0.5 were performed with the stan-
stateinthefollowingway.WeimposedFouriermodeswitha dardresolution.AtmorenegativevaluesofN2,however,we
randomphaseandanamplitudethatfollowsaKolmogorov observed that the more vigorous turbulence led to smaller
spectrum.Thespectrumisnormalizedsuchthatthekinetic dissipative scales. Therefore, we increased the resolution to
and magnetic energy of the perturbations amount to 1% of (n ,n ,n )=(384,192,96)forN2/Ω2 =−1(bothforP =
x y z e
the energy in the background magnetic field. The results 100andP =1000),andto(n ,n ,n )=(512,256,128)for
e x y z
in the non-linear phase are not very sensitive to the initial the simulation with N2/Ω2 =−1.5 and P =100. We note
e
conditions,whichinfluencethe durationof the initialexpo- that in all our simulations the resolution is twice lower in
nential growth phase (higher amplitude perturbations lead the azimuthal direction than in the radial or vertical direc-
to a shorter growth phase). tion. This is rather usual in numerical studies of the MRI,
Note that the simulations presented in this paper, because the presence of shear makes structures more elon-
meant to describe the flow in a fast rotating neutron star gated in the azimuthal direction. We have checked that our
at a radius of 20km, can be rescaled to other physical sit- results are not affected by the lower azimuthal resolution.
uations sharing the same values of the dimensionless num-
bers.Forexample,themoderatelyrotatingPNSmodelstud-
ied by Guilet et al. (2015) shares the same dimensionless 2.5 Turbulent energies and energy injection rates
numbers at a radius of 10km with the following parame-
We define the kinetic, magnetic, and thermal energy densi-
ters: Ω = 200s−1, ν = 2×109cm2s−1, ρ = 1014gcm−3,
ties as
B = 7 × 1012G, box size (L ,L ,L ) = (4,4,1)km,
0 x y z 1
Prandtl numbers Pr=0.1−1, and magnetic Prandtl num- Ekin ≡ 2ρ0u2, (21)
berP =4.TheBrunt-Va¨isa¨la¨frequenciesconsideredthen
m
range from N2 =−6×104s−2 to 4×105s−2 in that situa-
1
tion. E ≡ B2, (22)
mag 8π
and
2.3 Units
1
E ≡ ρ N2θ2, (23)
Throughout the paper, we normalise our results using the th 2 0
verticalsizeofthedomainLz =1km,theangularfrequency respectively, where thevelocity uisdefinedwith respectto
Ω = 103s−1, and the density ρ0 = 2×1013gcm−3. Time the equilibrium shear profile v0 according to u ≡ v−v0.
is therefore measured in units of 1ms, velocity in units of Forbrevity,wewillrefertouasthevelocity,andtoenergy
108cms−1, the magnetic field in units of 1.6×1015G, the densities as energies.
energydensityinunitsof2×1029ergcm−3,thespecificen- WefurtherdenotetheaverageofaquantityAoverthe
ergyinunitsof1016erg/g,theenergyinjectionratedensity wholecomputationalvolumeby(cid:104)A(cid:105).Thesumofthevolume
in units of 2×1032ergs−1cm−3, and the specific energy averaged kinetic and magnetic energies obeys the following
injection rate in units of 1019ergs−1g−1. conservation equation
d(cid:104)E (cid:105)+(cid:104)E (cid:105)
kin mag =qΩ[(cid:104)R (cid:105)+(cid:104)M (cid:105)]+(cid:104)W(cid:105)−(cid:104)(cid:15) (cid:105)−(cid:104)(cid:15) (cid:105),
2.4 Numerical methods dt xy xy ν η
(24)
InordertosolvetheBoussinesqMHDequations(1)–(5),we
used the pseudo-spectral code snoopy (Lesur & Longaretti where the Reynolds stress
2005, 2007). snoopy solves the 3D shearing box equations
R ≡ρ u u (25)
using a spectral Fourier method, where the shear is han- xy 0 x y
dledthroughtheuseofashearingwavedecomposition(with and the Maxwell stress
timevaryingradialwavevector)andaperiodicremapproce-
−B B
dure.Nonlineartermsarecomputedwithapseudo-spectral Mxy ≡ 4xπ y (26)
(cid:13)c 0000RAS,MNRAS000,000–000
MRI with buoyancy 5
aretheangularmomentumfluxesintheradialdirectiondue instability criterion becomes 0 < κ2 < 4Ω2, i.e. the MRI
totheflowandthemagneticfield,respectively.Thesumof unstable region now comprises a region that is stable ac-
qΩ(cid:104)R (cid:105)andqΩ(cid:104)M (cid:105)givestheenergydensityinjectionrate cordingtoequation(30)(turquoiseregionintherightpanel
xy xy
into turbulent flow extracted from the shear energy of the of Fig. 1).
flow by Reynolds and Maxwell stresses. Anotherconsequenceofthermaldiffusionistheappear-
ance of a new hydrodynamic instability (Klahr & Hubbard
W ≡−N2ρ θu (27)
0 x 2014;Lyra2014),whichwerefertoasoscillatoryconvection.
istheworkdonebythebuoyancyforceperunitoftimeand It occurs for N2 < 0 in a parameter regime that is stable
per volume, i.e. (cid:104)W(cid:105) in equation 24 represents the energy according to equation (29) (orange region in Fig. 1). Oscil-
exchangeratefromthermaltokineticenergythroughradial latoryconvectioncanbeconsideredtobeananalogofsemi-
transportofentropyandleptonnumbercausedbybuoyancy. convection(inwhichthermaldiffusioncausestheflowtobe
The remaining quantities in equation 24, (cid:15) ≡ρ ν|∇×u|2 unstable in the presence of an otherwise stabilising compo-
ν 0
and(cid:15) ≡η|∇×B|2/4π,aretheenergydissipationratesdue sitiongradient),rotationassumingtheroleofthestabilising
η
to the viscosity and the resistivity, respectively. composition gradient in the case of semi-convection.
In a time average sense, the total energy injection rate Theparameterregimeexploredwithoursimulationsis
into turbulence, qΩ[(cid:104)R (cid:105)+(cid:104)M (cid:105)]+(cid:104)W(cid:105), is equal to the shownwithplussignsinFig.1.WevariedtheBrunt-Va¨isa¨la¨
xy xy
sumoftheviscousandresistiveenergydissipationrates(cid:104)(cid:15) + frequency in the range −1.5 < N2/Ω2 < 10 and fixed the
ν
(cid:15) (cid:105). epicyclic frequency to a value κ2/Ω2 = 1.5. Hence, we ex-
η
plored three regimes: non-diffusive MRI for 0 < N2 < 2.5,
diffusive MRI (i.e. allowed by thermal diffusion) for N2 >
2.5,andamixedregimewhereMRIandoscillatoryconvec-
3 LINEAR GROWTH PHASE
tion can co-exist for −1.5 < N2 < 0. Note that in none
Letusstartthissectionbydiscussingstabilityconditionsin of our simulations the flow is unstable to axisymmetric
the presence of buoyancy and differential rotation both in classical (non-oscillatory) convection.
hydrodynamic and MHD flows. We restrict this discussion Figure2showsthetimeevolutionoftheMaxwellstress
tothespecialcaseconsideredinthispaper:alocalboxinthe averaged over the computational box. After a linear phase
equatorial plane of the PNS, where the angular frequency, of exponential growth, the Maxwell stress saturates in the
entropyandcompositiondependonlyonradius.Wefurther non-linear phase. In this section we analyse only the linear
restrict the discussion to axisymmetric perturbations1 and phase and postpone the analysis of the non-linear phase to
to flows which are stable against the hydrodynamic shear nextsection.Allbuttwosimulationsshowanon-oscillatory
instability, i.e. complying with the Rayleigh stability crite- growth,whichisduetotheMRImodifiedbybuoyancy.For
rion κ2 > 0, where κ is the epicyclic frequency defined by low thermal diffusion (P = 1000; right panel), the growth
e
rate varies by orders of magnitudes and becomes extremely
1 d(r4Ω2) smallforlargepositivevaluesofN2,uptothepointthatno
κ2 ≡ r3 dr . (28) growth occurs at all for N2/Ω2 =10. By contrast, for high
thermal diffusion (P =100; left panel), the effect of buoy-
We first consider stability conditions in the absence of e
ancy on the growth of the MRI is much less pronounced.
any diffusive process. If
TheMRIcanstillgrowinashorttimeevenforlargevalues
κ2+N2 <0 (29) of N2. A good agreement is obtained between the growth
rates obtained from our numerical simulations and the pre-
the flow is convectively unstable because the Solberg-
dictions from the linear analysis presented in Appendix A
Høiland stability criteria is violated (red region in Fig. 1).
(Fig. 3).
In the presence of magnetic field the condition reads
Our results are in qualitative agreement with those
κ2+N2 <4Ω2. (30) of Masada et al. (2007), who showed that the MRI can
grow unimpeded by stable buoyancy if the thermal diffu-
This criterion is less restrictive than the hydrodynamic one
sivity is sufficiently large compared to the viscosity, i.e. if
due to the appearance of the MRI for 0 < κ2+N2 < 4Ω2
χ > νN/Ω. This condition is indeed verified for high ther-
(dark blue region, and red and orange regions with blue
mal diffusion (P = 100 and χ/ν = 10), in which case
e
stripes in the right panel of Fig. 1). The classical MRI (i.e.
the growth of the MRI is only mildly affected by buoy-
in the absence of buoyancy effect) is obtained for N2 = 0
ancy.Bycontrast,thegrowthoftheMRIisdrasticallysup-
and 0<κ2 <4Ω2. pressedinthecaseoflowthermaldiffusionandlargeN2/Ω2
When thermal diffusion is present (but no resistivity
(P =1000andχ=ν),i.e.iftheaboveconditionisnotful-
e
and viscosity), the instability region is more extended for
filled.Notethat,evenforahighdiffusivity,theMRIgrowth
two reasons. First, thermal diffusion alleviates the stabilis-
rates are significantly smaller than those of the ideal MRI
ing effect of buoyancy on the MRI (Acheson 1978; Menou
(σ/Ω = q/2 = 0.625). This stands in contrast to the situ-
etal.2004;Masadaetal.2007).AsaconsequencetheMRI
ation considered by Masada et al. (2007) and is due to the
fact that we considered a magnetic field strength that lies
intheviscousregime,wheretheMRIgrowthissignificantly
1 Due to the presence of shear, non-axisymmetric perturbations
slowed down by neutrino viscosity (Guilet et al. 2015).
cannot lead to exponential growth but can nonetheless be tran-
siently growing. While this may not be considered as leading to Oneofoursimulationswithunstablebuoyancyandhigh
genuine linear instability, it can be important for the non-linear thermal diffusion (N2/Ω2 =−1.5 and Pe =100) showed an
turbulentdynamics. oscillatory exponential growth due to the oscillatory con-
(cid:13)c 0000RAS,MNRAS000,000–000
6 Guilet & Mu¨ller
Hydrodynamical stability diagram MHD stability diagram
4 4
oscillatory oscillatory
convection convection
3 3 diffusive MRI
/MRI
2 2
Ω Ω
2 stable 2
2/ 2/
κ κ
1 1
convection convection MRI
/MRI
0 0
−4 −2 0 2 4 6 −4 −2 0 2 4 6
N2/Ω2 N2/Ω2
Figure 1. Stability diagram as a function of the square of the Brunt-Va¨is¨al¨a frequency N2 (characterizing the influence of buoyancy)
and the square of the epicyclic frequency κ2 (characterizing the effect of rotation and shear). The left panel applies to hydrodynamic
flowsandshowstheregionsunstableagainstclassicalconvection(red)andoscillatoryconvectionenabledbythermaldiffusion(orange).
The right panel applies to MHD flows and shows the region unstable against the MRI (blue) including the also convectively unstable
regions (red and orange regions with blue stripes) and the diffusive MRI enabled by thermal diffusion (turquoise). The plus signs give
theparametersofoursimulations(thoseforN2/Ω2=8and10arenotshown).
high diffusion low diffusion
100 100
−2
ss 10−2 ss 10
e e
str str −4 2
well 10−4 N22=−1.5 well 10 NN2==−01
Max N2=0 Max 10−6 N2=4
N =4
10−6 2 N2=10
N =10 −8
10

0 100 200 300 0 200 400 600 800 1000
time (ms) time (ms)
Figure2.TimeevolutionoftheboxaveragedMaxwellstressforsimulationswithdifferentstratifications,i.e.differentvaluesofN2/Ω2,
for high (Pe = 100; left panel) and low (Pe = 1000; right panel) thermal diffusion. Note that the time interval is different in the two
panels.
vection discussed above. The transition of the fastest grow- Obergaulinger et al. 2009; Pessah 2010). This leads to the
ingmodefromMRItooscillatoryconvectionisclearlyseen disruption of the channel modes, thereby ending the expo-
in Fig. 3 for high thermal diffusion, where it happens at nentialgrowthphase.Fromthatpointonthedynamicscon-
N2 ∼ −Ω2 in agreement with the linear prediction. In sists, in most cases, of 3D statistically steady MHD turbu-
the regime of oscillatory convection, the growth rate pre- lence.Asomewhatdifferentdynamicstakesplace,however,
dictedbythehydrodynamiclinearanalysispresentedinAp- in the case of low thermal diffusion (P =1000) and a suf-
e
pendix B matches that of the MHD linear analysis and the ficiently stabilising buoyancy N2/Ω2 (cid:38) 3. In that case, we
simulations,clearlydemonstratingthehydrodynamicnature observed the recurrent emergence of exponentially growing
of the instability. channel flows and their subsequent disruption by parasitic
instabilities. The resulting dynamics is fundamentally non-
steadyduetothisquasi-periodicemergenceofthelargescale
channelflows.Forhighthermaldiffusionwedidnotobserve
4 NON-LINEAR PHASE
recurrent channel flows, but nonetheless noticed the devel-
When the exponentially growing axisymmetric (i.e. y- opment of long-lived coherent flow structures in the regime
independent) MRI channel modes reach non-linear ampli- where the flow is stable against buoyancy. Therefore, we
tudes, non-axisymmetric secondary (or ”parasitic”) insta- split the following discussion of the properties of the non-
bilities start to grow upon the primary modes (Goodman linearphaseintotwosubsections.Thefirstonewedevoteto
& Xu 1994; Pessah & Goodman 2009; Latter et al. 2009; the coherent flows (either recurrent or statistically steady),
(cid:13)c 0000RAS,MNRAS000,000–000
MRI with buoyancy 7
0.25 1.000
0.20
e e
t t 0.100
a a
r 0.15 r

h h
t t
w 0.10 w 0.010
o o
r r
G 0.05 G
0.001
0.00
-2 0 2 4 6 8 10 12 -2 0 2 4 6 8 10 12
N2 N2
Figure3.ComparisonofthegrowthratesasafunctionofthestrengthofbuoyancyasmeasuredbyN2/Ω2obtainedfromournumerical
simulations(diamonds)andourlinearanalysis(solidlines;seeAppendixAandB)forhigh(Pe=100;leftpanel)andlow(Pe=1000,
rightpanel)thermaldiffusion.Thedottedlinegivesthepredictionofthelineargrowthrateinahydrodynamicflow(AppendixB).Blue
colour is used when the growth rate is negative, i.e. when it gives the damping rate. We normalized both the growth rate σ and the
Brunt-V¨ais¨alafrequencyN bytheangularfrequencyΩ.
varying but horizontally uniform structure (i.e. flow quan-


tities only depend on coordinate z), and zonal flows with

a radially varying but azimuthally and vertically uniform
structure(i.e.flowquantitiesonlydependoncoordinatex).

Channel flows are related to the most unstable lin-
ear modes of the MRI in case of a vertical background
 magnetic field. These modes, which have a purely vertical
wavenumber, i.e. k =(0,0,k ), are known to be non-linear
z
growing solutions of the MHD equations (Goodman & Xu

 1994).TheyareunstabletoparasiticinstabilitiesofKelvin-
Helmholtz or resistive tearing mode type (Goodman & Xu

1994; Latter et al. 2009; Pessah & Goodman 2009; Pessah
 2010; Obergaulinger et al. 2009), the development of which

canstopthegrowthofthechannelmodeandtriggerturbu-
lence. The typical structure of MRI channel modes consists
of a sinusoidal profile that holds for the horizontal (i.e. x
and y) components of both velocity and magnetic field 2,
and in the presence of buoyancy for the buoyancy variable
θ (equation 6), too. As can be deduced from Eqs. (A16)–
(A18), in MRI channel modes the radial (x) and azimuthal
 (y) velocity components are typically in phase 3, while the
 radial and azimuthal magnetic field components are typi-
cally in anti-phase with each other 4. They are shifted also
by a quarter of a wavelength with respect to the horizontal
Figure 4. Schematic representation of the shearing box geom- (x,y)velocitycomponents,whereasthebuoyancyvariableθ
etry (upper panel) and of the geometry of the coherent flows
is in phase with these components (equation A19).
(lower panels): channel flows (left-hand lower panel) and zonal
Zonalflowshavebeenobservedinshearingboxsimula-
flows(right-handlowerpanel).
while we discuss in the second subsection the time and vol-
2 Inashearingboxwithshear-periodicboundaryconditions,the
ume averaged properties of the flow during the non-linear
horizontallyaveragedvertical(z)componentofbothvelocityand
phase. magnetic field is conserved and uniform due to the divergence
free condition (Eqs. 4 and 5), and has a value of zero and B0
respectively.
4.1 Coherent flows 3 Intheregimewheretheflowisunstabletobuoyancy,theycan
bealsoinanti-phase,i.e.theassociatedReynoldsstressbecomes
In the stable buoyancy regime (N2 > 0), we observe the
negative(seesection4.2andAppendixA).
appearance of either long lived (for Pe = 100) or recurrent 4 In the regime where the flow is unstable to buoyancy, they
(forPe =1000)coherentstructuresinthenon-linearstateof canbeinprincipleinphasetoo,i.e.theMaxwellstresswouldbe
theMRI.Twodifferenttypesofsuchcoherentstructurescan negative.Thiswas,however,notobservedfortheflowparameters
be distinguished (Figure 4): channel flows with a vertically consideredinthisstudy.
(cid:13)c 0000RAS,MNRAS000,000–000
8 Guilet & Mu¨ller
P =100 P =1000
e e
1 1
axisymmetric
y y
g 0.8 g 0.8
er channel er
en 0.6 zonal en 0.6
c c
eti 0.4 eti 0.4
n n
g g
a 0.2 a 0.2
m m
0 0
1 1
y 0.8 y 0.8
g g
er er
n 0.6 n 0.6
e e
c c
eti 0.4 eti 0.4
n n
ki 0.2 ki 0.2
0 0
1 1
y 0.8 y 0.8
g g
r r
e e
n 0.6 n 0.6
e e
al al
m 0.4 m 0.4
r r
e e
h 0.2 h 0.2
t t
0 0
−2 0 2 4 6 8 10 −2 0 2 4 6 8 10
2 2
N N
Figure 5.Fractionofmagnetic(toprow),kinetic(middlerow),andthermal(bottomrow)energydensitiescontainedinaxisymmetric
structures(black+symbols),channelflows(red×symbols),andzonalflows(greencircles)forsimulationsperformedwithhighthermal
diffusion(Pe=100;leftcolumn)andlowthermaldiffusion(Pe=1000;rightcolumn).
tions of MRI turbulent stratified accretion discs (Johansen Figure 5 demonstrates that the flow becomes more ax-
et al. 2009; Simon et al. 2012; Simon & Armitage 2014). isymmetric for increasing N2, i.e. for flows that are more
They mostly possess a radially varying azimuthal velocity, stable to buoyancy. This holds for the magnetic, kinetic,
which is associated with a radial variation of the pressure. andthermalenergiesatbothlowandhighthermaldiffusion
In the presence of a mean vertical magnetic field, a radial but is more pronounced at low thermal diffusion, in quali-
variation of the vertical magnetic field has been reported tative agreement with the idea that high thermal diffusion
recently too (Bai 2014; Bai & Stone 2014). somewhat alleviates the influence of a stabilising buoyancy.
The prominence of channel flows also increases monoton-
Figure5illustratestheoccurrenceofthedifferenttypes ically with N2 for all three energies and at both thermal
of coherent flows as a function of the Brunt-V¨aisa¨la¨ fre-
diffusivities. This behaviour is again more pronounced at
quency by showing the fractions of kinetic, magnetic, and
lowthermaldiffusion,wherethechannelflowscontainupto
thermal energy densities contained in axisymmetric struc-
about 50% of the energy. As already mentioned the chan-
tures,channelflows,andzonalflows,respectively.Toobtain
nel flows have a different time behaviour at P = 100 and
e
the energy densities contained in axisymmetric structures
P =1000.Athighthermaldiffusiontheyareapproximately
e
wefirstaveragedthemagneticfield,velocity,andbuoyancy
steady (see Section 4.1.1), while at low thermal diffusion
variable in azimuthal (y) direction. Then we used these av-
they show recurrent phases of exponential growth and dis-
eragestocomputetherespectiveenergydensitiesaveraging
ruption (see Section 4.1.3).
over the box and over time. Instead, we first averaged the
flow variables horizontally (i.e. both in azimuthal (y) and Zonalflowsstickoutmostclearlyregardingthekinetic
radial (x) direction) to compute the energy densities con- energy (remarkably, their fraction can reach up to 90%),
tainedinchannelflows.Finally,forthezonalflows,wefirst whiletheyamounttoalesserfractionofthethermalenergy,
averagedtheflowvariablesinvertical(z)andazimuthal(y) and an almost negligible fraction of the magnetic energy.
direction. Their dependence on N2 is more complex than that of the
(cid:13)c 0000RAS,MNRAS000,000–000
MRI with buoyancy 9
channelflows.WhiletheirimportancefirstincreaseswithN2
(up to N2 =4−6 at P =100, and N2 =2 at P =1000),
e e
it starts to decrease at still higher N2. As for the channel
flows, zonal flows have a different time behaviour at high
andlowdiffusion.Theyareapproximatelysteadyandlong-
livedatP =100(Section4.1.2),butrecurrentatP =1000
e e
(Section 4.1.3).
4.1.1 Quasi-stationary channel flows
Figure 6 shows snapshots of the azimuthal magnetic field,
the azimuthal velocity, and the buoyancy variable for two
simulations at high thermal diffusion for a flow unstable
(N2/Ω2 = −1; left column) and stable (N2/Ω2 = 10; right
column) to buoyancy. In the latter case, the flow is much
more axisymmetric, and we observe a channel flow struc-
ture superposed to a disordered flow.
The channel flow structure and time-evolution can be
made more apparent by calculating horizontal averages of
the flow variables. We show space-time diagrams of these
averages in Fig. 7 for the simulation with N2/Ω2 =10 and
P =100. A first transient phase of strong channel flow oc-
e
curs between t ∼ 100ms and 200ms corresponding to the
MRIchannelsof theexponentialgrowthphasedescribedin
Section 3. After the disruption of these channels by para-
sitic instabilities and the onset of turbulence, one observes
no clear channel flow until t ∼ 300ms. At that moment,
anotherchannelflowappearswhoseamplitudeandphaseis
remarkably steady until the end of the simulation, except
foraphaseshiftthattakesplacearoundt∼750ms.Toour
knowledge,thisisthefirsttimethatsuchastationarychan-
nel flow has been seen to appear in MRI turbulence (only
recurrent channel flows were sometimes observed as will be
discussed in Section 4.1.3).
Let us now compare the structure of this stationary
channel flow to that of growing linear MRI modes. From time averaged profiles
thespace-timediagrams(Fig.7)itisalreadyclearthatthe 0.2
wavelength of the stationary flow (one wavelength fits into v *10
x
the box) is longer than that of the exponentially growing
0.1
mode (a mixture of modes with n = 2 and n = 3 wave- v
y
lengthsfitintothebox),thelatteronecorrespondingtothe
linear prediction 5. Obviously, the quasi-stationary chan- 0 Bx*10
nel flow does not correspond to the most unstable MRI
B
mode. The time-averaged profiles (Fig. 7; bottom panel) −0.1 y
shownonethelessinterestingqualitativesimilaritieswithlin-
θ*10
ear MRI modes. The profiles resemble quite closely that of
a sinusoidal function whose wavelength is equal to the ver- −0.2
tical size of the box, except for the profiles of the radial −0.5 0 0.5
velocity and the magnetic field which both show some in- z
dication of smaller scale structure on top of the dominant
n = 1 sinusoid. Interestingly, the radial velocity, azimuthal Figure 7. Quasi-stationary channel flow for N2/Ω2 = 10 and
velocity,andbuoyancyvariableareinphase,whiletheradial high thermal diffusion Pe = 100. The upper five panels show
space time diagrams (from top to bottom) of the horizontally
and azimuthal magnetic fields are in anti-phase and shifted
averagedradialvelocity,azimuthalvelocity,radialmagneticfield,
by a quarter of a wavelength with respect to the velocity
azimuthal magnetic field, and buoyancy variable θ. The bottom
profiles. This result is in accordance with the structure of
panelgivestheverticalprofilesofthesequantitiesaveragedover
thetimeinterval400ms(cid:54)t(cid:54)700ms.
5 Thelinearanalysispredictsthatthen=3modeisthefastest
growing mode (σ/Ω = 0.066), while the n = 2 mode is the sec-
ondfastest growingone witha growthrate onlyslightlysmaller
(σ/Ω=0.061).Twofurthermodes(n=1,4)areunstablewitha
roughlytwicesmallergrowthrate.
(cid:13)c 0000RAS,MNRAS000,000–000
10 Guilet & Mu¨ller
Figure 6. Three-dimensional snapshots at time t = 553ms of two simulations with Pe = 100 for a buoyantly unstable flow with
N2/Ω2 =−1 (left column) and buoyantly stable flow with N2/Ω2 =10 (right column). The three rows show the azimuthal magnetic
field (top), azimuthal velocity (middle), and buoyancy variable θ (bottom). Note that we adjusted the colour scales to the level of
turbulence,i.e.theyaredifferentforthetwocolumns.
a linear MRI channel mode, and it suggests a link between 4.1.2 Quasi-stationary zonal flows
the stationary channel flow and linear MRI modes.
Zonalflowscanconvenientlybevisualisedwithaspace-time
diagramshowingtheradialprofileofazimuthallyandverti-
callyaveragedquantities.Fig.8showssuchadiagramfora
simulationwithN2/Ω2 =4andP =100,whereaverysta-
e
In view of the qualitative similarity but quantitative ble quasi-stationary zonal flow is clearly visible. It shows
differencebetweenthestationarychannelflowandthelinear that azimuthal velocity has by far the largest amplitude
MRI modes, we interpret the stationary channel flow as an amountingto∼90%ofthekineticenergy,butamodulation
MRIchannelmodesignificantlymodifiedbytheconcurrent withawavelengthequaltotheradialsizeoftheboxisalso
turbulence,whichpresumablybringsittomarginalstability. recognisableintheverticalmagneticfieldandthebuoyancy
Furtherstudiesofthephysicsofthestationarychannelflow variable. The modulation of the vertical magnetic field is
and its interaction with turbulence may shed light on the shiftedbyaquarterofawavelengthwithrespecttothatof
mechanismdeterminingthelevelofMHDturbulence.Using azimuthal velocity (i.e. the maximum of the magnetic field
enlargeddiffusioncoefficientstoroughlydescribetheimpact correspondstoazeroofthevelocity),inagreementwiththe
of turbulence on the channel flow may provide the basis results of Bai & Stone (2014). We note that the magnetic
for modelling MRI saturation and could allow for a deeper fieldzonalmodulationrepresentsanegligiblefractionofthe
interpretation of the numerical simulations. One could, for total magnetic energy (see Fig. 5), but it amounts to a sig-
example,speculatethatthelevelofturbulenceadjustssoas nificant fraction of the mean vertical magnetic field. Thus,
to bring the largest scale MRI mode to marginal stability it may play a significant role in the flow dynamics.
(Ogilvie & Livio 2001; Guilet & Ogilvie 2012, 2013). ThezonalflowshowninFig.8isrepresentativeofflows
(cid:13)c 0000RAS,MNRAS000,000–000