9512.net
甜梦文库
当前位置:首页 >> >>

The Excitation, Propagation and Dissipation of Waves in Accretion Discs The Non-linear Axis



Mon. Not. R. Astron. Soc. 000, 000–000 (0000)

Printed 1 February 2008

a (MN LTEX style ?le v2.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 1
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA of Physics, University of Exeter, Stocker Road, Exeter EX4 4QL 3 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
2 School

arXiv:astro-ph/0201149v1 10 Jan 2002

1 February 2008

ABSTRACT

We analyse the non-linear propagation and dissipation of axisymmetric waves in accretion discs using the ZEUS-2D hydrodynamics code. The waves are numerically resolved in the vertical and radial directions. Both vertically isothermal and thermally strati?ed accretion discs are considered. The waves are generated by means of resonant forcing and several forms of forcing are considered. Compressional motions are taken to be locally adiabatic (γ = 5/3). Prior to non-linear dissipation, the numerical results are in excellent agreement with the linear theory of wave channelling in predicting the types of modes that are excited, the energy ?ux by carried by each mode, and the vertical wave energy distribution as a function of radius. In all cases, waves are excited that propagate on both sides of the resonance (inwards and outwards). For vertically isothermal discs, non-linear dissipation occurs primarily through shocks that result from the classical steepening of acoustic waves. For discs that are substantially thermally strati?ed, wave channelling is the primary mechanism for shock generation. Wave channelling boosts the Mach number of the wave by vertically con?ning the wave to a small cool region at the base of the disc atmosphere. In general, outwardly propagating waves with Mach numbers near resonance Mr > 0.01 undergo ? shocks within a distance of order the resonance radius. Key words: accretion, accretion discs – binaries: general – hydrodynamics – waves.

1

INTRODUCTION

Gaseous discs play a signi?cant role in the evolution of binary star and planetary systems. Generally, the tidal forces from stellar or planetary companions distort the discs from an axisymmetric form. However, where resonances occur in the discs, the tidal forces also generate waves that transport energy and angular momentum. The resulting resonant torques cause orbital evolution of the perturbing objects (e.g. Goldreich & Tremaine 1980; Lin & Papaloizou 1993; Lubow & Artymowicz 1996). The waves also cause the discs to evolve since they transfer their energy and angular momentum to the disc in the locations where they damp. Damping may be provided by shocks, radiative damping (Cassen & Woolum 1996), or other viscous damping mechanisms. Waves in discs have also been considered as possible explanations of quasi-periodic variability in the X-ray emission of systems involving accretion discs around black holes (e.g. Kato & Fukue 1980; Nowak & Wagoner 1991; Kato 2001).

Initial studies of wave propagation approximated the disc as two-dimensional and ignored the e?ects of vertical structure (perpendicular to the disc plane). Goldreich & Tremaine (1979) developed a two-dimensional linear theory for resonant tidal torques and the associated wave propagation. Later studies of the non-linear case found that the torques were within a few percent of those predicted by the two-dimensional linear formula (Shu, Yuan & Lissauer 1985; Yuan & Cassen 1994). For a thin disc that is vertically isothermal and has an isothermal thermodynamic response (γ = 1), the only wave excited in the disc has a twodimensional structure. The wavefronts are perpendicular to the disc plane and the motion is purely horizontal (parallel to the disc plane). Thus, a two-dimensional treatment is valid in a thin disc if both the vertical structure of the disc and its thermodynamic response are locally isothermal. However, this is not realistic for most discs, such as those in cataclysmic variables, circumstellar and circumbinary discs around pre-main-sequence stars, and protoplanetary discs. In such discs, the waves are no longer two-dimensional and the vertical structure of the disc must be taken into account. A study of three-dimensional wave propagation was

?

E-mail: mbate@astro.ex.ac.uk

2

M. R. Bate et al.
structure of the equilibrium discs that we model and the grid layout for the ZEUS-2D calculations. Section 4 brie?y reviews the semi-analytical linear method for solving the three-dimensional wave propagation problem, describes our method of wave excitation, and discusses the requirements for convergence of the numerical calculations. Sections 5 and 6 contain the results for a vertically isothermal disc, and a polytropic disc with a vertically isothermal atmosphere, respectively. Intermediate discs with optical depths not much larger than unity are considered in Section 7. A summary of the results is contained in Section 8.

made using linearised numerical simulations (Lin, Papaloizou & Savonije 1990a,b). Unfortunately, the limited range of physical parameter space that was covered led to the interpretation that waves in a vertically thermally strati?ed disc are refracted upwards into the atmosphere where they dissipate via shocks. Later, it was realised that the linear problem could be solved semi-analytically (Lubow & Pringle 1993; Korycansky & Pringle 1995; Lubow & Ogilvie 1998; Ogilvie & Lubow 1999). These studies found that the propagation of waves in a thin disc is similar to the propagation of electromagnetic waves along a waveguide (Jackson 1962; Landau & Lifshitz 1960). The motions in the disc can be described in terms of modes. For each mode, the vertical wave structure is determined in detail at each disc radius. The horizontal variations of the mode are treated by means of a WKB approximation in which the horizontal (radial) wavenumber is determined as a function of wave frequency. A ?ux conservation law determines the mode amplitude as a function of radius. Depending on the vertical structure of the disc and the thermodynamic response of the gas, the wave energy may be channelled towards the surface of the disc, towards the mid-plane of the disc, or occupy the entire vertical extent as it propagates in radius. Ideally, one would like to determine the non-linear, nonaxisymmetric response of a variety of disc models (vertically isothermal and thermally strati?ed) to resonant forcing. Unfortunately, this is a formidable numerical problem at present. Instead we begin, in this paper, with an analysis of the non-linear axisymmetric response. The axisymmetric case o?ers the major advantage of being expressed as a two-dimensional problem, which is currently accessible numerically. Although axisymmetric waves only transport energy and not angular momentum, they should resemble the response of the disc to low-azimuthal number tidal forcing, as arises in binary star systems. The local physics of the waveguide is determined within a region whose size is a few times the semi-thickness of the disc. On this scale a non-axisymmetric wave of low azimuthal wavenumber is almost indistinguishable from an axisymmetric wave. This explains why the local physics of the waveguide can be studied within the shearing-sheet approximation (Lubow & Pringle 1993) and why the azimuthal wavenumber appears in the dispersion relation only through a Doppler shift of the wave frequency. Consequently, in this paper, we study the non-linear propagation and dissipation of axisymmetric, but fully resolved, waves in accretion discs. The purpose of the paper is threefold. First, we wish to determine, using non-linear hydrodynamic calculations, whether or not the semi-analytical linear theory provides an accurate description of wave propagation in accretion discs. Secondly, although linear theory can predict where non-linearity is likely to become important, shocks and the process of energy deposition (and angular momentum deposition with non-axisymmetric waves) cannot be handled. We aim to determine how and where dissipation occurs through shocks and how accurately this can be predicted from the linear analysis. Finally, we wish to determine how well non-linear hydrodynamic calculations can model the problem and what resolution is required. The outline of this paper is as follows. In Section 2, we brie?y review the ZEUS-2D hydrodynamic code that is used to obtain the non-linear results. In Section 3, we discuss the

2

NUMERICAL METHOD

The non-linear hydrodynamic calculations presented here were performed with the ZEUS-2D code developed by Stone & Norman (1992). No modi?cations to the code were required, except the addition of damping boundary conditions for test purposes only (see below). The ZEUS-2D code comes with various options for the interpolation scheme and arti?cial viscosity. We use the van Leer (second-order) interpolation scheme. Linear and quadratic arti?cial viscosities are provided in the code, with two forms of the quadratic viscosity: the standard von Neumann & Richtmyer (1950) form, and a tensor form. We use only a quadratic arti?cial viscosity to handle shocks. Calculations were performed both with the von Neumann & Richtmyer form and with the tensor form. We ?nd no signi?cant di?erence in the results. All the results presented here use the standard von Neumann & Richtmyer form. The strength of the quadratic viscosity is controlled by a parameter qvisc . Generally, we use qvisc = 4, but we present some results from calculations that were performed with qvisc = 0.5. No other type of viscosity (e.g. one that includes a shear stress) is used. All of the calculations presented here use re?ecting boundaries. In order to verify that our results were not affected by re?ections from the boundaries, some calculations were performed with damping boundaries at the inner and outer radii of the numerical grid. Damping is provided by adding an acceleration ?v = ?2ωv (1) ?t to the inner or outer 10?20 radial zones of the grid, where v is the velocity of the gas and ω is the frequency of the forcing that is used to generate the waves (see Section 4). There is no signi?cant di?erence between the results obtained using re?ecting and damping boundaries.

3

MODELLING THE ACCRETION DISCS

We consider the propagation and dissipation of axisymmetric (m = 0) waves in two types of accretion disc: a locally vertically isothermal disc, and a polytropic disc (polytropic index n = 1.5) that joins smoothly on to an isothermal atmosphere. We ignore all e?ects of disc turbulence. A full description of the equilibrium structures is given in Appendix A; we only brie?y list their properties here. The disc is described in cylindrical coordinates (R, φ, z). As described in Appendix A, we shall work in dimensional units, such that

Waves in accretion discs
the disc angular velocity is unity when the dimensionless disc radius R = 1. The discs are Keplerian, with angular velocity ? = R?3/2 , apart from small corrections. In both discs, the value of the polytropic exponent that describes the adiabatic pressure-density relation for the waves is γ = 5/3. For convenience in modelling, we choose discs with ?nite inner and outer radii, R1 and R2 . We also choose a scale-height that increases linearly with radius, which requires that the sound speed varies as c ∝ R?1/2 . The scale-height of the vertically isothermal disc is H ∝ R. The polytropic disc has a well-de?ned polytropic layer with a semi-thickness Hp , that smoothly changes into an isothermal atmosphere with its own scale-height Hi , above z ≈ Hp . Both Hp and Hi are proportional to R. From this point on, we will refer to Hp simply as H for the polytropic disc. We use discs with H/R = 0.05, 0.1, and 0.2. The polytropic disc is a fair representation of an optically thick disc in which the temperature is greatest at the mid-plane and declines with increasing height until the photosphere is reached. Above this height, the atmosphere is nearly isothermal. We quote an approximate e?ective optical depth at the mid-plane, τ , through the relation (e.g. Bell et al. 1997) τ = 8 3 Tm Ti
4

3

,

(2)

where Tm and Ti are the temperatures at the mid-plane and in the isothermal atmosphere, respectively. For our standard parameters (Appendix A and Section 6), the e?ective optical depth is τ ≈ 25000. The mid-plane density pro?les are power laws, ρ(R, φ, z = 0) ∝ R?3/2 , throughout the main body of the disc. The density drops smoothly at the edges of the discs over a distance equal to the local value of H. For the vertically isothermal discs, we adopt R1 = 0.5 and R2 = 10.0. The polytropic discs have outer radii of R2 = 3.0. The value of the inner radius is R1 = 0.2 for all but the highestresolution calculations, for which R1 = 0.6. To model axisymmetric waves in these discs in ZEUS2D, we use spherical polar coordinates (r, θ). For the vertically isothermal discs, we adopt inner and outer radial grid boundaries rin = 0.35, rout = 13, and the grid encompasses 8 vertical scale-heights above and below the mid-plane (i.e. θmin = π/2 ? 8H/r and θmax = π/2 + 8H/r). For the polytropic discs, the outer grid radius is chosen as rout = 4, and 3 vertical scale-heights are modelled above and below the mid-plane (except in Section 7, where 6 vertical scale-heights are modelled). The inner grid radius rin depends on the resolution. We adopt rin = 0.15 for all but the highest resolution, for which rin = 0.6. For all cases, we adopt a logarithmic radial grid in which (?r)i+1 /(?r)i =
Nr

zones per vertical scale-height and the same radial spacing, but each zone is either stretched or compressed vertically. In practice, calculations of vertically isothermal discs were performed with Nr × Nθ zones of: 250 × 111, 500 × 221, or 1000 × 441. Calculations of the standard polytropic discs were computed with: 250 × 45, 500 × 91, 1000 × 183, or 1156 × 365 zones (the last of which has a resolution equivalent to 2000 × 365 zones but models a smaller range of radii than the lower-resolution calculations). The highest resolution calculations took up to 6000 CPU hours on 440 MHz Sun Workstations with each factor of 2 increase in resolution taking a factor of ≈ 8 longer (a factor of 4 from the increase in the number of zones and a factor of 2 from the decrease in the time step). In some cases even higher resolution would be desirable, but doubling the resolution again is impractical. The models are run until any transient structure disappears. Typically, for discs with H/r = 0.1, this takes ≈ 30 orbital periods at the resonant radius (usually r = 1, see below). The number of periods required is inversely proportional to H/r. We calculate the wave energy, mean Mach number, radial energy ?ux, and dissipation rate by averaging these quantities over N ≈ 10 wave periods, P , after the disc has completed 40 rotations at the resonant radius. Let (vr , vθ ) denote the velocity deviations from Keplerian. The mean (time-averaged) local wave kinetic energy density is calculated as E(r, θ) = 1 NP
NP 0

1 2 2 ρ(vr + vθ ) dt. 2

(5)

The mean local Mach number is calculated as M(r, θ) = 1 NP
NP 2 2 ρ(vr + vθ )/(γp) dt. 0

(6)

The mean vertically integrated radial energy ?ux is calculated as
θmax

fr (r) = 2π
θmin

1 NP

NP 0

vr (p ? p0 ) dt r 2 sin θ dθ

(7)

where p0 is the mean pressure, which is calculated over the previous wave period. The dissipation rate is calculated as the average over several wave periods of the rate, per unit volume, at which thermal energy is deposited by the arti?cial viscosity.

4

AXISYMMETRIC WAVE EXCITATION AND PROPAGATION

rout /rin ,

(3)

where Nr is the number of radial grid zones. The θ grid is uniform with Nθ = θmax ? θmin 0.1 (H/r) (?r)i+1 /(?r)i ? 1 (4)

The radial propagation of linear axisymmetric waves in a vertically isothermal accretion disc was analysed by Lubow & Pringle (1993). Subsequently, this work was extended to study the radial propagation of waves in vertically thermally strati?ed discs (Korycansky & Pringle 1995; Lubow & Ogilvie 1998; Ogilvie & Lubow 1999) and magnetized discs (Ogilvie 1998). 4.1 The two-dimensional wave

zones. This choice means that discs with H/r = 0.1 have grid zones that have nearly equal radial and vertical lengths. Discs with other values of H/r have the same number of

The simplest wave structure is that of a two-dimensional wave that propagates radially through a vertically isothermal disc and has no vertical motion (e.g. Goldreich &

4

M. R. Bate et al.
the dominant forces involved in their propagation. The f (fundamental), p (pressure-dominated), and g (buoyancydominated) modes propagate where ω > ?(r), while r modes (rotationally dominated modes) propagate where ω < ?(r). There are only two f modes, f e and f o . Each of the other classes contains an in?nite sequence of modes with increasing numbers of nodes in the vertical structure. These may be classi?ed as pe , pe , pe , . . . , po , po , po , . . . , etc., where 3 1 2 3 1 2 the number refers to the order of the mode, and ‘e’ or ‘o’ refers to even or odd symmetry about the mid-plane. In this paper, we focus our attention on the following modes. ? the f e mode, which propagates where ω > ?(r). In a vertically isothermal disc, this mode is the two-dimensional mode described in the previous subsection. For other types of discs, the f e mode behaves like the two-dimensional wave close to the resonance, but behaves like a surface gravity wave away from the resonance (Ogilvie 1998; Lubow & Ogilvie 1998). √ ? the pe mode, which propagates where ω > ?(r) γ + 1. 1 Its vertical velocity at resonance is proportional to z (or cos θ). ? the f o /ro mode, which propagates at all radii, but has 1 a resonance where ω = ?. Its vertical velocity eigenfunction at resonance is independent of z, while its radial velocity is proportional to z. The f o /ro mode is a corrugation wave or 1 disc warp, where the mid-plane of the disc oscillates up and down. To excite these modes e?ciently, we apply an acceleration (ar , aθ ) (force per unit mass) to our equilibrium discs of the following forms. (i) ar = A(r) sin ωt and aθ = 0 for the f e mode, with A(r) = A0 r ?2 ; (ii) ar = 0 and aθ = A(r)r cos θ sin ωt for the pe mode, 1 with A(r) = A0 r ?3 ; (iii) ar = 0 and aθ = A(r) sin ωt for the f o /ro mode, with 1 A(r) = A0 r ?2 . In each case, the power-law dependence of A on r is chosen so that the amplitude of the non-resonant response can be reasonably small (linear) at both small and large radii. We choose a wave frequency ω = 1 so that the resonance in our √ dimensionless units (ω = ? or ω = ? γ + 1) lies at r = 1 (for the f e and corrugation modes) or r = (γ + 1)1/3 ≈ 1.39 (for the pe mode). The simulated disc contains these reso1 nance locations. The acceleration is applied from the beginning of the calculations and its strength is parameterised by A0 . Each form of the acceleration listed above generally results in the excitation of a variety of modes. The total energy ?ux generated at the resonance that is carried by all modes for an axisymmetric disc can be determined by adapting the torque formula of Goldreich & Tremaine (1979). In Appendix B, we derive the total energy ?ux for each form of acceleration listed above. In practice, it is more convenient to express the strength of the forcing in terms of the spatial peak of the timeaveraged wave Mach number at the disc mid-plane that occurs near the resonance, which we denote as Mr . In terms of our notation, we have

Tremaine 1979). The dispersion relation for the axisymmetric two-dimensional wave is γp 2 k (8) ω 2 = κ2 + ρ where ω > 0 is the wave frequency, κ is the local epicyclic frequency of the disc, and k is the radial wavenumber. For a Keplerian disc, we have that κ = ? and the radius where ω = κ is an outer Lindblad resonance. Its behaviour in launching waves is essentially like that of an outer Lindblad resonance encountered in non-axisymmetric cases (e.g. Goldreich & Tremaine 1979). The two-dimensional axisymmetric wave can propagate only outside this radius (i.e. where ω > ?) for a Keplerian disc. If such a wave behaves purely isothermally (i.e. γ = 1), it will propagate without loss to the outer radius of the disc. However, as we will see in this paper, adiabatic waves with γ = 5/3 steepen into shocks after propagating a ?nite distance. Shocks occur because the sound speed at the crest of the wave is greater than that in the trough. Thus, a wave that is launched with a sinusoidal pro?le eventually steepens into a saw-tooth form and its energy is dissipated in a shock. The number of wavelengths required for the shock formation to occur can be ? c/?c, where ?c is the di?erence in sound speed between the crest and the trough of the wave. This argument ignores various geometrical corrections that occur in a multi-dimensional disc. In addition, it ignores the possibility that dispersive e?ects could lessen the wave steepening (Larson 1990; Papaloizou & Lin 1995). This estimate suggests that the propagation distance before shocks set in depends on the initial amplitude of the wave, i.e. the wave amplitude at the resonance.

4.2

General wave properties

Along with the simple two-dimensional wave, a vertically isothermal disc admits series of other waves, all of which possess vertical motion (Lubow & Pringle 1993). These wave modes can be categorised by their vertical velocity structure. In a thermally strati?ed disc, as pointed out by Lin, Papaloizou & Savonije (1990), the two-dimensional wave does not exist because γp/ρ is a function of z (cf. equation 8). In such a disc, all modes involve vertical motions, as follows from Korycansky & Pringle (1995). In the spherical geometry of the simulations, z = r cos θ. The linear wave modes then have the form Q(r, θ, t) = Re ? Q(r, θ) exp ?iωt + i k(r) dr , (9)

where Q is any perturbation quantity. The frequency ω is de?ned to be positive. The wavenumber k(r) satis?es |kr| ? 1 except near resonances, where this WKB form breaks down. The wavenumber can be positive or negative, depending on the radial direction (inward or outward) of the phase ve? locity. The amplitude Q varies slowly with r. The vertical structure of the wave, and the relation between ω and k, are determined by an eigenvalue problem, local in r, which admits a set of discrete modes. Each mode has a di?erent dispersion relation ω = ω(k, r) which must be satis?ed at every r. Using the notation of Ogilvie (1998), these modes can be categorised into four classes which are determined by

Waves in accretion discs

5

Figure 1. The integrated radial energy ?ux, |fr |, (top) and mean (time-averaged) Mach number at the mid-plane of the disc, M(r, θ = π/2), (bottom) are plotted as functions of radius for the f e mode in a vertically isothermal disc with H/r = 0.1. The results are shown for four di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre-left), Mr = 0.03 (centre-right), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions and/or viscosity to test for convergence: 250 × 111 and qvisc = 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). The horizontal dotted line indicates 75 percent of the ?ux predicted by the Goldreich & 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 deposited in the f e mode. As can be seen from the ?gure, the hydrodynamic results are in excellent agreement with linear theory for the low-amplitude cases. In the high-amplitude cases, the radial energy ?ux is attenuated by dissipation.

