Instability and warping in vertically oscillating accretion discs

ABSTRACT

Many accretion discs have been found to be distorted: either warped due to a misalignment in the system, or non-circular as a result of orbital eccentricity or tidal deformation by a binary companion. Warped, eccentric, and tidally distorted discs are not in vertical hydrostatic equilibrium, and thus exhibit vertical oscillations in the direction perpendicular to the disc, a phenomenon that is absent in circular and flat discs. In extreme cases, this vertical motion is manifested as a vertical ‘bouncing’ of the gas, potentially leading to shocks and heating, as observed in recent global numerical simulations. In this paper, we isolate the mechanics of vertical disc oscillations by means of quasi-2D and fully 3D hydrodynamic local (shearing-box) models. To determine the numerical and physical dissipation mechanisms at work during an oscillation we start by investigating unforced oscillations, examining the effect of initial oscillation amplitude, as well as resolution, boundary conditions, and vertical box size on the dissipation and energetics of the oscillations. We then drive the oscillations by introducing a time-dependent gravitational potential. A key result is that even a purely vertically oscillating disc is (parametrically) unstable to developing inertial waves, as we confirm through a linear stability analysis. The most important of these has the character of a bending wave, whose radial wavelength depends on the frequency of the vertical oscillation. The nonlinear phase of the instability exhibits shocks, which dampen the oscillations, although energy can also flow from the bending wave back to the vertical oscillation.

1 INTRODUCTION

The classical picture of accretion discs assumes they are circular, flat (not warped), and not tilted or misaligned with respect to the rotation axis of their central object (or binary orbital plane in the case of discs in binaries). In fact, many accretion discs in nature have been observed or deduced to be distorted from this simple form, in that they are warped, eccentric or tidally deformed. ALMA and VLTI observations of discs around young stars have inferred warps from the shadows cast by a misaligned inner disc on the outer part of the disc (Kluska et al. 2020; Bohn et al. 2022). Other examples include the misaligned (with the binary orbit) discs in the young binary HK Tau (Jensen & Akeson 2014; Manara et al. 2019), the triple star system GW Ori (Bi et al. 2020; Kraus et al. 2020), the tilted circumbinary disc of KH 15D (Chiang & Murray-Clay 2004; Poon, Zanazzi & Zhu 2021), and the warped disc around the young star IRAS 04368 + 2557 (Sakai et al. 2019; Villenave et al. 2024); the X-ray binaries Her X-1, SS 433, SMC MAXI J1820 + 070 (among others) are believed to have precessing tilted discs (Katz 1973; Gerend & Boynton 1976; Clarkson et al. 2003; Begelman, King & Pringle 2006; Kotze & Charles 2012; Thomas et al. 2022), and warped discs have also been observed in the active galactic nuclei NGC 4258, NGC 1068, and the Circinus galaxy (Miyoshi et al. 1995; Herrnstein, Greenhill & Moran 1996; Greenhill et al. 2003). Precessing eccentric discs have been deduced, for example, from observations of superhumps in close binary stars (Patterson et al. 2005) and from variable asymmetric emission lines from discs around Be stars (Okazaki 2016) and white dwarfs (Manser et al. 2021); eccentric discs have also been inferred from emission lines from some tidal disruption events (e.g. Wevers et al. 2022). Non-axisymmetric structure in tidally distorted discs in close binary stars has been measured using Doppler tomography (Marsh 2005).

Distorted discs exhibit a variety of dynamics not found in their idealized circular and planar counterparts. In particular, they are not in hydrostatic equilibrium in the direction perpendicular to the (local) orbital plane, resulting in periodic vertical contraction and expansion of the disc. This type of vertical oscillation, which we will refer to as bouncing, can result in heating of the gas through adiabatic or dissipative effects, leading to non-axisymmetric emission (Ogilvie 2002a; Liska et al. 2021; Kaaz et al. 2023).

1.1 Phenomena in distorted discs

In addition to the aforementioned vertical oscillations (bouncing), distorted accretion discs can also exhibit behaviour such as nodal and/or apsidal precession of the disc (Bardeen & Petterson 1975; Lubow 1991; Ogilvie 2001; Deng & Ogilvie 2022; Tremaine 2023), disc breaking/tearing when the disc is strongly tilted or warped (Lodato & Price 2010; Nixon & King 2012; Tremaine & Davis 2014; Drewes & Nixon 2021; Liska et al. 2021; Raj, Nixon & Doğan 2021; Nealon et al. 2022; Kaaz et al. 2023; Young et al. 2023), and hydrodynamic instability (known as parametric instability). In this paper, we will isolate the dynamics of vertical oscillations, though a new and interesting result is that these oscillations are parametrically unstable. We review both in turn below.

1.1.1 Vertical oscillations or bouncing

Radial resonances (usually known as Lindblad resonances) – and, more generally, radial oscillations and waves – in tidally forced discs have been studied in great detail (Goldreich & Tremaine 1979; Meyer-Vernet & Sicardy 1987; Artymowicz & Lubow 1994; Lubow, Martin & Nixon 2015; Nixon & Lubow 2015; Ju, Stone & Zhu 2016; Armitage 2022). However, vertical oscillations and resonances, involving symmetric motion¹ perpendicular to the plane of the disc, are also of interest, although they have been more sparsely studied in the literature (Lubow 1981; Stehle 1999; Stehle & Spruit 1999; Ogilvie 2002b; Fairbairn & Ogilvie 2021). We will refer to the latter interchangeably as vertical oscillations, vertical breathing modes, or bouncing, though some care is needed with the terminology, which we detail below.

The specific vertical oscillations we are interested in here involves uniform expansion and contraction of the disc in the |$z$|-direction perpendicular to the local orbital plane, i.e. the displacement and velocity of the gas are proportional to |$z$|⁠. See Fig. 1 for an illustrative example in a tidally distorted circumstellar disc in a circular binary. These oscillations can range from mild vertical breathing motion (Ogilvie 2002b, 2008; Ryu, Krolik & Piran 2021) to a more extreme bouncing motion in which pressure support is unimportant most of the time and the disc expands and contracts in a ballistic manner in the vertical direction, with pressure acting impulsively to prevent the disc from collapsing (Fairbairn & Ogilvie 2021; Lynch & Ogilvie 2021).² It is worth emphasizing that these vertical oscillations of the gas occur around an orbit, in the sense that the disc (as viewed from a global reference frame) will appear ‘pinched’ at certain locations. This can occur even in very thin discs and does not, in general, involve coupling between vertical and radial motion. It is therefore distinct from the homogeneous contraction and expansion of the entire disc in both the radial and vertical directions. The latter oscillations can arise in slim or thick discs (tori), and are sometimes referred to as |$+$|- (plus) and breathing-modes (Blaes, Arras & Fragile 2006).

Bouncing can be excited in all types of distorted disc. It has been studied in the context of warped discs (Ogilvie & Latter 2013; Fairbairn & Ogilvie 2021), eccentric discs (Ogilvie & Barker 2014; Ryu et al. 2021; Chan, Piran & Krolik 2023), and also in tidally distorted discs (Stehle & Spruit 1999; Ogilvie 2002a), though mechanisms driving the bouncing differ in each case. In an eccentric disc, the vertical component of gravity due to the central star varies along each eccentric orbit (Ogilvie 2001; Barker & Ogilvie 2014), thus the vertical oscillation frequency is equal to the local orbital frequency. In a tidally distorted disc, the frequency at which the non-axisymmetric potential of the binary drives vertical oscillations range from about once to twice per orbit (in the frame of the fluid) depending on radial location in the disc (see Fig. 1), and vertical resonances can occur at specific locations (Lubow 1981). In a warped disc, the mechanism driving vertical oscillations is more subtle: an observer moving along a fixed orbit sees neighbouring orbits tilt downwards and upwards along the trajectory. The combination of vertical motion (due to the tilting) and radial motion (due to radial pressure imbalances between neighbouring vertical columns of fluid) induces vertical compression and rarefaction of the disc along an orbit at a frequency of twice per orbit, which can be enhanced by a non-linear resonance (Ogilvie & Latter 2013; Fairbairn & Ogilvie 2021).

1.1.2 Parametric instability

Distorted discs are unstable to a hydrodynamic instability involving the resonant amplification of inertial waves known as parametric instability.³ It was first discovered (theoretically) in the context of tidally distorted discs (Goodman 1993; Lubow, Pringle & Kerswell 1993; Ryu & Goodman 1994), and later also in warped discs (Papaloizou & Terquem 1995; Gammie, Goodman & Ogilvie 2000; Ogilvie & Latter 2013) and eccentric discs (Papaloizou 2005a, b; Barker & Ogilvie 2014; Ogilvie & Barker 2014; Wienkers & Ogilvie 2018).

1.2 Modelling distorted discs using local models

Distorted discs have been modelled in global simulations using Lagrangian (smoothed-particle hydrodynamics [SPH]) (Lodato & Price 2010; Deng & Ogilvie 2022; Nealon et al. 2022), Eulerian (grid-based) (Kimmig & Dullemond 2024), and meshless approaches (Deng, Ogilvie & Mayer 2021; Deng & Ogilvie 2022). However, global simulations struggle to capture some of the key physics of distorted discs, in particular the parametric instability, although this has been achieved with grid-based (Pierens, McNally & Nelson 2020) and meshless (Deng et al. 2021) codes. For example, SPH simulations employ a smoothing length (typically some fraction of the disc thickness), which dampens the instability and also artificially damps the bouncing motion (confusing it with a shock-forming flow), while grid-based simulations introduce an advection error when the flow is not aligned with the grid.

These difficulties have motivated local (shearing-box) studies of distorted discs, which have led to warped (Ogilvie & Latter 2013; Paardekooper & Ogilvie 2019), as well as eccentric (Ogilvie & Barker 2014; Wienkers & Ogilvie 2018) extensions to the standard shearing-box formalism. Perhaps surprisingly, much of the physics of distorted discs can actually also be captured in a standard shearing box (Ogilvie 2022), which is the approach we take in this work. There are limitations to this approach, however. First, the standard shearing box follows circular streamlines in a disc. Thus, we do not take into account any geometrical effects due to eccentric or tidally distorted orbits, in order to focus solely on the dynamics of the vertical oscillations. Secondly, the standard shearing box cannot pick up a large-scale warp, only a radially periodic warp (also known as a corrugation or bending wave/mode), although a non-periodic ring or torus with a large-scale warp can be studied in the local model (Fairbairn & Ogilvie 2021).

