Wavelike nature of the vertical shear instability in global protoplanetary discs

ABSTRACT

The vertical shear instability (VSI) is a robust phenomenon in irradiated protoplanetary discs (PPDs). The majority of previous numerical simulations have focused on the turbulent properties of its saturated state. However, the saturation of the VSI manifests as large-scale coherent radially travelling inertial waves. In this paper, we study inertial-wave-disc interactions and their impact on VSI saturation. Inertial-wave linear theory is developed and applied to a representative global 2D simulation using the athena++ code. It is found that the VSI saturates by separating the disc into several radial wave zones roughly demarcated by Lindblad resonances (turning points); this structure also manifests in a modest radial variation in the vertical turbulence strength. Future numerical work should employ large radial domains to accommodate this radial structure of the VSI, while concurrently adopting sufficiently fine resolutions to resolve the parametric instability that attacks the saturated VSI inertial waves.

hydrodynamics, instabilities, waves, protoplanetary discs

1 INTRODUCTION

The vertical shear instability (VSI) is a robust phenomenon in irradiated protoplanetray discs. Originally discovered in the context of differentially rotating stars, it is related to the Goldreich–Schubert–Fricke instability (Goldreich & Schubert 1967; Fricke 1968) and hence is centrifugal in nature, but with double diffusive aspects (Barker & Latter 2015; Latter & Papaloizou 2018). Its application to PPDs was first successfully demonstrated in Nelson, Gressel & Umurhan (2013). Subsequent analytical studies and numerical simulations have illustrated its linear behaviour (e.g. Lin & Youdin 2015; Cui & Lin 2021; Latter & Kunz 2022) and non-linear evolution (e.g. Nelson et al. 2013; Stoll & Kley 2014).

The saturation of the VSI in current numerical simulations is dominated by remarkably coherent ‘corrugation’ oscillations (Nelson et al. 2013). However, the majority of these simulations have treated the saturation as mere turbulence and focused on including more physics (e.g. Flock et al. 2017; Cui & Bai 2020; Schäfer, Johansen & Banerjee 2020) and exploring numerical convergence (e.g. Richard, Nelson & Umurhan 2016; Flores-Rivera, Flock & Nakatani 2020; Manger et al. 2020). On the other hand, linear theory reveals that (in vertically global models) the VSI is an overstability, a destabilized inertial wave (Barker & Latter 2015); indeed, the simulated corrugation modes correspond to n = 1, vertically standing and radially travelling inertial waves.

Upon radial propagation, VSI corrugation modes must respond to variations in the background disc structure. For instance, (1) as they travel, their radial wavenumbers will evolve (though not their wave frequencies), possibly reaching zero at special resonant radii (turning points); (2) the conditions for maximum VSI growth will vary with radius for a given wave, thus an unstable travelling wave might move to a region where it does not grow as fast (or at all). It follows that the disc may divide itself into several radial wave zones. At the start of a zone, a wave is excited with the frequency that maximizes its VSI growth rate, and the end of a zone may correspond to the turning point. As it approaches its turning point, a wave is supplanted by a faster growing wave, and a new wave zone starts. It follows that each zone may possess distinct turbulence strengths, leading to the large-scale radial structuring of the disc. To accommodate the above mentioned disc–wave interaction, a wide radial simulation domain is required. Very narrow domains likely impede the radial propagation of waves, accentuate boundary effects, and thus potentially misrepresent VSI saturation.

The paper is organized as follows. In Section 2, we introduce the linear theory for inertial waves. In Section 3, we conduct global numerical simulation and compare the theory to the simulation data. We summarize and discuss the main findings in Section 4.

2 LINEAR WAVE THEORY

In this section, we introduce the linear theory of inertial waves so as to facilitate the simulation data analysis in Section 3. The inertial waves’ dispersion relation and eigenfunctions in a vertical stratified shearing box model are derived in Section 2.1. Their propagation over, and interactions with, a global disc are described in Section 2.2. The linear theory of the VSI from a purely local model is recapped in Section 2.3 and applied to the travelling waves.

2.1 Dispersion relation and eigenfunctions

Working in a vertically stratified, radially localized Keplerian shearing sheet model, co-orbiting with the gas with frequency Ω, the Cartesian coordinates x, y, z represent the radial, azimuthal, and vertical directions, respectively. The governing equations are (Latter & Papaloizou 2017)

