The Excitation, Propagation and Dissipation of Waves in Accretion Discs: The Non-linear Axisymmetric Case

36 Pages·1.1 MB·English

by  M. R. Bate

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed1February2008 (MNLaTEXstylefilev2.2)
The excitation, propagation and dissipation of waves in
accretion discs: the non-linear axisymmetric case
M. R. Bate,1,2⋆ G. I. Ogilvie,1 S. H. Lubow,1,3 and J. E. Pringle1,3
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA
2 2School of Physics, Universityof Exeter, StockerRoad, ExeterEX4 4QL
0 3Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
0
2

n 1February2008
a
J

0 ABSTRACT
1 We analyse the non-linear propagation and dissipation of axisymmetric waves in ac-
cretion discs using the ZEUS-2D hydrodynamics code. The waves are numerically
1
resolvedintheverticalandradialdirections.Bothverticallyisothermalandthermally
v
stratified accretion discs are considered. The waves are generated by means of reso-
9
nant forcing and several forms of forcing are considered. Compressional motions are
4
takento be locally adiabatic (γ =5/3). Priorto non-linear dissipation, the numerical
1
results are in excellent agreement with the linear theory of wave channelling in pre-
1
0 dicting the types of modes that are excited, the energy flux by carried by each mode,
2 and the vertical wave energy distribution as a function of radius. In all cases, waves
0 are excited that propagate on both sides of the resonance (inwards and outwards).
/ Forverticallyisothermaldiscs,non-lineardissipationoccursprimarilythroughshocks
h
that result from the classical steepening of acoustic waves. For discs that are sub-
p
- stantially thermally stratified, wave channelling is the primary mechanism for shock
o generation. Wave channelling boosts the Mach number of the wave by vertically con-
r fining the wave to a small cool region at the base of the disc atmosphere. In general,
t
s outwardly propagating waves with Mach numbers near resonance Mr>∼0.01 undergo
a
shocks within a distance of order the resonance radius.
:
v
i Key words: accretion, accretiondiscs – binaries: general– hydrodynamics – waves.
X
r
a
1 INTRODUCTION Initial studies of wave propagation approximated the
disc as two-dimensional and ignored the effects of vertical
Gaseous discs play a significant role in the evolution of bi-
structure (perpendicular to the disc plane). Goldreich &
nary starand planetary systems.Generally, thetidal forces
Tremaine (1979) developed atwo-dimensional linear theory
from stellar or planetary companions distort the discs from
for resonant tidal torques and the associated wave propa-
an axisymmetric form. However, where resonances occur in
gation. Later studies of the non-linear case found that the
the discs, the tidal forces also generate waves that trans-
torques were within a few percent of those predicted by
port energy and angular momentum. The resulting reso-
the two-dimensional linear formula (Shu, Yuan & Lissauer
nant torques cause orbital evolution of the perturbing ob-
1985; Yuan & Cassen 1994). For a thin disc that is verti-
jects (e.g. Goldreich & Tremaine 1980; Lin & Papaloizou
cally isothermal and has an isothermal thermodynamic re-
1993;Lubow&Artymowicz1996).Thewavesalsocausethe
sponse(γ =1),theonlywaveexcitedinthedischasatwo-
discs to evolve since they transfer their energy and angular
dimensional structure.The wavefronts are perpendicularto
momentum to the disc in the locations where they damp.
the disc plane and the motion is purely horizontal (paral-
Damping may be provided by shocks, radiative damping
lel to the disc plane). Thus, a two-dimensional treatment
(Cassen & Woolum 1996), or other viscous damping mech-
is valid in a thin disc if both the vertical structure of the
anisms.
discanditsthermodynamicresponsearelocallyisothermal.
Wavesindiscshavealsobeenconsideredaspossibleex-
However,thisisnotrealisticformostdiscs,suchasthosein
planationsofquasi-periodicvariabilityintheX-rayemission
cataclysmic variables, circumstellar and circumbinary discs
ofsystemsinvolvingaccretion discsaroundblackholes(e.g.
around pre-main-sequence stars, and protoplanetary discs.
Kato & Fukue1980; Nowak & Wagoner 1991; Kato 2001).
In such discs, the waves are no longer two-dimensional and
theverticalstructureofthediscmustbetakenintoaccount.
⋆ E-mail:[email protected] A study of three-dimensional wave propagation was
2 M. R. Bate et al.
made using linearised numerical simulations (Lin, Pa- structure of the equilibrium discs that we model and the
paloizou & Savonije 1990a,b). Unfortunately, the limited grid layout for the ZEUS-2D calculations. Section 4 briefly
range of physical parameter space that was covered led to reviews the semi-analytical linear method for solving the
theinterpretationthatwavesinaverticallythermallystrat- three-dimensional wave propagation problem, describes our
ified disc are refracted upwards into the atmosphere where method of wave excitation, and discusses the requirements
they dissipate via shocks. Later, it was realised that the forconvergenceofthenumericalcalculations.Sections5and
linear problem could be solved semi-analytically (Lubow & 6 contain the results for a vertically isothermal disc, and a
Pringle1993;Korycansky&Pringle1995;Lubow&Ogilvie polytropic disc with a vertically isothermal atmosphere, re-
1998; Ogilvie& Lubow1999). Thesestudiesfound that the spectively.Intermediatediscswith optical depthsnot much
propagation of waves in a thin disc is similar to the propa- largerthanunityareconsideredinSection7.Asummaryof
gationofelectromagneticwavesalongawaveguide(Jackson theresults is contained in Section 8.
1962; Landau&Lifshitz1960).Themotionsinthedisccan
bedescribedintermsofmodes.Foreachmode,thevertical
wave structure is determined in detail at each disc radius. 2 NUMERICAL METHOD
Thehorizontalvariationsofthemodearetreatedbymeans
of a WKB approximation in which the horizontal (radial) The non-linear hydrodynamic calculations presented here
wavenumberis determined as a function of wave frequency. wereperformedwiththeZEUS-2DcodedevelopedbyStone
A flux conservation law determines the mode amplitude as & Norman (1992). No modifications to the code were re-
a function of radius. Dependingon thevertical structure of quired,excepttheadditionofdampingboundaryconditions
the disc and the thermodynamic response of the gas, the for test purposes only (see below).
wave energy may be channelled towards the surface of the The ZEUS-2D code comes with various options for the
disc,towardsthemid-planeofthedisc,oroccupytheentire interpolation scheme and artificial viscosity. We use the
vertical extent as it propagates in radius. van Leer (second-order) interpolation scheme. Linear and
quadraticartificialviscositiesareprovidedinthecode,with
Ideally,onewouldliketodeterminethenon-linear,non-
twoformsofthequadraticviscosity:thestandardvonNeu-
axisymmetric response of a variety of disc models (verti-
mann &Richtmyer(1950) form, and atensor form. Weuse
cally isothermal and thermally stratified) to resonant forc-
onlyaquadraticartificialviscositytohandleshocks.Calcu-
ing. Unfortunately, this is a formidable numerical problem
lationswereperformedbothwiththevonNeumann&Richt-
at present. Instead we begin, in this paper, with an anal-
myerform and with thetensor form. Wefindno significant
ysis of the non-linear axisymmetric response. The axisym-
difference in the results. All the results presented here use
metric case offers the major advantage of being expressed
thestandardvonNeumann&Richtmyerform.Thestrength
as a two-dimensional problem, which is currently accessible
numerically. Although axisymmetric waves only transport ofthequadraticviscosityiscontrolled byaparameterqvisc.
energy and not angular momentum, they should resemble Generally,weuseqvisc =4,butwepresentsomeresultsfrom
theresponseofthedisctolow-azimuthalnumbertidalforc- calculations that were performed with qvisc =0.5. No other
type of viscosity (e.g. one that includes a shear stress) is
ing, as arises in binary star systems. The local physics of
used.
the waveguide is determined within a region whose size is
All of the calculations presented here use reflecting
a few times the semi-thickness of the disc. On this scale a
boundaries. In order to verify that our results were not af-
non-axisymmetricwaveoflow azimuthal wavenumberisal-
fectedbyreflectionsfrom theboundaries,somecalculations
mostindistinguishablefromanaxisymmetricwave.Thisex-
were performed with damping boundaries at the inner and
plainswhythelocalphysicsofthewaveguidecanbestudied
outer radii of the numerical grid. Damping is provided by
within theshearing-sheet approximation (Lubow & Pringle
addingan acceleration
1993) and why the azimuthal wavenumber appears in the
dispersionrelationonlythroughaDopplershiftofthewave ∂v
= 2ωv (1)
frequency. ∂t −
Consequently, in this paper, we study the non-linear totheinnerorouter10 20radialzonesofthegrid,wherev
propagation and dissipation of axisymmetric, but fully re- −
isthevelocityofthegasandωisthefrequencyoftheforcing
solved,wavesinaccretiondiscs.Thepurposeofthepaperis
that is used to generate the waves (see Section 4). There is
threefold. First, we wish to determine, using non-linear hy-
no significant difference between the results obtained using
drodynamiccalculations,whetherornotthesemi-analytical
reflecting and damping boundaries.
lineartheoryprovidesanaccuratedescriptionofwaveprop-
agation in accretion discs. Secondly, although linear theory
can predict where non-linearity is likely to become impor-
3 MODELLING THE ACCRETION DISCS
tant, shocks and the process of energy deposition (and an-
gular momentum deposition with non-axisymmetric waves) Weconsiderthepropagation anddissipation of axisymmet-
cannot be handled. We aim to determine how and where ric (m = 0) waves in two types of accretion disc: a locally
dissipation occurs through shocks and how accurately this vertically isothermal disc, and a polytropic disc (polytropic
can be predicted from the linear analysis. Finally, we wish index n=1.5) that joins smoothly on to an isothermal at-
todeterminehowwellnon-linearhydrodynamiccalculations mosphere.Weignorealleffectsofdiscturbulence.Afullde-
can model theproblem and what resolution is required. scription of theequilibrium structuresis given in Appendix
The outlineof thispaperis as follows. InSection 2, we A; we only briefly list their properties here. The disc is de-
brieflyreviewtheZEUS-2Dhydrodynamiccodethatisused scribed in cylindrical coordinates (R,φ,z). As described in
toobtainthenon-linearresults.InSection3,wediscussthe Appendix A, we shall work in dimensional units, such that
Waves in accretion discs 3
the disc angular velocity is unity when the dimensionless zones per vertical scale-height and thesame radial spacing,
disc radius R = 1. The discs are Keplerian, with angular but each zone is either stretched or compressed vertically.
velocity Ω=R−3/2, apart from small corrections. In practice, calculations of vertically isothermal discs
In both discs, the value of the polytropic exponent wereperformedwithNr Nθ zonesof:250 111,500 221,
× × ×
thatdescribestheadiabaticpressure-densityrelationforthe or1000 441. Calculations ofthestandardpolytropicdiscs
×
waves is γ = 5/3. For convenience in modelling, we choose were computed with: 250 45, 500 91, 1000 183, or
× × ×
discs with finite inner and outer radii, R1 and R2. We also 1156 365 zones(thelast of which hasaresolution equiva-
×
choose a scale-height that increases linearly with radius, lentto2000 365zonesbutmodelsasmallerrangeofradii
which requires that the sound speed varies as c R−1/2. than the low×er-resolution calculations). The highest resolu-
∝
The scale-height of the vertically isothermal disc is H R. tion calculations took up to 6000 CPU hours on 440 MHz
∝
Thepolytropicdischasawell-definedpolytropiclayerwith Sun Workstations with each factor of 2 increase in resolu-
asemi-thicknessHp,thatsmoothlychangesintoanisother- tion taking a factor of 8 longer (a factor of 4 from the
≈
malatmospherewithitsownscale-heightHi,abovez Hp. increase in the number of zones and a factor of 2 from the
≈
Both Hp and Hi are proportional toR.From thispoint on, decrease in the time step). In some cases even higher reso-
we will refer to Hp simply as H for thepolytropic disc. We lutionwouldbedesirable,butdoublingtheresolutionagain
use discs with H/R=0.05,0.1, and 0.2. is impractical.
The polytropic disc is a fair representation of an opti- Themodelsarerununtilanytransientstructuredisap-
cally thick disc in which the temperature is greatest at the pears. Typically, for discs with H/r = 0.1, this takes 30
≈
mid-planeanddeclineswithincreasingheightuntilthepho- orbital periods at the resonant radius (usually r = 1, see
tosphere is reached. Above this height, the atmosphere is below).Thenumberofperiodsrequiredisinverselypropor-
nearly isothermal. We quote an approximate effective opti- tional to H/r. We calculate the wave energy, mean Mach
caldepthatthemid-plane,τ,throughtherelation(e.g.Bell number, radial energy flux, and dissipation rate by averag-
et al. 1997) ing these quantitiesover N 10 waveperiods, P, after the
≈
disc has completed 40 rotations at the resonant radius. Let
τ = 8 Tm 4, (2) (vr,vθ) denote the velocity deviations from Keplerian. The
3 Ti mean (time-averaged) local wave kinetic energy density is
(cid:16) (cid:17)
calculated as
whereTm andTi arethetemperaturesatthemid-planeand
intheisothermalatmosphere,respectively.Forourstandard 1 NP 1
E(r,θ)= ρ(v2+v2) dt. (5)
parameters(AppendixAandSection6),theeffectiveoptical NP 2 r θ
depth is τ 25000. Z0
≈ The mean local Mach numberis calculated as
The mid-plane density profiles are power laws,
ρ(R,φ,z = 0) R−3/2, throughout the main body of the 1 NP
∝ (r,θ)= ρ(v2+v2)/(γp) dt. (6)
disc. The density drops smoothly at the edges of the discs M NP r θ
over a distance equal to the local value of H. For the ver- Z0
p
tically isothermal discs, we adopt R1 =0.5 and R2 = 10.0. The mean vertically integrated radial energy flux is calcu-
ThepolytropicdiscshaveouterradiiofR2 =3.0.Thevalue lated as
orefsothluetioinnncearlcrualdaituisoniss, fRo1r w=hic0h.2Rf1or=a0ll.6.but the highest- fr(r)=2π θmax N1P NP vr(p−p0) dt r2sinθ dθ (7)
To model axisymmetric waves in these discs in ZEUS- Zθmin Z0
2D, we use spherical polar coordinates (r,θ). For the verti- where p0 is the mean pressure, which is calculated over the
cally isothermal discs, we adopt inner and outer radial grid previous wave period. The dissipation rate is calculated as
boundaries rin =0.35, rout =13, and the grid encompasses the average over several wave periods of the rate, per unit
8 vertical scale-heights above and below the mid-plane (i.e. volume,atwhichthermalenergyisdepositedbytheartificial
θmin=π/2 8H/r and θmax =π/2+8H/r). viscosity.
−
Forthepolytropicdiscs,theoutergridradiusischosen
as rout =4, and 3 vertical scale-heights are modelled above
and below the mid-plane (except in Section 7, where 6 ver- 4 AXISYMMETRIC WAVE EXCITATION AND
tical scale-heights are modelled). The inner grid radius rin PROPAGATION
depends on the resolution. We adopt rin = 0.15 for all but
The radial propagation of linear axisymmetric waves in a
thehighest resolution, for which rin =0.6.
vertically isothermal accretion disc was analysed by Lubow
Forallcases,weadoptalogarithmicradialgridinwhich
& Pringle (1993). Subsequently, this work was extended to
(∆r)i+1/(∆r)i = Nr rout/rin, (3) study the radial propagation of waves in vertically ther-
mallystratifieddiscs(Korycansky&Pringle1995;Lubow&
where N is the nupmber of radial grid zones. The θ grid is Ogilvie 1998; Ogilvie & Lubow 1999) and magnetized discs
r
uniform with (Ogilvie 1998).
Nθ = (H0./1r) (∆rθ)mi+a1x/−(∆θrm)iin 1 (4) 4.1 The two-dimensional wave
−
zones. This choice means that discs with H/r = 0.1 have The simplest wave structure is that of a two-dimensional
gridzonesthathavenearlyequalradialandverticallengths. wave that propagates radially through a vertically isother-
Discs with other values of H/r have the same number of mal disc and has no vertical motion (e.g. Goldreich &
4 M. R. Bate et al.
Tremaine1979).Thedispersionrelation fortheaxisymmet- the dominant forces involved in their propagation. The f
ric two-dimensional wave is (fundamental), p (pressure-dominated), and g (buoyancy-
γp dominated)modespropagatewhereω>Ω(r),whilermodes
ω2=κ2+ k2 (8)
ρ (rotationally dominated modes)propagatewhereω<Ω(r).
There are only two f modes, fe and fo . Each of the other
where ω > 0 is the wave frequency, κ is the local epicyclic
classes contains an infinite sequence of modes with increas-
frequency of the disc, and k is the radial wavenumber. For
ing numbers of nodes in the vertical structure. These may
a Keplerian disc, we have that κ=Ω and the radius where
be classified as pe, pe, pe, ..., po, po, po, ..., etc., where
ω = κ is an outer Lindblad resonance. Its behaviour in 1 2 3 1 2 3
the number refers to the order of the mode, and ‘e’ or ‘o’
launchingwavesisessentially likethatofanouterLindblad
refers to even or odd symmetry about themid-plane.
resonanceencounteredinnon-axisymmetriccases (e.g.Gol-
In this paper, we focus our attention on the following
dreich& Tremaine1979). Thetwo-dimensional axisymmet-
modes.
ric wave can propagate only outside this radius (i.e. where
ω>Ω) for a Keplerian disc. the fe mode, which propagates where ω > Ω(r). In a
Ifsuchawavebehavespurelyisothermally (i.e. γ =1), ver•tically isothermal disc, this modeis thetwo-dimensional
itwillpropagatewithoutlosstotheouterradiusofthedisc. mode described in the previous subsection. For other types
However, as we will see in this paper, adiabatic waves with ofdiscs,thefe modebehaveslikethetwo-dimensionalwave
γ = 5/3 steepen into shocks after propagating a finite dis- close to the resonance, but behaves like a surface grav-
tance.Shocksoccurbecausethesoundspeedatthecrest of ity wave away from the resonance (Ogilvie 1998; Lubow &
the wave is greater than that in the trough. Thus, a wave Ogilvie 1998).
that is launched with a sinusoidal profile eventually steep- thepe mode,whichpropagateswhereω>Ω(r)√γ+1.
ens into a saw-tooth form and its energy is dissipated in a Its•vertica1l velocity at resonance is proportional to z (or
shock. The number of wavelengths required for the shock cosθ).
formation to occur can be c/∆c, where ∆c is the differ- the fo/ro mode, which propagates at all radii, but has
ence in sound speed betwe∼en the crest and the trough of • 1
aresonancewhereω=Ω.Itsverticalvelocityeigenfunction
the wave. This argument ignores various geometrical cor- at resonance is independentof z, while its radial velocity is
rectionsthatoccurinamulti-dimensionaldisc.Inaddition, proportional to z. The fo/ro mode is a corrugation wave or
1
it ignores the possibility that dispersive effects could lessen
discwarp,wherethemid-planeofthediscoscillates upand
thewavesteepening(Larson 1990; Papaloizou &Lin 1995). down.
Thisestimatesuggeststhatthepropagationdistancebefore
shocks set in depends on the initial amplitude of the wave, To excitethesemodesefficiently,weapply an accelera-
i.e. thewave amplitudeat the resonance. tion (a ,a ) (force per unit mass) to our equilibrium discs
r θ
of thefollowing forms.
4.2 General wave properties (i) ar = (r)sinωt and aθ = 0 for the fe mode, with
(r)= 0r−A2;
Along with the simple two-dimensional wave, a vertically A(ii) aA= 0 and a = (r)rcosθsinωt for the pe mode,
isothermal disc admits series of other waves, all of which with (rr)= 0r−3;θ A 1
possessverticalmotion(Lubow&Pringle1993).Thesewave (iii)Aa =0Aanda = (r)sinωtforthefo/ro mode,with
modescanbecategorisedbytheirverticalvelocitystructure. (r)= r0r−2. θ A 1
In a thermally stratified disc, as pointed out by Lin, Pa- A A
paloizou & Savonije (1990), the two-dimensional wave does Ineach case, thepower-law dependenceof on r ischosen
A
not exist because γp/ρ is a function of z (cf. equation 8). so that the amplitude of the non-resonant response can be
Insuchadisc,allmodesinvolveverticalmotions,asfollows reasonably small (linear) at both small and large radii. We
from Korycansky & Pringle (1995). chooseawavefrequencyω=1sothattheresonance inour
Inthesphericalgeometryofthesimulations,z=rcosθ. dimensionless units (ω = Ω or ω = Ω√γ+1) lies at r = 1
The linear wave modes then havethe form (for thefe and corrugation modes) or r=(γ+1)1/3 1.39
(for the pe mode). The simulated disc contains these≈reso-
1
Q(r,θ,t)=Re Q˜(r,θ)exp iωt+i k(r)dr , (9) nancelocations. Theacceleration isappliedfromthebegin-
−
(cid:26) (cid:20) Z (cid:21)(cid:27) ningofthecalculationsanditsstrengthisparameterisedby
whereQisanyperturbationquantity.Thefrequencyωisde- 0.
A
finedtobepositive.Thewavenumberk(r)satisfies kr 1 Each form oftheacceleration listed abovegenerally re-
| |≫
exceptnearresonances,wherethisWKBformbreaksdown. sultsintheexcitationofavarietyofmodes.Thetotalenergy
Thewavenumbercan bepositive ornegative, dependingon fluxgenerated at theresonance thatiscarried byall modes
the radial direction (inward or outward) of the phase ve- for an axisymmetric disc can be determined by adapting
locity. The amplitude Q˜ varies slowly with r. The vertical thetorqueformula of Goldreich & Tremaine(1979). InAp-
structure of the wave, and the relation between ω and k, pendix B, we derive the total energy flux for each form of
are determined by an eigenvalue problem, local in r, which acceleration listed above.
admits a set of discrete modes. Each mode has a different Inpractice,itismoreconvenienttoexpressthestrength
dispersion relation ω = ω(k,r) which must be satisfied at of the forcing in terms of the spatial peak of the time-
every r. averaged waveMach numberat thediscmid-planethat oc-
Using the notation of Ogilvie (1998), these modes can curs near the resonance, which we denote as r. In terms
M
be categorised into four classes which are determined by of our notation, we have
Waves in accretion discs 5
Figure1. Theintegratedradialenergyflux, fr,(top)andmean(time-averaged)Machnumberatthemid-planeofthedisc, (r,θ=
π/2),(bottom)areplottedasfunctionsofradiu|s|forthefemodeinaverticallyisothermaldiscwithH/r=0.1.TheresultsaresMhownfor
four different wave amplitudes: r =0.001 (left), r =0.01(centre-left), r =0.03 (centre-right), r =0.1 (right). Each case was
M M M M
performedwithdifferentnumericalresolutionsand/orviscositytotestforconvergence:250 111andqvisc=4(long-dashed),500 221
× ×
and qvisc = 4 (triple-dot-dashed), 500 221 and qvisc = 0.5 (short-dashed), 1000 441 and qvisc = 4 (dot-dashed), 1000 441 and
× × ×
qvisc=0.5(solid).Thehorizontal dottedlineindicates75percentofthefluxpredictedbytheGoldreich&Tremaineformula(equation
B1) and the associated Mach number of the wave. The value of 75 percent is that predicted by linear theory to be deposited in the fe
mode. As can be seen from the figure, the hydrodynamic results are in excellent agreement with linear theory for the low-amplitude
cases.Inthehigh-amplitudecases,theradialenergyfluxisattenuated bydissipation.
r =max[ (r,θ=π/2)], (10) where A and λ are the amplitude and wavelength of the
M M
wave, respectively.
where is defined by equation (6) and the maximum
M There are several ways to determine whether a wave
is taken over a suitable neighbourhood of the resonance
dissipates numerically or physically. First, the resolution of
(0.8 < r < 1.2). This measure of forcing strength is used
the calculations can be increased. If the dissipation is nu-
in comparing the simulations involving different disc mod-
merical,thewavewillpropagatefurtherbeforeitdissipates.
els.IncarryingoutasimulationwithachosenvalueofMr, Secondly, keeping the resolution fixed, the strength of the
we vary A0 until thedesired valueof Mr is achieved. artificial viscosity, qvisc, can be varied. If the dissipation is
numerical, the wave will propagate further before it dissi-
pates with lower qvisc. Note, however, that if qvisc 1,
≪
4.3 Numerical convergence energy will simply be lost from the calculation instead of
being converted to thermal energy (e.g. calculations with
As described above, in a vertically isothermal disc, the fe qvisc = 0.1 and qvisc = 0.01 give almost identical results).
mode should propagate outwards from the resonance with- The third way is to examine the form of the wave to see
out a loss of flux until the wave steepens to form shocks. if wave steepening occurs just before the wave dissipates.
Oneofthemaindifficultiesinobtainingnon-linearresultsis If the wave maintains a sinusoidal profile, the dissipation
spatially resolving the wave until this wave steepening and is numerical. In order to determine which of the non-linear
shock formation occur. If the flow becomes unresolved (i.e. results are due to physical dissipation, we use all three of
itvariesonascalesimilartothezonespacing),thewavewill theseindicators.
dissipateenergy.Thisdissipationoccursinoneoftwoways.
Intheabsenceofanartificialviscosity,theenergycarriedby
thewaveis lost from thecalculation. Ifthereisan artificial
viscosity, as is employed in the calculations presented here,
5 THE VERTICALLY ISOTHERMAL DISC
the wave energy is converted into thermal energy. Artificial
viscosityprovidesameansofsimulatingshocks(flowdiscon- Inthissection,westudytheexcitation,propagationanddis-
tinuities)inthecode.NumericaldissipationintheZEUS-2D sipationofadiabatic(γ =5/3)wavesinaverticallyisother-
code is proportional to the square of the velocity difference mal accretion disc. First, we consider the fe mode (a plane
across a grid zone. Thus, for an approximately linear wave, wave). Subsequently,we discuss the pe mode and the fo/ro
1 1
the numerical dissipation is proportional to qviscA2/(λNr), (corrugation) mode.
6 M. R. Bate et al.
mid-plane,therewillbesomefiniteheightatwhichthemo-
tions are supersonic and we would expect shocks to form.
However,themajority of thewaveenergy iscontained near
themid-planeandthereforethisdoesnotdiminishthewave
fluxappreciably.
The r modes will be discussed in detail in Section 5.3.
Briefly,however,there modeshouldpropagateinwardsand
1
bechannelled towards the mid-planeof thedisc.
5.1.2 Low-amplitude, non-linear results
We turn now to the non-linear hydrodynamic calculations
performed with ZEUS-2D. Figure 1 plots the magnitude of
the vertically integrated radial energy flux, fr (equation
| |
7), and the mean Mach number (equation 6) in the mid-
planeofthedisc, (r,θ=π/2),asfunctionsofr.Theyare
M
plotted for four values of the Mach number near resonance
r.Noticethattheradialenergyfluxisinfactnegativefor
rM<0.7.Thiseffectisduetheinwardlypropagatingre mode
1
that is excited along with the fe mode.
The dotted lines in the upper panels of Figure 1 show
thecontributionfromthefe modetothetotalradialenergy
flux that is predicted by linear theory. In the lower panels,
weplotthecorrespondingMachnumbersthatarepredicted
inthemid-planebylineartheory.Thenon-linearresultsare
in excellent agreement with the linear theory.
Consider the weakest fe mode (left panels, Figure 1).
Figure 2. The Mach number of the radial motion of the wave
Excitation of the fe mode occurs over a finite radial extent
inthemid-planeofthediscisplottedataparticularinstantasa
functionofradiusforthefemodeintheverticallyisothermaldisc (a few H) around r = 1. For r > 1, the flux reaches a
withH/r=0.1.Variouswaveamplitudes areshown: r=0.01 peakvaluejust beyondtheresonance.Thedropinfluxjust
M
(top), r=0.03(middle)and r=0.1(bottom). Forthelow- beyond the peak is probably due to the rapid dissipation
M M
amplitude case, r = 0.01, the wave maintains its sinusoidal of g modes, which account for about 20% of the total flux.
M
profile until numerical dissipation occurs. For the higher ampli- Thefe modecontinuestopropagatetogreater radiiandits
tudes, Mr =0.03and Mr =0.1, the wave steepens, shocks and fluxremainsnearlyconstantfor2<r<8.Thevalueofthe
dissipates.Thisindicatesthatthephysicaldissipationisnumeri- flux there agrees well with the predicted fe mode flux. The
callyresolved.
dissipation that occurs at large radii is numerical as seen
bythelackofconvergencewithincreasingresolution.Thus,
thewaveshouldcontinuetopropagateoutwards,beyondthe
5.1 fe mode grid boundary.
Forr<1,there modepropagatesinwardsandoverlaps
5.1.1 Linear theory 1
withtheevanescenttailofthefe mode.There modeisonly
1
As discussed in Section 4, we adopt a purely radial driving modelledforr >0.5becausetheedgeoftheequilibriumdisc
force perunitmass of ar(r,t)= (r)sinωt,in ordertoeffi- is at R1 =0.5.∼The two modes havefluxesof opposite sign,
ciently excite the fe mode. The fAe mode corresponds to the and thenetflux changessign near r=0.7. The peak radial
two-dimensionalwaveinaverticallyisothermaldisc.Linear energy flux from the simulation in this region is about 3%
theory predicts that the radial forcing of Section 4 should ofthetotalflux.This result isconsistent with theexpected
lead to theoutwardly propagating fe mode receiving 74.5% re mode negative flux of about 6% of the total overlapping
1
of the total radial energy flux (see Appendix B1). The in- with the tail of the positive fe mode flux. The r modes will
wardly propagating re mode should receive about 6.0% of bediscussed in moredetailin section 5.3 wherean ro mode
1 1
thetotal flux.The remainder is goes intog modes (see Ap- receives 50% of thetotal flux.
pendix B1). These results are reinforced by Figure 3, which gives
Accordingtolineartheory(Lubow&Pringle1993),the contour plots of the wave energy density (equation 5) and
fe mode should occupy the entire vertical thickness of the the energy dissipation rate for the low-amplitude calcula-
disc and propagate outwards in radius. The mean Mach tions (left panels). There is a lot of dissipation at high al-
number at the mid-plane of the disc should increase as titude in the disc in the range 1 < r < 2, as expected for
(r,θ = π/2) r1/2 for r 1. This weak dependence the g modes. In the bulk of the disc, there is no significant
MoftheMach num∝beronr mean≫sthatalow-amplitudewave dissipation and the fe mode propagates outwards with vir-
is expected to propagate a large distance before it steep- tuallynolossofflux(asmallamountislostathighaltitudes
ens into a shock and dissipates. The vertical structure of where themotions become supersonic).
the wave can also be predicted by linear theory. The mean Thebehaviouroftheinwardlypropagatingre modecan
1
Mach number varies vertically as exp(z2/5H2), for alsobeseenintheleft-handpanelsofFigure3.Aspredicted
M ∝
γ = 5/3. Thus, although a wave may be linear near the bylinear theory(Lubow&Pringle 1993), thewaveischan-
Waves in accretion discs 7
Figure 3. The case considered here is the fe mode in a vertically isothermal disc with H/r =0.1. The kinetic energy density in the
wave (logE, equation 5) is shown in the upper panels. The energy dissipation rate per unit volume is shown in the lower panels as a
logarithmic grey-scale. The grey-scales cover three orders of magnitude. The left panels are for amplitude r = 0.001 and the right
panelsareforamplitude r=0.1.Atlowamplitudethefe modepropagatestolargeradiusthroughoutthebModyofthediscunaffected
M
by the small amount of dissipation that occurs at high z. In the high-amplitude case dissipation occurs in the body of the disc and
| |
significantlyattenuates thewave(seeFigure1).InbothcasesthegridresolutionisNr×Nθ =1000×441andqvisc=4.0.
nelledtowardsthemid-plane.Thiswillbediscussedfurther gate a shorter distance before numerical dissipation occurs
in Section 5.3. (cf. Figure 1, results with r =0.001, and r =0.01).
M M
Physicaldissipationoccursifthesinusoidalwavesteep-
ens into a shock, which takes place in c/∆c wavelengths
5.1.3 Dependence on wave amplitude ∼
fromtheresonance(Section4.1).Forwaveswithverysmall
As the driving amplitude is increased, both numerical and values of r, steepening requires many wavelengths and
M
physical dissipation occur more quickly. Numerical dissipa- numerical dissipation sets in first, as we saw in thecase for
tion increases as the square of the amplitude of the wave r = 0.001. Physical dissipation in our simulations was
M
(Section 4.3). Thus, waves with larger values of r propa- obtained for cases with r = 0.03 and 0.1. Consider the
M M
8 M. R. Bate et al.
Figure4. Theintegratedradialenergyflux, fr,(top)andmean(time-averaged)Machnumberatthemid-planeofthedisc, (r,θ=
π/2), (bottom) are plotted as functions of ra|diu|s for the fe mode in a vertically isothermal disc with H/r = 0.05. The reMsults are
shown for three different wave amplitudes: r = 0.001 (left), r = 0.01 (centre), r = 0.1 (right). Each case was performed with
M M M
differentnumericalresolutionstotestforconvergence:250 111andqvisc=4(long-dashed),500 221andqvisc=4(triple-dot-dashed),
× ×
1000 441andqvisc=4(dot-dashed).Thehorizontaldottedlineindicates75percentofthefluxpredictedbytheGoldreich&Tremaine
×
formula (equation B1)and the associated Mach number of the wave. The value of 75 percent is that predicted by linear theory to be
depositedinthefe mode.Ascanbeseenfromthefigure,thehydrodynamicresultsareinexcellentagreementwithlineartheoryforthe
lowestamplitudecase.Inthehigh-amplitudecases,theradialenergyfluxisattenuated bydissipation.
panels for r = 0.03 in Figure 1. Convergence of flux and Figures1,4and5).Thus,forexample,waveswith r =0.1
M M
(r,θ = π/2) is obtained for r < 4 with resolutions of are well resolved in the hottest disc, H/r = 0.2, even with
M500 221 and 1000 441. With∼an initial amplitude of N N =250 111, and dissipate dueto shocks, whereas
r θ
× × × ×
r =0.1, evena resolution of 250 111 issufficient to ob- inthecoldestdisc,H/r=0.05,eventhehighestresolutions
M ×
tain physical dissipation. In the r = 0.03 case, the wave do not give convincing convergence (numerical dissipation
M
flux has been reduced by an order of magnitude by r 10, still plays a significant role).
≈
whereas in the r =0.1 case this level of reduction occurs Strictly,lineartheoryassumesH r.However,wefind
M ≪
at r 4. that,evenwithH/r=0.2,lineartheoryprovidesagoodde-
≈
To emphasise that this dissipation is physical rather scription of thewave excitation and propagation. A further
than numerical, we plot in Figure 2 the radial velocity prediction that can be derived from linear theory is that
profile of the wave in the mid-plane for the cases with theratiooftheMachnumberattheresonancetotheMach
r = 0.01,0.03, and 0.1. Steepening of the wave towards number at r rres should be nearly independent of H/r.
M ≫
the profile of a shock-wave is clearly visible in those cases This is confirmed by the non-linear calculations presented
with r =0.03 and0.1, whereasinthe r =0.01calcula- here.
M M
tiontheprofileremainssinusoidalindicatingonlynumerical,
rather than physical, dissipation. This result is consistent
5.2 pe mode
with Figure 3 which shows that increasing r decreases 1
thevolume in which thefe mode propagates.M 5.2.1 Linear theory
Asdiscussed inSection 4,weadoptadrivingforceperunit
massthatispurelyintheθ-directionoftheforma (r,θ,t)=
5.1.4 Dependence on disc thickness θ
(r)rcosθsinωt,in order to efficiently excite thepe mode.
InFigures4and 5,weplot theenergy fluxandmean Mach TAhe pe mode is launched at a vertical resonance r1 1.39
1 ≈
number (r,θ = π/2) as functions of r for discs with and propagates outwards. We predict its radial energy flux
M
H/r = 0.05 and H/r = 0.2, respectively. The main dif- in AppendixB2. No otherwave is expected.
ference between the waves in these discs and those in the Unlikethefe mode, thepe mode hasvertical structure
1
H/r=0.1discisthatthewavelengthisproportionaltoH/r. (seeLubow&Pringle1993).Inthevicinityoftheresonance,
Since the numerical dissipation is 1/(λN ) (Section 4.3), the wave energy has a minimum on the mid-plane of the
r
∝
awaveinathickdisccanbefollowed tolargerradiithanin discandtwomaximaatz/H 2.Thisisduetotheformof
≈
a thin disc for the same radial resolution (cf. r =0.01 in the excitation which depends linearly on z (Section 4). For
M
Waves in accretion discs 9
Figure 6. Theintegratedradialenergyflux,|fr|,isplottedasafunctionofradiusforthepe1 modeinaverticallyisothermaldiscwith
H/r =0.1. The results areshown for three different wave amplitudes: r =0.001 (left), r =0.01 (centre), r =0.1 (right). Each
M M M
case was performed with different numerical resolutions to test for convergence: 250 111 and qvisc =4 (long-dashed), 500 221 and
× ×
qvisc = 4 (triple-dot-dashed). The horizontal dotted line indicates the flux predicted by linear theory (equation B20). As can be seen
fromthefigure,thehydrodynamicresultsareinexcellentagreementwithlineartheoryforthelowamplitudecases.Athighamplitude,
theradialenergyfluxisattenuated bydissipation.
r >1.7,thischangestothreemaxima,oneonthemid-plane
an∼d two at z/H 3, with two minima at z/H 2. Lubow
≈ ≈
& Pringle (1993) interpreted this as the wave energy rising
upwithinthediscasthewavepropagatesoutwards.Infact,
however, the vertical distribution of the wave energy does
not change beyond r 2. The wave propagates outwards
≈
maintainingtheabovedistributionofenergyindefinitely.As
withthefemode,thewaveiseventuallyexpectedtoundergo
steepening into a shock and dissipate, but a low-amplitude
wave will propagate for a large distance before this takes
place.
5.2.2 Low-amplitude, non-linear results
Figure6plotsthemagnitudeoftheverticallyintegratedra-
dialenergyflux, fr (equation7)asafunctionofr.Wecon-
| |
siderthreevaluesofthemeanMachnumbernearresonance
Mr.Aswiththecaseforthefe mode,thelow-amplitudepe1
modepropagatesoutwardsfromtheresonanceuntilmostof
theenergy is dissipated numerically.
A small negative energy flux is seen for r < 0.9. We
attribute this to a low-amplitude r mode that is∼excited at
theLindbladresonanceatr=1.Inaverythindisc,vertical
Figure 5. The integrated radial energy flux, fr, (top) and
| | forcingofthetypeappliedhereshouldnotexciteanymodes
mean(time-averaged)Machnumberatthemid-planeofthedisc,
(r,θ = π/2), (bottom) are plotted as functions of radius for attheLindbladresonance.However,whenH/r=0.1,there
Mthe fe mode ina vertically isothermal discwith H/r =0.2. The issomeambiguitybetweenthemeaningsof‘horizontal’ and
results are shown for two different wave amplitudes: r =0.01 ‘vertical’awayfromthemid-plane.Asaresult,somelaunch-
M
(left), r = 0.1 (right). Each case was performed with differ- ing at the Lindblad resonance is possible at the level of
M
ent numerical resolutions to test for convergence: 250 111 and (H/r)2.
×
qvisc = 4 (long-dashed), 500 221 and qvisc = 4 (triple-dot- InFigure7weplotthewaveenergyandenergydissipa-
×
tdaalshdeodt)t,ed10li0n0e×in4d4ic1ataensd75qvpisecrc=ent4o(fdotht-edaflsuhxedp)r.edTichteedhobryiztohne- tion rate for the low-amplitude pe1 mode (left panels). The
vertical distribution of the wave energy density is in excel-
Goldreich&Tremaineformula(equation B1)andtheassociated
lentagreementwithlineartheoryinthebulkofthediscwith
Mach number of the wave. The value of 75 percent is that pre-
dicted by linear theory to be deposited in the fe mode. As can onlyasmallamountofdissipationathighaltitudewherethe
beseenfromthefigure,thehydrodynamicresultsfortheenergy motions become sonic.
fluxareingoodagreement withlineartheoryforthelowestam-
plitudecase,butforthisrelativelythickdiscthepredictedMach
numberunderestimatesthecomputedvaluebyabout30percent.
5.2.3 Wave amplitude and disc thickness
In the high-amplitude case, the radial energy flux is attenuated
bydissipation. As with the fe mode, the dissipation occurs more rapidly
with increased driving. Calculations with r =0.1 resolve
M
the wave until physical dissipation due to wave steepening
10 M. R. Bate et al.
Figure 7. The case considered here is the pe mode in a vertically isothermal disc with H/r =0.1. The kinetic energy density in the
1
wave(logE,equation5)isshownintheupperpanelscoveringfourordersofmagnitude.Theenergydissipationrateperunitvolumeis
shown in the lower panels as a logarithmic grey-scale covering six orders of magnitude. The left panels are for amplitude r = 0.001
athneddtihsecuringahfftepctaendelbsyartheefosrmaamllpalmituoudnetMofrd=issi0p.1a.tiAonttlhowataomccpulritsuadtehtihgehpze1.mTohdeepperompoadgeatdeissptloaylasrtgheervaedrituicsatlhsrtoruugchtuoruetpMtrheedibcoteddyboyf
| | 1
lineartheory.Inthehigh-amplitudecasedissipationoccursinthebodyofthediscandsignificantlyattenuates thewave(seeFigure6).
InbothcasesthegridresolutionisNr×Nθ =500×221andqvisc=4.0.
occurs,butthelowerMr casesarestilllimitedbynumerical 5.3 fo/ro1 mode
dissipation.
5.3.1 Linear theory
We have not performed calculations of the pe mode
1
for discs with different thicknesses since, as we found for Asdiscussed inSection 4,weadoptadrivingforceperunit
fe mode,wedonotexpectanysignificantdependenceofthe mass that purely in the θ-direction of the form a (r,t) =
θ
propagation of the pe mode on H/r. As with the fe mode, (r)sinωt, in order to efficiently excite the fo/ro mode.
the longer wavelengt1h in hotter discs would make it easier TAhe fo/ro mode, or corrugation mode, is also laun1ched at
1
to resolve thewave to larger radii. a vertical resonance although, in a Keplerian disc, this co-