Mr = max[M(r, θ = π/2)],

(10)

where M is de?ned by equation (6) and the maximum is taken over a suitable neighbourhood of the resonance (0.8 < r < 1.2). This measure of forcing strength is used in comparing the simulations involving di?erent disc models. In carrying out a simulation with a chosen value of Mr , we vary A0 until the desired value of Mr is achieved.

4.3

Numerical convergence

As described above, in a vertically isothermal disc, the f e mode should propagate outwards from the resonance without a loss of ?ux until the wave steepens to form shocks. One of the main di?culties in obtaining non-linear results is spatially resolving the wave until this wave steepening and shock formation occur. If the ?ow becomes unresolved (i.e. it varies on a scale similar to the zone spacing), the wave will dissipate energy. This dissipation occurs in one of two ways. In the absence of an arti?cial viscosity, the energy carried by the wave is lost from the calculation. If there is an arti?cial viscosity, as is employed in the calculations presented here, the wave energy is converted into thermal energy. Arti?cial viscosity provides a means of simulating shocks (?ow discontinuities) in the code. Numerical dissipation in the ZEUS-2D code is proportional to the square of the velocity di?erence across a grid zone. Thus, for an approximately linear wave, the numerical dissipation is proportional to qvisc A2 /(λNr ),

where A and λ are the amplitude and wavelength of the wave, respectively. There are several ways to determine whether a wave dissipates numerically or physically. First, the resolution of the calculations can be increased. If the dissipation is numerical, the wave will propagate further before it dissipates. Secondly, keeping the resolution ?xed, the strength of the arti?cial viscosity, qvisc , can be varied. If the dissipation is numerical, the wave will propagate further before it dissipates with lower qvisc . Note, however, that if qvisc ? 1, energy will simply be lost from the calculation instead of being converted to thermal energy (e.g. calculations with qvisc = 0.1 and qvisc = 0.01 give almost identical results). The third way is to examine the form of the wave to see if wave steepening occurs just before the wave dissipates. If the wave maintains a sinusoidal pro?le, the dissipation is numerical. In order to determine which of the non-linear results are due to physical dissipation, we use all three of these indicators.

5

THE VERTICALLY ISOTHERMAL DISC

In this section, we study the excitation, propagation and dissipation of adiabatic (γ = 5/3) waves in a vertically isothermal accretion disc. First, we consider the f e mode (a plane wave). Subsequently, we discuss the pe mode and the f o /ro 1 1 (corrugation) mode.

6

M. R. Bate et al.
mid-plane, there will be some ?nite height at which the motions are supersonic and we would expect shocks to form. However, the majority of the wave energy is contained near the mid-plane and therefore this does not diminish the wave ?ux appreciably. The r modes will be discussed in detail in Section 5.3. Brie?y, however, the re mode should propagate inwards and 1 be channelled towards the mid-plane of the disc. 5.1.2 Low-amplitude, non-linear results

Figure 2. The Mach number of the radial motion of the wave in the mid-plane of the disc is plotted at a particular instant as a function of radius for the fe mode in the vertically isothermal disc with H/r = 0.1. Various wave amplitudes are shown: Mr = 0.01 (top), Mr = 0.03 (middle) and Mr = 0.1 (bottom). For the lowamplitude case, Mr = 0.01, the wave maintains its sinusoidal pro?le until numerical dissipation occurs. For the higher amplitudes, Mr = 0.03 and Mr = 0.1, the wave steepens, shocks and dissipates. This indicates that the physical dissipation is numerically resolved.

5.1 5.1.1

f e mode Linear theory

As discussed in Section 4, we adopt a purely radial driving force per unit mass of ar (r, t) = A(r) sin ωt, in order to e?ciently excite the f e mode. The f e mode corresponds to the two-dimensional wave in a vertically isothermal disc. Linear theory predicts that the radial forcing of Section 4 should lead to the outwardly propagating f e mode receiving 74.5% of the total radial energy ?ux (see Appendix B1). The inwardly propagating re mode should receive about 6.0% of 1 the total ?ux. The remainder is goes into g modes (see Appendix B1). According to linear theory (Lubow & Pringle 1993), the f e mode should occupy the entire vertical thickness of the disc and propagate outwards in radius. The mean Mach number at the mid-plane of the disc should increase as M(r, θ = π/2) ∝ r 1/2 for r ? 1. This weak dependence of the Mach number on r means that a low-amplitude wave is expected to propagate a large distance before it steepens into a shock and dissipates. The vertical structure of the wave can also be predicted by linear theory. The mean Mach number varies vertically as M ∝ exp(z 2 /5H 2 ), for γ = 5/3. Thus, although a wave may be linear near the