$$\begin{eqnarray*} \frac{\partial \rho }{\partial t}+\mathbf {u} \cdot \nabla \rho =-\rho \nabla \cdot \mathbf {u}, \end{eqnarray*}$$

(1)

$$\begin{eqnarray*} \frac{\partial \mathbf {u}}{\partial t}+\mathbf {u} \cdot \nabla \mathbf {u}=-\frac{1}{\rho } \nabla P-2 \Omega \boldsymbol{e}_{z} \times \mathbf {u}-\nabla \Phi , \end{eqnarray*}$$

(2)

where |$\Phi = -\frac{3}{2}\Omega ^2 x^2+\frac{1}{2}\Omega ^2 z^2$| is the tidal potential. The gas has an isothermal equation of state, |$P=c_s^{2} \rho$|⁠, where c_s is a constant sound speed.

The above equations admit the equilibrium solution

$$\begin{eqnarray*} \mathbf {u}_{0} =-\frac{3}{2} \Omega x \boldsymbol{e}_{y}, \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} \rho _0 = \rho _\text{mid}\,\, \text{exp}[-z^2/(2H^2)], \end{eqnarray*}$$

(4)

where ρ_mid denotes the constant mid-plane density, and H = c_s/Ω is the pressure scale height. This steady state is perturbed with small axisymmetric perturbations |$\tilde{\rho }$|⁠, |$\tilde{\mathbf {u}}$| of the form |$\tilde{\rho }= \rho ^{\prime }(z) \mathrm{e}^{\mathrm{i} k_x x-\mathrm{i} \omega t}$|⁠, etc. We assume, without loss of generality, that ω > 0; thus positive k_x corresponds to radially outward travelling waves, and negative k_x to inward travelling waves. The linearized perturbation equations are now

$$\begin{eqnarray*} \mathrm{i} \omega u_x^{\prime } = -2\Omega u_y^{\prime }+\mathrm{i} k_xh^{\prime }, \end{eqnarray*}$$

(5)

$$\begin{eqnarray*} \mathrm{i} \omega u_y^{\prime } =\frac{1}{2} \Omega u_x^{\prime }, \end{eqnarray*}$$

(6)

$$\begin{eqnarray*} \mathrm{i} \omega u_z^{\prime } =\frac{\mathrm{d} h^{\prime }}{\mathrm{d} z}, \end{eqnarray*}$$

(7)

$$\begin{eqnarray*} \mathrm{i} \omega \frac{h^{\prime }}{c_s^{2}} =\frac{\mathrm{d} u_z^{\prime }}{\mathrm{d} z}+\frac{\mathrm{d} \ln \rho _{0}}{\mathrm{d} z} u_z^{\prime }+ \mathrm{i} k_x u_x^{\prime }, \end{eqnarray*}$$

(8)

where h′ = P′/ρ₀, and |$P^{\prime }=c_s^2\rho ^{\prime }$|⁠. Eliminating the velocity variables, we arrive upon a form of Hermite’s equation

$$\begin{eqnarray*} \frac{\mathrm{d}^{2} h^{\prime }}{\mathrm{d} z^{2}}-\frac{z}{H^{2}} \frac{\mathrm{d}h^{\prime }}{\mathrm{d} z}+\frac{\xi }{H^{2}} h^{\prime }=0, \end{eqnarray*}$$

(9)

where

$$\begin{eqnarray*} \xi =\frac{\omega ^{2}}{\Omega ^{2}}-H^{2} k_x^{2}\left(\frac{\omega ^{2}}{\omega ^{2}-\Omega ^{2}}\right). \end{eqnarray*}$$

(10)

Solutions that have vanishing momentum as |z| → ∞ require ξ to be an integer n, and are thus Hermite polynomials of degree n (Lubow & Pringle 1993). The condition ξ = n gives the dispersion relation

$$\begin{eqnarray*} \omega ^{4}-\left[\Omega ^{2}(n+1)+c_s^{2} k^{2}\right] \omega ^{2}+\Omega ^{4} n=0. \end{eqnarray*}$$

(11)

It describes three groups of waves: the low-frequency inertial waves (ω < Ω), the high-frequency acoustic waves (ω > Ω), and density waves, for which n = 0.