1.3 Motivation and aims of paper

Our main aim is to isolate the dynamics of vertical (i.e. perpendicular to the local plane of the disc) oscillations or bouncing expected to occur in tidally distorted, eccentric, and warped accretion discs. We do so by means of quasi-2D and fully 3D hydrodynamic shearing box simulations. For simplicity, we employ an isothermal equation of state. Note that in astrophysical contexts we expect the oscillations to be forced (driven). However, to determine the main mechanisms damping the oscillations, we first consider free (unforced) oscillations.

The structure of the paper is as follows. In Section 2.1, we discuss our set-up, initial conditions, parameters, and key diagnostics. In Section 3, we present a 1D theory (including a stability analysis) of vertical disc oscillations. We present simulations of free oscillations in Section 4. In Section 5, we discuss forced oscillations (including astrophysically relevant forcing frequencies and amplitudes). We conclude in Section 6.

2 METHODS

2.1 Governing equations

We use the shearing-box approximation (Goldreich & Lynden-Bell 1965; Hawley, Gammie & Balbus 1995; Latter & Papaloizou 2017). Rather than modelling the entire disc, this formalism considers only a small patch (a Cartesian box) in the disc. The box is centred around some fiducial radius |$r = r_0$| and orbits with the angular frequency of the disc at that radius, |$\Omega _0 \equiv \Omega (r_0)$|⁠. A point in the box has Cartesian coordinates |$(x, y, z)$| along the radial, azimuthal, and vertical directions, respectively. In this rotating frame, the equations of gas dynamics are as follows:

with the symbols taking their usual meanings. We close the system with the equation of state for an isothermal gas, |$P = c_\text{s}^2 \rho$|⁠, where |$c_\text{s}$| is the constant sound speed.

All our simulations are stratified in the vertical direction. In our simulations of a freely vertically oscillating disc, the effective gravitational potential is embodied in the usual (constant-in-time) tidal acceleration |$\boldsymbol {g_\text{eff}}=2q\Omega _0^2x\, \boldsymbol {e}_x-\Omega _0^2 z\, \boldsymbol {e}_z$| [third term on the right-hand side of equation (2)], where |$q$| is the dimensionless shear parameter |$q \equiv -\left.d\ln {\Omega }/d\ln {r}\right|_{r=r_0}$|⁠. For Keplerian discs |$q=3/2$|⁠, a value we adopt throughout this paper.

In order to drive vertical oscillations/bouncing (Section 5), we introduce a time-dependent gravitational acceleration in the vertical direction, by modulating the usual gravity of the shearing box according to |$g_{\text{eff},z} = -\Omega _0^2 z (1+a\cos {\omega t})$|⁠, where |$a$| is the dimensionless forcing amplitude and |$\omega$| the forcing frequency. With these two parameters we can mimic the driving of the vertical component of gravity found, e.g. in a binary system in which a companion external to the disc modulates the vertical gravity experienced by gas in the disc. We determine realistic values of these parameters to use in our forced bounce simulations in Section 5.1.

Finally, the viscous stress tensor is given by |$\boldsymbol {T} \equiv 2\rho \nu \boldsymbol {S}$|⁠, where |$\nu$| is the kinematic viscosity, and |$\boldsymbol {S} \equiv (1/2)[\nabla \boldsymbol {u} + (\nabla \boldsymbol {u})^\text{T}] - (1/3)(\nabla \cdot \boldsymbol {u})\boldsymbol {I}$| is the traceless shear tensor (Landau & Lifshitz 1987). We keep the viscosity fixed in space and time in any given simulation. We have also run select simulations without an explicit viscosity.

2.2 Numerical set-up

2.2.1 Code

Our simulations are run using the conservative, finite-volume code pluto (Mignone et al. 2007). In terms of numerical algorithms, we use the HLLC Riemann solver, second-order-in-space linear interpolation, and the second-order-in-time Runge–Kutta algorithm. To allow for longer time-steps, we take advantage of the fargo scheme (Mignone et al. 2012). When we include explicit viscosity |$\nu$|⁠, we also use the Super-Time-Stepping (STS) scheme (Alexiades, Amiez & Gremaud 1996), which further reduces computational time. Ghost zones are used to implement the boundary conditions (BCs).

We use the built-in shearing box module in pluto (Mignone et al. 2012). Rather than solving equations (1)–(2) (primitive form), pluto solves the governing equations in conservative form.

2.2.2 Units

Note that from this point onwards, all quantities are given in terms of dimensionless (code) units. Time units are selected so that |$\Omega _0 = 1$|⁠. The length unit is chosen so that the initial sound speed |$c_{\text{s}0} = 1$|⁠, which in turn defines a reference equilibrium scale height |$H_0\equiv c_{\text{s}0}/\Omega _0=1$|⁠, which can be thought of as the disc thickness in vertical hydrostatic equilibrium. Finally the mass unit is set by the initial mid-plane density |$\rho _0 = 1$|⁠. Pressure, stresses, and energy densities are expressed in units of |$c_{\text{s}0}^2 \rho _0$|⁠. From here onwards, we will drop the subscripts on |$\Omega _0$| and |$c_{\text{s}0}$|⁠, since these quantities are constant. However, we will retain the subscript on |$H_0$| (the equilibrium scale height) to distinguish it from the time-dependent dynamical scale height |$H(t)$| (formally defined below).

2.2.3 Initial conditions

We employ different initial conditions depending on whether the vertical oscillations are free or forced. In the former case, we start from a disc that has been uniformly stretched in the vertical direction, with a Gaussian density profile of the form

where |$H(0)$| is the dynamical scale height (defined below) at initialization (⁠|$t=0$|⁠). In other words, |$H(0)\lt 1$| (⁠|$H(0)\gt 1$|⁠) corresponds to a disc that is compressed (stretched) in the vertical direction compared to its equilibrium thickness |$H_0$|⁠. Note that this initial condition is not in equilibrium by design, in order to induce the disc to oscillate. In our freely bouncing disc simulations, we select |$H(0)$| to study bounces of small and large amplitude, respectively. In our forced bounce simulations, we set |$H(0)$| based on the theoretical predicted value for the disc thickness at a given forcing frequency (see Section 5.2).

The background velocity is given by |$\boldsymbol {u}_0 = -q \Omega x \, \boldsymbol {e}_y$|⁠. At initialization, we usually perturb all the velocity components with random noise exhibiting a flat power spectrum (different random noise is used for each component). The perturbations |$\delta \boldsymbol {u}$| have maximum amplitude of about |$5\times 10^{-2}\, c_\text{s}$|⁠.

2.2.4 Box size and resolution

The majority of our simulations are fully three-dimensional (3D) and run at a resolution of 32 cells per scale height |$H_0$| (denoted |$32/H_0$|⁠), in a box of size |$L_y = 4H_0$| in the azimuthal direction and |$L_z = 12H_0$| in the vertical direction. (We also repeated our fiducial 3D run of free [unforced] vertical oscillations at |$64/H_0$|⁠.) As the radial wavelength of the fastest growing mode of the parametric instability is sensitive to the forcing frequency |$\omega$|⁠, we vary the radial box size at different forcing frequencies to ensure we always capture the fastest growing mode, from |$L_x = 4H_0$| (at |$\omega /\Omega = 1$|⁠) to |$L_x = 30H_0$| (at |$\omega /\Omega = 2$|⁠). To investigate the dependence on the vertical box size, we have also repeated select unforced runs in boxes of size |$[4,4,8]H_0$| and |$[4,4,16]H_0$| (keeping the resolution per scale height fixed at 32 cells per |$H_0$|⁠). Finally, we have also run a series of quasi-2D simulations in boxes of size |$[8,0.008, 12]H_0$|⁠. These enable us to reach resolutions up to |$512/H_0$|⁠. All the simulations described in this paper are listed in Tables B1–B5 in Appendix  B.

2.2.5 Boundary conditions

We use standard shear-periodic BCs in the |$x$|-direction (see Hawley et al. 1995), and periodic BCs in the |$y$|-direction. In the vertical direction, we mostly employ reflective BCs (⁠|$\rho \rightarrow \rho$|⁠, |$u_{x,y} \rightarrow u_{x,y}$|⁠, and |$u_z \rightarrow -u_z$|⁠). However, we have also checked our fiducial 3D simulation with outflow BCs, for which we set |$\partial _z u_i = 0$| (⁠|$i \in {x,y,z}$|⁠), described in greater detail in Held, Mamatsashvili & Pessah (2024). Note that when we use outflow BCs mass (and energy) can be advected across the vertical boundaries, and are not added back into the domain using a mass source term as was done in Held et al. (2024).

2.3 Diagnostics

2.3.1 Averaged quantities

The volume-average of a quantity |$X$| is denoted |$\langle X \rangle$| and is defined as

where |$V=L_x L_y L_z$| is the volume of the box.

We denote the horizontal average of a quantity |$X$| as |$\langle {X}\rangle _{xy}$|⁠, defined as

where |$A = L_x L_y$|⁠.

2.3.2 Reynolds stress

In purely hydrodynamic accretion discs, the radial transport of angular momentum is related to the |$xy$|-component of the Reynolds stress |$R_{xy} \equiv \rho u_x \delta u_y$|⁠, where |$\delta u_y \equiv u_y + q\Omega x$| is the perturbation of the y-component of the total velocity |$u_y$| about the background Keplerian flow |$u_{0y} = -q \Omega x$|⁠.⁴ (From now on we will drop the |$\delta$| and simply refer to the perturbed velocity as |${\bf u}$|⁠.)

2.3.3 Dynamical scale height

A key diagnostic in our simulations is the dynamical scale height|$H(t)$| which measures the time-variability of the disc thickness as it oscillates in the vertical direction. We define the square of this quantity as follows:

where |$\Sigma \equiv \int _{-L_z/2}^{L_z/2} \langle \rho \rangle _{xy} dz$| is the surface density of the disc. Thus, |$H$| is related to the second moment of the density and can be interpreted as the standard deviation of |$z$| (assuming zero mean) of the mass elements of the disc.

Finally, we define the amplitude of the bounce as

where |$H_{\text{max},n}$| and |$H_{\text{min},n}$| correspond to the maximum and minimum scale heights, respectively, in the |$n$|th oscillation cycle (we define a cycle as starting when the disc is at maximum thickness, so |$H_{\text{min},n}$| occurs half a cycle later).

2.3.4 Energetics

Another useful diagnostic is the total energy density, which is given by