We turn now to the non-linear hydrodynamic calculations performed with ZEUS-2D. Figure 1 plots the magnitude of the vertically integrated radial energy ?ux, |fr | (equation 7), and the mean Mach number (equation 6) in the midplane of the disc, M(r, θ = π/2), as functions of r. They are plotted for four values of the Mach number near resonance Mr . Notice that the radial energy ?ux is in fact negative for r < 0.7. This e?ect is due the inwardly propagating re mode 1 that is excited along with the f e mode. The dotted lines in the upper panels of Figure 1 show the contribution from the f e mode to the total radial energy ?ux that is predicted by linear theory. In the lower panels, we plot the corresponding Mach numbers that are predicted in the mid-plane by linear theory. The non-linear results are in excellent agreement with the linear theory. Consider the weakest f e mode (left panels, Figure 1). Excitation of the f e mode occurs over a ?nite radial extent (a few H) around r = 1. For r > 1, the ?ux reaches a peak value just beyond the resonance. The drop in ?ux just beyond the peak is probably due to the rapid dissipation of g modes, which account for about 20% of the total ?ux. The f e mode continues to propagate to greater radii and its ?ux remains nearly constant for 2 < r < 8. The value of the ?ux there agrees well with the predicted f e mode ?ux. The dissipation that occurs at large radii is numerical as seen by the lack of convergence with increasing resolution. Thus, the wave should continue to propagate outwards, beyond the grid boundary. For r < 1, the re mode propagates inwards and overlaps 1 with the evanescent tail of the f e mode. The re mode is only 1 > modelled for r ? 0.5 because the edge of the equilibrium disc is at R1 = 0.5. The two modes have ?uxes of opposite sign, and the net ?ux changes sign near r = 0.7. The peak radial energy ?ux from the simulation in this region is about 3% of the total ?ux. This result is consistent with the expected re mode negative ?ux of about 6% of the total overlapping 1 with the tail of the positive f e mode ?ux. The r modes will be discussed in more detail in section 5.3 where an ro mode 1 receives 50% of the total ?ux. These results are reinforced by Figure 3, which gives contour plots of the wave energy density (equation 5) and the energy dissipation rate for the low-amplitude calculations (left panels). There is a lot of dissipation at high altitude in the disc in the range 1 < r < 2, as expected for the g modes. In the bulk of the disc, there is no signi?cant dissipation and the f e mode propagates outwards with virtually no loss of ?ux (a small amount is lost at high altitudes where the motions become supersonic). The behaviour of the inwardly propagating re mode can 1 also be seen in the left-hand panels of Figure 3. As predicted by linear theory (Lubow & Pringle 1993), the wave is chan-

Waves in accretion discs

7

Figure 3. The case considered here is the f e mode in a vertically isothermal disc with H/r = 0.1. The kinetic energy density in the wave (log E, 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 Mr = 0.001 and the right panels are for amplitude Mr = 0.1. At low amplitude the f e mode propagates to large radius throughout the body of the disc una?ected 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 signi?cantly attenuates the wave (see Figure 1). In both cases the grid resolution is Nr × Nθ = 1000 × 441 and qvisc = 4.0.

nelled towards the mid-plane. This will be discussed further in Section 5.3. 5.1.3 Dependence on wave amplitude

As the driving amplitude is increased, both numerical and physical dissipation occur more quickly. Numerical dissipation increases as the square of the amplitude of the wave (Section 4.3). Thus, waves with larger values of Mr propa-

gate a shorter distance before numerical dissipation occurs (cf. Figure 1, results with Mr = 0.001, and Mr = 0.01). Physical dissipation occurs if the sinusoidal wave steepens into a shock, which takes place in ? c/?c wavelengths from the resonance (Section 4.1). For waves with very small values of Mr , steepening requires many wavelengths and numerical dissipation sets in ?rst, as we saw in the case for Mr = 0.001. Physical dissipation in our simulations was obtained for cases with Mr = 0.03 and 0.1. Consider the

8

M. R. Bate et al.

Figure 4. The integrated radial energy ?ux, |fr |, (top) and mean (time-averaged) Mach number at the mid-plane of the disc, M(r, θ = π/2), (bottom) are plotted as functions of radius for the f e mode in a vertically isothermal disc with H/r = 0.05. The results are shown for three di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 111 and qvisc = 4 (long-dashed), 500 × 221 and qvisc = 4 (triple-dot-dashed), 1000 × 441 and qvisc = 4 (dot-dashed). The horizontal dotted line indicates 75 percent of the ?ux predicted by the Goldreich & 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 deposited in the f e mode. As can be seen from the ?gure, the hydrodynamic results are in excellent agreement with linear theory for the lowest amplitude case. In the high-amplitude cases, the radial energy ?ux is attenuated by dissipation.

panels for Mr = 0.03 in Figure 1. Convergence of ?ux and < M(r, θ = π/2) is obtained for r ? 4 with resolutions of 500 × 221 and 1000 × 441. With an initial amplitude of Mr = 0.1, even a resolution of 250 × 111 is su?cient to obtain physical dissipation. In the Mr = 0.03 case, the wave ?ux has been reduced by an order of magnitude by r ≈ 10, whereas in the Mr = 0.1 case this level of reduction occurs at r ≈ 4. To emphasise that this dissipation is physical rather than numerical, we plot in Figure 2 the radial velocity pro?le of the wave in the mid-plane for the cases with Mr = 0.01, 0.03, and 0.1. Steepening of the wave towards the pro?le of a shock-wave is clearly visible in those cases with Mr = 0.03 and 0.1, whereas in the Mr = 0.01 calculation the pro?le remains sinusoidal indicating only numerical, rather than physical, dissipation. This result is consistent with Figure 3 which shows that increasing Mr decreases the volume in which the f e mode propagates.

Figures 1, 4 and 5). Thus, for example, waves with Mr = 0.1 are well resolved in the hottest disc, H/r = 0.2, even with Nr × Nθ = 250 × 111, and dissipate due to shocks, whereas in the coldest disc, H/r = 0.05, even the highest resolutions do not give convincing convergence (numerical dissipation still plays a signi?cant role). Strictly, linear theory assumes H ? r. However, we ?nd that, even with H/r = 0.2, linear theory provides a good description of the wave excitation and propagation. A further prediction that can be derived from linear theory is that the ratio of the Mach number at the resonance to the Mach number at r ? rres should be nearly independent of H/r. This is con?rmed by the non-linear calculations presented here. 5.2 5.2.1 pe mode 1 Linear theory

5.1.4

Dependence on disc thickness

In Figures 4 and 5, we plot the energy ?ux and mean Mach number M(r, θ = π/2) as functions of r for discs with H/r = 0.05 and H/r = 0.2, respectively. The main difference between the waves in these discs and those in the H/r = 0.1 disc is that the wavelength is proportional to H/r. Since the numerical dissipation is ∝ 1/(λNr ) (Section 4.3), a wave in a thick disc can be followed to larger radii than in a thin disc for the same radial resolution (cf. Mr = 0.01 in

As discussed in Section 4, we adopt a driving force per unit mass that is purely in the θ-direction of the form aθ (r, θ, t) = A(r)r cos θ sin ωt, in order to e?ciently excite the pe mode. 1 The pe mode is launched at a vertical resonance r ≈ 1.39 1 and propagates outwards. We predict its radial energy ?ux in Appendix B2. No other wave is expected. Unlike the f e mode, the pe mode has vertical structure 1 (see Lubow & Pringle 1993). In the vicinity of the resonance, the wave energy has a minimum on the mid-plane of the disc and two maxima at z/H ≈ 2. This is due to the form of the excitation which depends linearly on z (Section 4). For

Waves in accretion discs

9

Figure 6. The integrated radial energy ?ux, |fr |, is plotted as a function of radius for the pe mode in a vertically isothermal disc with 1 H/r = 0.1. The results are shown for three di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre), Mr = 0.1 (right). Each case was performed with di?erent 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 ?ux predicted by linear theory (equation B20). As can be seen from the ?gure, the hydrodynamic results are in excellent agreement with linear theory for the low amplitude cases. At high amplitude, the radial energy ?ux is attenuated by dissipation.

> r ? 1.7, this changes to three maxima, one on the mid-plane and two at z/H ≈ 3, with two minima at z/H ≈ 2. Lubow & Pringle (1993) interpreted this as the wave energy rising up within the disc as the wave propagates outwards. In fact, however, the vertical distribution of the wave energy does not change beyond r ≈ 2. The wave propagates outwards maintaining the above distribution of energy inde?nitely. As with the f e mode, the wave is eventually expected to undergo 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

Figure 5. The integrated radial energy ?ux, |fr |, (top) and mean (time-averaged) Mach number at the mid-plane of the disc, M(r, θ = π/2), (bottom) are plotted as functions of radius for the f e mode in a vertically isothermal disc with H/r = 0.2. The results are shown for two di?erent wave amplitudes: Mr = 0.01 (left), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 111 and qvisc = 4 (long-dashed), 500 × 221 and qvisc = 4 (triple-dotdashed), 1000 × 441 and qvisc = 4 (dot-dashed). The horizontal dotted line indicates 75 percent of the ?ux predicted by the Goldreich & 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 deposited in the f e mode. As can be seen from the ?gure, the hydrodynamic results for the energy ?ux are in good agreement with linear theory for the lowest amplitude case, but for this relatively thick disc the predicted Mach number underestimates the computed value by about 30 percent. In the high-amplitude case, the radial energy ?ux is attenuated by dissipation.

Figure 6 plots the magnitude of the vertically integrated radial energy ?ux, |fr | (equation 7) as a function of r. We consider three values of the mean Mach number near resonance Mr . As with the case for the f e mode, the low-amplitude pe 1 mode propagates outwards from the resonance until most of the energy is dissipated numerically. < A small negative energy ?ux is seen for r ? 0.9. We attribute this to a low-amplitude r mode that is excited at the Lindblad resonance at r = 1. In a very thin disc, vertical forcing of the type applied here should not excite any modes at the Lindblad resonance. However, when H/r = 0.1, there is some ambiguity between the meanings of ‘horizontal’ and ‘vertical’ away from the mid-plane. As a result, some launching at the Lindblad resonance is possible at the level of (H/r)2 . In Figure 7 we plot the wave energy and energy dissipation rate for the low-amplitude pe mode (left panels). The 1 vertical distribution of the wave energy density is in excellent agreement with linear theory in the bulk of the disc with only a small amount of dissipation at high altitude where the motions become sonic.

5.2.3

Wave amplitude and disc thickness

As with the f e mode, the dissipation occurs more rapidly with increased driving. Calculations with Mr = 0.1 resolve 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 (log E, equation 5) is shown in the upper panels covering four orders of magnitude. The energy dissipation rate per unit volume is shown in the lower panels as a logarithmic grey-scale covering six orders of magnitude. The left panels are for amplitude Mr = 0.001 and the right panels are for amplitude Mr = 0.1. At low amplitude the pe mode propagates to large radius throughout the body of 1 the disc una?ected by the small amount of dissipation that occurs at high |z|. The pe mode displays the vertical structure predicted by 1 linear theory. In the high-amplitude case dissipation occurs in the body of the disc and signi?cantly attenuates the wave (see Figure 6). In both cases the grid resolution is Nr × Nθ = 500 × 221 and qvisc = 4.0.

occurs, but the lower Mr cases are still limited by numerical dissipation. We have not performed calculations of the pe mode 1 for discs with di?erent thicknesses since, as we found for f e mode, we do not expect any signi?cant dependence of the propagation of the pe mode on H/r. As with the f e mode, 1 the longer wavelength in hotter discs would make it easier to resolve the wave to larger radii.

5.3 5.3.1

fo /ro mode 1 Linear theory

As discussed in Section 4, we adopt a driving force per unit mass that purely in the θ-direction of the form aθ (r, t) = A(r) sin ωt, in order to e?ciently excite the f o /ro mode. 1 The f o /ro mode, or corrugation mode, is also launched at 1 a vertical resonance although, in a Keplerian disc, this co-

Waves in accretion discs

11

Figure 8. The integrated radial energy ?ux, |fr |, is plotted as a function of radius for the f o /ro (corrugation) mode in a vertically 1 isothermal disc with H/r = 0.1. The results are shown for three di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 111 and qvisc = 4 (longdashed), 500 × 221 and qvisc = 4 (triple-dot-dashed). The horizontal dotted line indicates the ?ux predicted by linear theory (equation B31). For this case, the inwardly propagating ro mode (r < 1) and the outwardly propagating f o mode (r > 1) carry equal amounts of 1 energy. As can be seen from the ?gure, the hydrodynamic results are in reasonable agreement with linear theory for the low-amplitude high-resolution case. At high amplitude, the radial energy ?ux is attenuated by dissipation.

incides with the Lindblad resonance making it a hybrid vertical/Lindblad resonance. In the discs studied here, this hybrid resonance occurs at r = 1. The f o /ro mode is able to 1 propagate throughout the entire disc with equal amounts of ?ux in the inwardly propagating ro mode and the outwardly 1 propagating f o mode. In Appendix B3, we determine the total energy ?ux produced at this resonance. The ro mode is rapidly channelled towards the mid1 plane of the disc with an associated rapid rise in the Mach number of the wave that might be expected to lead to dissipation of the wave in shocks (Lubow & Pringle 1993). The f o mode has two peaks in the vertical distribution of its wave energy with a minimum on the disc mid-plane.