The waves that dominate VSI simulations are the low frequency n = 1 inertial waves, i.e. the corrugation modes identified by Nelson et al. (2013). Corrugation modes possess |$u_x^{\prime }\propto z$| and |$u_y^{\prime }\propto z$|⁠, but |$u_z^{\prime }$| does not depend on z, which leads to the radial sequence of upward and downward motion familiar from simulations (all velocity components vary sinusoidally with x). The meridional velocity vector |$(u_x^{\prime }, u_z^{\prime })$| of a corrugation mode is plotted in the meridional plane in the top panel of Fig. 1. Notice its characteristic vertical ‘roll’ structure, which crosses the mid-plane. These roll features have indeed been observed in the meridional plane of recent VSI simulations (see bottom panel of fig. 7 in Flores-Rivera et al. 2020). In addition, we plot the higher order n = 6 inertial wave in the lower panel of Fig. 1, which exhibits a vertical pattern of smaller rolls. These structures have only been witnessed in the highest resolution simulations of Flores-Rivera et al. (2020) (top panel of their fig. 1), and in that instance have probably been excited by parametric instability attacking the primary VSI corrugation mode, as detailed in Cui & Latter (2022).

$Arrows: meridional velocity vectors $(u_x^{\prime }, u_z^{\prime })$ of n = 1 (top) and n = 6 (bottom) inertial wave modes. Background colours: equilibrium density profile ρ0/ρmid by equation (4).$

Figure 1.

Arrows: meridional velocity vectors |$(u_x^{\prime }, u_z^{\prime })$| of n = 1 (top) and n = 6 (bottom) inertial wave modes. Background colours: equilibrium density profile ρ₀/ρ_mid by equation (4).

Open in new tab Download slide

2.2 Global wave propagation

The waves described above are radially local, in that their radial wavelength is assumed to be much shorter than the length-scale of any background variation in the disc. They may, however, travel over significant radial distances and thus sample the large-scale structure of the disc. As a result, their properties will evolve over radius, in particular their radial wavelength, though we expect their wave frequencies to stay fixed on the time-scale of propagation (Whitham 1974; Lighthill 1978). We now describe this evolution.

2.2.1 Resonances or ‘turning points’

As waves travel large distances, they may encounter special radii where they interact strongly with the background flow, i.e. resonances (Kato 2001). The most important radii for our axisymmetric waves are where the radial wavelength approaches infinity (k_x = 0), called turning points, which also correspond to Lindblad resonances. These turning points split the disc into radial regions where a general inertial wave can propagate and where it is prohibited. Waves that are incident on a turning point usually reflect. Rearranging equation (11), we obtain

$$\begin{eqnarray*} c_s^2 k_x^2 = \frac{(\Omega ^2-\omega ^2)(n\Omega ^2-\omega ^2)}{\omega ^2}. \end{eqnarray*}$$

(12)

Turning points can be obtained by setting the left side equal to zero. Corrugation modes of n= 1 are unique for several reasons. In a Keplerian disc, their turning points occur when ω = Ω and correspond to a hybrid vertical/Lindblad resonance. Moreover, and quite unusually, corrugation modes can propagate on either side of the turning point, and potentially can even travel through the point, despite the wave possessing an infinite wavelength there (Bate et al. 2002). If we assume a Keplerian disc and set Ω = Ω ₀(R/ R₀) ^−3/2 (for fixed Ω ₀ and R₀), we find that the radial location R_tp of the turning points are

$$\begin{eqnarray*} R_{\rm tp} = R_0 (\Omega _0/\omega)^{2/3}. \end{eqnarray*}$$

(13)

2.2.2 Radial wavenumber

Furthermore, as a wave packet of given frequency ω travels radially from its point of excitation, its radial wavenumber will change on account of the background disc structure. This change is described by equation (12), which gives us k_x as a function of R (through c_s and Ω).

Far from the turning point, so that ω ≪ Ω, a low-frequency n = 1 inertial wave’s wavenumber will vary as |k_x| ≈ Ω²/(c_sω). If, in addition, the disc possesses a background temperature structure so that |$T\propto R^{-q_T}$|⁠, the radial wavenumber scales as |$k_x\sim R^{q_T/2-3}$|⁠. Early VSI simulations by Stoll & Kley (2014) used q_T = 1 and they indeed reported a radial wavenumber that varied as R^−5/2. Note that the more accurate equation (12) should be used more generally because it properly accounts for behaviour near the resonance.

2.2.3 Phase and group velocities