where the first two terms on the right-hand side correspond to the kinetic |$E_{\text{kin}}$| and gravitational potential energy density, respectively. The term |$c_\text{s}^2 \rho \ln {\rho }$| embodies the variable part of the Helmholtz-free energy, which replaces the internal energy in an isothermal gas.⁵ In our freely bouncing disc simulations, |$\Phi =-\frac{3}{2}\Omega ^2x^2+\frac{1}{2}\Omega ^2z^2$|⁠, which is the standard (time-independent) effective gravitational potential in the shearing-box approximation for a Keplerian disc. In our forced bouncing disc simulations, the vertical component of the gravitational potential is time-dependent (to model the time-variation in vertical gravity along an orbit in a distorted disc). In this case, a term of the form |$\rho \partial \Phi /\partial t$| appears as a source term on the right-hand side of the total energy equation.

3 THEORY

In this section, we derive some basic theoretical predictions related to disc bouncing in the absence of damping or forcing. In Section 3.1, we write down the governing equations for 1D vertical bouncing, and, by considering first hydrostatic solutions, use these to help derive a governing equation for the dynamical scale height |$H(t)$| in the non-hydrostatic case. In Section 3.2, we examine the stability of the basic bouncing flow to radially dependent perturbations, and discover that the bouncing flow is parametrically unstable. We defer a discussion of forced bouncing to Section 5.2.

3.1 1D theory of a vertically bouncing disc

3.1.1 Governing equations

Using the shearing-box model of Section 2.1, we consider a vertical motion of the disc, independent of |$x$| and |$y$|⁠, in which the vertical velocity |$u_z(z,t)$| and the density |$\rho (z,t)$| satisfy the continuity equation and the |$z$|-component of the momentum equation,

together with |$p=c_\text{s}^2\rho$|⁠, where the material derivative is given by |$D/Dt \equiv \partial _t + u_z \partial _z$|⁠. As for the horizontal components of the momentum equation, the orbital shear flow is unaffected by the vertical motion.

3.1.2 Hydrostatic equilibrium solution

In hydrostatic equilibrium, we take |$u_z = 0$|⁠, |$\rho = \rho (z)$|⁠, and |$p=p(z)$| and the governing equations reduce to the structural equation

Together with the isothermal condition |$p=c_\text{s}^2\rho$|⁠, this equation determines the usual Gaussian structure of a hydrostatic disc, which we write in the form

where

is a normalized Gaussian function and |$\zeta =z/H$| is a dimensionless vertical coordinate. For hydrostatic equilibrium, the scale height is given as usual by |$H=c_\text{s}/\Omega$|⁠.

3.1.3 Non-hydrostatic solutions

To study the bouncing motion, we now aim to derive an equation governing the time-evolution of the disc thickness, or scale height. We look for non-hydrostatic solutions in which the scale height |$H(t)$| used in the non-dimensionalization above depends on time. The solutions we are interested in have the same dimensionless density (and pressure) profile as the hydrostatic solution; they are achieved by a uniform, time-dependent expansion or contraction of the model. The dimensionless variable |$\zeta$| becomes a Lagrangian coordinate, which means that the vertical velocity is given by |$u_z=(\mathrm{d}H/\mathrm{d}t)\zeta =(\mathrm{d}H/\mathrm{d}t)(z/H)$|⁠. Mass conservation implies that the surface density |$\Sigma$| is independent of time. Substituting the density (12), with a time-dependent scale height, into the momentum equation (10), we find that all terms are proportional to |$\zeta$|⁠, and thereby obtain an evolution equation for the disc thickness or dynamical scale height |$H(t)$|⁠:

The hydrostatic disc corresponds to the equilibrium solution |$H=H_0=c_\text{s}/\Omega$|⁠.

In the units defined in Section 2.2.2, we have

where a dot denotes a time-derivative.

We note the derivation above can be extended to allow for adiabatic thermodynamics and different vertical structure profiles, in which case the dimensionless scale height satisfies |$\ddot{H} = -H + H^{-\gamma }$|⁠, where |$\gamma$| is the adiabatic exponent (cf. Ogilvie & Barker 2014).

3.1.4 Periodic solutions

We now derive predictions for the oscillation period in the small- and large-oscillation amplitude limits. Let us first consider the case of small-amplitude oscillations |$H=1+\epsilon (t)$|⁠. Substituting into equation (15) and linearizing, we obtain |$\ddot{\epsilon }=-2\epsilon$|⁠, which describes small-amplitude oscillations of angular frequency |$\sqrt{2}$| (or |$\sqrt{\gamma +1}$| when adiabatic thermodynamics is taken into account).

To consider large-amplitude oscillations, we first rewrite the right-hand side of equation (15) in terms of the gradient of a potential |$V$|⁠,

where |$V(H) = (1/2)H^2-\ln H$|⁠. Integrating, we obtain

where the sum of kinetic and potential energies |$E$| is a constant. Rearranging this equation and integrating between the minimum and maximum disc thickness over half an oscillation, we obtain an expression for the period of the oscillation for arbitrary oscillation amplitude:

In the limit of large amplitude, we find |$2\pi /T \rightarrow 2$| (Fairbairn & Ogilvie 2021).

The angular frequency increases from |$\sqrt{2}$| to 2 as the amplitude |$H_\text{max}-H_\text{min}$| increases from 0 to |$\infty$| (see Fig. 2).

3.2 Parametric instability of solutions

Motivated by the simulation results to be described below, we analyse the stability of a periodically bouncing disc to axisymmetric perturbations that are independent of |$y$| but depend on |$x$| through a common factor |$\mathrm{e}^{\mathrm{i}kx}$|⁠, where |$k$| is the radial wavenumber. The linearized equations for small perturbations can be formulated in terms of the Lagrangian displacement |$\boldsymbol{\xi }$|⁠. We find