cupies the bulk of the disc as it propagates to large radii, in a similar manner to the f e mode except that it has vertical structure. As expected from linear theory (Lubow & Pringle 1993), the inwardly propagating ro mode is strongly focused to1 wards the mid-plane of the disc (Figure 9). Both linear theory and the non-linear results also show that the wavelength of the ro mode decreases towards the centre. 1

5.3.3

Wave amplitude

5.3.2

Low-amplitude, non-linear results

Figure 8 plots the magnitude of the vertically integrated radial energy ?ux, |fr | (equation 7) as a function of r. We consider three values of the mean Mach number near resonance Mr . As predicted by linear theory, this mode propagates throughout the entire disc. Outward of the resonance (r > 1), the wave propagates as an f o mode. Inward of the resonance, the wave propagates as an ro mode. An equal 1 amount of ?ux goes into each of the modes, and the magnitude of the ?ux is in good agreement with linear theory, particularly for the calculations with the highest resolution. Once again, for low-amplitude driving, the distance that the waves propagate is limited by numerical dissipation rather than physical dissipation (shocks). The energy of the f o mode is concentrated away from the mid-plane of the disc (Figure 9) and its vertical distribution is in excellent agreement with linear theory for < |z/H| ? 3. Above this region and particularly within the < range 1 ≤ r ? 2.5, shocks form and a small amount of energy is dissipated. However, the majority of the wave en< ergy remains within |z/H| ? 3 and propagates outwards inde?nitely. Lubow & Pringle (1993) noted the rise of wave energy away from the mid-plane that occurs from near the resonance to r = 3 and interpreted this as the wave energy rising into the atmosphere. However, we see here that this ‘rising’ does not continue inde?nitely; the wave simply oc-

With strong driving (e.g. Mr = 0.1), physical dissipation through shocks occurs in the simulations and those with di?erent resolution converge, at least for the f o mode. By the time the f o mode reaches r = 10, the ?ux has been reduced by a factor of ? 103 . The wave is con?ned to a region very near the mid-plane of the disc; shocks are present at large |z| (Figure 9, right panels). We have not performed calculations of the corrugation mode with other disc thicknesses.

6

THE POLYTROPIC DISC

We now turn our attention to polytropic discs. The disc studied in this section has an optical depth τ = 25000, and only a tenuous isothermal atmosphere. Discs that have τ = 5 ? 100 (i.e. intermediate between a vertically isothermal disc and a very optically thick disc) are brie?y described in Section 7.

6.1 6.1.1

fe mode Linear theory

As discussed in Section 4, we adopt a purely radial driving force per unit mass of ar (r, t) = A(r) sin ωt, in order to e?ciently excite the f e mode. The f e mode corresponds to the two-dimensional wave in a vertically isothermal disc. Linear theory predicts that this radial forcing should lead to the outwardly propagating f e mode receiving 98.3% of the total radial energy ?ux (see Appendix B1). A greater fraction

12

M. R. Bate et al.

Figure 9. The case considered here is the f o /ro (corrugation) mode in a vertically isothermal disc with H/r = 0.1. The kinetic energy 1 density in the wave (log E, equation 5) is shown in the upper panels covering ?ve orders of magnitude. The energy dissipation rate per unit volume is shown in the lower panels as a logarithmic grey-scale covering six orders of magnitude. The left panels are for amplitude Mr = 0.001 and the right panels are for amplitude Mr = 0.1. At low amplitude the f o mode propagates outwards from r = 1 to large radius throughout the body of the disc una?ected by the small amount of dissipation that occurs at high |z|. The ro mode propagates 1 inwards from r = 1 and is channelled towards the mid-plane of the disc. Both the f o mode and the ro mode display the vertical structure 1 predicted by linear theory. In the high-amplitude case dissipation occurs in the body of the disc and signi?cantly attenuates the wave (see Figure 8). In both cases the grid resolution is Nr × Nθ = 500 × 221 and qvisc = 4.0.

goes into the f e mode in the polytropic disc than in the vertically isothermal disc. The remainder goes into an inwardly propagating re mode. 1 The linear behaviour of the f e mode in a polytropic disc has been studied by Korycansky & Pringle (1995), Lubow & Ogilvie (1998) and Ogilvie & Lubow (1999). Although the f e mode is excited throughout the entire vertical extent of the

disc at the resonance, it is channelled towards the surface of the polytropic disc, unlike in the vertically isothermal disc. The rate of channelling (the distance from the resonance at which the wave energy becomes concentrated away from the disc mid-plane) depends only on the degree of vertical strati?cation, which is determined by the index of the polytrope, and the azimuthal wavenumber m, which is ?xed at m = 0

Waves in accretion discs
for the axisymmetric case. The rate is independent of the disc thickness H, in the limit that H ? r. This concentration of the wave energy into a small volume of cool, lowdensity gas near the disc surface (or atmosphere) leads to a rapid increase of the Mach number of the wave with radius. Thus we expect that waves in polytropic discs will undergo shocks more quickly (closer to the resonance) than waves in the isothermal case described earlier. Lubow & Ogilvie (1998) speculated that once the Mach number of the wave became of order unity, the wave would shock and thus, such waves would not propagate over large distances. One of the aims of the non-linear hydrodynamic calculations is to determine where this dissipation occurs. 6.1.2 Low-amplitude, non-linear results
H/r Mr =0.01 Mr =0.03 Mr =0.1 0.05 1.7 1.5 1.4 0.1 1.8 1.4 1.4 0.2 1.9+ 1.5 1.4

13

Table 1. The radius r at which the f e mode dissipates is tabulated for polytropic discs in which the calculations resolve the shock dissipation. Values are given for di?erent disc thicknesses, H/r, and for di?erent initial wave amplitudes Mr . The ‘+’ indicates that this value is a lower limit. Note that the radius of dissipation is almost independent of the disc thickness as predicted by the theory of wave channelling.

We turn now to the non-linear hydrodynamic calculations performed with ZEUS-2D. Figure 10 plots the magnitude of the vertically integrated radial energy ?ux, |fr | (equation 7) as a function of r. We consider four values of the Mach number near resonance Mr . The calculations with Mr = 0.001 do not have su?cient resolution to follow the f e mode (r > 1) until physical dissipation occurs. This can be seen by the absence of convergence of the ?ux in Figure 10. Because of the higher concentration of wave energy compared with the vertically isothermal case, it is more di?cult to provide adequate resolution. However, it is clear that the wave propagates to at least r = 2.3 with little dissipation, over a distance that is many times greater than the disc thickness H. Notice that the radial energy ?ux is negative for r < 0.7. This e?ect is due the inwardly propagating re mode that is 1 excited along with the f e mode. As predicted by linear theory, most of the ?ux goes into an outwardly propagating f e mode. The magnitude of this ?ux is indistinguishable from the linear prediction. We ?nd ≈ 1.2% goes into an inwardly propagating re mode (Figure 10, Mr = 0.001), similar to 1 the linear prediction. The vertical thermal strati?cation of the polytropic disc has a strong e?ect on the distribution of wave energy of the f e and re modes compared to the isothermal disc (Figure 1 11, top left). As predicted by linear theory, instead of occupying the entire thickness of the disc, as in the isothermal case, the energy in the f e mode is channelled towards the surface of the disc. It e?ectively becomes a wave that travels along the surface of the disc. For the re mode, rather 1 than being channelled towards the mid-plane, it continues to occupy the entire thickness of the disc as it propagates towards the central object. This behaviour di?ers from the vertically isothermal case, where the re mode is channelled 1 towards the mid-plane. The reason for this di?erence is that the channelling of the re mode is due to vertical buoyancy 1 (Lubow & Pringle 1993). The di?erence in behaviours is due to the vanishing buoyancy frequency throughout the vertically polytropic disc, unlike the vertically isothermal case. The wave energy and dissipation for our weakest driving calculation, Mr = 0.001, are shown in Figure 11 (left panels). The dissipation is only suggestive, since numerical convergence has not occurred at larger radii. However, the channelling of the f e mode to the surface of the disc is obvious. The wave dissipates at r ≈ 2.3 in the highest resolution calculation (Nr × Nθ = 1156 × 365) and at r ≈ 2.2 with half of this resolution (Nr × Nθ = 1000 × 183).

The low-amplitude re mode propagates towards the cen1 tre of the disc (Figure 11 (left panels). The zig-zag pattern that can be seen in the distribution of wave energy is due to the superposition of the re mode (for which the wave 1 energy has a minimum on the mid-plane) and the evanescent tail of the f e mode. The wavelength of the re mode is 1 equal to the wavelength of this zig-zag pattern and decreases rapidly with decreasing radius. Linear theory predicts that the wavelength should be ∝ r 5/2 and this is con?rmed by the non-linear calculations. As the wave propagates towards the centre of the disc it dissipates numerically because the wave becomes unresolved (the size of the grid zones decreases logarithmically, while the wavelength decreases more rapidly). 6.1.3 Dependence on wave amplitude

As described above, the resolution of the calculations with Mr = 0.001 is not su?cient to resolve physical dissipation of the f e mode due to shocks. However, calculations with Mr ≥ 0.01 are able to resolve the f e mode until physical dissipation occurs. Considering the convergence of |fr | in Figure 10, we ?nd that physical dissipation is marginally resolved in the Mr = 0.01 calculations with Nr × Nθ = 1156 × 365 (there is very little di?erence in the radius at which |fr | rapidly decreases between this calculation and that with Nr × Nθ = 1000 × 183). Calculations with Mr ≥ 0.03 are easily resolved with Nr ×Nθ = 1000×183. The radii to which these waves propagate before physical dissipation occurs are given in Table 1. Compared with the case of the vertically isothermal disc, the wave propagates a shorter distance before dissipating, although still much greater than H. The increased non-linearity is due to the increase in the magnitude of velocity perturbations and the decrease in local sound speed that occurs as the wave is channelled into the small region near the disc surface. In Figure 12, we plot the height above the mid-plane of the peak of the wave energy density as a function of radius and compare it with the prediction from linear theory. We also plot the mean Mach number of the wave along this peak in wave energy M(r, θ) and compare this to linear theory. > Away from the resonance (r ? 1.3), where the linear theory should give an accurate description, the agreement between the linear theory and the non-linear hydrodynamic calculations is excellent for both the location of the peak in the wave energy density and the Mach number of the wave along this maximum. It appears that shock dissipation occurs once the mean Mach number M exceeds approximately 0.15 (i.e.

14

M. R. Bate et al.

Figure 10. The integrated radial energy ?ux, |fr |, is plotted as a function of radius for the f e mode in a polytropic disc with H/r = 0.1 and τ = 25000. The results are shown for four di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre-left), Mr = 0.03 (centre-right), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 45 and qvisc = 4 (long-dashed), 500 × 91 and qvisc = 4 (triple-dot-dashed), 1000 × 183 and qvisc = 4 (dot-dashed), 1156 × 365 and qvisc = 4 (solid). The horizontal dotted line indicates 98 percent of the ?ux predicted by the Goldreich & Tremaine formula (equation B1). The value of 98 percent is that predicted by linear theory to be deposited in the f e mode. As can be seen from the ?gure, the hydrodynamic results are in excellent agreement with linear theory for the low-amplitude cases. In the high-amplitude cases, the radial energy ?ux is attenuated by dissipation.

Figure 11. The case considered here is the f e mode in a polytropic disc with H/r = 0.1 and τ = 25000. The kinetic energy density in the wave (log E, 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 four orders of magnitude. The left panels are for amplitude Mr = 0.001 and the right panels are for amplitude Mr = 0.1. At low amplitude the f e mode propagates outwards from r = 1 and becomes increasingly channelled to z ≈ H. The dissipation occurs at r ≈ 2.2 as the wave becomes non-linear (upper panel, Figure 13). In the high-amplitude case, wave channelling of the f e mode still occurs but the high amplitude results in dissipation closer to the resonance. In both cases, the pattern visible at r < 1 is generated by a superposition of the evanescent tail of the f e mode and a low-amplitude inwardly propagating r mode. In both cases the grid resolution is Nr × Nθ = 1000 × 183 and qvisc = 4.0.

Waves in accretion discs

15

Figure 12. Comparison of the non-linear fe mode with linear theory for the polytropic disc with H/r = 0.1 and τ = 25000. Top panel: We plot as a function of radius, the height, z, in the disc where the wave energy density reaches a maximum. The linear prediction is given by the dotted line and the low-amplitude Mr = 0.001 computation is given by the solid line. For higher amplitude waves the maxima lie roughly along the same curve, but attempting to extract them from the calculations results in curves that are so noisy that plotting them makes the ?gure unintelligible. Bottom panel: The time-averaged wave Mach number at the location of the peak in wave energy density (as de?ned in the upper panel) is plotted as a function of radius. Four cases are shown: non-linear calculations with (lower to upper solid lines) Mr = 0.001, Mr = 0.01, Mr = 0.03, Mr = 0.1 and the corresponding linear predictions. The numerical results display excellent agreement with linear theory both for the location of the maximum wave energy and for the Mach number of the wave at these maxima, until the wave dissipates in shocks (Mr = 0.01, Mr = 0.03, Mr = 0.1) or numerical dissipation (Mr = 0.001).

Figure 13. The Mach number of the radial motion of the wave at the height in the disc at which the wave energy density peaks (see Figure 12) is plotted at a particular instant as a function of radius for the fe mode in the polytropic disc with H/r = 0.1 and τ = 25000. Various wave amplitudes are shown: Mr = 0.001 (top), Mr = 0.01 (middle) and Mr = 0.03 (bottom). For the lowest amplitude case, Mr = 0.001, the wave maintains its sinusoidal pro?le, but with increasing amplitude as it is channelled towards the surface, until numerical dissipation occurs. For the higher amplitudes, Mr = 0.01 and Mr = 0.03 the wave becomes nonlinear very rapidly, shocks and dissipates.

6.1.4

Dependence on disc thickness

the peak-to-trough mean Mach number exceeds ≈ 0.4). Note that the Mach number reached by the Mr = 0.001 calculation at r = 2.5 is close to this value, implying that even the lowest-amplitude wave may be marginally resolved with the resolution of Nr × Nθ = 1156 × 365. Finally, we note that the e?ciency of the dissipation depends on Mr , in that stronger waves result in signi?cant ?ux beyond the dissipation radius (cf. Mr = 0.01 and Mr = 0.1 in Figure 10). This might be because the strong wave dissipation modi?es the equilibrium disc structure in a manner that makes non-linear dissipation more di?cult.

We have seen that in a polytropic disc the f e mode is channelled towards the surface of the disc and dissipates due to wave steepening when the mean Mach number at the location where the wave energy is concentrated exceeds M ≈ 0.15. We examine here how this behaviour depends on the thickness of the disc. Recall also that that the linear theory has been developed under a thin-disc approximation. Figures 14 and 15 correspond to Figure 10, but for values of H/r equal to 0.05 and 0.2 respectively. In Table 1, we give the radius of dissipation of the f e mode for those calculations in which the wave dissipates due to shocks (rather than numerically). Although H/r is varied by a factor of 4, the radius at which the wave dissipates is almost independent of H/r. This result is consistent with the expectations of wave channelling, since channelling occurs over a distance that is independent of H/r (Lubow & Ogilvie 1998). The main di?erence between discs of di?erent thickness is that thinner discs seem to result in more e?cient dissipation of the wave. For example, consider the calculations with Mr = 0.01. The radial energy ?ux drops by nearly a factor of ? 1000 between r = 1.7 and r = 2.5 for the H/r = 0.05 disc, whereas it only drops by a factor of ≈ 10 for H/r = 0.1

16

M. R. Bate et al.

Figure 14. The integrated radial energy ?ux, |fr |, is plotted as a function of radius for the f e mode in a polytropic disc with H/r = 0.05 and τ = 25000. The results are shown for four di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre-left), Mr = 0.03 (centre-right), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 45 and qvisc = 4 (long-dashed), 500 × 91 and qvisc = 4 (triple-dot-dashed), 1000 × 183 and qvisc = 4 (dot-dashed), 1156 × 365 and qvisc = 4 (solid). The horizontal dotted line indicates 98 percent of the ?ux predicted by the Goldreich & Tremaine formula (equation B1). The value of 98 percent is that predicted by linear theory to be deposited in the f e mode. As can be seen from the ?gure, the hydrodynamic results are in excellent agreement with linear theory for the low-amplitude cases. In the high-amplitude cases, the radial energy ?ux is attenuated by dissipation.