As discussed by Lubow & Pringle (1993), wave fronts travel radially, with the disc acting as a wave-guide. The phase speed is defined as c_p = ω/k_x, which is approximately ±(ω/Ω)²c_s for n = 1 inertial waves far from the turning point (ω ≪ Ω), with the ± indicating inward or outward propagation. The radial location of any given wave crest R_p(t) can then be determined from dR_p/dt = c_p. To facilitate comparison with our later simulations (Section 3), we set q_T = 1 for the rest of this section, and thus c_s = c₀(R/R₀)^−1/2; taking the wave crests to move inwards, we find

$$\begin{eqnarray*} R_\text{p}\approx R_i\left[1+\alpha (t-t_i)\right]^{-2/3}, \end{eqnarray*}$$

(14)

where α = 3hω²(R_i/R₀)^3/2/(2Ω₀), h = c₀/(ΩR₀), and the crest begins at R = R_i at t = t_i. Thus the wavecrests slow down as they move inwards.

In a vertically unstratified local model, the group velocity of a freely travelling inertial wave is perpendicular to the phase velocity. But in a vertically stratified shearing box the group and phase velocities are oppositely directed (Lubow & Pringle 1993). Taking the derivative of equation (12) with respect to k_x, provides the following expression for the group velocity of n = 1 inertial waves,

$$\begin{eqnarray*} c_\text{g} = \frac{\mathrm{d}\omega }{\mathrm{d}k_x}=\mp \frac{c_s}{1+\Omega ^2/\omega ^2} , \end{eqnarray*}$$

(15)

where the minus (plus) sign is taken if the phase speed is positive (negative). Far from the resonance we have ω ≪ Ω, and so c_g ≈ ∓(ω/Ω)²c_s ≈ −c_p. The group and phase velocities share the same magnitude but have opposite sign. A packet of inertial waves travels at the group speed, but wavetrains that extend over many individual wavelengths are often exposed to modulations, jumps, or other inhomogeneities in their phases and amplitudes; these ‘defects’ also travel with speed c_g. If a wave defect is sufficiently localized in radius, its trajectory can be calculated by solving dR_d/dt = c_g, where R_d(t) is the radial location of the defect. Far from the resonance we can use the approximation above to obtain

$$\begin{eqnarray*} R_\text{d}(t) \approx R_i [1-\alpha (t-t_i)]^{-2/3}, \end{eqnarray*}$$

(16)

assuming the defect begins at R = R_i at t = t_i and c_g points outwards.

2.3 Linear growth rate of the VSI

So far we have dealt with free inertial waves, and not with the circumstances of their excitation. As discussed in early studies, the VSI provides a mechanism to drive the growth of inertial waves with certain wavevectors and frequencies (Nelson et al. 2013; Barker & Latter 2015). Of special interest is how the VSI chooses the frequencies of the inertial waves that ultimately dominate the disc.

At a given radius we might expect that the fastest growing VSI mode will control the dynamics. To characterize this mode, we adopt the purely local incompressible model of Latter & Papaloizou (2018).¹ We define the dimensionless vertical shear rate as q_z = −R∂ln Ω/∂z, which is of order h from the thermal wind equation (Knobloch & Spruit 1986; Latter & Kunz 2022). Denoting the vertical wavenumber of the mode by k_z, the maximum VSI growth occurs when

$$\begin{eqnarray*} \frac{k_x}{k_z} = -\frac{1}{q_z}, \end{eqnarray*}$$

(17)

and growth ceases altogether for sufficiently small | k _x / k _z |. To obtain the associated wave frequency of fastest growth, we substitute k_x into equation (12) and solve the quadratic equation for ω. The difficulty is to relate n to k_z. A corrugation mode has n = 1 and thus its vertical ‘wavelength’ will be of order the disc thickness, i.e. some multiple of H. Anticipating our global simulation (Section 3), we set k_z ∼ 1/(5H) and q_z ∼ h = 0.05, which gives us

$$\begin{eqnarray*} \omega _{\rm VSI} = (\sqrt{5}-2)\Omega . \end{eqnarray*}$$

(18)

Because ω_VSI scales as Ω, the VSI favours lower frequency waves further out in the disc. Note also that ω_VSI is always less than Ω, where the resonance occurs, and thus it always lies in the inertial wave zone.