Transforming from Eulerian |$(z,t)$| to Lagrangian |$(\zeta ,t)$| variables [recalling (Section 3.1.3) that |$D/Dt \equiv \partial _t + u_z \partial _z$| and |$u_z = (dH/dt)(z/H) \equiv \dot{H} \zeta$|⁠, and using the Gaussian dimensionless density profile (equation (13) to evaluate |$(1/\rho)\partial _z \rho$| in the last equation], we find

By analogy with the linear modes of a disc in hydrostatic equilibrium (Okazaki, Kato & Fukue 1987), we seek a solution in Hermite polynomials:

where the vertical mode number |$n$| is a non-negative integer. The linearized equations are satisfied provided that

For an equilibrium disc (⁠|$H=H_0$|⁠), these equations have constant coefficients and therefore have harmonic solutions with angular frequency |$\omega$| given by the dispersion relation

in agreement with Okazaki et al. (1987). (Note that here, and here only, |$\omega$| refers to the mode frequency, whereas below and elsewhere in the paper it refers to the forcing frequency at which the vertical oscillations are driven.) For a bouncing disc, we have instead a fourth-order system of ordinary differential equations with time-periodic coefficients [because |$H(t)$| is periodic] and there is the possibility of instability. Floquet’s method can be used to analyse the solutions and their growth rates. We solve the dimensionless version of equations (29)–(31), in which the factors of |$\Omega$| and |$c_\text{s}$| drop out and |$H(t)$| satisfies |$\ddot{H}=-H+H^{-1}$| (see equation 15). Starting at |$t=0$| with an initial value |$H(0)\gt 1$| and with |$\dot{H}(0)=0$|⁠, meaning that the disc is released from rest above its equilibrium scale height. By integrating four different solutions of the systems (29)–(31) with four linearly independent initial conditions such as |$(X,Z,\dot{X},\dot{Z})=(1,0,0,0)$|⁠, (0,1,0,0), (0,0,1,0), and (0,0,0,1) and integrating for one oscillation period |$T$| of |$H(t)$|⁠, we construct the monodromy matrix that maps the initial values to the final values one period later (Teschl 2024). We calculate its eigenvalues, which are the Floquet multipliers |$\lambda$|⁠, and convert these into exponential growth rates |$s$| via |$\exp (sT)=|\lambda |$|⁠.

The linearized equations in the form written above are also valid when the bouncing motion is forced, provided that |$H(t)$| now satisfies the appropriate forced equation such as |$\ddot{H}=-(1+a\cos \omega t)H+H^{-1}$|⁠, where |$\omega$| now refers to the forcing frequency (not the mode frequency).

Typical results for the instability of freely bouncing discs are shown in Fig. 3. The nature of the instability is very similar to that analysed for eccentric discs analysed by Barker & Ogilvie (2014). An inertial wave with a natural frequency close to half the bounce frequency is destabilized by the oscillating geometry of the disc. The instability is an example of three-wave coupling in which a parent mode (the bouncing motion) couples to two daughter modes (inertial waves travelling radially inwards and outwards). For the coupling to be exactly resonant, the frequencies and wavenumbers of the three modes must sum to zero for some choice of sign combinations. Since the parent mode has radial wavenumber |$k=0$|⁠, the two daughter modes have the same value of |$k$|⁠. Selection rules regarding the vertical mode number |$n$| imply that the daughter modes are identical except for the direction of propagation; the combined growing mode has the character of a standing wave. For each |$n$| there is a single |$k$| at which the inertial-wave frequency is half the bounce frequency and the resonance condition is satisfied. The growth rate at exact resonance increases with the bounce amplitude. Each resonant band has a width that also increases with the bounce amplitude, because parametric resonance can tolerate some detuning.

Analytical expressions can be obtained for the radial wavenumber and growth rate of the parametric instability in the limit of small bounce amplitude. We omit the details of this calculation. If |$\omega$| is the angular frequency of the bouncing motion (whether free or forced), then for vertical mode |$n$| the resonance is centred at

and has a maximum growth rate of

which is proportional to the semi-amplitude |$H_\text{max}-1$|⁠. The fastest-growing mode is |$n=1$|⁠, with growth rate

In the case of free bouncing (⁠|$\omega =\sqrt{2}$|⁠), this simplifies to |$(\sqrt{2}/6)(H_\text{max}-1)\approx 0.236(H_\text{max}-1)$| and the resonant wavenumber is |$k=1/\sqrt{2}$|⁠, corresponding to a wavelength of nearly nine scale heights.

It is interesting to compare these results with the case of eccentric discs. Papaloizou (2005a) first analysed the parametric instability of inertial waves in eccentric discs, using a model that neglects vertical gravity. In this case, the instability is due to the oscillatory horizontal geometry of the disc associated with the non-circular orbital motion. Barker & Ogilvie (2014) revisited the problem and found that the growth rate was significantly enhanced by the vertical bouncing of the disc driven by the variation of vertical gravity around an elliptical orbit. In this work, we are finding a related instability that derives solely from the bouncing motion and not from any non-circular orbital motion.

4 SIMULATIONS OF UNFORCED OSCILLATIONS

Having presented our theoretical model in the previous section, we now carry out simulations of free (unforced) vertical oscillations. We initialize each simulation with an initial condition in which the Gaussian disc profile (see equation 3) has been stretched in the vertical direction. We perturb each velocity component at initialization with random perturbations. We begin with a description of our fiducial simulation in Section 4.1, and identify the main dissipative mechanisms at work during a bounce by means of an energy analysis. We then discuss a series of quasi-2D simulations, which enable us to check the key results at much higher resolutions (Section 4.2). Finally, in Section 4.3 we systematically investigate the effects of changing various physical and numerical parameters such as the initial oscillation amplitude, vertical box size, vertical BCs, and explicit versus numerical viscosity.

4.1 Fiducial simulation

In this section, we describe our fiducial simulation. The disc has been stretched in the vertical direction at initialization, so that the initial density profile is |$\rho (z,t=0) = \rho _0 \exp (-z^2/2H(0)^2)$|⁠, where we take |$H(0) = 2H_0$|⁠.⁶ At initialization, we perturb all three velocity components with (different) random noise of amplitude about |$0.05 c_\text{s}$|⁠. The box size is |$[L_x,L_y,L_z] = [8,4,12]H_0$| at a resolution of 32 cells per |$H_0$| in all three directions. We use reflective BCs in the vertical direction, and run the simulation for 50 orbits (⁠|$\sim 314\, \Omega ^{-1}$|⁠). Note that we repeated this fiducial simulation at double the resolution (⁠|$64/H_0$|⁠) and obtained very similar results. We also carried out a much more comprehensive resolution study, but in 2D, which we will discuss in Section 4.2.

4.1.1 Flow field

In Fig. 4, we plot the density in the |$xz$|-plane from snapshots at different times to illustrate the evolution of our initially out-of-equilibrium disc. After initialization (the top left-hand panel of Fig. 4), the disc rapidly collapses towards the mid-plane, reaching a thickness of around |$0.3H_0$| after around 0.33 orbits (top-middle panel), before expanding again due to pressure support. For the first six orbits, the disc continues oscillating in this manner. We observe two shocks propagating towards the disc mid-plane from the reflective vertical boundaries shortly after the disc reaches its maximum thickness (we check the effects of more realistic outflow BCs below). The bounce amplitude decreases over the course of each bounce, with the dissipation dominated by shocks as discussed below.

However, around orbit 7 the disc begins to exhibit a noticeable warp or bending wave (top right-hand panel), which develops into a standing wave as the disc continues to bounce. This wave grows as a result of the parametric instability described in Section 3.2, where it corresponds to the |$n=1$| mode. As the amplitude of the warp grows we observe counter-propagating radial fluid motions (which arise as a result of the radial pressure mismatch between neighbouring fluid columns in the bending wave). These eventually become so strong that they compress the disc in the radial direction by a factor of nearly 3.5 around orbit 10 (the bottom left-hand panel),⁷ corresponding to non-linear saturation of the parametric instability (see below). After non-linear saturation the disc continues to oscillate in the vertical direction (albeit with reduced amplitude), while maintaining a time-dependent radially periodic warp (i.e. a corrugation/bending wave) in the form of a standing wave. The mid-plane and atmosphere of the disc during this period are characterized by radially propagating shocks and waves (the bottom middle and bottom right-hand panels in Fig. 4).

To investigate the structures, we observe in the 2D flow field in greater detail; in Fig. 5, we show radial slices through |$z = 2H_0$| and |$z = 4H_0$| taken from a snapshot shortly after non-linear saturation (the bottom-middle panel of Fig. 4). At |$z = 2H_0$| (bottom two panels), we find four discontinuities in the density at |$x/H_0 = -2.64, 0.23, 2.03$| and 3.51, respectively. These coincide with discontinuities in the radial velocity. The profiles resemble N-waves resulting from wave steepening and are likely shocks.⁸ In the disc atmosphere (top two panels), we find clear evidence of strong shocks at the radial locations |$x = -3.68H_0$| and |$x=1.46H_0$|⁠, respectively. Here, the motion is supersonic with the radial velocity at the discontinuities reaching around |$2 c_\text{s}$|⁠. (Note that we have repeated this analysis in one of our quasi-2D simulations at a much higher resolution of |$512/H_0$|⁠. The wave- and shock-fronts are sharper at higher resolution, but the disc structure is the same.)

4.1.2 Time-series

To quantify the behaviour we observe in the flow field, we next turn to the top panel of Fig. 6 where show the time-series of the dynamical scale height (see equation 6). After initialization the disc rapidly collapses, reaching a thickness of around |$0.3H_0$| half-way through the first bounce. Pressure support then causes the disc to expand. The oscillation or bounce period |$\Delta T$| is around 0.66 orbits or around 1.5 bounces per orbit, which is consistent with the theoretical prediction (Fig. 2) when |$H_\text{max} = 2$|⁠. The oscillation amplitude (cf. equation 7) decays gradually over the first few orbits of the simulation from an initial amplitude of |$\Delta H \sim 1.67H_0$| at initialization to around |$1.55 H_0$| at orbit 8. Thereafter the amplitude of the oscillation drops rapidly over several orbits, reaching a minimum of |$\sim 0.55 H_0$| around orbit 10. As we will describe shortly, this occurs during the nonlinear saturation of the parametric instability during which a warp (bending mode) is seeded in the disc, which feeds off the bouncing motion. Surprisingly, during the non-linear phase of the instability (orbits 10–20) we observe periods during which the oscillation is ‘re-invigorated’ or increased over several orbits. The amplitude of the oscillation is smaller in each successful re-invigoration phase, indicating that the bouncing motion is decaying with time, as expected since there is nothing to drive the oscillation in this simulation. In Appendix  A, we provide a mathematical model of coupled and damped oscillators that can explain the observed re-invigoration.

Next, we turn to the time-evolution of the volume-averaged radial kinetic energy density, shown for this simulation in the bottom panel of Fig. 6, which is a good diagnostic for tracking the development of the parametric instability of the background vertically oscillating flow. On short (orbital) time-scales this oscillates with the oscillation of the disc scale height |$H(t)$| (top panel). More markedly, however, it exhibits exponential growth between around orbits 4 and 10, before saturating at around |$7\times 10^{-2}$| (in units of |$c_\text{s}^2 \rho _0$|⁠). Note that by the time the parametric instability starts to grow, the maximum thickness has been damped from the initial value of |$H(0) = 1.97H_0$| to |$H_\text{max} \sim 1.8H_0$|⁠. The growth rate of the amplitude (i.e. half the growth rate of the radial kinetic energy) during the linear phase is around |$0.169\Omega ^{-1}$|⁠, in very good agreement with value predicted from theory of |$0.17\Omega ^{-1}$| for an oscillation with |$H _\text{max}= 1.8H_0$| (see the right-hand panel of Fig. 3). After nonlinear saturation, the radial kinetic energy density also exhibits phases of decay and re-invigoration, similar to (but out of phase with) the dynamical scale height. This predator–prey-type behaviour can be understood as follows. First, the vertical oscillation (bouncing) excites a warp/bending mode by means of parametric instability, leading to a decrease in the oscillation amplitude, but an increase in radial kinetic energy. The warp/bending mode then excites vertical oscillations due to the combination of tilt and shear, as explained in the last paragraph of Section 1.1.1. This cyclic behaviour is damped, however, because of the dissipation associated with shocks.

Finally, we have also looked at the time-evolution of the various Reynolds stress tensor components (not shown). During most of the linear phase the Reynolds stress is zero. Before and after non-linear saturation, however, the stress becomes highly oscillatory, peaking in amplitude at 0.04 (in units of |$c_\text{s}^2 \rho _0$|⁠) around non-linear saturation of the instability at orbit 10. Overall, the stress oscillates around zero indicating that the instability does not lead to a net transport of angular momentum; averaged from orbits 5 to 20, we find |$\langle \langle R_{xy} \rangle \rangle _t \sim (-5\times 10^{-5}) \pm (7\times 10^{-3})$|⁠. The |$xy$|-component of the stress dominates, but we also observe oscillations of non-negligible amplitude in the |$xz$|- and |$yz$|-components of the stress.

4.1.3 Energetics: first few bounces

To determine the mechanisms damping the oscillations before the parametric instability sets in, we conduct an energy analysis of the first few bounces (an oscillation or bounce being defined as the interval between successive peaks in the dynamical scale height |$H$|⁠). During this interval the disc undergoes six full oscillations. The bounce amplitude decreases by 6.8 per cent from |$\Delta H/H_0 \sim 1.673$| to |$\sim 1.559$|⁠, while the total energy decreases by around 12.5 per cent. The main energy terms of interest at this stage are the vertical kinetic energy density |$E_{\text{kin},z}$| (which tracks the vertical oscillatory motion of the disc as it collapses and expands periodically), the Helmholtz free energy |$c_\text{s}^2 \rho \ln {\rho }$| (which takes into account losses/gains in total energy due to compression/expansion), and the total energy itself |$E_\text{tot}$| (see equation 8), so we will focus on these terms.

During the very first half oscillation (as the disc collapses from its initial thickness of |$H \sim 1.97H_0$| at orbit 0 down to |$H \sim 0.297H_0$| at orbit 0.33, the total energy actually increases very slightly over the first 0.16 orbits. Total energy can only increase due to a positive Reynolds stress extracting energy from the radial shear, and indeed during these very early stages of the simulation we observe a sharp spike in the Reynolds stress (not shown) soon after initialization, reaching |$\sim 3\times 10^{-6}$|⁠. This probably results from the transient amplification of non-axisymmetric perturbations present in the initial noise. The stress then drops rapidly, reaching around |$\sim 5\times 10^{-7}$| at orbit 0.16. Over the course of the next 4 orbits, the Reynolds stress oscillates but does not exceed |$\sim 5\times 10^{-8}$|⁠. Therefore, the stress plays a negligible role in the energy budget during this interval, so we will restrict our subsequent energy budget analysis to the second and third oscillations, in order to avoid the transient behaviour immediately after initialization.

The second oscillation/bounce (orbit 0.66–1.33) begins with the disc at a maximum thickness of |$H \sim 1.95H_0$| at orbit 0.66 and ends with the disc at maximum thickness of |$H \sim 1.93H_0$| one cycle later. (Subsequent oscillations exhibit similar behaviour.) We plot the total energy over the course of this oscillation in the top panel of Fig. 7. Other energy components are plotted in the bottom panel. The overall behaviour of the key energy terms over the course of the bounce is as follows: (i) the total energy decreases monotonically by around 2.5 per cent over the course of the entire oscillation (orbit 0.66–1.32), from |$\sim 0.391$| to |$\sim 0.381$| (in units of |$c_\text{s}^2 \rho _0$|⁠); (ii) the vertical kinetic energy density (red curve, bottom panel) increases from zero to a maximum value and decreases back to zero again during each half cycle; (iii) the Helmholtz free energy (dashed orange curve, bottom panel) increases as the disc collapses, peaks when the disc is most compressed, and decreases again as the disc expands.

Although the total energy decreases monotonically during the course of the bounce, it does not do so at a constant rate: in the top panel of Fig. 7, we have identified four different stages. These four stages can be understood by combining information about the flow field (density |$\rho$| and vertical velocity |$u_z$|⁠) and the total energy. In Fig. 8, we show 2D (in the |$xz$|-plane; left-hand panel) and 1D vertical slices (through |$x=0$|⁠; right-hand panel) of the density logarithm (top) and |$u_z$| (bottom), taken at orbit 1.09 (as the disc is expanding). The middle panels show the evolution of the disc thickness (top) and total energy density (bottom) over the course of the entire bounce (orbit 0.66–1.33). During the initial stage of disc contraction (orbit 0.66–0.875), a shock propagates down from the vertical boundaries through the atmosphere. This is associated with a mild decrease in total energy due to shock dissipation. When the shock hits the bulk of the disc (around orbit 0.875) the dissipation rate increases as it passes through denser gas. The disc starts to expand after orbit 1, during which the total energy is nearly constant. The rate of energy dissipation then increases again after orbit 1.15. This coincides with time at which highly supersonic fluid (Mach number |$\text{M}\sim 8$|⁠) hits the rigid outer boundary, which launches a new shock from the boundary. Overall, the total energy density decreases by around 3 per cent during the course of the bounce, all of which is due to shock dissipation and ultimately related to the interaction of the bouncing motion with the unphysical vertical boundaries.

4.1.4 Energetics: non-linear saturation

Next, we analyse a bounce during nonlinear saturation of the parametric instability. We focus on a vertical oscillation between around orbit 9.35 and orbit 10.01 (see Fig. 9). At this point, the disc is noticeably warped and strong compression and rarefaction can be seen in the disc due to radial motion induced by the instability. The dynamics is complicated: in addition to the periodic vertical oscillation of the disc and the behaviour and dissipation associated with that motion (described in the previous subsection), the warped/corrugated disc is now also oscillating as a standing wave and this corrugation drives radial flows.

The bounce amplitude decreases by around 22 per cent during this cycle, while the total energy decreases by around 25 per cent, a marked increase in dissipation compared to when the disc was undergoing pure vertical oscillations (orbits 0–4). Interestingly, the total energy (Fig. 9, bottom panel, black curve) does not decrease monotonically over the course of the oscillation, but actually increases as the disc compresses (roughly orbits 9.35– 9.68). Since the only source term in the total energy equation is flux of energy through radial boundaries due to the Reynolds stress, this must be due to a non-zero Reynolds stress. Indeed, we find the Reynolds stress (not shown) is positive over this short time interval. The implication is that the unstable bending wave gains energy from the Keplerian shear (via the Reynolds stress) during this half cycle over which the total energy is increasing. [A related phenomenon can occur in the vertical shear instability, in which an axisymmetric inertial wave gains energy both from the vertical shear and from the Keplerian shear (Ogilvie, Latter & Lesur, in preparation).]

The vertical component of the kinetic energy density (solid red curve) is non-zero even when the disc is maximally compressed, on account of vertical motion associated with the corrugation rather than the background oscillation of the disc. This can be seen more clearly by inspecting the perturbed part of the vertical kinetic energy density, defined as |$E_{\text{kin},z,\text{pert}} \equiv \langle (1/2)\rho (u_z-\langle u_z \rangle _{xy})^2 \rangle$| (dashed red line), which peaks at around orbit 9.68. At this point in the oscillation cycle, the disc has actually reached its minimum thickness and has not begun to expand yet, so there is no vertical motion associated with expansion/contraction of the disc itself. Instead, vertical motion during this phase is driven by the corrugation as the disc evolves from planar (at orbit 9.60) to warped (easily visible at orbit 9.68). Finally, we also show the (perturbed) azimuthal kinetic energy (green dashed curve). This is out of phase with the radial kinetic energy (solid red curve) due to epicyclic motions (conservation of angular momentum).

4.2 Quasi-2D simulations

To determine whether non-axisymmetry is important to the problem, we have run a series of quasi-2D simulations. See Table B4 in Appendix B. The reduced dimensionality also enabled us to probe much higher resolutions, from |$32/H_0$| (the same resolution used in our fiducial 3D simulation described in the previous section) up to |$512/H_0$|⁠. The box size across all simulations was fixed at |$[L_x, L_y, L_z] = [8,0.008,12]H_0$|⁠. Furthermore, we employed reflective BCs in the vertical direction, and included a small explicit viscosity (⁠|$\text{Re} = 4687$|⁠), as in our fiducial 3D simulation.

In Fig. 10, we show the time-evolution of the dynamical scale height |$H(t)$| (top) and volume-averaged radial kinetic energy density |$E_{\text{kin},x}$| (bottom). Reassuringly, we find almost no difference in the bounce period and amplitude between our lowest resolution run (32 cells per equilibrium scale height |$H_0$|⁠) and highest resolution run (⁠|$512/H_0$|⁠). During the first half-cycle, the minimum disc thickness is around |$H \sim 0.26$|⁠. This corresponds to around 133 cells across the disc at the highest resolution, so the disc remains well resolved even when it is most compressed. Parametric instability sets at slightly later times as the resolution is increased, on account of the lower level of initial noise in higher-resolution simulations.⁹ The growth rate of the instability is very similar regardless of the resolution, e.g. |$s/\Omega \sim 0.177$| and 0.175 at resolutions of |$32/H_0$| and |$512/H_0$|⁠, respectively. The level of non-linear saturation and the behaviour after non-linear saturation are very similar at all resolutions.

We compared the flow field (density in the |$xz$|-plane) at different resolutions in these quasi-2D simulations (not shown). Regardless of the resolution, after nonlinear saturation the disc exhibits corrugation/warping and radially propagating waves and shocks. At higher resolutions we can see finer structure where the shocks meet (at |$x \sim -3H_0$| and |$x \sim 1.5H_0$|⁠), but the overall structure of the flow is the same.

In conclusion, our high-resolution simulations indicate that the resolution of |$32/H_0$| used in our fiducial 3D simulation is sufficient to capture all the salient physics associated with the vertical oscillations, including damping of the oscillations by shocks, the development of parametric instability, and the subsequent corrugation of the disc. The fact that the quasi-2D simulations reproduce all the physics of the fiducial fully 3D simulation confirms that the problem is axisymmetric (which we also determined directly through visual inspection of density snapshots in the |$yz$|-plane in our fiducial 3D simulation).

4.3 Effect of numerical and physical parameters on the bouncing motion

4.3.1 Effect of initial bounce amplitude

In Fig. 11, we compare three different initial disc thicknesses: |$H(0) = 1.10, 1.22$|⁠, and 1.97, which we refer to as ‘small’, ‘intermediate’ and ‘large’-amplitude oscillations. We carried out the small-amplitude run in a slightly larger box (⁠|$L_x = 9H_0$| compared to |$L_x=8H_0$|⁠) so that the box size fits comfortably within the peak of the |$n=1$| mode predicted by theory (see the left-hand panel of Fig. 3). As predicted by the theory (see Fig. 2), the initial oscillation amplitude determines the subsequent oscillation period. We find that the disc oscillates |$\sim 1.41$| times per orbit in the small-amplitude case, exactly as predicted by theory in the limit of small oscillations. For the intermediate-amplitude case the disc oscillates |$\sim 1.43$| times per orbit, while in the large-amplitude case the disc oscillates |$\sim 1.51$| times per orbit.

All three simulations exhibit parametric instability, with growth rates |$s/\Omega \sim 0.024$|⁠, 0.044 and 0.169, respectively, in good agreement with the theoretically predicted values (see Table B1). The onset of non-linear saturation is delayed as the initial bounce amplitude is decreased, and does not occur until 100 orbits after initialization in the small amplitude case (⁠|$H(0) \sim 1.10H_0$|⁠). For this very weakly oscillating case, the radial kinetic energy density saturates at around |$10^{-3}$| (in units of |$c_\text{s}^2 \rho _0$|⁠), about 10 times smaller than the saturation amplitude in the large-amplitude oscillation simulation (⁠|$H(0) \sim 1.97H_0$|⁠).

4.3.2 Effect of vertical box size

Since we employ reflective BCs, the periodic expansion of the disc results in reflection of material at the vertical boundaries. Ideally we want the box to be as tall as possible to avoid significant boundary effects such as this from interfering with the dynamics closer to the bulk of the disc, which lies within |$|z|\lt 2H_0$|⁠. Thus, we have repeated our fiducial simulation in shorter (⁠|$L_z = 8H_0$|⁠) and taller (⁠|$L_z = 16H_0$|⁠) boxes (keeping the resolution the same). See Table B2 in Appendix B. In Fig. 12, we show the results of this study with the radial kinetic energy density in the top panel and the dynamical scale height |$H(t)$| in the lower panel. Note that in the shortest box a significant fraction of the disc atmosphere is cut off by the vertical boundaries and the initial dynamical scale height is therefore |$H(0) \sim 1.67H_0$|⁠, compared to |$H(0) \sim 1.97$| and |$\sim 1.99$| in the |$L_z=12H_0$| and |$16H_0$| boxes, respectively. We see that the results are converged for boxes of size |$L_z \gt 8H_0$|⁠, with the growth rate of the parametric instability and evolution of the bounce amplitude nearly identical in the |$L_z = 12 H_0$| (fiducial run) and |$L_z = 16H_0$| boxes. In the shortest box |$L_z = 8H_0$| (Fig. 12, black curve), however, first the oscillation amplitude is smaller to begin with because a significant fraction of the initial disc profile is cut off by the boundaries, and second the oscillation amplitude is strongly damped even before the parametric instability starts to grow because the boundaries are lower. Thus, the growth rate of the instability is significantly reduced in the shortest box.

4.3.3 Effect of vertical BCs

To compare vertical boundary condition effects (zBCs), we have repeated our fiducial simulation (Section 4.1) using outflow rather than reflective BCs in the vertical direction. See Table B3 in Appendix B. The dynamical scale height and radial kinetic energy density are shown in .

Overall the solution in the simulation with outflow BCs is very similar to that with reflective BCs (see Fig. 13), although we do find small qualitative differences between the two. The oscillation is damped somewhat faster in the outflow case during the linear phase, although by the end of the simulation the oscillation amplitude is very similar in the two cases. The linear growth rate is somewhat larger in the reflective case (⁠|$0.169\Omega$|⁠, compared to |$0.154\Omega$| in the outflow case), which means the reflective zBC simulation reaches nonlinear saturation earlier. This results in a slight phase difference (of a few orbits) in the oscillation of the disc between the two simulations by orbit 20.

Although the overall result over 20 orbits is very similar in the two cases, we find more marked differences between the two runs when we focus on the first few bounces (before parametric instability sets in). We have carried out an energy budget analysis of the outflow simulation during the first four orbits (during which the disc bounces around six times), similar to the one we presented for the fiducial simulation with reflective zBC (see Section 4.1.4). This interval is before the onset of parametric instability, so the disc is undergoing only vertical oscillations. The total energy in the outflow run drops by around 25 per cent between orbits 0 and 4 (6 oscillations) compared to a drop of around 13 per cent in the reflective case.

In the outflow zBC simulation, we did not use a mass source term (to avoid interfering with the vertically oscillating flow), and thus mass and energy were able to leave the box through the vertical boundaries. Indeed, over the entire duration of the simulation (50 orbits) around 2.2 per cent of the initial mass in the box was lost through the vertical boundaries due to the periodic expansion of the disc, although most of this mass-loss occurs during the first 10 orbits (linear phase), with over half of the total mass-loss (1.2 per cent of the initial mass) taking place during the first four orbits alone. The rapid drop in mass-loss after orbit 10 coincides with the onset of non-linear saturation of the parametric instability, which significantly damps the oscillations, and thus the resultant mass-loss through the vertical boundaries.

The behaviour of the total energy over the second bounce (orbits 0.66–1.33) in the outflow simulation is shown in Fig. 14. The total energy is constant during the initial stages of disc contraction since there are no shocks from the boundaries to dissipate the energy. A shock eventually does form as the disc contracts supersonically, which results in a rapid drop in the total energy during the final moments of contraction and the subsequent initial moments of expansion (orbits 0.907–1.05). When this shock (moving at Mach number |${\textit M}\sim 8$|⁠) hits the outer boundary, there is another sharp drop in the total energy, this time because mass and energy are advected across the vertical boundary. Overall, the total energy density decreases by around 5 per cent during the course of the oscillation (most of that due to loss of energy through the vertical boundaries).

4.3.4 Effect of explicit viscosity

In our fiducial 3D run and in most of our quasi-2D runs, we included a small explicit viscosity, corresponding to a Reynolds number of |$\text{Re} = 4687$|⁠. To ensure the resolution we used is large enough to capture the viscous length scale in our non-ideal simulations, and if so, to determine the effect of the explicit viscosity on the dissipation of oscillations, we have tested the effects of explicit and numerical viscosity. We do so primarily using quasi-2D simulations since these enable us to achieve much higher resolutions.¹⁰ At |$32/H_0$|⁠, we find that the oscillation amplitude and growth rates are nearly identical whether we include explicit viscosity or not. However, the viscous scale [which is given by |$l_{\nu } \sim 2\pi /\sqrt{\text{Re}}$|⁠; see discussion in Held & Mamatsashvili (2022)] in the non-ideal run is poorly resolved at |$\text{Re}=4687$| (⁠|$\sim 3$| cells per viscous length scale). Thus, we have repeated our ideal and non-ideal runs at |$256/H_0$|⁠. In this case, the viscous length scale is well resolved (⁠|$\sim 25$| cells). Again we find very little difference between the evolution of the dynamical scale height, the parametric instability growth rate (⁠|$\sim 0.17\Omega$| in both cases), or any of the dynamics during the nonlinear phase between the ideal and non-ideal runs. Thus, we conclude that the explicit viscosity we have included in our simulations plays a negligible role in dissipating the vertical oscillations.

5 SIMULATIONS OF FORCED OSCILLATIONS

In the previous section, we analysed the behaviour of a disc undergoing free oscillations in the vertical direction in order to isolate the physical and numerical mechanisms that dissipate the oscillation. We now turn to the case in which the vertical motion is forced. This can occur when the vertical component of gravity varies along a streamline, for example in an eccentric disc, or in a disc that is tidally distorted by a binary companion. In a warped disc, the forcing arises due to the combination of tilt and shear even in the absence of an external perturber (Ogilvie & Latter 2013; Fairbairn & Ogilvie 2021).

In order to drive vertical oscillations, we introduce a time-dependent gravitational acceleration in the vertical direction of the form |$g_{\text{eff},z} = -\Omega ^2 z (1+a\cos {\omega t})$| (see Section 2.1). This allows us to set the forcing frequency |$\omega$| (in units of the angular frequency of the disc |$\Omega$|⁠) and the dimensionless forcing amplitude |$a$|⁠. Thus we have two parameters which we can use to control the forcing. We first discuss astrophysically relevant forcing frequencies and forcing amplitudes in Section 5.1. Then, in Section 5.2, we extend our 1D model for vertical disc oscillations to include forcing, and calculate the expected parametric growth rates for six select forcing frequencies. We also calculate the amplitude of oscillations of the disc scale height |$H(t)$| for each forcing frequency (see Fig. 15), which we subsequently use to determine the appropriate initial condition for each simulation. Finally, we discuss fully 3D simulations with forcing in Section 5.3.

5.1 Astrophysical applications: forcing frequencies and amplitudes

5.1.1 Forcing frequencies

Let us first consider which forcing frequencies are relevant for each type of distorted disc (eccentric, tidally distorted, and warped).

In an eccentric disc, the distance to the central star varies along each eccentric streamline. The vertical component of gravity due to a central point source is strongest at pericentre and weakest at apocentre, so the forcing frequency experienced by the gas in an eccentric disc is simply equal to the orbital frequency (⁠|$\omega = \Omega$|⁠), if we neglect any precession of the eccentricity. A sinusoidal variation of the vertical gravity is a good model for small eccentricities.

In a tidally distorted disc (i.e. a circumstellar disc orbited by a binary companion, or a circumbinary disc), the distance – and thus the vertical component of gravity – also varies along any particular orbit, but the forcing frequency is a function of cylindrical radius |$r$|⁠. Consider a circumstellar disc in a circular binary system: in the frame rotating with the binary angular frequency |$\Omega _\text{b}$|⁠, the disc can have a stationary, non-circular shape as a result of the tidal force due to the companion (see Fig. 1). We can decompose the non-axisymmetric part of the gravitational potential into Fourier modes with azimuthal mode number |$m$| (Lubow 1981; Ogilvie 2002a, b). At any given cylindrical radius |$r$| in the disc, the angular frequency of the disc relative to the tidal deformation is |$\Omega (r)-\Omega _\text{b}$|⁠, where |$\Omega (r) \equiv (GM/r^3)^{1/2}$| is the Keplerian angular frequency of the disc (ignoring effects due to the companion). The forcing frequency (in units of the disc angular frequency) associated with mode number |$m$| is given by |$\omega /\Omega = m(1-\Omega _\text{b}/\Omega)$|⁠. Lubow (1981) showed that vertical resonances occur where |$m(1-\Omega _\text{b}/\Omega) = \sqrt{\gamma +1}$| and the vertical oscillation mode is excited. In practice |$m = 2$| is the most important mode, both because of the need to fit the resonance inside the disc and because |$m = 2$| is the strongest component of the tidal deformation. Note that the forcing frequency decreases as we move radially outwards from the primary. Closest to the central star, the forcing frequency (for the |$m=2$| mode) is |$\omega \sim 2\Omega$| (since |$\Omega \gg \Omega _\text{b}$|⁠), although here the deformation of the orbit is smallest and so the forcing actually vanishes. Farthest from the central star (approaching the 2:1 resonance, which actually lies outside the disc), |$\Omega \gtrsim 2\Omega _\text{b}$| and so |$\omega \gtrsim \Omega$|⁠. The resonance should lie somewhere in between these two locations, but within the disc.¹¹

Finally, a warped disc can be thought of as a series of nested circular orbits, each with a different tilt, so the vertical component of gravity does not vary along any given streamline. However, the combination of tilt and shear experienced by the gas along any given orbit results in periodic vertical compression and rarefaction as if the vertical component of gravity were periodic (Ogilvie & Latter 2013; Ogilvie 2022). The forcing frequency is equal to twice the orbital frequency (⁠|$\omega = 2\Omega$|⁠), if we neglect any precession of the warp.

5.1.2 Forcing amplitudes

In the case of an eccentric disc, the dimensionless forcing amplitude |$a$| can be directly related to the eccentricity |$e$| through |$a=3e$|⁠. The factor of 3 comes from the fact that the vertical gravity scales with |$r^{-3}$|⁠.

Next, we obtain an estimate of the appropriate forcing amplitude |$a$| for a tidally distorted disc. The gravitational (plus indirect) potential due to a primary star of mass |$M_1$| and secondary star of mass |$M_2$| in a circular binary of separation |$r_\text{b}$| is given by (Lubow 1981)

where |$\Omega _\text{b} = (GM/r_\text{b}^3)^{1/2}$| is the binary angular frequency, |$M=M_1+M_2$|⁠, and the relative azimuthal coordinate is |$\phi ^{\prime } = \phi -\Omega _\text{b} t$|⁠. Expanding the non-axisymmetric part for |$r \ll r_\text{b}$| (at |$z=0$|⁠), we obtain |$\Phi _2 \approx -(3/4)G M_2 (r^2/r_\text{b}^3) \cos {2\phi ^{\prime }}$|⁠. This gives a horizontal force that distorts the orbital motion from a circular form, causing the vertical gravity due to |$M_1$| to vary around the non-circular orbit. (There is also a weaker component, beyond the quadrupolar approximation, that gives a direct contribution to the vertical gravity from |$M_2$|⁠. Both effects scale linearly with the mass ratio |$q=M_2/M_1$|⁠.)

The radial displacement of the tidal deformation can be calculated (neglecting pressure) from

where |$m=2$| and |$\omega =2(\Omega -\Omega _\text{b})$|⁠. This gives

where

with |$\tilde{r}\equiv r/r_\text{b}$|⁠. Since the vertical gravity due to |$M_1$| is |$GM_1z/r^3\propto r^{-3}$|⁠, the amplitude of its fractional variation (assumed small) is three times the amplitude of |$\xi _r/r$|⁠, i.e.

For |$q=1$| and |$\tilde{r}=r/r_\text{b}$| between 0.2 and 0.3, |$a$| is approximately between 0.04 and 0.2, so |$a \sim 0.1$| is a plausible forcing amplitude.

5.2 Solutions from theory

In Section 3, we presented a 1D theory of vertical oscillations, resulting in an evolution equation for the disc thickness or dimensionless scale height |$H(t)$| (equation 15). We now modify this equation to include the forcing term |$a \cos {\omega t}$|⁠:

(From here onwards, we will consider only the case |$\gamma = 1$| since our simulations are isothermal.) To determine appropriate initial conditions for the simulations we solve this equation for a given forcing amplitude |$a$| and frequency |$\omega$|⁠. The response curve for the scale height |$H(0)$| (⁠|$t=0$| being the time when forcing is maximal) is shown in Fig. 15 for forcing frequencies between |$\omega /\Omega = 1.0$| and |$\omega /\Omega = 2.0$|⁠. We plot the response curve for two cases: small-amplitude forcing (⁠|$a=0.01$|⁠, blue curve) and large-amplitude forcing (⁠|$a=0.1$|⁠, red curve). Linear resonance (Lubow 1981) occurs at |$\omega = \sqrt{2} \approx 1.41$| (indicated by the dashed vertical line). Below the resonant frequency, only one (stable) solution branch exists. At frequencies higher than the resonant one, multiple solutions may exist: there is a weakly oscillating upper branch, and a strongly oscillating lower branch, and hysteresis may be possible. (There are also unstable solutions denoted by the dashed lines.) The non-linearity of the resonance can be seen in the fact that the endpoint of the upper branch (which in fact folds over and becomes unstable at this point) moves to the right as the forcing amplitude is increased.

Next we analyse the stability of the forced bouncing solutions to radially dependent perturbations, in the same manner as we did for the freely bouncing case (see Section 3.2). In Fig. 16, we plot the growth rates for six different forcing frequencies (at a constant forcing amplitude of |$a=0.1$|⁠). We show the growth rate as a function of radial wavenumber for the first ten modes (highlighting |$n=1$| [blue], |$n=2$| [green] and |$n=3$| [red]). At forcing frequencies below resonance (⁠|$\omega /\Omega = \sqrt{2}$|⁠), increasing the forcing frequency results in an increase in the growth rate of the fastest growing mode, while the peak simultaneously broadens and moves to the left (i.e. to longer wavelengths). On the other side of the resonance, the growth rate of the fastest growing mode decreases with increasing forcing frequency (as does the oscillation amplitude on the upper branch), and the peak narrows. At |$\omega /\Omega = 2$|⁠, we find a very narrow peak in the limit of very large wavelength, sitting on top of a broader peak. We will analyse this behaviour in more detail when we discuss the corresponding simulation (Section 5.3.6).

5.3 Simulations at different forcing frequencies

Having presented the theory of forced (driven) vertical oscillations, we turn to the simulations. We discuss, in turn, six runs with forcing frequencies ranging from |$\omega /\Omega = 1$| (one oscillation or bounce per orbit), which is relevant to an eccentric disc, or the outer part of a tidally distorted single disc, to |$\omega /\Omega = 2$| (two bounces per orbit), which is relevant to a warped disc, or the inner part of a tidally distorted single disc). Forcing frequencies between these two limits correspond to different radial locations in a tidally distorted disc. We used a forcing amplitude of |$a = 0.1$| for all runs. Note that we adjust the radial box size in each simulation to ensure that we can capture the |$n=1$| mode of the parametric instability for that forcing frequency (see Fig. 16). We initialize each run with a disc thickness |$H(0)$| that is equal to (or very close to) the initial thickness predicted from theory for that forcing frequency and amplitude (blue curve in Fig. 15), because otherwise an additional free oscillation with a different period is superimposed on the forced oscillation. In the azimuthal and vertical directions, we keep the box size fixed at |$L_y = 4H_0$| and |$L_z = 12H_0$|⁠, respectively. Finally, we use a resolution of |$32/H_0$| for all runs, and run each simulation for 100–200 orbits (628–|$1257\, \Omega ^{-1}$|⁠). We show the results (time-evolution of the dynamical scale-height at different forcing frequencies) in Fig. 17.

5.3.1 The case |$\omega /\Omega = 1$|

The lowest forcing frequency we investigated is the case where the forcing frequency is equal to the local orbital frequency, i.e. |$\omega /\Omega = 1$|⁠, which is far from resonance (top left-hand panel of Fig. 17). This could be the case for an eccentric disc around a single star (see e.g. Wienkers & Ogilvie 2018). The growth rates of the parametric instability for this forcing frequency and amplitude are shown in the top left-hand panel of Fig. 16. After initialization the disc oscillates once per orbit, and its thickness varies between |$H \sim 1.1H_0$| and |$\sim 0.9H_0$|⁠. We observe very little change in the amplitude of the bounce over the first 100 orbits. The radial kinetic energy density grows exponentially after a few orbits, but with a relatively weak growth rate of |$s/\Omega = 0.012$| (the predicted value from theory for a wavelength equal to the radial box size |$L_x=4H_0$| is |$s/\Omega = 0.014$|⁠). Non-linear saturation of the parametric instability does not occur until around orbit 110 at |$\langle E_\text{kin,x} \rangle \sim 7\times 10^{-4}$|⁠. Around this time, we can begin to see the disc develop a bending mode (or radially periodic warp/corrugation); however, even in the non-linear phase the behaviour is quite gentle compared to our fiducial freely bouncing simulation, consisting only of small oscillations together with periodic corrugation even after 200 orbits.

5.3.2 The case |$\omega /\Omega = 1.3$|

Next, we consider a slightly larger forcing frequency of |$\omega /\Omega = 1.3$|⁠, corresponding to a radius larger than the radius at which linear resonance occurs in a tidally distorted single disc in a binary (top middle panel of Fig. 17). During the linear phase, the disc thickness |$H(t)$| varies between |$0.7H_0$| and |$1.35H_0$|⁠. In addition, the scale height exhibits a mild beating on a timescale of around 8–9 orbits. This arises due to the mismatch between the forcing frequency at which the oscillation is driven (⁠|$\omega /\Omega = 1.3$|⁠), and the natural or resonant frequency (⁠|$\sqrt{2} \approx 1.41\Omega$|⁠).¹²

The oscillating disc undergoes parametric instability with a measured growth rate of |$0.041 \Omega$| in the linear phase, compared to the predicted value of |$\sim 0.055 \Omega$| at |$kH = 2\pi /8 \approx 0.78$|⁠, and non-linear saturation occurs around orbit 35. Around this time, the oscillation is considerably damped by shocks resulting from radially converging flows around |$x = \pm 2H_0$| due to the pressure misalignment between neighbouring fluid columns in the corrugated disc. Note that after nonlinear saturation the dynamical scale height oscillates around a slightly larger value (⁠|$1.06H_0$|⁠) than before non-linear saturation, which is due to this diagnostic picking up the standing-wave corrugation of the disc. During the nonlinear phase we observe many of the features that we observed in our fiducial unforced simulation (Section 4.1): the disc continues to expand and contract in the vertical direction (due to the forcing) and simultaneously oscillates as a standing wave. On top of these large-scale motions, we observe many radially propagating waves and shocks.

5.3.3 The case |$\omega /\Omega = \sqrt{2}$| (resonance)

Next, we consider a disc oscillating at the resonant frequency |$\omega /\Omega = \sqrt{2} \approx 1.41$|(top right-hand panel of Fig. 17). The resonant forcing causes the disc to oscillate with a large amplitude |$\Delta H \sim 1.3H_0$|⁠, most closely resembling that in our fiducial unforced simulation (Section 4.1). Between initialization and around orbit 5 the time-averaged disc thickness appears to increase, possibly due to a very long beat phenomenon in this case (as the natural frequency of nonlinear oscillations differs from the forcing frequency by an amount that depends on the amplitude). The linear phase growth rate measured between orbits 5 and 10 is |$0.15 \Omega$|⁠, compared to the value of |$0.162 \Omega$| predicted from theory at |$k H_0 = 2\pi /8 \approx 0.78$|⁠. Around orbit 7, we see a corrugation begin to develop in the disc, and non-linear saturation sets in around orbit 10, significantly earlier than at other forcing frequencies, owing to the larger growth rate of the instability at resonance. After non-linear saturation, the disc exhibits strong shocks and radial compression, which significantly damp the oscillation amplitude (orbits 10–15).

Nevertheless, we observe at least five phases during which the oscillation is re-invigorated after nonlinear saturation, most noticeably between orbits 15 and 18, similar to what we observed in our fiducial freely bouncing simulation (Section 4.1.2). During each of these phases we see the disc become increasingly corrugated until the resultant radial shocks in the disc atmosphere are so strong that the subsequent bouncing is damped; this can be clearly seen from a snapshot (not shown) taken at orbit 19.257 which is taken half a cycle after the bounce amplitude peaks during the first of these re-invigoration phases. The disc is highly corrugated, and two intersecting waves near the mid-plane of the disc at |$x \sim -1H_0$| and |$\sim +3H_0$| have steepened into shocks (with Mach number |$M \sim 4$| in the disc atmosphere). The peak amplitude in the vertical oscillation during successive phases of re-invigoration is always weaker than the last. From orbit 50 onwards, the disc settles down into a quasi-steady equilibrium consisting of a complicated flow field with many intersecting waves and shocks.

5.3.4 The case |$\omega /\Omega = 1.6$|

Next, we turn to a simulation with a forcing frequency of |$\omega /\Omega = 1.6$|⁠, i.e. a forcing frequency that is greater than the resonant frequency (bottom left-hand panel of Fig. 17). In a tidally distorted single disc in a binary this is expected to occur at a radius inside the resonant radius. The disc undergoes parametric instability with a measured growth rate in the linear phase of |$0.033 \Omega$|⁠, compared to the theoretical prediction of |$\sim 0.041\Omega$| for the |$n = 1$| mode at this wavelength. Like at forcing frequencies less than the resonant frequency (⁠|$\omega /\Omega = \sqrt{2}$|⁠), we observe a beating of the dynamical scale height, but the envelope is more complicated than in the |$\omega /\Omega = 1.3$| simulation. We measure an average beat period (interval between successive peaks in the envelope) of around 5.8 orbits, compared to the expected value of around |$1/|1.6-\sqrt{2}| \sim 5.4$| orbits for this forcing frequency. Non-linear saturation sets in at around orbit 55, after which the time-averaged dynamical scale height increases to |$\langle H \rangle _t \sim 1.08H_0$|⁠. During the non-linear phase the behaviour of the dynamical scale height is complicated, exhibiting alternating phases during which the amplitude of the oscillation increases (e.g. between orbits 55 and 65 the oscillation amplitude increases from around |$\Delta H \sim 0.1H_0$| to |$\sim 0.25 H_0$|⁠), and then decreases again (orbits 65–75).

5.3.5 The case |$\omega /\Omega = 1.8$|

At a forcing frequency of |$\omega = 1.8 \Omega _0$| the disc first oscillates with an amplitude of around |$\Delta H \sim 0.2 H_0$| and a beat of |$\sim 2.6$| orbits (compared to the expected beat of |$1/(1.8-\sqrt{2}) \sim 2.6$| orbits). See bottom middle panel of Fig. 17. We measured a growth rate of |$\sim 0.019 \Omega _0$|⁠, in excellent agreement with that predicted by theory for this radial box size (⁠|$L_x = 30H_0$|⁠, i.e. |$kH_0 \sim 0.21$|⁠) of |$\sim 0.20 \Omega _0$| (see the middle panel of Fig. 16). The disc begins to corrugate around orbit 80 – this is development of the bending mode, as seen in our other simulations.

5.3.6 The case |$\omega /\Omega = 2$|

Finally we consider the case where the disc is forced to oscillate twice per orbit, relevant to a warped disc around a single star (bottom right-hand panel of Fig. 17). The shape of the dispersion relation for the |$n=1$| mode (we remind the reader that here |$n$| is the vertical mode number) is markedly different in this case compared to the other forcing frequencies (see the bottom right-hand panel of Fig. 16). Only the |$n=1$| mode grows, and only at long (⁠|$\gtrsim 21H_0$|⁠) wavelengths. Thus, we need a large radial domain to capture this mode: we chose a domain of size |$L_x = 30H_0$| in the radial direction, corresponding to a wavenumber of |$kH_0 \sim 0.2$|⁠. Note that this box size is rather extreme,¹³ but we include this run for completeness.

The simulation at this forcing frequency differs markedly from the earlier cases. The disc initially exhibits small-amplitude oscillations about the equilibrium value (⁠|$\Delta H \sim 0.1H_0$|⁠), in keeping with the theoretical prediction at this forcing frequency (Fig. 15). The radial kinetic energy density grows exponentially, but with a small growth rate of around |$s \sim 0.007\Omega$| (the theoretical growth rate at this box size/radial wavenumber is |$s_{\text{th}} \sim 0.01\Omega$|⁠). Unlike at smaller forcing frequencies, the disc does not begin to develop a warp as the mode grows. Instead, starting at around orbit 40, the entire disc begins to be oscillate uniformly about the mid-plane (⁠|$z=0$|⁠) twice per orbit even as the disc itself continues to expand and contract at the same rate. The amplitude of this oscillation continues to grow until, at around orbit 85, the disc hits the reflective vertical boundaries at |$z = \pm 6H_0$| and is subsequently disrupted.¹⁴

What we are seeing is the excitation of an |$n=1$| mode with radial wavenumber |$k=0$|⁠. The up-and-down motion (about |$z=0$|⁠) of the entire disc is how a tilted disc would appear to an observer moving along an untilted circular orbit, and can be related to the tilt instability at the 3:1 resonance¹⁵ in a tidally distorted disc in binary (Lubow 1992). This can understood quantitatively by deriving an evolution equation (we omit the details of this calculation) for the height of the centre of mass of the disc (defined as |$Z(t)=(1/M)\iint \rho z\, dx\, dz$|⁠, where |$M=\iint \rho \, dx\, dz$| is the total mass) with periodic forcing

which is a form of the Mathieu equation. In particular, parametric resonance occurs for |$\omega$| sufficiently close to |$2\Omega$|⁠. The expected growth rate is

where the quantity inside the square root is positive, i.e. where the detuning is sufficiently small, i.e. for |$\omega /\Omega$| between |$2\pm (a/2)$|⁠. The maximum growth rate of |$a\Omega /4$| is expected at |$\omega =2\Omega$|⁠. For a forcing amplitude of |$a=0.1$| this works out to be |$\sim 0.025\Omega$|⁠, corresponding to the spike near |$k=0$| in Fig. 16, and in good agreement with the value of |$\sim 0.022\Omega$| measured directly in the simulation from the time-series of vertical kinetic energy density (not shown).

If the case |$\omega =2\Omega$| is identified with the forcing of vertical oscillations in a warped disc, then to say that such oscillations can parametrically excite a large-scale warp illustrates the fact that there are nonlinear couplings between such a warp and vertical oscillations that can exchange energy is either direction. This is consistent with the picture derived in the oscillating torus model of Fairbairn & Ogilvie (2021).

6 CONCLUSIONS

We have studied vertical oscillations in accretion discs, involving periodic expansion and contraction of the gas perpendicular to the plane of the disc. These breathing or bouncing motions are expected to occur in all types of distorted accretion discs (warped, tidally distorted and eccentric discs) on account of the lack of vertical hydrostatic equilibrium in those objects. In theory, periodic nonlinear oscillations are exact solutions of the equations of ideal gas dynamics in a local model of the disc. To investigate the dynamics of these oscillations we carried out both quasi-2D and fully 3D local (shearing-box) hydrodynamic simulations with an isothermal equation of state. To complement our numerical simulations, we have also presented a 1D model of vertical oscillations and analysed their stability.

To isolate the mechanisms damping these vertical oscillations, we first considered the case of a freely oscillating disc (i.e. in the absence of any forcing). When the oscillation amplitude is sufficiently large, the disc collapses and expands supersonically, with pressure providing an impulsive restoring force. The resultant shocks dominate the damping of the vertical oscillations, at least in the very early stages of the simulation when the disc is undergoing only vertical oscillations, but their existence and properties are dependent on the artificial BCs at the top and bottom of the box.

A significant new result is that the vertical oscillations can excite a corrugation or bending mode by means of the parametric instability. This instability, which is an example of the parametric instability of inertial waves in distorted discs, further dampens the vertical oscillations as their kinetic energy is converted into that of the bending mode. The disc develops a radially periodic corrugation or warp that oscillates like a standing wave. The resultant radial pressure gradients in the corrugated disc lead to radially converging flows. After non-linear saturation the mid-plane of the corrugated disc consists of many intersecting waves which steepen into shocks in the atmosphere (⁠|$z \gtrsim 2H_0$|⁠). The dissipation due to these shocks ultimately dampens the vertical oscillations in the absence of forcing. However, on short time intervals (⁠|$\sim 2{\!-\!}5$| orbits), the oscillations can actually be re-invigorated by means of a three-mode coupling between the vertical oscillation and the oppositely directed travelling waves that form the warp or standing wave. These results are robust to changes in resolution, vertical box size, vertical BCs, and numerical and physical (explicit) viscosity.

We then investigated the case in which the oscillations were driven by means of periodic variation in the vertical component of gravity for various forcing frequencies |$\omega$|⁠, ranging from |$\omega =\Omega$| (relevant to an eccentric disc) to |$\omega = 2\Omega$| (relevant to a warped disc). Various values between these two limits correspond to the forcing frequency at different radial locations in a tidally distorted disc. Here, too, a disc that is initially oscillating only in the vertical direction becomes unstable, leading to the development of a bending mode or warp. However, for forcing frequencies greater than the resonant frequency (⁠|$\sqrt{2}\Omega$|⁠) we observe, in addition to the bending mode, vertical oscillations in the centre of mass of the disc. Viewed in the plane perpendicular to the disc in a local (shearing-box) frame this is manifested as uniform vertical oscillation of the entire disc, which becomes more extreme as the forcing frequency approaches twice the orbital frequency. These oscillations are related to the tilt instability expected to occur at the 3:1 resonance in tidally distorted binary discs (Lubow 1992).

ACKNOWLEDGEMENTS

This research was funded by STFC (Science and Technology Facilities Council) through grant ST/X001113/1. The authors would like to thank Hongping Deng, Henrik Latter, Callum Fairbairn, and Jim Stone for helpful discussions. Loren E. Held thanks Masaru Shibata for providing access to the Cobra and Raven clusters at the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany, and to the Yamazaki workstations at the Max Planck Institute for Gravitational Physics (Albert Einstein Institute) in Potsdam, Germany.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: MODEL OF COUPLED AND DAMPED OSCILLATORS

To try to explain the behaviour seen in several of the simulations, in which there are phases of re-invigoration of the bouncing motion, we consider the following model:

where |$g=1+a\cos \omega t$| is the strength of vertical gravity. This model is based on equations (29)–(31) in the case |$n=1$|⁠, with |$Y=\mathrm{i}X$|⁠, but also includes the conservative non-linear feedback of the daughter modes on the parent oscillation; the term involving |$Y^{*}Z$| in the first equation can be derived from a Lagrangian formulation of the three-mode system. It also includes non-linear damping terms (the final term in each equation) to model the dissipative effects of shocks. We have not made an extensive exploration of the model parameters |$\alpha _i$| and |$\beta _i$|⁠, but in Fig. A1 present an illustrative solution for |$H(t)$| in a set-up corresponding to the fiducial model of Section 4.1; this can be compared with Fig. 6 and shows phases in which the vertical oscillation alternately loses and gains energy from the daughter modes, combined with an overall damping. Note that when we set the damping parameters |$\alpha _i$| to zero, there is no dissipation anymore; instead we see quasi-periodic variation in |$H$| for small-amplitude oscillations. For larger amplitudes, it appears to be chaotic.

APPENDIX B: TABLES OF SIMULATIONS