Figure 15. The integrated radial energy ?ux, |fr |, is plotted as a function of radius for the f e mode in a polytropic disc with H/r = 0.2 and τ = 25000. The results are shown for three di?erent wave amplitudes: Mr = 0.01 (left), Mr = 0.03 (centre), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 45 and qvisc = 4 (long-dashed), 500 × 91 and qvisc = 4 (triple-dot-dashed), 1000 × 183 and qvisc = 4 (dot-dashed), 1156 × 365 and qvisc = 4 (solid). The horizontal dotted line indicates 98 percent of the ?ux predicted by the Goldreich & Tremaine formula (equation B1). The value of 98 percent is that predicted by linear theory to be deposited in the f e mode. As can be seen from the ?gure, the hydrodynamic results are in reasonable agreement with linear theory for the most highly resolved low-amplitude case. In the high-amplitude cases, the radial energy ?ux is attenuated by dissipation.

and only a factor of ≈ 3 for H/r = 0.2. The reason for this dependence of the dissipation e?ciency on the disc thickness is not clear. In summary, linear theory, although derived for thin discs, gives a good description of the wave propagation even in discs as hot as H/r = 0.2, although the dissipation of the waves may be less e?cient for thicker discs.

6.2.2

Low-amplitude, non-linear results

6.2 6.2.1

pe mode 1 Linear theory

As discussed in Section 4, we adopt a driving force per unit mass that is purely in the θ-direction of the form aθ (r, θ, t) = A(r)r cos θ sin ωt, in order to e?ciently excite the pe mode. 1 As described earlier, the pe mode is launched at a vertical 1 resonance at r = 1.39 and propagates outwards. No other wave is expected and the radial energy ?ux of the pe mode is 1 predicted in Appendix B2. In a vertically polytropic disc, the pe mode energy is expected to channel to the disc surface. 1

Figure 16 plots the magnitude of the vertically integrated radial energy ?ux, |fr | (equation 7) as a function of r. We consider three values of the mean Mach number near resonance Mr . The wave excitation is similar to that in the isothermal disc with most of the energy going into an outwardly propagating pe mode which is excited at r = 1.39 1 and a small fraction producing an inwardly propagating re 1 mode which is excited at r = 1. We ?nd that the ratio of the two ?uxes is ≈ 1.1% (Figure 16, Mr = 0.01). As noted previously, this e?ect is of order (H/r)2 and should vanish in the limit of a thin disc. In the isothermal disc, the peak of wave energy of the pe mode in the vicinity of the resonance (r = 1.39) occurs 1 < < o? the mid-plane for 1.0 ? r ? 1.7. (There is one maximum above, and one below the mid-plane). This structure changes as the wave propagates to larger radii to a peak on > the mid-plane and two o?-mid-plane peaks for r ? 1.7. In the polytropic disc (Figure 17, top left), this general struc-

Waves in accretion discs

17

Figure 16. The integrated radial energy ?ux, |fr |, is plotted as a function of radius for the pe mode in a polytropic disc with H/r = 0.1 1 and τ = 25000. The results are shown for three di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 45 and qvisc = 4 (long-dashed), 500 × 91 and qvisc = 4 (triple-dot-dashed), 1000 × 183 and qvisc = 4 (dot-dashed). The horizontal dotted line indicates the ?ux predicted by linear theory (equation B20). As can be seen from the ?gure, the hydrodynamic results are in excellent agreement with linear theory for the low-amplitude cases. At high amplitude, the radial energy ?ux is attenuated by dissipation.

Figure 17. The case considered here is the pe mode in a polytropic disc with H/r = 0.1 and τ = 25000. The kinetic energy density 1 in the wave (log E, equation 5) is shown in the upper panels covering three orders of magnitude. The energy dissipation rate per unit volume is shown in the lower panels as a logarithmic grey-scale covering four orders of magnitude. The left panels are for amplitude Mr = 0.001 and the right panels are for amplitude Mr = 0.1. At low amplitude the pe mode propagates outwards from r ≈ 1.39 and 1 becomes increasingly channelled to z ≈ H. The mode displays vertical structure as predicted by linear theory (cf. the vertically isothermal disc, Figure 7). The dissipation occurs at r ≈ 2.5 as the wave becomes non-linear. In the high-amplitude case, wave channelling of the pe mode still occurs but the high amplitude results in dissipation closer to the resonance. In both cases, the pattern visible at r < 1 is 1 generated by a superposition of the evanescent tail of the f e mode and a low-amplitude inwardly propagating r mode. In both cases the grid resolution is Nr × Nθ = 1000 × 183 and qvisc = 4.0.

18

M. R. Bate et al.

Figure 18. The integrated radial energy ?ux, |fr |, (top) and mean (time-averaged) Mach number at the mid-plane of the disc, M(r, θ = π/2), (bottom) are plotted as functions of radius for the f o /ro (corrugation) mode in a polytropic disc with H/r = 0.1 and 1 τ = 25000. The results are shown for three di?erent wave amplitudes: Mr = 0.001 (left), Mr = 0.01 (centre), Mr = 0.1 (right). Each case was performed with di?erent numerical resolutions to test for convergence: 250 × 45 and qvisc = 4 (long-dashed), 500 × 91 and qvisc = 4 (triple-dot-dashed), 1000 × 183 and qvisc = 4 (dot-dashed). The horizontal dotted line indicates the ?ux predicted by linear theory (equation B31). For this case, the inwardly propagating ro mode (r < 1) and the outwardly propagating f o mode (r > 1) carry 1 equal amounts of energy. As can be seen from the ?gure, the hydrodynamic results are in reasonable agreement with linear theory for the low-amplitude high-resolution case. At high amplitude, the radial energy ?ux is attenuated by dissipation. It can be seen that the Mach number of the ro mode increases as it propagates inwards caused by the rapid decrease of the group velocity (or radial wavelength, 1 Figure 18).

ture is repeated, but with the added complexity of the wave simultaneously being channelled to the disc surface in a similar manner to the f e mode. Thus, in the vicinity of the res< < onance (1.0 ? r ? 1.7) the wave energy has a minimum > on the mid-plane as in the isothermal disc. At r ? 1.7, the two o?-mid-plane peaks which were present in the isothermal disc are near the surface of the polytropic disc. The peak which was on the mid-plane in the isothermal disc is < < also present in the polytropic disc for 1.7 ? r ? 2.1, but is > rapidly channelled towards the surface of the disc for r ? 2.1 in a manner which is similar to the channelling of the f e > mode. For r ? 2.5, there is very little wave energy near the mid-plane of the disc. The re mode, which is excited at r = 1 and propagates 1 inwards, behaves in the same way as the re mode which 1 was excited by the forcing for the f e mode. It occupies the entire thickness of the disc and shows no sign of physical dissipation as it propagates to the centre. It is dissipated numerically because of the rapid decrease of wavelength.

the f e mode (c.f. Figure 10) because the wave is excited at r = 1.39 instead of r = 1. 6.3 6.3.1 fo /ro mode 1 Linear theory

As discussed in Section 4, we adopt a driving force per unit mass that is purely in the θ-direction of the form aθ (r, t) = A(r) sin ωt, in order to e?ciently excite the f o /ro mode. The 1 general properties of the waves were described in Section 5.3 for the vertically isothermal case. In the polytropic case, the primary new feature is that wave channelling occurs. 6.3.2 Low-amplitude, non-linear results

6.2.3

Wave amplitude

We only have su?cient resolution to ?nd the radius at which physical dissipation occurs for the case with Mr = 0.1 (Figure 16). However, in cases with the same radial resolutions, the radius of dissipation is larger for the pe mode than for 1

Figure 18 plots the magnitude of the vertically integrated radial energy ?ux, |fr | (equation 7), and the mean Mach number (equation 6) in the mid-plane of the disc, M(r, θ = π/2), as functions of r. They are plotted for three values of the mean Mach number near resonance Mr . The wave excitation is similar to that in the isothermal disc with 50% of the ?ux going into each of the inwardly and outwardly propagating modes. In the isothermal disc, the wave energy of the f o mode peaked o? the mid-plane. In the polytropic disc (Figure 19, < < top left), this is apparent near the resonance (1 ? r ? 1.3), > 1.3, channelling of the wave to the disc surface but for r ?

Waves in accretion discs

19