The above calculation gives a rough indication of which inertial wave frequency might be selected at any given radius. However, in reality, we do not expect to observe the ω of the waves to smoothly vary with radius. Indeed, in our simulation (Fig. 4) we see distinct wave zones emerge, each extending over a range of radii. We expect a zone that starts at a given radius to possess an ω that yields maximized VSI growth, via equation (18). If the wave travels outwards towards its turning point, its radial wavenumber will decrease towards 0 and its driving by the VSI will end; it may then be overtaken by a different growing VSI mode, and a new wave zone will begin.

Left: snapshot of vertical velocity normalised by local sound speed vz/cs at t = 500P0, displayed in logarithmic scale. Right: meridional velocity vector (vx, vz) for R ∈ [1.7, 2].

Figure 2.

Left: snapshot of vertical velocity normalised by local sound speed v_z/c_s at t = 500P₀, displayed in logarithmic scale. Right: meridional velocity vector (v_x, v_z) for R ∈ [1.7, 2].

Open in new tab Download slide

Figure 3.

Space–time diagram of corrugation modes. Colours: strength of vertical velocity at the mid-plane. Yellow dashed curve: fitting to phase velocity with R_i = 1.65, t_i = 2100 by equation (14). White dashed curve: fitting to group velocity with R_i = 3.97, t_i = 2100 by equation (16). Notice that white curve is slightly shifted upwards in order to display wave defect pattern from colours. Vertical dotted lines: locations of turning points (Table 1).

Open in new tab Download slide

Figure 4.

Unnormalized Fourier amplitudes of the corrugation modes at different radii. Dotted line denotes frequency for maximum VSI growth (equation 18). Solid line denotes the Keplerian frequency (where resonance occurs, ω = Ω).

Open in new tab Download slide

3 NUMERICAL SIMULATIONS

Global hydrodynamic simulations of the VSI are performed to verify and extend the linear inertial wave theory developed in the preceding section. We now detail our numerical methods (Section 3.1), disc model (Section 3.2), simulation setup (Section 3.3), and simulation results (Section 3.4).

3.1 Dynamical equations

We employ the grid-based high-order Godunov code Athena++ (Stone et al. 2020) to conduct global 2D hydrodynamic simulations of the VSI. The continuity, momentum, and energy equations in their conservative form are

$$\begin{eqnarray*} \frac{\partial\rho} {\partial t} + {\nabla} \cdot {(\rho {\rm v})} = 0, \end{eqnarray*}$$

(19)

$$\begin{eqnarray*} \frac{\partial(\rho {\rm v})} {\partial t}+\nabla \cdot { (\rho {\rm v} {\rm v}+P {\bf I})} = - \rho \nabla {\Phi }, \end{eqnarray*}$$

(20)

$$\begin{eqnarray*} \frac{\partial E} {\partial t}+\nabla \cdot {(E+P) {\rm v}} =-\rho ( {\rm v}\cdot \nabla \Phi) -\Lambda _\mathrm{c}. \end{eqnarray*}$$

(21)

Here, |${\bf v}$|⁠, ρ, and P denote gas velocity, density, and pressure, and |${\bf I}$| is the identity tensor. The total energy density is E = ϵ + ρv²/2, where ϵ denotes the internal energy density. This is related to the gas pressure by the ideal gas equation of state P = (γ − 1)ϵ. An adiabatic index of γ = 7/5 is adopted to match the molecular gas in a protoplanetary disc. The gravitational potential of the central star is given by Φ = −GM/r with stellar mass M, and Λ_c is a cooling term. The simulation is conducted in spherical polar coordinates (r, θ, ϕ), though cylindrical coordinates (R, z, ϕ) are occasionally used to improve presentation.

3.2 Disk model

An equilibrium disc model with constant vertical temperature is adopted, i.e. the gas is vertically isothermal. The temperature and mid-plane density are prescribed by radial power laws

$$\begin{eqnarray*} \rho _{\mathrm{mid}}=\rho _{0}\left(R/R_{0}\right)^{-q_\mathrm{D}}, \end{eqnarray*}$$

(22)

$$\begin{eqnarray*} c_s^2=c_0^2\left(R/R_{0}\right)^{-q_\mathrm{T}}, \end{eqnarray*}$$

(23)

where ρ₀ and c₀ denote density and sound speed at the inner radial boundary R₀. Following Nelson et al. (2013), the equilibrium solutions for density and angular velocity are

$$\begin{eqnarray*} \rho (R,z) = \rho _0 (R/R_0)^{-q_\mathrm{D}}\exp \left[\frac{GM}{c_s^2}\bigg (\frac{1}{r}- \frac{1}{R}\bigg)\right], \end{eqnarray*}$$

(24)

$$\begin{eqnarray*} \Omega (R,z) = \Omega _\text{K}[-(q_\mathrm{D}+q_\mathrm{T})h^2 + (1-q_\mathrm{T})+q_\mathrm{T}R/r]^{1/2}. \end{eqnarray*}$$

(25)

We choose q_T = 1, and since |$H \propto R^{(-q_\mathrm{T}+3)/2}$|⁠, the aspect ratio h = H/R is a constant, which we set to 0.05. Thus the vertical shear is

$$\begin{eqnarray*} \frac{\partial\Omega} {\partial z} \approx -\frac{1}{2}\Omega _\text{K}\frac{z}{R^2}, \end{eqnarray*}$$

(26)

and hence q_z ∼ h. Lastly, the power-law index for density is q_D = 2.

Rapid Newtonian cooling is employed to help excite the VSI (Nelson et al. 2013). In the code, we adjust the temperature at each time-step Δt by

$$\begin{eqnarray*} \Delta T =(T_0-T)[1-\exp (-\Delta t/\tau)]. \end{eqnarray*}$$

(27)

The thermal relaxation time-scale τ is taken to be extremely small, so as to enforce locally isothermality; it is 10⁻¹⁰ innermost orbital periods.

3.3 Simulation setup

The radial domain of the simulation spans r ∈ [1, 30] and we use logarithmic gridding, in which the spacing increases by a constant factor of 1.005 per radial grid cell. The uniform spacing of the θ domain ranges as θ ∈ [π/2 − 5h, π/2 + 5h]. Our 2D simulation is equipped with 704 × 448 cells in r × θ. This allows us to reach a spatial resolution of 10 cells per H in r and 12 cells per H in θ.

In code units, we set GM = R₀ = 1. The simulation is run up to 3000P₀, where P₀ = 2π/Ω₀ is the orbital period at the inner boundary. To initialize the VSI, background noise is introduced to the velocities and is set to 1 per cent of the local sound speed. Variables are set strictly to the equilibrium solution at the boundaries.

3.4 Simulation results

Here, we present our numerical simulation results and compare with the linear theory developed in Section 2. Analyses are performed within a time interval from 2000 to 3000P₀ and at the disc mid-plane.

We first give an overview of the simulation results, with a focus on wave zones, wave dynamics, and numerical boundary effects. The left-hand panel of Fig. 2 shows the vertical velocity normalized by the local sound speed in the R–z plane. The large-scale vertical oscillations witnessed here (as in all extant published VSI simulations) correspond to the inertial waves that our linear theory describes in Section 2. These inertial waves possess a vertical node number n = 1 (Section 2.1) because of the near absence of vertical structure and are thus corrugation modes (e.g. Nelson et al. 2013). The right-hand panel of Fig. 2 shows the meridional velocity vector (v_x, v_z) but in a narrower radial domain. These circulations, obtained from the simulation, though slightly bent and possessing a steeper radial variation, have the same morphology as in the top panel of Fig. 1, which shows the same quantity predicted by linear theory.

Fig. 3 shows the space (radial) time diagram of the VSI vertical velocities. First, it is clearly seen that our simulation domain is divided into four radial wave zones, roughly demarcated by the vertical dotted lines (see below). Next, the most discernible pattern in Fig. 3 are the coloured bands: red and blue colours correspond to the crests and troughs of the radially travelling inertial waves (corrugation modes) depicted in Fig. 2. They travel inwards and slow down as they move towards smaller radii. Narrow grey strips, especially prominent in the innermost wave zone, correspond to wave defects that propagate outwards. Furthermore, in between the wave zones, there appears an interference pattern as different waves interact. The boundaries between zones appear relatively stable over time, and thus permit a quasi-steady large-scale structure to be sustained as part of the VSI saturation. Lastly, the wave propagation is cut off by the outer radial boundary, possibly resulting in numerical artefacts in the outermost zone; our large radial domain, however, ensures the rest of the wave zone dynamics is unaffected and is physically realistic.