Figure 19. The case considered here is the f o /ro (corrugation) mode in a polytropic disc with H/r = 0.1 and τ = 25000. The kinetic 1 energy density in the wave (log E, 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 four orders of magnitude. The left panels are for amplitude Mr = 0.001 and the right panels are for amplitude Mr = 0.1. At low amplitude the f o mode propagates outwards from r = 1 and becomes increasingly channelled to z ≈ H. Dissipation of the f o mode occurs at r ≈ 2.2 as the wave becomes non-linear. In the high-amplitude case, wave channelling of the f o mode still occurs but the high amplitude results in dissipation closer to the resonance. In both cases, the pattern visible at r < 1 is generated by a superposition of the evanescent tail of the f o mode and the ro mode that propagates inwards from 1 r = 1. Note that as the ro mode propagates inwards its radial wavelength (and hence its group velocity) decreases. In both cases the 1 grid resolution is Nr × Nθ = 1000 × 183 and qvisc = 4.0.

dominates the distribution of wave energy. In fact, apart from the minimum in the wave energy on the mid-plane < < for 1 ? r ? 1.3, the distribution of wave energy is almost identical to that of the f e mode (cf. Figure 11, top left). The ro mode, which is excited at r = 1 and propagates 1 inwards, occupies the entire thickness of the disc as did the re mode which was excited by the forcing for the f e and 1 pe modes. The mean Mach number at mid-plane increases 1 rapidly with decreasing radius for the ro mode (Figure 18, 1 lower panels). In the lowest amplitude, highest resolution calculation, it increases by a factor of ≈ 4 in going from the resonant radius r = 1 to radius r = 0.3. In reality, this would result in the formation of shocks and dissipation of the wave as it approaches the centre of the disc. However, in the calculations presented here, we do not have su?cient resolution to resolve physical dissipation as can be seen by the lack of convergence of |fr | with increasing resolution (Figure 18). Instead, the ro mode is dissipated numerically because of 1 the rapid decrease of wavelength.

6.3.3

Wave amplitude

For a given wave amplitude, it appears that the f o mode dissipates at about the same radius as the f e mode. This con?rms what was found in comparing the distance of propagation of the f e and pe modes, namely that the channelling 1 of the waves towards the disc surface and their resulting radius of dissipation occurs at about the same distance from resonance. This result is probably a consequence of the fact that wave channelling occurs when kH is greater than unity (Lubow & Ogilvie 1998). This condition is determined by the wave dispersion relations and occurs at roughly the same distance from resonance (namely ? r) for the modes considered here.

7

DISCS WITH τ = 5 ? 100

Thus far, we have considered vertically isothermal discs, with optical depth τ ? 0, or polytropic discs with tenuous

20

M. R. Bate et al.

Figure 20. We consider the f e mode in polytropic discs with H/r = 0.1 and varying optical depth. From top to bottom: τ = 5, τ = 10, τ = 30, and τ = 100. In each case, the amplitude is Mr = 0.01. The kinetic energy density in the wave (log E, equation 5) is shown in the left-hand panels. The energy dissipation rate per unit volume is shown in the middle panels as a logarithmic grey-scale. The grey-scales cover four orders of magnitude. The right panels give the integrated radial energy ?ux as functions of radius. The horizontal dotted lines indicate the ?ux predicted by linear theory for the f e mode. As can be seen from the ?gure, the predicted ?uxes for discs of intermediate optical depth are in excellent agreement with linear theory. As the optical depth increases, there is a gradual change from the vertically isothermal behaviour (i.e. the wave occupying the entire thickness of the disc, top) to wave channelling (bottom). This is predicted by linear theory (Ogilvie & Lubow 1999). As wave channelling becomes more important, the attenuation of the wave increases. In each case the grid resolution is Nr × Nθ = 500 × 183 and the disc is modelled over 6 vertical scale-heights.

isothermal atmospheres and τ ≈ 25000. Wave propagation and dissipation are quite di?erent for these two types of discs. In this section we brie?y examine cases of intermediate optical depth, as de?ned by equation 2. This has been studied in linear theory by Ogilvie & Lubow (1999). They found that as the optical depth of the disc is decreased, outwardly propagating waves are channelled into thicker layers. This results in a milder increase of the Mach number as the wave propagates outwards. In this section, we consider only the f e mode. We perform calculations with discs that have τ = 5, 10, 30 and 100.

Their individual parameters were chosen so that the thickness of the polytropic layer relative to r is the same for each disc (H/r ≈ 0.1). Each calculation had Mr = 0.01. The results are provided in Figure 20. As predicted by linear theory, there is a gradual change from the vertically isothermal behaviour (i.e. the wave occupying the entire thickness of the disc, top) to wave channelling (bottom) as the optical depth is increased. For discs of increasing optical depth, a greater fraction of the f e mode’s energy is channelled to the surface of the disc where it dissipates. Thus, the fraction of the f e mode ?ux that propagates past r ≈ 2

Waves in accretion discs
decreases (right panels). The ?ux in the τ = 10 disc is only reduced by a factor of ≈ 2, while the f e mode loses more > than 80% of its ?ux if τ ? 100. These results demonstrate that wave channelling is important even in discs of quite low optical depth.

21

8 8.1

CONCLUSIONS AND DISCUSSION Comparison with linear theory

We have performed non-linear numerical simulations of the excitation, propagation and dissipation of axisymmetric hydrodynamic waves in an accretion disc. The disc is Keplerian, vertically resolved and inviscid apart from an arti?cial numerical bulk viscosity that is required in order to resolve shocks. The waves are excited by applying a temporally periodic acceleration throughout the computational domain, and are launched resonantly at radii where the forcing frequency coincides with a natural oscillation frequency of the disc. We have examined Lindblad resonances, vertical resonances and hybrid vertical/Lindblad resonances. In each case the waves propagate radially through the disc away from the resonance, changing continuously in vertical structure as they do so, and ultimately dissipate. This mechanism of generating the waves is directly comparable with the way that non-axisymmetric waves are excited in discs by tidal forcing from an orbiting companion in a binary or protoplanetary system. A low-mass companion such as a sub-Jovian planet orbiting a star exerts its tidal in?uence on the disc mainly through launching f e -mode waves at the Lindblad resonances (Goldreich & Tremaine 1980). Lindblad resonances are also important in systems of larger mass ratio, such as circular orbit binary star systems of extreme mass ratio (Lin & Papaloizou 1979) and eccentric orbit binaries (Artymowicz & Lubow 1994). Related resonances play a role in the growth of the eccentric disc mode in superhump binaries (Lubow 1991) and possibly binaries with brown dwarf secondaries (Papaloizou, Nelson & Masset 2001). Vertical resonances are an important aspect of the tidal interaction in close binary stars (Lubow 1981; Ogilvie 2002), where the Lindblad resonances are all excluded from the disc. The corrugation-mode resonances we have studied are the axisymmetric equivalent of the bending-mode resonances observed in Saturn’s rings (Shu, Cuzzi & Lissauer 1983). Wherever the waves are linear, our results are in excellent and detailed agreement with the linear theory developed by Lubow & Pringle (1993), Korycansky & Pringle (1995), Ogilvie (1998), Lubow & Ogilvie (1998) and Ogilvie & Lubow (1999). Elements of this theory can also be found in the work of Loska (1986) and Okazaki, Kato & Fukue (1987). The disc acts as a waveguide that allows a variety of wave modes to propagate in the radial direction. The modes can be classi?ed similarly to stellar oscillations, and we have discussed all four classes (f, p, r and g) to some extent in this paper. Each mode has a dispersion relation that connects the wave frequency and the radial wavenumber at every radius. In general, waves of a given frequency can propagate only in restricted radial intervals, bounded by turning points where the wavenumber vanishes. The turning points correspond to the resonant radii where waves are launched by an applied periodic forcing.

The linear theory is con?rmed by the simulations in a number of respects. First, the concept of the disc as a waveguide that supports discrete propagating modes is validated. The waves are vertically con?ned by the continuous strati?cation of the disc, not by zero-density surfaces, arti?cial boundaries or other discontinuities. The solutions exhibit the same ‘separation of variables’ that also features in the linear analysis. That is, the waveform is an oscillatory function of radius (and time) multiplied by a vertical pro?le that changes slowly with radius. The fact that the velocity amplitudes of the linear solutions may formally diverge high in the atmosphere does not invalidate the rest of the solution when the energy density of the wave is vertically con?ned in the bulk of the disc. Secondly, the energy ?uxes imparted to the di?erent wave modes through launching at di?erent types of resonance are in excellent agreement with the predictions of linear theory presented in Appendix B. This analysis of resonant wave launching is closely related to the well-known resonant torque formula of Goldreich & Tremaine (1979) and to the analysis of vertical resonances by Lubow (1981). Finally, the detailed waveforms predicted by linear theory are con?rmed in the simulations (see, e.g., Figure 21 and the next subsection), even when the disc is not especially thin. 8.2 Wave channelling versus refraction

The continuous change of the vertical pro?le of a wave as it propagates radially is a process we have termed ‘wave channelling’ (Lubow & Ogilvie 1998). The energy density of the wave can be channelled either towards the surfaces of the disc or towards the mid-plane as it propagates away from its site of launching. The results of this paper provide ample evidence for both types of wave channelling. In particular, the f and p modes in a disc with a vertical temperature gradient are channelled towards the surfaces (e.g. Figures 11, 17 and 19), while the r modes in a disc with a vertical entropy gradient are channelled towards the mid-plane (e.g. Figures 3 and 9) as predicted by Lubow & Pringle (1993). The wave energy of the f e mode, which is the principal mode launched at a Lindblad resonance, rises towards the surfaces of a thermally strati?ed disc. Close to the resonance, the f e mode occupies the full vertical extent of the disc and behaves compressibly. As the mode propagates away from the resonance, its energy within the disc midplane region rises and becomes con?ned (or channelled) to a layer at the base of the disc atmosphere. The extent of channelling depends on the degree of thermal strati?cation. In this regime, the wave becomes a surface gravity wave and behaves incompressibly (Ogilvie 1998; Lubow & Ogilvie 1998). The process by which this occurs is wave channelling and not acoustic refraction. This is demonstrated as follows. If refraction were responsible (e.g. Lin et al. 1990a,b), the wave would be directed upwards into the isothermal atmosphere after travelling a radial distance of order H. This is not observed in the simulations. Plots of the wave energy for low-amplitude waves in polytropic discs (Figures 11, 17 and 19) clearly show that the waves propagate along the base of the atmosphere before dissipating in shocks. Propagation does not switch from horizontal to vertical, as would be expected if simple refraction were involved. In fact, the f e mode is launched even in an incompressible ?uid with a

22

M. R. Bate et al.

Figure 21. Comparison of the simulations with linear theory, and comparison of wave channelling with acoustic refraction. In the upper √ grey-scale ?gure, we plot a snapshot of the radial velocity multiplied by the square root of the density, vr ρ, from a simulation involving e mode in a polytropic disc that joins smoothly on to a low-mass isothermal atmosphere (H/r = 0.1, τ = 25000). (The square of an f this quantity is the radial part of the instantaneous wave energy density, c.f. equation 5). The amplitude of the wave is Mr = 0.001 and the grid resolution is Nr × Nθ = 1000 × 183. In the lower grey-scale ?gure, we plot the corresponding semi-analytical prediction of the linear theory of wave channelling, ignoring the r mode and using the WKB approximation in the radial direction (which results in a mild singularity at the Lindblad resonance). The agreement is excellent until dissipation occurs. Over the upper grey-scale, we trace the rays that would result from the refraction of an initially vertical wavefront at the resonant radius. The refracted wavefronts, which are orthogonal to the rays, would rapidly become severely tilted and the wave would be predicted to propagate almost vertically into the isothermal atmosphere. However, the simulation shows that the wavefronts in fact remain nearly vertical as the wave propagates horizontally along the base of the disc atmosphere (shown by the dashed line, z = H). The f e mode is generated throughout the thickness of the disc at r = 1 but its energy rises from within the polytropic layer owing to wave channelling and becomes con?ned near the base of the atmosphere. The wave propagates without loss of ?ux to r ≈ 2.2 (Figure 10, dot-dashed line, left panel) before most of its energy is dissipated in two regions at the base of the atmosphere (white contour lines at r ≈ 2.2). The dissipation in this example is partly numerical, as suggested by Figure 10. If we were able to increase the resolution further, the wave would propagate a greater distance (c.f. di?erent resolutions in the left panel of Figure 10) and might agree even more closely with the linear prediction.

Waves in accretion discs
vertical energy distribution that is similar to the compressible case. Consequently, the process of energy concentration cannot be due to acoustic refraction. Figure 21 shows an explicit comparison between the results of the simulations and the predictions of acoustic refraction. In order to show both the wavefronts and the vertical extent of the wave, we plot a snapshot of the radial velocity, normalized by the RMS velocity, and multiplied by the mean wave energy density, from a simulation involving an f e mode in a polytropic disc. On top of the grey-scale, we trace the rays that would result from the refraction of an initially vertical wavefront at the resonant radius. The refracted wavefronts, which are orthogonal to the rays, would rapidly become severely tilted and the wave would be predicted to propagate almost vertically into the isothermal atmosphere. However, the simulations show that the wavefronts in fact remain nearly vertical as the wave propagates horizontally along the base of the disc atmosphere. The wave propagates without loss of ?ux to r ≈ 2.2 before most of its energy is dissipated in two regions at the base of the atmosphere (see Figure 10 for the radial energy ?ux). Figure 21 demonstrates that the linear theory of wave channelling accurately predicts the propagation of the wave until nonlinear dissipation occurs. In summary, our non-linear hydrodynamic calculations clearly show that for axisymmetric waves the disc acts as a waveguide and that the behaviour of the waves is determined primarily by wave channelling and not refraction. Explanations based on simple refraction are incorrect. Ray tracing, if applied correctly, does o?er a valid description of modes with a short vertical wavelength, i.e. those of large vertical mode number (see Figure 11 of Lubow & Pringle 1993). However, such modes are unlikely to be excited by tidal resonant forcing and the f e mode, in particular, is vertically evanescent and has a vertical mode number of zero. As mentioned in Section 1, non-axisymmetric waves of low < azimuthal wavenumber m ? r/H are almost indistinguishable from axisymmetric waves on the scale of a few H. We therefore expect that such waves will also propagate as if through a waveguide and be subject to similar wave channelling in agreement with linear theory. Several papers (e.g. Lin et al. 1990a,b; Terquem 2001) have interpreted the behaviour of low-m (e.g. m = 2) non-axisymmetric waves as being due to refraction. Given our results for m = 0 waves, we conclude that the dominant behaviour of low-m nonaxisymmetric waves is almost certainly determined by wave channelling and not refraction. 8.3 Dissipation of waves

23

The crest of the wave travels faster than the trough and the wave steepens as it propagates until shocks form. This is the principal e?ect leading to the dissipation of the f e mode in a vertically isothermal disc with γ = 5/3 (Figure 2). The situation is not entirely straightforward because the wave Mach number is not uniform across the wavefront. In addition, the steepening may be accelerated by the gradual increase of wave amplitude as the wave propagates outwards into material of lower surface density while conserving its energy ?ux. On the other hand, the steepening may be temporarily postponed by the dispersive character of the wave close to the resonance. The r modes in the vertically isothermal disc are channelled into an increasingly thin layer near the mid-plane and simultaneously develop very short radial wavelengths. These modes are therefore susceptible to viscous damping even in a nearly inviscid disc. In our simulations, these modes develop scales shorter than the grid spacing and are ultimately lost. In the polytropic disc, classical wave steepening does not occur in any simple sense because none of the waves behaves like a plane acoustic wave. The f e mode, in particular, is highly dispersive, which tends to resist non-linear steepening. It undergoes rapid wave channelling to the base of the isothermal atmosphere and behaves like a surface gravity wave. The wave channelling enhances the amplitude of the wave by concentrating its energy into a smaller region of space. In earlier work (Lubow & Ogilvie 1998; Ogilvie & Lubow 1999) we showed that this e?ect is typically su?cient to amplify the wave to sonic velocities where steepening and dissipation are unavoidable. At the same time, the radial and vertical length-scales of the wave also decrease as it propagates outwards (Figure 21), making linear viscous damping a possible source of dissipation. Because wave channelling occurs over a distance that is almost independent of the thickness of the disc, the dissipation of waves with the same initial Mach number near the reasonance occurs at the same radius regardless of the disc’s thickness. We have performed calculations with discs whose thicknesses vary by a factor of 4 to demonstrate this e?ect (Section 6.1.4 and Table 1). The results of the simulations suggest that the f e mode in a polytropic disc dissipates energy at the base of the isothermal layer at a Mach number of approximately 0.2 (Figure 12), which it acquires through the e?ect of wave channelling. The wave energy does not propagate vertically within the isothermal atmosphere, contrary to the suggestion of Papaloizou & Lin (1995). 8.4 Outlook

In situations where waves in discs are excited by tidal forcing from an orbiting companion, the site at which a wave ultimately dissipates is of some importance, because the energy and angular momentum carried by the wave are transferred to the disc there. The simulations provide valuable information on the location and means of wave dissipation, matters about which we were previously able only to speculate (e.g. Lubow & Ogilvie 1998). In a nearly inviscid disc, dissipation can occur only if the wave develops very large velocity gradients. One way to achieve this is the classical non-linear steepening of an acoustic wave in a gas with γ > 1 (e.g. Lighthill 1978).

In future we hope to address some questions that remain unanswered in this paper. In particular, a more detailed study of the competition between linear dispersion and nonlinear steepening of waves in discs would be valuable. Some important wave phenomena will require non-axisymmetric simulations. For example, a Keplerian disc supports slowly varying m = 1 density and bending waves that e?ect an eccentric distortion and a warping of the disc, respectively. Unlike those studied in this paper, these waves can propagate over long distances without experiencing wave channelling. Finally, attention should be given to studying wave propagation in possibly more realistic models such as fully turbu-

24

M. R. Bate et al.
Yuan C., Cassen P., 1994, ApJ, 437, 338

lent discs and layered discs (Gammie 1996). These challenges provide ample opportunities for further developments.

ACKNOWLEDGMENTS We are grateful to Jim Stone for providing the ZEUS-2D code and for many useful discussions. We also gratefully acknowledge support from NASA grants NAG5-4310 and NAG5-10732, the STScI visitor program, and the Institute of Astronomy visitor programme. Some of the calculations discussed in this paper were performed on the GRAND computer. GIO acknowledges the support of the Royal Society through a University Research Fellowship.

APPENDIX A: DISC MODELS A1 A polytropic disc with an isothermal atmosphere

We consider a di?erentially rotating ?uid with density ρ(R, z), pressure p(R, z) and angular velocity ?(R, z), referred to cylindrical polar coordinates (R, φ, z). The equations of equilibrium are ? R?2 0 = = ?Φ 1 ?p ? , ρ ?R ?R 1 ?p ?Φ ? ? , ρ ?z ?z ? (A1) (A2)

REFERENCES
Abramowitz M., Stegun I. A., 1965, Handbook of Mathematical Functions, Dover, New York Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651 Bell K. R., Cassen P. M., Klahr H. H., Henning Th., 1997, ApJ, 486, 372 Cassen P., Woolum D. S. 1996, ApJ, 472, 789 Gammie C. F., 1996, ApJ, 457, 355 Goldreich P., Tremaine S. 1979, ApJ, 233, 857 Goldreich P., Tremaine S. 1980, ApJ, 241, 425 Jackson J. D., 1962, Classical Electrodynamics, Wiley, New York Kato S., Fukue J., 1980, PASJ, 32, 377 Kato S., 2001, PASJ, 53, 1 Korycansky D. G., Pringle J. E., 1995, MNRAS, 272, 618 Landau L. D., Lifshitz E. M., 1960, Electrodynamics of Continuous Media, Pergamon, New York Larson R. G., 1990, MNRAS, 243, 588 Lighthill M. J., 1978, Waves in Fluids, Cambridge Univ. Press, Cambridge Lin D. N. C., Papaloizou J. C. B., 1979, MNRAS, 186, 799 Lin D. N. C., Papaloizou J. C. B., 1993, in Levy E. H., Lunine J. I., eds, Protostars and Planets III, Univ. Arizona Press, Tucson, p. 749 Lin D. N. C., Papaloizou J. C. B., Savonije G. J., 1990a, ApJ, 364, 326 Lin D. N. C., Papaloizou J. C. B., Savonije G. J., 1990b, ApJ, 365, 748 Loska Z., 1986, Acta Astron., 36, 43 Lubow S. H., 1981, ApJ, 245, 274 Lubow S. H., 1991, ApJ, 381, 259 Lubow S. H., Artymowicz P., 1996, in Wijers R. A. M. J., Davies M. B., Tout C. A., eds, Evolutionary Process in Binary Stars, Kluwer, Dordrecht, p. 53 Lubow S. H., Ogilvie G. I., 1998, ApJ, 504, 983 Lubow S. H., Pringle J. E., 1993, ApJ, 409, 360 Nowak M. A., Wagoner, R. V., 1991, ApJ, 378, 656 Ogilvie G. I., 1998, MNRAS, 297, 291 Ogilvie G. I., 2002, MNRAS, in press Ogilvie G. I., Lubow S. H., 1999, ApJ, 515, 767 Okazaki A. T., Kato S., Fukue J., 1987, PASJ, 39, 457 Papaloizou J. C. B., Lin D. N. C., 1995, ARA&A, 33, 505 Papaloizou J. C. B., Nelson R. P., Masset F., 2001, A&A, 366, 263 Shu F. H., Cuzzi J. N., Lissauer J. J., 1983, Icarus, 53, 185 Shu F. H., Yuan C., Lissauer J. J., 1985, ApJ, 291, 356 Stone J. M., Norman M. L., 1992, ApJS, 80, 753 Terquem C. E. J. M. L. J., 2001, in Mathieu R. D., Zinnecker H., eds, The Formation of Binary Stars, ASP Conf. Ser., vol. ???, p.??? von Neumann J., Richtmyer R. D., 1950, J. Appl. Phys., 21, 232

where the gravitational potential is Φ(R, z) = ?GM (R2 + z 2 )?1/2 . (A3)

Let the pressure and density be related by (cf. Lin, Papaloizou & Savonije 1990) p = c2 ρ 1 n+1 ρ ρs
1/n

+1 ,

(A4)

where n is a positive constant, while c(R) and ρs (R) are positive functions. At high densities ρ ? ρs , the gas becomes vertically polytropic, with polytropic index n. At low densities ρ ? ρs , the gas becomes vertically isothermal, with sound speed c. The transitional density ρs de?nes the approximate surface of the disc. De?ne the pseudo-enthalpy h = c2 ρ ρs
1/n

+ ln

ρ ρs

?1 .

(A5)

Then it follows that dp = dh + F dR, ρ where F = d(c2 ) n ? dR n+1 + c2 d ln ρs dR ρ ρs
1/n

(A6)

? ln ρ ρs
1/n

ρ ρs +1 .

+2

1 n+1

(A7)

Consider ?rst the vertical equilibrium at each R separately. Equation (A2) becomes 0=? ?h ?Φ ? , ?z ?z (A8)

with solution h = ?Φ + f (R). By considering the mid-plane z = 0 we determine f = c2 ρc ρs
1/n

(A9)

+ ln

ρc ρs

?1 ?

GM , R

(A10)

Waves in accretion discs
where ρc (R) is the mid-plane density. In general, therefore, c2 ρ ρs
1/n

25

? √

ρc ρs

1/n

+ ln

ρ ρc

= (A11)

R1 ?, and w2 = R2 ?. The value of R1 = 0.2 for all but the highest resolution calculations, for which R1 = 0.6. Results are given for three di?erent values of the disc thickness, ? = Hp /R = 0.05, 0.1, 0.2. A2 A vertically isothermal disc

GM GM ? . R R2 + z 2

Given the functions ρc (R), ρs (R) and c(R), this equation can be solved numerically to determine ρ at any given point, and p follows from equation (A4). Now consider the radial equilibrium. Equation (A1) becomes ?h ?Φ R?2 = +F + . (A12) ?R ?R Thus, from equation (A9), R?2 = = df +F dR d(c2 ) GM + 2 R dR d(c2 ) + 1 ? ln dR + c2 d ln ρs dR 1 n ρc ρs ρ ρc ρ ρs
1/n 1/n 1/n

In the case of a vertically isothermal disc, we take p = c2 ρ, h = c2 ln ρ ρc , (A19) (A20)

with c = c(R), ρ = ρc (R). Then the analysis proceeds as before, but now with F = d(c2 ) ? ln dR ρ ρc + 1 + c2 d ln ρc dR (A21)

?

n n+1

ρ ρs

1/n

and f=? GM . R The density is obtained from c2 ln i.e. ρ ρc = √ GM GM ? , R R2 + z 2 GM GM ? R R2 + z 2 (A22)

1 n+1 ρc ρs

?

1 n

ρc ρs

1/n

(A23)

+c

2 d ln ρc

dR

+1 ,

(A13)

ρ = ρc exp

1 c2



.

(A24)

which determines the angular velocity at every point. Let us assume that the disc is thin and that ρc ? ρs so that there is a well-de?ned polytropic layer with an isothermal atmosphere outside. The scale-height of the isothermal layer is c Hi ≈ , (A14) ? while the thickness of the polytropic layer is Hp ≈ 21/2 ρc ρs
1/2n

The angular velocity is given by R?2 = GM d(c2 ) + ? ln R2 dR ρ ρc + 1 + c2 d ln ρc . dR (A25)

To obtain a scale-height that increases linearly with radius we take c=? GM R
1/2

,

(A26)

Hi .

(A15)

If we take c ∝ R?1/2 and ρc /ρs = constant then Hi ∝ Hp ∝ R so that the scale-height of the disc increases linearly with radius. We de?ne β = ρs /ρc and ? to be the nominal value of Hp /R. Thus c = 2?1/2 β 1/2n ? and ρs = βρc . (A17) GM R
1/2

where, as before, ? is a small constant, so that the isothermal scale-height H = c/? satis?es H/R ≈ ?. As in the polytropic case, the pro?le of ρc (R) is given by equation A18. We use the same values of G, M , σ, δ for the vertically isothermal discs as we use for the polytropic discs, but take R1 = 0.5, R2 = 10.0, w1 = R1 ? and w2 = R2 ?, while n and β are not used.

(A16) APPENDIX B: RESONANT LAUNCHING OF AXISYMMETRIC WAVES IN LINEAR THEORY B1 The f e mode

Finally, we take the pro?le of ρc (R) to be a power law but tapered to give inner and outer edges, ρc = R?σ δ + R ? R1 1 tanh 2 w1 + 1 R2 ? R tanh 2 w2 ,(A18)

where δ ? 1 is the factor by which the density is reduced outside the interval R1 < R < R2 occupied by the disc, and ? ? w1 and w2 are the widths of the inner and outer tapers. The parameters of the model are then n, ?, β, σ, δ, R1 , R2 , w1 , w2 . For our standard polytropic discs, we choose G = 1, M = 1, n = 1.5, β = 0.01, σ = 1.5, δ = 0.01, R2 = 3.0, w1 =

The launching of non-axisymmetric waves at a Lindblad resonance in a three-dimensional disc has been analysed by Lubow & Ogilvie (1998). The total torque exerted at the resonance is identical to that determined by Goldreich & Tremaine (1979) using a two-dimensional model that neglects the vertical structure of the disc. However, in a three-dimensional disc the torque is partitioned into di?erent wave modes. In an inviscid disc, the total torque is equal to the sum of the radial ?uxes of angular momentum of the launched waves, averaged over time, integrated azimuthally

26

M. R. Bate et al.
√ γ = 5/3, as in the present paper, f2D = 5/3 ≈ 74.5%. A further 6.0% goes into the re mode, and yet smaller frac1 tions into higher-order r modes. Interestingly, 18.9% of the expected ?ux is not accounted for by previously documented modes, which evidently do not form a complete set. When we calculate the eigenfunctions of the modes at a Lindblad resonance in a vertically isothermal disc using the vertical boundary conditions imposed in the simulations, the ?ux fractions imparted to the two-dimensional and r modes agree with the numbers given above, indicating that the boundaries are su?ciently distant not to a?ect the modes. At the same time, we discover a further sequence of modes that fully account for the ‘missing’ 18.9% of the expected ?ux. These modes propagate on the same side of the resonance as the f e mode but are strongly a?ected by the boundaries. They become spectrally dense in the limit that the boundaries are removed to a great height, and disappear in the limit γ → 1. Most of the wave energy of these modes is contained away from the mid-plane, near the boundaries. We have identi?ed these additional modes as g modes that are arti?cially con?ned by the boundaries. In a vertically isothermal disc with γ > 1, the buoyancy frequency is proportional to z and therefore increases inde?nitely with height. In the absence of vertical boundaries, internal gravity waves launched at a Lindblad resonance may propagate with a positive but ever diminishing vertical group velocity and never reach a turning point. As a result, they can never satisfy a standing wave condition and form a vertically con?ned g mode. However, they constitute a continuous spectrum of outgoing waves and, as such, receive a non-zero share of the total energy ?ux. High in the atmosphere of a real disc, the adiabatic exponent γ = 1, since the compressional energy of the wave is radiated away in less than a wave period. In that case, we expect that some g modes will be excited and will be vertically con?ned. The individual g modes that are excited should take up the required energy ?ux in a manner that is independent of the location of any vertical boundaries. If we start from our standard polytropic model and continuously reduce the e?ective optical thickness towards zero by increasing the parameter β so that the isothermal atmosphere acquires more mass, the ?ux fraction for the f e mode decreases continuously from 98.3% towards 74.5%. B2 The pe mode 1

and vertically, and reckoned in the direction of propagation at some distance from the resonance. In the case of axisymmetric forcing no torque is exerted. However, the work done at the resonance can still be determined from the formula of Goldreich & Tremaine (1979) in an appropriate limit, using the fact that the energy ?ux of each wave is equal to the angular momentum ?ux multiplied by the angular pattern frequency. For a Keplerian disc the rate at which work is done is equal to F= π 2 R2 A2 Σ , 3ω (B1)

where Σ = ρ dz is the surface density, and the applied radial acceleration is aR = A(R) exp(?iωt). All quantities are evaluated at the location of the resonance. The Lindblad resonance condition ω = κ corresponds to ω = ? in a Keplerian disc. In an inviscid disc, F is equal to the total energy ?ux in the launched waves, i.e. the sum of their radial energy ?uxes, averaged over time, integrated azimuthally and vertically, and reckoned in the direction of propagation at some distance from the resonance. The fraction of the energy ?ux imparted to any wave mode is given by equation (45) of Lubow & Ogilvie (1998),
2

f=

ρu′ R

dz

Σ

ρ|u′ |2 dz, R

(B2)

where u′ (z) is the eigenfunction of the radial velocity perR turbation for the mode considered, and the integrals are over the full vertical extent of the disc. Again, all these quantities are evaluated at the resonance. The imparted fraction therefore depends on an overlap integral between the applied forcing and the modes of the disc. Lubow & Ogilvie (1998) found that, in a polytropic disc, more than 95% of the ?ux is taken up by the f e mode. The remainder goes mainly into the re mode. We have repeated 1 the calculation for the standard polytropic model adopted in the present paper, including the e?ects of the isothermal atmosphere and the vertical boundary conditions imposed in the simulations. We ?nd that 98.3% and 1.6% should go into the f e and re modes, respectively. We recall that the 1 r modes propagate away from the resonance in the opposite direction to the f e mode. The vertical boundaries are su?ciently distant that they have no e?ect on the modes or on the ?ux fractions. In a vertically isothermal disc the eigenfunctions are known analytically from the work of Lubow & Pringle (1993). The radial velocity perturbation for the twodimensional (or f e ) mode is of the form1 u′ ∝ exp R γ?1 γ z2 , 2H 2 (B3)

while, of course, ρ ∝ exp(?z 2 /2H 2 ). We then ?nd for the fraction imparted to the two-dimensional (or f e ) mode f2D = [γ(2 ? γ)]1/2 . (B4)

Only in the case of isothermal perturbations (γ = 1) does the forcing match perfectly to the two-dimensional mode. When
1

There is a sign error in the argument of the exponential in the ?rst and second lines of equation (37) of Lubow & Pringle (1993).

The pe mode is launched at a vertical resonance of the type 1 analysed by Lubow (1981). In the case of axisymmetric forcing the resonance condition is ω = (γ + 1)1/2 ?z , where ?z is the vertical oscillation frequency, equal to ? is a Keplerian disc. Near the resonance, the pe mode involves a predomi1 nantly vertical velocity perturbation u′ ∝ z, corresponding z to a ‘breathing’ motion of the disc, and is excited e?ciently by an applied vertical acceleration az = A(R)z exp(?iωt). We give below an informal derivation of the wave equation satis?ed by the pe mode near the resonance, in order to 1 predict the energy ?ux in the launched wave. A rigorous derivation can be given, if desired, using asymptotic expansions. The exact linearized equations governing axisymmetric Eulerian perturbations of the disc in the presence of vertical forcing of the above type are

Waves in accretion discs
? iωu′ ? 2?u′ = ? R φ ? iωu′ + φ ρ′ ?p 1 ?p′ + 2 , ρ ?R ρ ?R (B5) (B6) (B7)

27

equation (B14) is transformed into the inhomogeneous Airy equation, d2 y = 1, (B17) dx2 as also appears in the analysis of Lindblad resonances. The desired solution, representing a wave emitted from the resonance and travelling into X > 0 with no incoming component, is xy + y = ?π [Gi(?x) + i Ai(?x)] , (B18)

? u′ ? R (R2 ?) + u′ (R?) = 0, z R ?R ?z 1 ?p′ ρ′ ?p ? iωu′ = ? + 2 + Az, z ρ ?z ρ ?z ? iωρ′ + u′ R ? iωp′ + u′ R

1 ? ?u′ ?ρ ?ρ z + u′ = ?ρ (Ru′ ) + , (B8) z R ?R ?z R ?R ?z 1 ? ?u′ ?p ?p z + u′ = ?γp (Ru′ ) + .(B9) z R ?R ?z R ?R ?z

In a thin, Keplerian disc we may take ? = (GM/R3 )1/2 , neglecting its vertical variation, and set ?p/?z = ?ρ?2 z for vertical hydrostatic equilibrium. For a disturbance with radial wavelength short compared to R, we may also neglect the terms involving radial derivatives of p in equation (B5), of ρ and of R in equation (B8), and of p and of R in equation (B9). With these approximations we obtain, after an integration by parts, iωA ρz dz = ω ? (γ + 1)? ?
2 2 2

where Ai and Gi are homogeneous and inhomogeneous Airy functions (e.g. Abramowitz & Stegun 1965). From the asymptotic form for large positive x, y ? ?π 1/2 x?1/4 exp i 2 3/2 π x + 3 4 , (B19)

we determine the radial energy ?ux (averaged in time, and integrated azimuthally and vertically) in the wave at some distance from the resonance, F = = πR u′ p′ dz, R ρz 2 dz. (B20)
?

ρu′ z z ?u′ R ?R

dz dz. (B10)

?p γp + z ?z

π 2 R2 A2 3ω

The quantity in square brackets vanishes at the resonance. If X = R ? Rres is the radial distance from the resonance, we may apply the linear approximation ?2 (B11) ω 2 ? (γ + 1)?2 ≈ 3(γ + 1) X. R e The leading approximation to the p1 mode near the resonance is u′ z iωzY (X), ?ρ ρ′ = ρ+z Y (X), ?z ?p p′ = γp + z Y (X), (B12) ?z where Y (X) is a function to be determined. The horizontal components of the equation of motion then determine u′ = R iω ?2 ? ω 2 1 ?p γp + z ρ ?z dY . dX (B13) =

This is directly analogous to equation (B1) for a Lindblad resonance. B3 The f o /ro mode 1

The integral relation (B10) then yields the ordinary di?erential equation d2 Y = C3 , dX 2 where the coe?cients C1 XY + C2 C1 C2 C3 = = = ?2 3(γ + 1) R 1 γ?2 A ρz 2 dz,
2

(B14)

Generically, the f o mode, which e?ects a corrugation or tilt of the disc, is launched at a vertical resonance. In the case of axisymmetric forcing the resonance condition is ω = ?z . In a Keplerian disc this coincides with the Lindblad resonance and the f o /ro mode is launched at a hybrid vertical/Lindblad 1 resonance that has not been discussed in the literature. Near the resonance, the f o /ro mode involves a vertical 1 velocity perturbation u′ independent of z and horizontal z velocity perturbations u′ , u′ ∝ z of comparable magnitude. R φ It is excited e?ciently by an applied vertical acceleration az = A(R) exp(?iωt). We give below an informal derivation of the wave equation satis?ed by the f o /ro mode near the 1 resonance. Again, a rigorous derivation can be given using asymptotic expansions. The equations governing the perturbations are the same as in the previous subsection except that the vertical acceleration no longer contains the factor z. We now obtain the integral relation iωA ρ dz = (ω 2 ? ?2 ) ρu′ dz ? z ?p ?u′ R dz, ?z ?R (B21)

1 ?p γp + z ρ ?z ρz dz
2

dz, (B15)

in which the quantity in brackets vanishes at the resonance. If X = R ? Rres is the radial distance from the resonance, we may apply the linear approximation 3?2 X. (B22) R The leading approximation to the f o /ro mode near the res1 onance is ω 2 ? ?2 ≈ u′ z (B16) ρ′ = = iωY (X), ?ρ Y (X), ?z

are all to be evaluated at the location of the resonance. By means of the change of variables Y (X) X = = C1
?2/3 ?1/3

C2

?1/3 1/3

C3 y(x),

C1

C2

x,

28
p′ =

M. R. Bate et al.

?p Y (X). (B23) ?z The horizontal components of the equation of motion then determine 1 ?p dY iω . (B24) u′ = R ?2 ? ω 2 ρ ?z dX The denominator of this expression may be represented by the same linear approximation, equation (B22). The integral relation (B21) then yields the ordinary di?erential equation C1 XY + C2 d dX 1 dY X dX = C3 , (B25)

where the coe?cients C1 C2 C3 = = = 3?2 R R 3?2 A ρ dz, 1 ρ ρ dz ?p ?z
2

dz, (B26)

are all to be evaluated at the location of the resonance. By means of the change of variables Y (X) X = = C1
?3/4 ?1/4

C2

?1/4 1/4

C3 y(x), (B27)

C1

C2

x,

equation (B25) is transformed into xy + d dx 1 dy x dx = 1. (B28)

Unlike equation (B17), this equation allows wave propagation on both sides of the resonance, x > 0 and x < 0. The complementary functions of equation (B28) are cos( 1 x2 ) 2 1 and sin( 2 x2 ), and the general solution of the inhomogeneous equation can therefore be expressed in terms of Fresnel integrals. The desired solution, representing one wave travelling into X > 0 and another wave travelling into X < 0, with no incoming component, is y = 1 ? iπ 1/2 ? 2
x 1 sin( 1 t2 ) dt cos( 2 x2 ) 2 0 x 1 cos( 2 t2 ) dt sin( 1 x2 ). 2 0

1 + ? iπ 1/2 + 2

(B29)

From the asymptotic form for large positive x, y?? π 2
1/2

exp i

1 2 π x + 2 4

,

(B30)

we determine the energy ?ux in the wave at some distance from the resonance on the positive side, F = = πR u′ p′ dz, R (B31)
?

π 2 R2 A2 Σ . 6ω

This is exactly half the value given by equation (B1) for a Lindblad resonance. From the symmetry property y(?x) = ?y ? (x), it follows that the wave emitted into x < 0 carries an equal and oppositely directed energy ?ux. Therefore the total energy ?ux is identical to that for a Lindblad resonance, although of course A is the amplitude of the vertical, not horizontal, acceleration.

This figure "fig03.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig07.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig09.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig11.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig17.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig19.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig20.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1

This figure "fig21.gif" is available in "gif" format from: http://arXiv.org/ps/astro-ph/0201149v1



更多相关文章:
更多相关标签:

All rights reserved Powered by 甜梦文库 9512.net

copyright ©right 2010-2021。
甜梦文库内容来自网络,如有侵犯请联系客服。zhit325@126.com|网站地图