Fig. 4 shows the Fourier amplitudes of the VSI corrugation modes at different radii, by performing temporal Fourier transforms of the mid-plane vertical velocities. It can be immediately seen that four distinct bands of frequencies appear, which correspond to the four wave zones in Fig. 3. We note that the frequency zones partly overlap, though Fig. 3 indicates that normally only one wave dominates throughout its wave zone (also see Fig. 5 later), presumably because it possesses the strongest Fourier amplitude. We have gone through by eye, selected one frequency associated with each band, and tabulated them in Table 1.

The radial wavelengths λ of corrugation modes versus radius over a time interval from 2000 to 3000P0. Colours, in logarithmic scale, are the probability that a specific wavelength can occur at a fixed radius. Overlaid curves are theoretical predictions by equation (12). Verical dashes lines are locations of turning points (Table 1).

Figure 5.

The radial wavelengths λ of corrugation modes versus radius over a time interval from 2000 to 3000P₀. Colours, in logarithmic scale, are the probability that a specific wavelength can occur at a fixed radius. Overlaid curves are theoretical predictions by equation (12). Verical dashes lines are locations of turning points (Table 1).

Open in new tab Download slide

Table 1.

Open in new tab

A list of frequencies and turning points.

	1	2	3	4
ω/Ω₀	0.075	0.025	0.0094	0.0031
R_tp	5.6	11.7	22.4	47.6

Table 1.

Open in new tab

A list of frequencies and turning points.

	1	2	3	4
ω/Ω₀	0.075	0.025	0.0094	0.0031
R_tp	5.6	11.7	22.4	47.6

Now, we compare the numerical results with the analytical theory outlined in Section 2.2.3 and explain features seen in Fig. 3. Wave crests and troughs are expected to travel at the phase speed, and hence their trajectory is delineated by equation (14). Indeed, the yellow dashed curve in Fig. 3, denoting the theoretical prediction for a selected crest, shows good agreement with the simulation data. In the theory, the phase velocity can either point radially inwards or outwards, with the corresponding group velocity possessing the same amplitude but having opposite sign (equation 15). Nevertheless, our numerical simulations, as well as Stoll & Kley (2014) and Pfeil & Klahr (2021), always find the group velocity pointing outwards. One explanation is that VSI waves generated in inner regions emerge first, where the time-scales are faster, and then colonize the outer regions before the VSI there has time to saturate. Finally, wave defects are expected to travel at the inertial wave group velocity. Similarly, we anticipate their trajectories to be described by equation (16) and so overlay the theoretical prediction (white dashed curve), which perfectly agrees with the simulation data.

In Fig. 4, the solid line denotes the Keplerian frequency. The intersections between the Keplerian frequency and the yellow bands yield the locations of turning points (equation 12). These turning points are located roughly at where each wave zone terminates, consistent with Fig. 3. The dotted line corresponds to the frequency for the maximum VSI growth (equation 18). The radial domain that each wave zone spans seems to be confined by these two characteristic frequencies, as predicted in Section 2.3.

Fig. 5 shows the radial wavelengths λ of the corrugation modes versus radius over a time interval between 2000 and 3000P₀. The wavelengths are estimated by setting the distance between two consecutive sign changes in the mid-plane vertical velocity to equal one full wavelength (Stoll & Kley 2014). The colours, shown in logarithmic scale, represent the probability that a specific wavelength λ occurs at a radius R. More specifically, at each R, 1000 wavelengths are obtained over the selected time interval. We bin the wavelengths, and at each R, calculate the probabilities by taking the ratio of data points collected in each λ bin to the total data points collected (=1000). It can be clearly seen that four wave zones are present. Each of them terminates near the locations of their turning points (vertical dashed lines). However, before reaching exactly the turning point, each wavetrain is supplanted by another wavetrain with frequency closer to that giving fastest VSI growth. As is quite striking, within each wave zone the wavelengths increase as the radius increase. This was indeed anticipated in Section 2.2.2, with the result k_x ∼ R^−5/2 away from the turning points. We use the more accurate equation (12) to obtain the wavelength λ(R), with frequencies gained from Fig. 4. These predictions are overlaid in Fig. 5, which match perfectly with the coloured bands.

So far, we have seen that the disc is divided into four wave zones in the simulation domain, which appear relatively stable over time (Fig. 3). Whether mean disc properties vary from wave zone to wave zone is thereby of interest. Fig. 6 shows the root mean square vertical velocity δv_z normalized by the local sound speed. It is clearly seen that the mean vertical turbulence level increases with radius. Superimposed on this mean increase is a sawtooth pattern, separated into four zones, roughly coincident with the radial turning points. The level of variation in δv_z/c_s is about a factor of two between the highest and lowest values in our domain. Therefore, a modest correlation between the saturated vertical level of turbulence and wave zones is found. However, other disc properties are found to be very weakly affected, consistent with previous numerical simulations with extended radial domains (e.g. Flock et al. 2017, 2020; Cui & Bai 2020).

Figure 6.

Level of vertical turbulence normalized by local sound speed δv_z/c_s as a function of radius at the mid-plane. The equation used to calculate this quantity can be found by equation (31) in Cui & Bai (2021).

Open in new tab Download slide

4 CONCLUSIONS AND DISCUSSION

In this paper, we revisited the global VSI hydrodynamic simulation and studied its wave-like properties upon saturation. As witnessed in previous simulations, the non-linear saturation of the VSI is dominated by radially travelling corrugation modes (n = 1 inertial waves). In the first part of this paper, a linear theory for these waves is developed in a vertically global but radially local model, and their response to the background disc structure upon radial propagation over a global disc is determined. In the second part of this paper, a representative 2D hydrodynamic simulation was conducted with the athena++ code and compared with the linear theory.

The comparison between the theory and the simulation results is satisfactory. It is found that the simulation domain is separated into four radial wave zones. Each zone roughly starts at locations that maximize the linear VSI growth rate and ends near the locations of the turning points (resonances), where the modes' wavelengths go to infinity. In each wave zone, inertial wave crests travel radially inwards at the phase velocity. The wave defects travel radially outwards at the group velocity. Wave zones possess vertical turbulence levels that are modestly distinct from each other.

In future work, it is necessary for numerical simulations to employ large radial domain to accommodate the VSI wave propagation, in particular to describe distinct wave zones. This paper studies the large-scale coherent corrugation modes that dominate in low resolution global simulations. This coherence could, however, be disrupted by parametric instability, though this requires finer resolutions to resolve (Cui & Latter 2022). It should be a high priority to verify that the saturation properties described here are qualitatively correct when the parametric instability can emerge.

ACKNOWLEDGEMENTS

We thank the anonymous referee, whose comments helped improve and clarify this manuscript, and Gordon Ogilvie who pointed out some errors in a previous version of the manuscript. We thank Lizxandra Flores-Rivera for useful discussions. ES acknowledges support from the Philippa Fawcett Internship Programme at the Centre for Mathematical Sciences, University of Cambridge. CC and HNL acknowledge funding from STFC grant ST/T00049X/1. Numerical simulations are conducted on the FAWCETT cluster at the Department of Applied Mathematics and Theoretical Physics, University of Cambridge.

DATA AVAILABILITY

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

Footnotes

Ideally, we would use a vertically stratified model to determine the k_x of maximum VSI growth for corrugation modes of n = 1. Unfortunately, this is not possible as from the outset these models assume that k_x is very large (Nelson et al. 2013; Barker & Latter 2015) and prove erroneous when k_x takes smaller values.

REFERENCES

Bate

Ogilvie

Lubow

Pringle

2002

MNRAS

332

575

Month:	Total Views:
June 2022	138
July 2022	508
August 2022	220
September 2022	269
October 2022	270
November 2022	192
December 2022	131
January 2023	580
February 2023	222
March 2023	199
April 2023	152
May 2023	267
June 2023	110
July 2023	122
August 2023	125
September 2023	475
October 2023	114
November 2023	137
December 2023	344
January 2024	142
February 2024	320
March 2024	245
April 2024	297
May 2024	498
June 2024	214
July 2024	280
August 2024	170
September 2024	101
October 2024	263
November 2024	257
December 2024	247
January 2025	163
February 2025	164
March 2025	175
April 2025	217

Article Contents

Wavelike nature of the vertical shear instability in global protoplanetary discs

ABSTRACT

1 INTRODUCTION

2 LINEAR WAVE THEORY

2.1 Dispersion relation and eigenfunctions

2.2 Global wave propagation

2.2.1 Resonances or ‘turning points’

2.2.2 Radial wavenumber

2.2.3 Phase and group velocities

2.3 Linear growth rate of the VSI

3 NUMERICAL SIMULATIONS

3.1 Dynamical equations

3.2 Disk model

3.3 Simulation setup

3.4 Simulation results

4 CONCLUSIONS AND DISCUSSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only