Free Access
Issue
A&A
Volume 650, June 2021
Article Number A49
Number of page(s) 20
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202040094
Published online 07 June 2021

© ESO 2021

1. Introduction

When sufficiently massive, the dynamics of astrophysical disks is controlled by their own gravity in addition to the object they orbit. A number of observed class 0/I protostellar disks could fit in this category and thus exhibit self-gravitating phenomena (e.g., Eisner et al. 2005; ALMA Partnership 2015; Mann et al. 2015; Liu et al. 2016; Forgan et al. 2016). These phenomena can be global as they result from the long-range interaction of distant disk parcels. Even so, they may be described in terms of the local properties near a given disk parcel, such as the stability of the disk to gravitational collapse (Safronov 1960; Toomre 1964; Goldreich & Lynden-Bell 1965a). In this paper, we consider gravitationally unstable gaseous disks whose mass is a significant fraction of that of the central object. These include the outer regions of circumstellar disks or active galactic nuclei (Shlosman & Begelman 1987; Durisen et al. 2007; Kratter & Lodato 2016), but this excludes weakly collisional systems such as galactic disks or planetary rings (e.g., Salo 1992; Richardson 1994).

Gravitational instability (GI) can trigger the fragmentation of the disk into stellar and substellar companions (Kolykhalov & Syunyaev 1980; Podolak et al. 1993; Boss 1998, 2000), provided that radiative cooling is fast enough (Gammie 2001; Rice et al. 2003; Lin & Kratter 2016; Booth & Clarke 2019). We consider cooling times longer than orbital times for which the instability sustains gravito-turbulence instead. The gas orbital energy is converted into heat by the dissipation of supersonic turbulence, balancing radiative cooling while allowing matter to accrete onto the central object (Paczynski 1978; Lin & Pringle 1987). Gravito-turbulence appears to organize itself into large-scale spiral arms amidst small-scale fluctuations (Laughlin & Bodenheimer 1994; Riols et al. 2017). These spiral arms are believed to induce variability in the mass accretion rate (Lodato & Rice 2005; Vorobyov & Basu 2005, 2015), and they could play important roles in planet formation via dust grain accumulation (Rice et al. 2004; Clarke & Lodato 2009; Dipierro et al. 2015), thermal processing, and collisional fragmentation (Boss & Durisen 2005; Boley & Durisen 2008; Ida et al. 2008; Booth & Clarke 2016). The morphology of observed spirals may also be used to identify self-gravitating disks and to constrain their mass and temperature distributions (Dong et al. 2015, 2016, 2018; Hall et al. 2016; Meru et al. 2017; Forgan et al. 2018; Cadman et al. 2020).

The “shearing box” approximation (Hill 1878; Goldreich & Lynden-Bell 1965b) has often been used to study the saturation of gravitational disk instabilities at high spatial resolution in a controlled environment (Gammie 2001; Shi & Chiang 2014; Riols et al. 2017; Hirose & Shi 2017; Baehr et al. 2017; Booth & Clarke 2019). However, massive disks may also support global wave modes which the shearing box cannot capture (Adams et al. 1989; Papaloizou & Savonije 1991). Although tightly wound waves of short wavelengths can fit in a shearing box, the assumed periodicity forces them to interact with the same fluid elements without ever leaving the box. Should a global wave train travel through the disk, its coherence would be tied to keeping a well-defined angular pattern speed (Lin & Shu 1964). In this case, the wave could extract orbital energy at one location and transport it radially before thermalization, in contrast to the local energy conversion of viscous α disk models (Shakura & Sunyaev 1973; Balbus & Papaloizou 1999) and of local boxes. As a consequence, the angular momentum flux and the resulting mass accretion rate would no longer be determined by a condition of local thermal balance (Gammie 2001).

Several authors have addressed the issue of nonlocal energy transport in global simulations of gravito-turbulent disks. Using cylindrical grid-based simulations and focusing on low-order spirals, Mejía et al. (2005), Boley et al. (2006), and Steiman-Cameron et al. (2013) identified traveling waves by their distinct pattern speed. However, these waves could only travel over relatively narrow radial intervals. In addition, their nonlinear coupling with smaller-scale modes may have led to local energy dissipation and interfered with wave-induced transport if present (Michael et al. 2012). Using smoothed particle hydrodynamics (SPH), Lodato & Rice (2004, 2005) measured a stress compatible with local energy conversion in disks as massive as the central object. Later, Cossins et al. (2009) and Forgan et al. (2011) argued for the preponderance of nonlocal energy transport for disk mass approaching the central object’s, which is in apparent conflict with Lodato & Rice (2005). However, their simulation diagnostics relied on the linear dispersion relation for tightly wound two-dimensional (2D) isothermal density waves, with a simple correction for the disk stratification. The use of such a dispersion relation is questionable given the nonlinear nature of these waves (Shu et al. 1985; Laughlin et al. 1997, 1998) and the nonisothermal structure of the disk (Lubow & Ogilvie 1998).

In this paper, we study the spiral structures emerging in global three-dimensional (3D) simulations of self-gravitating disks. In particular, we examine their morphology, coherence, and propagation through the disk while paying special attention to their deviations from locality. Unlike most previously published results, we used a finite-volume code in spherical geometry to capture discontinuities and to resolve the vertical motions of the flow from the disk midplane to its low-density corona. In Sect. 2 we present our disk model, provide some theoretical expectations, and explain the numerical method designed to confront them. We chose a simulation as a reference to present in Sect. 3 before considering different disk masses and cooling times in Sect. 4. We summarize our findings in Sect. 5 and enumerate several issues left open by this work.

2. Model, theory, and method

Our goal is to characterize steady gravito-turbulent disks within a global framework. In this section we describe our physical modeling of the problem, some background theory including our a priori expectations, and the numerical method we used.

2.1. Disk model

We consider a rotationally supported gaseous disk orbiting a central mass m. We assume that the disk is geometrically thin and that its mass mdisk is sufficiently large for its own gravity to be dynamically important. Finally, we assume that the disk radiates heat away at a prescribed rate, thereby evolving toward a cold and gravitationally unstable state.

We place the central object at the origin of the spherical coordinate system (r,θ,φ), with θ = 0 along the rotation axis of the disk. Let also (R,φ,z) denote cylindrical coordinates, with z = rcos(θ) measuring the height relative to the disk midplane. We isolate the disk by considering only a finite radial interval r ∈ [rin,rout], excluding both the central mass and the environment of the disk on larger scales.

2.2. Governing equations

The gas density ρ, velocity v, and pressure P evolve according to

t ρ = · ( ρ v ) , $$ \begin{aligned}&\partial _t \rho = - \nabla \cdot \left(\rho \boldsymbol{v}\right), \end{aligned} $$(1)

ρ t v = ρ v · v P ρ Φ , $$ \begin{aligned}&\rho \partial _t \boldsymbol{v} = -\rho \boldsymbol{v}\cdot \nabla \boldsymbol{v} - \nabla P - \rho \nabla \Phi , \end{aligned} $$(2)

t P = v · P γ P · v ( γ 1 ) q , $$ \begin{aligned}&\partial _t P = -\boldsymbol{v}\cdot \nabla P - \gamma P \nabla \cdot \boldsymbol{v} - \left(\gamma -1\right) q^-, \end{aligned} $$(3)

assuming an ideal equation of state P = (γ−1)ρϵ with ϵ the specific internal energy and γ = 5/3 a constant adiabatic exponent. The gravitational potential Φ is the sum of two contributions: Φ = −Gm/r from the central mass and Φdisk from the disk itself, which satisfies Poisson’s equation

Δ Φ disk = 4 π G ρ , $$ \begin{aligned} \Delta \Phi _{\mathrm{disk} } = 4\uppi G \rho , \end{aligned} $$(4)

and is assumed to vanish at infinity.

We took the reference frame to be inertial, hence we omitted the net force induced by the disk on the central object. Doing so precludes the development of a strong dipolar moment of the gas distribution, as manifested by eccentric or one-armed spirals (Adams et al. 1989; Heemskerk et al. 1992). One motivation for taking this approximation issues from the mass and momentum flows through the boundaries of the domain, which make it impossible to track the position of a true barycenter over time. In Appendix B we confirm that this acceleration is unimportant as far as our diagnostics on the disk dynamics are concerned (see also Michael & Durisen 2010).

The rightmost term in the internal energy Eq. (3) represents radiative cooling, which is modeled as a local heat sink:

( γ 1 ) q = P ρ c floor 2 ( R ) β Ω 1 , $$ \begin{aligned} \left(\gamma -1\right) q^- = \frac{P-\rho c^2_{\mathrm{floor} }\left(R\right)}{\beta \Omega _{\star }^{-1}}, \end{aligned} $$(5)

where β is a prescribed constant, Ω = G m / R 3 $ \Omega_{\star} = \sqrt{G m_{\star} / R^3} $ is the local Keplerian frequency, and cfloor(R) ≡ 10−3ΩR serves as a precaution to keep the disk aspect ratio h/R ≳ 10−3.

2.3. Units and dimensionless parameters

The above set of equations can be transposed to arbitrary scales. We measured masses relative to the central mass with G = m = 1. We took the disk inner radius rin = 1 as length unit and the inner Keplerian frequency Ω(rin) = 1 as frequency unit. The inner Keplerian period is denoted tin ≡ 2π/Ωin.

A global disk model admits two parameters: its total mass1 relative to the central mass M ≡ mdisk/m and its dimensionless cooling time β. Additionally, the disk thickness can be characterized by its aspect ratio h/R, where h denotes the (local) hydrostatic pressure scale height. In the following, the aspect ratio h/R is always measured as a density-weighted average of cs, iso/vφ, where c s , iso P / ρ $ c_{s,\mathrm{iso}} \equiv \sqrt{P/\rho} $ is the isothermal sound speed – the adiabatic sound speed is simply noted c s γ c s , iso $ c_s\equiv \sqrt{\gamma}c_{s,\mathrm{iso}} $. From the gas surface density Σ + ρ d z $ \Sigma \equiv \smallint_{-\infty}^{+\infty} \rho \mathrm{d} z $ and epicyclic frequency κ one can evaluate the stability of a rotationally supported disk to self-gravitating disturbances via the Toomre number:

Q κ c s , iso π G Σ · $$ \begin{aligned} Q\equiv \frac{\kappa c_{s,\mathrm{iso} }}{\uppi G \Sigma }\cdot \end{aligned} $$(6)

Typically, when Q ≲ 1 the GI extracts orbital energy and converts it into heat via turbulence and shocks. We therefore expect gravito-turbulence to regulate itself and keep the disk close to marginal stability, that is, Q ∼ 1.

2.4. Theoretical background

The Toomre number Eq. (6) can decide on local linear stability: when Q < 1 axisymmetric density disturbances grow exponentially in razor-thin disk models. For slightly larger Q there exist no unstable normal modes; nonetheless the disk may be nonlinearly unstable to sufficiently large perturbations (Riols et al. 2017). On the other hand, there also exist global linear instabilities able to amplify eccentric (Adams et al. 1989) or spiral structures (Noh et al. 1991; Laughlin & Rozyczka 1996). These global modes are standing waves in radius and typically rely on radial features such as boundaries or vortensity extrema (Papaloizou & Savonije 1991). In between the strictly local and global approaches is the study of tightly wound (WKBJ) density waves. WKBJ wave modes are useful to interpret the morphology and transport properties of the large-scale structures emerging in global models of GI turbulence. In this section, we adopt this intermediate approach and run through the ideas most pertinent to our simulations.

2.4.1. Dispersion relations

In the WKBJ approximation, the spiral modes of the disk may be described locally by radial and azimuthal wavenumbers (k,m) such that

d R d φ = m k = R tan ( i ) $$ \begin{aligned} \frac{\mathrm{d} R}{\mathrm{d} \varphi } = -\frac{m}{k} = - R \tan \left(i\right) \end{aligned} $$(7)

describes the morphology of a wave front, with i the pitch angle between the spiral front and a circle around the central object. In this local approximation, linear waves in razor-thin 2D disks obey the dispersion relation

m 2 ( Ω s Ω ) 2 = κ 2 2 π G Σ | k | + c s 2 k 2 = κ 2 ( 1 Q 2 ) + c s 2 ( k k 0 ) 2 , $$ \begin{aligned} m^2 \left(\Omega _s - \Omega \right)^2&= \kappa ^2 - 2\uppi G \Sigma \left|k \right|+ c_s^2 k^2 \nonumber \\&= \kappa ^2 \left(1 - Q^{-2}\right) + c_s^2 \left(k - k_0\right)^2, \end{aligned} $$(8)

where Ω ≡ vφ/R is the gas angular velocity, Ωs ≡ ω/m is the angular pattern speed of a mode with (constant) frequency ω, and k 0  πG Σ/ c s 2 $ k_0 \equiv \uppi G \Sigma / c_s^2 $ is the wavenumber most strongly influenced by self-gravity (see, for example, Binney & Tremaine 2008). When Q < 1, there exist unstable modes whose growth rate is maximal for |k| = k0. When Q = 1, the modes with |k| = k0 are marginally stable and their pattern is corotating with the gas (Ωs = Ω). The linear theory makes no prediction for the preferred azimuthal wavenumber m, however.

It is customary to adopt a quasi-linear theory of GI saturation and attribute the maintenance of a marginal (Q ≈ 1) turbulent state to the action of the marginally stable modes. Their properties are thus of considerable importance and care should be taken so that they are not mischaracterized. A critical feature of the marginal modes at Q = 1 is that their radial wavelength equals that of the disk scale height (k0 = 1/h), and is thus relatively short. This property brings into question the validity of the 2D disk model itself, which requires k ≪ 1/h. It is possible to amend the dispersion relation Eq. (8) so as to account for the “dilution” of the disk’s gravity by its finite thickness:

m 2 ( Ω s Ω ) 2 κ 2 2 π G Σ | k | 1 + | k | / k 0 + c s 2 k 2 $$ \begin{aligned} m^2 \left(\Omega _s - \Omega \right)^2 \simeq \kappa ^2 - \frac{2\uppi G \Sigma \left|k \right|}{1 + \left|k \right|/ k_0} + c_s^2 k^2 \end{aligned} $$(9)

(see Bertin 2014, Eq. (15.18)). The stability threshold is then shifted to Q ≲ 0.647 but marginally stable modes are still corotating with the gas.

However, the disk’s gravity is only one part of the picture. The “quasi-3D” relation Eq. (9) omits the modes’ vertical structure and velocities, which become increasingly important as k approaches 1/h. In nonisothermal disks especially, 3D effects can radically alter the nature of the density waves and their dispersion relations (Lubow & Ogilvie 1998). Moreover, when including self-gravity, Mamatsashvili & Rice (2010) showed that instability first appears on the (usually incompressible) inertial wave branch, a complication that may impact on its radial propagation.

Last but not least, spiral waves excited to significant amplitudes (and consequently exhibiting steep gradients) are nonlinear and thereby couple a range of k and m (Michael et al. 2012). If a nonlinear dispersion relation can be found, it may relate ω and the preponderant (k,m) differently to that predicted by linear theory (e.g., Shu et al. 1985; Fromang & Papaloizou 2007). If the waves shock and thus lose energy, the linear and nonlinear relations will deviate further again. With these problems in mind, one must be wary when constructing a saturation theory, and in particular simulation diagnostics, that rely too heavily on the details of the dispersion relation Eq. (9).

2.4.2. Local turbulence and large-scale spirals

One contentious aspect of gravito-turbulence is whether it can be treated as a local process, coupling the local gas accretion and heating rates as in an α disk (Balbus & Papaloizou 1999). Spiral density waves that carry energy over large distances before damping are central to this problem. From the 2D linear theory of WKBJ modes, the ratio of the wave to turbulent energy fluxes is simply ξ = |(Ωs−Ω)/Ω| (Cossins et al. 2009). Hence the importance of nonlocal energy transport associated with a given wave depends on its proximity to its corotation radius – at which it is marginally stable. But in a quasi-linear theory of GI saturation, gravito-turbulence at any given radius is controlled by the local marginally stable modes. If these modes travel radially they will leave the region that can sustain them and then scatter and dissipate on more vigorous modes at their new location. One might then expect that WKBJ waves never travel far from their corotation radii and thus never contribute appreciably to nonlocal energy transport (Gammie 2001).

We are then left with a picture of gravito-turbulence that comprises an incoherent patchwork of modes with short radial wavelengths and localized to their corotation radii. This picture contrasts with the large-scale spirals reported by numerical simulations and often taken as a starting point to test the GI hypothesis in observed circumstellar disks (e.g., Forgan et al. 2018). One way to resolve this tension is to posit that at any given radius there exist several spiral wave packets distributed across the 2π in azimuth, each corotating with the background orbital flow and each exhibiting the same m and k = k0 ≈ 1/h. Over the course of several orbital periods, localized spirals on neighboring radii will intermittently coalesce and form loose composite structures on a larger radial scale. The apparent coherence of these larger structures will issue from the similar pitch angles of their constituent wakes. If there are a sufficient number of localized wave packets per radius, any snapshot of a simulation will likely show up a sequence of transient agglomerations over several radii of this type. Later in this paper we test this idea.

2.5. Numerical method

We used the version 4.3 of the PLUTO code (Mignone et al. 2007) to integrate Eqs. (1)–(3) in time. We calculated solutions of the self-gravitating disk dynamics for the initial and boundary conditions described below.

2.5.1. Computational domain

The computational domain was defined by intervals in spherical coordinates. The radial interval r ∈ [rin,32rin] was meshed with 518 grid cells with a uniform log(r) sampling. The interval in colatitude θ ∈ [π/2±0.35] was meshed with two different grid spacings. We uniformly meshed the midplane region θ ∈ [π/2±0.05] with 32 cells. We meshed the extreme θ intervals with 32 “stretched” cells on both sides of the midplane, whose angular increment Δθ increases geometrically from θ = π/2 ± 0.05 to the boundaries θ = π/2 ± 0.35. The azimuthal interval φ ∈ [0,2π] was uniformly meshed with 512 grid cells, corresponding to the highest value considered in the resolution studies of Michael et al. (2012) and Steiman-Cameron et al. (2013).

2.5.2. Integration scheme

We used PLUTO to integrate Eqs. (1)–(3) in conservative form via a finite-volume spatial discretization and a second-order Runge–Kutta time stepping. We used a piecewise linear reconstruction (Mignone 2014) with Van Leer’s slope limiter (van Leer 1974) to determine the state of the primitive variables (ρ,v,P) on both sides of each cell interface. These states are converted into conservative variables (ρ,ρv,E), where E ≡ ρ(ϵ+v2/2) is the gas energy density. The energy equation as integrated by PLUTO can be written

t E = · [ ( E + P ) v ] · [ Φ ρ v ] Φ t ρ q , $$ \begin{aligned} \partial _t E = -\nabla \cdot \left[\left(E + P\right)\boldsymbol{v} \right] - \nabla \cdot \left[\Phi \rho \boldsymbol{v} \right] - \Phi \partial _t \rho - q^-, \end{aligned} $$(10)

rightfully excluding the time derivative ∂tΦ of the gravitational potential. We used the HLLc Riemann solver of Toro et al. (1994) to compute the interface fluxes from which the volume-averaged conservative variables are updated. The gas pressure is finally recovered as P/(γ−1) = E − ρv2/2, independent of the gravitational potential energy. The Rankine–Hugoniot jump conditions are built into the Riemann solver, so shock heating is handled without explicit dissipation terms in Eq. (10). This integration scheme is stable with the chosen Courant–Friedrichs–Lewy coefficient of 0.3.

Because solving Poisson’s Eq. (4) can be computationally costly, we updated the gravitational potential of the gas Φdisk only once per time step – not at every Runge–Kutta substep. While the integration scheme becomes formally first-order accurate in time, it is still able to properly resolve the growth of the GI. We present in Appendix A how Poisson’s equation was numerically solved and the related tests.

2.5.3. Initial conditions

We wish to resolve the vertical structure of the disk and to make the dimensionless properties of turbulence independent of radius. The first point above requires that the disk thickness h should be resolved by several grid cells at every radius, that is, cs, iso/vφ ≳ 10−2 near the disk midplane. The second point requires the dimensionless parameters to be independent of radius, including the disk aspect ratio h/R and the Toomre number Q. The angular frequency of a Keplerian disk is Ω ∝ R−3/2, so the aspect ratio is constant when cs, iso ∝ R−1/2. Assuming that gravito-turbulence will maintain Q ≈ 1 in quasi-steady state, we started near Q = 1 at every radius to avoid a long initial cooling phase and to keep the disk geometrically thick. This choice implies Σ ∝ R−2 for the initial gas surface density.

We set the initial Σ(R) = Σin(R/rin)−2 such that m disk ( x ) r in x 2 π Σ R d R $ m_{\mathrm{disk}}\left(x\right) \equiv \smallint_{r_{\mathrm{in}}}^{x} 2\uppi \Sigma R \mathrm{d} R $ equals the prescribed disk mass at the outer radius x = rout. The initial velocity was set to vr = vθ = 0 and

v φ ( R ) = G [ m + m disk ( R ) ] R · $$ \begin{aligned} { v}_{\varphi }\left(R\right) = \sqrt{\frac{G \left[m_{\star } + m_{\mathrm{disk} }\left(R\right)\right]}{R}}\cdot \end{aligned} $$(11)

The isothermal sound speed cs, iso was initialized such that Ωcs, isoGΣ = 1 at every radius. We then distributed the mass density as a simple Gaussian in the vertical direction:

ρ 0 ( R , z ) Σ ( R ) h ( R ) 2 π exp [ 1 2 ( z h ( R ) ) 2 ] . $$ \begin{aligned} \rho _0\left(R,z\right) \equiv \frac{\Sigma \left(R\right)}{h\left(R\right)\sqrt{2\uppi }} \exp \left[-\frac{1}{2}\left(\frac{z}{h\left(R\right)}\right)^2\right]. \end{aligned} $$(12)

These initial conditions are not exactly in quasi-static equilibrium due to omitting the radial pressure gradient and the vertical component of the gas self-gravity in the momentum equation. As we show below, the disk quickly relaxes to a quasi-steady state with no violent redistribution of mass or momentum.

2.5.4. Boundary conditions

Let f(xboundx) = ±f(xbound+x) represent symmetric (+) or antisymmetric (−) conditions for the flow variable f near the boundary of coordinate xbound. We used the same boundary conditions at the inner and outer radial boundaries: we symmetrized ρ and vφ/r, while we antisymmetrized r2vr and vθ. These conditions favor solid-body rotation and disfavor radial in and outflows. On the colatitude θ boundaries we symmetrized ρ and antisymmetrized all three velocity components. These conditions also prevent flows through the domain boundaries, and they disfavor rotational support. Finally, we imposed a small but nonzero pressure in the ghost cells of every boundary, corresponding to the pressure floor given in the section below. This choice prevents internal energy inputs from the boundaries.

These boundary conditions are meant to be simple and yet allow the disk to find a quasi-steady state inside the computational domain. We note that the (anti-)symmetry cannot be imposed exactly due to the nonuniform grid spacings and to the reconstruction and conversion steps involved in the finite-volume scheme. As a consequence, the total mass inside the domain does evolve by a few per cent per hundred inner orbits (see Appendix B), not affecting our conclusions.

2.5.5. Additional precautions

The vertical density stratification and the presence of sonic turbulence are challenging to resolve with a stable and accurate scheme. We took additional precautions to maintain numerical stability while retaining second-order spatial accuracy over most of the computational domain. We imposed a density floor ρfloor = 10−6ρ0(rout,0) from the initial conditions Eq. (12). Since the initial midplane density decreases as R−3, this floor only intermittently affects a few grid cells at the edges of the computational domain. We imposed a pressure floor P floor = ρ floor c floor 2 ( r out ) $ P_{\mathrm{floor}} = \rho_{\mathrm{floor}} c_{\mathrm{floor}}^2\left(r_{\mathrm{out}}\right) $, where cfloor(rout) is the threshold sound speed used in the cooling function Eq. (5). This pressure floor also affects a small fraction of grid cells located near the radial boundaries of the domain and in the vicinity of strong discontinuities. The density and pressure floors being applied independently, the gas temperature may artificially change in the concerned grid cells. As a precaution, we prevented the isothermal sound speed from going below cfloor(R). These floor values are chosen to have a negligible effect on the disk mass and internal energy. Finally, should the relative pressure difference between neighboring cells exceed a factor five, we switched to the more dissipative MINMOD slope limiter (Roe 1986) and to a two-state HLL Riemann solver (Harten et al. 1983) at the concerned cell interfaces. This switch affects between one and six per cent of all grid cells, decreasing with increasing disk thickness. Being restricted to the vicinity of shocks, it does not enhance artificial heating in smooth regions of the flow. To minimize the influence of such procedures we favor thick and vertically well-resolved disks.

2.6. Data reduction

We used various procedures to analyze the simulation data. We only present the main ones here, leaving more elaborate measurement techniques to be introduced in later sections when needed.

Let X be a flow quantity varying in time and over the simulation grid. The midplane distribution of X is denoted as Xmid(r,φ). The time average of X is denoted as X ¯ $ \overline{X} $ and taken over one outer Keplerian period 2π/Ω(rout) ≈ 181tin. Spatial averages are denoted with brackets and defined as a volume-weighted sum over grid cells in the corresponding dimensions. For example, affecting the indices (i,j,k) to cell attributes at coordinates (ri,θj,φk), the azimuthal average of X is:

X φ ( r i , θ j ) = φ k = 0 φ k < 2 π [ X δ V ] i , j , k φ k = 0 φ k < 2 π δ V i , j , k , $$ \begin{aligned} \langle X\rangle _{\varphi }\left(r_i,\theta _j\right) = \frac{\sum _{\varphi _k=0}^{\varphi _k<2\uppi } \left[ X \,\delta V\right]_{i,j,k}}{\sum _{\varphi _k=0}^{\varphi _k < 2\uppi } \,\delta V_{i,j,k}}, \end{aligned} $$(13)

where δV denotes the volume of a grid cell. Looking at geometrically thin disks on a spherical grid, we opted to use ⟨Xθ in place of true vertical averages for simplicity. Radial averages are taken over r/rin ∈ [2,16] and evaluated with a d​log(r) measure matching the logarithmic grid spacing. We define a density-weighted vertical average as

X ρ ( r i , φ k ) = θ j [ X ρ δ V ] i , j , k θ j [ ρ δ V ] i , j , k · $$ \begin{aligned} \langle X\rangle _{\rho }\left(r_i,\varphi _k\right) = \frac{\sum _{\theta _j} \left[ X \rho \,\delta V\right]_{i,j,k}}{\sum _{\theta _j} \left[\rho \,\delta V\right]_{i,j,k}}\cdot \end{aligned} $$(14)

From the vertically and azimuthally averaged profile ⟨Xθ, φ, we define the relative fluctuations

X ( r , φ ) = X ( r , θ , φ ) X θ , φ ( r ) θ X θ , φ ( r ) . $$ \begin{aligned} \tilde{X}\left(r,\varphi \right) = \frac{ \langle X\left(r,\theta ,\varphi \right) - \langle X\rangle _{\theta ,\varphi }\left(r\right)\rangle _{\theta }}{\langle X\rangle _{\theta ,\varphi }\left(r\right) }. \end{aligned} $$(15)

3. Reference simulation

We start by examining one simulation in detail, with a disk mass M = 1/3 of the central mass and a local cooling time β = 10. We label this simulation M3B10 and take it as our reference case for later comparison. We sample different disk masses and cooling times in Sect. 4.

3.1. Initial transient and disk spreading

We show in Fig. 1 the evolution in time of the azimuthally averaged surface density ⟨Σ⟩φ relative to the initial profile Σ(R,t = 0). Because the initial conditions are not in quasi-static equilibrium, the disk first contracts vertically to restore momentum balance. From a local perspective, the time needed to restore a hydrostatic equilibrium is h / c s Ω 1 $ {\sim} h/c_s \simeq \Omega_{\star}^{-1} $, shorter than the imposed cooling time β Ω 1 $ \beta \Omega_{\star}^{-1} $. The radius at which the product Ωt = 1 increases as t2/3, so the contraction front appears to propagate outward at a velocity ∝r−1/2. From a global perspective, this differential contraction sweeps material radially and induces an inertial-acoustic response visible in the gas surface density, as traced by the dashed curves in Fig. 1.

thumbnail Fig. 1.

Evolution in time (horizontal axis) of the radial profile (vertical axis) of the gas surface density relative to its initial profile in run M3B10. The dashed curves on the left fit the initial transient fronts with r(t) ∝ t2/3, i.e., a radial velocity proportional to r−1/2.

After the initial transient, we can distinguish narrow stripes of axisymmetric density perturbations traveling radially at roughly Keplerian velocities. These fluctuations evolve on local orbital timescales while the global redistribution of mass throughout the disk takes several 102tin. The disk has therefore reached a quasi-steady state thermally, at least in the inner r/rin ≲ 16, and it is the dissipation of turbulent motions that balances cooling losses. During this phase, the surface density decreases at small radii and increases at large radii, varying by less than a factor two relative to the initial profile on most of the radial interval. At 500tin the disk has lost fourteen per cent of its initial mass, equivalent to the mass initially located inside 1.6rin. Most of this mass has left the computational domain through the inner radial boundary, while the density increase in r/rin ∈ [4,23] can only result from an outward transport of mass. We explain both trends in Sect. 3.6 by measuring the radial flux of angular momentum.

3.2. Vertical structure of the disk

In the midplane of run M3B10 we measure cs, iso/vφ ≈ 0.04 upon radial and azimuthal averaging, so we resolve one midplane pressure scale height by roughly thirteen grid cells. If the disk was vertically isothermal, the computational domain would cover eight pressure scale heights on both sides of the midplane. As a result, we are well placed to probe the vertical structure and time-dependent dynamics of the disk.

For every grid cell at time 500tin, we divide the local density and sound speed by their azimuthally averaged midplane value at the same spherical radius. We then average these ratios over φ ∈ [0,2π] and r/rin ∈ [2,16]. The resulting profiles represent the mean vertical structure of the disk and are drawn in Fig. 2.

thumbnail Fig. 2.

Radial and azimuthal averages of the density (blue, left axis) and sound speed (red, right axis) relative to their azimuthally averaged midplane value at time 500tin in run M3B10.

The gas density decreases by roughly five orders of magnitude with height at any radius. The density profile is not Gaussian – as initially prescribed – but is instead narrower in the midplane with flat tails at high latitudes. In the midplane region |π/2−θ| ≲ 0.05, the gas density is primarily compressed by the disk’s own gravity. Indeed, the vertical component of the gravitational acceleration is approximately five times larger than the central Φ(r) contribution at small heights.

Away from the midplane, the sound speed increases with height by up to a factor eight, allowing vertical hydrostatic balance with a density profile flatter than Gaussian. This temperature increase is much steeper than what Shi & Chiang (2014, Fig. 1) and Riols et al. (2017, Fig. 2) reported from shearing box experiments with the same cooling prescription. We believe that these differences issue from our imposed boundary conditions. By antisymmetrizing vφ at the θ boundaries we hinder rotational support at high latitudes, and to reach a steady state the gas can either travel radially or become pressure supported. By antisymmetrizing r2vr at the radial boundaries we prevent the development of Bondi–Parker flows, leaving pressure to support the gas against gravity in the radial direction. The combination of a decreasing density and increasing temperature implies that the entropy quickly increases with height, stabilizing the disk against convective motions (Boley et al. 2006).

After azimuthal averaging, we measure radial infall velocities at high latitudes. We find ⟨vrφ ≈ −0.13⟨csφ at |π/2 − θ| = 0.25 and reaching ⟨vrφ ≈ −0.35⟨csφ near the θ boundaries. This mean inflow remains much slower than free fall, confirming that pressure makes up for the lack of rotational support in the disk corona. The turbulent motions being subsonic above the disk, the heat source is presumably localized near the midplane and heat is transported by the vertical velocity fluctuations ≳0.1cs to the corona. Despite the temperature increase with height, the internal energy stored above |π/2 − θ|> 0.1 represents only two per cent of the total internal energy of the disk. Since the corona extracts a negligible amount of heat from the background rotation, we do not expect it to significantly influence the saturation level of the GI.

3.3. Instantaneous equatorial flow

We now look at the instantaneous equatorial flow in run M3B10 at time 500tin, starting traditionally with the gas surface density in Fig. 3. The disk is clearly nonaxisymmetric, which Fig. 1 could not capture because of the azimuthal averaging. The gas density is organized in spiral arms with a relatively tight winding and order-unity density fluctuations along φ at any radius. The largest spiral arms span more than π in azimuth, and some spirals extend to several times their footpoint in radius. However, we are not able to identify a clear-cut global spiral pattern by visually inspecting Fig. 3.

thumbnail Fig. 3.

Gas surface density Σ at time 500tin in run M3B10; the color scale spans three orders of magnitude from black to white.

To examine these structures quantitatively, we move on to the (r,φ) plane and show the midplane distribution of the radial velocity in Fig. 4. This figure reveals transsonic motions at all radii. The velocity fluctuations are oriented both inward and outward in roughly equal proportions. We identify spiral wakes as the oblique lines where the velocity vr = 0 changes sign in a compressive fashion. Bearing in mind the logarithmic vertical scale of Fig. 4, these wakes appear to have a roughly logarithmic spiral morphology. Across each wake, the velocity jumps by more than twice the mean sound speed. We deduce that the flow must shock in the vicinity of the wakes, making them a primary location of entropy generation.

thumbnail Fig. 4.

Radial velocity relative to the azimuthally averaged sound speed in the midplane of run M3B10 at time 500tin. The velocity fluctuations reach supersonic amplitudes both outward (red) and inward (blue).

Because Fig. 4 only represents the instantaneous state of the flow, it is too early to describe the spiral structures as global waves propagating coherently through the disk. To address this issue, we measure their pattern speed and temporal coherence in Sects. 3.8.2 and 3.8.3. Because Fig. 4 only represents a slice of the disk midplane, we also cannot assume that the flow shocks over the entire height of the disk. We show in Sect. 3.8.4 that the compressive midplane motions in fact transit to vertical outflows near the wakes.

3.4. Toomre instability

Figure 4 shows both small-scale fluctuations and large-scale spirals with transsonic amplitudes. To confirm the predominant role of the GI in driving this turbulence, we measure the Toomre number in the quasi-steady state of run M3B10. We used the density-weighted mean orbital velocity Ωρ ≡ ⟨vφ/Rρ to estimate the epicyclic frequency

κ 2 Ω ρ R d R 2 Ω ρ d R $$ \begin{aligned} \kappa \simeq \sqrt{ \frac{2\Omega _{\rho }}{R} \frac{\mathrm{d} R^2\Omega _{\rho }}{\mathrm{d} R}} \end{aligned} $$(16)

in the definition Eq. (6) of the Toomre number. To distinguish the cold and dense disk midplane from its warm corona, we compute the mean sound speed with and without density weighting, denoting the corresponding Toomre numbers as ⟨Qρ and ⟨Qθ, φ respectively. The profiles of Q are averaged over one outer Keplerian period and drawn in Fig. 5.

thumbnail Fig. 5.

Radial profiles of the time-averaged Toomre number Q in run M3B10, using either the density-weighted sound speed ⟨cs, isoρ (solid black) or a simple vertical averaging (dashed red) in Eq. (6). The dotted horizontal line marks the threshold Q = 1 of marginal stability.

With a simple vertical averaging, the Toomre number ⟨Qθ, φ ≈ 5 is significantly larger than unity. However, this value is dominated by the warm corona above the disk midplane, see Fig. 2, and is thus unsuitable. If we instead focus on the midplane by taking the density-weighted mean Toomre number, we find ⟨Qρ ≈ 1 for almost all radii, with fluctuations of less than thirty per cent after time-averaging (Rafikov 2015). Since Q = 1 corresponds to marginal gravitational stability (Toomre 1964), the measured profile ⟨Qρ ≈ 1 verifies that the GI is the main driver of turbulence: if a region of the disk cools below Q < 1, the instability converts orbital energy into heat to restore Q ≳ 1. We therefore expect the orbital energy extraction rate, and a fortiori the angular momentum flux, to depend on the prescribed β cooling timescale (Gammie 2001).

3.5. Absence of other instabilities

The saturation of the GI could be complicated by other processes such as hydrodynamic instabilities unrelated to the self-gravitating nature of the disk, or even by unsuitable boundary conditions (Shu et al. 1990; Lin & Papaloizou 2011). We attempt to rule out several common possibilities by applying simple stability criteria to the quasi-steady state of run M3B10.

We can discard the vertical shear instability (Urpin & Brandenburg 1998) given the relatively long cooling times β > 1 and the stabilizing effect of the steep vertical entropy stratification (Nelson et al. 2013). On the other hand, while there is a vertical shear in the mean radial infall, it is only found high above the disk and is unlikely to affect the flow near the disk midplane.

On the top panel of Fig. 6 we draw the dimensionless shear rate q ≡ ⟨d​logΩ/d​log rρ. Apart from the inner radius where the boundary conditions enforces solid body rotation (q = 0), the shear profile remains Keplerian througout the disk (q = −3/2) and never approaches Rayleigh’s criterion q < −2, allowing us to discount the centrifugal instability.

thumbnail Fig. 6.

Time-averaged radial profiles of three key quantities to assess the disk stability: the dimensionless shear-rate (upper panel), the function ℒ introduced by Lovelace et al. (1999) and defined by Eq. (17), and the vertical component of the vortensity.

To test for the presence of the Rossby wave instability, we followed Lovelace et al. (1999) and computed the function

L ( R ) Σ ω z ( P Σ γ ) 2 / γ θ , φ , $$ \begin{aligned} \mathcal{L} \left(R\right) \equiv \langle \frac{\Sigma }{\omega _z} \left(\frac{P}{\Sigma ^{\gamma }}\right)^{2/\gamma }\rangle _{\theta ,\varphi }, \end{aligned} $$(17)

using ωz ≃ (1/R)d(R2Ωρ)/dR to estimate the mean vertical vorticity. As shown on the middle panel of Fig. 6, this function is maximal at the radial boundaries of the domain and has no extrema in the main body of the disk, thereby dismissing the Rossby wave instability as a source of waves and turbulence, except possibly near the inner radius.

Finally, we show the radial profile of the vortensity ωz/Σ on the bottom panel of Fig. 6. Though the flow is nonbarotropic and the vortensity is not conserved, it may nonetheless suggest nonaxisymmetric instability (Papaloizou & Savonije 1991). This vortensity profile features a small bump near the inner radius and a relatively smooth increase toward the outer radius. The bump could contribute to exciting nonaxisymmetric waves inside r/rin < 1.4, though its influence would likely remain localized near the inner boundary. Over the rest of the disk, the absence of a clear vortensity extremum discounts instabilities of nonaxisymmetric waves corotating with the extremum.

We acknowledge that our application of stability criteria to a fully nonlinear turbulent state is neither rigorous nor exhaustive. It nevertheless supports the idea that our simulations exhibit a relatively pure form of global gravito-turbulence.

3.6. Angular momentum transport

Mass accretion is governed by the equation of angular momentum transport, which can be put in conservative form even for self-gravitating disks (Lynden-Bell & Kalnajs 1972). The radial flux of angular momentum through the disk consists of two terms, one hydrodynamic and the other gravitational. The hydrodynamic (Reynolds) stress arises from correlated velocity fluctuations and can be defined as

R R φ = ρ ( v R v R ρ ) ( v φ v φ ρ ) . $$ \begin{aligned} {R}_{\rm R\varphi } = \rho \left( { v}_{\rm R} - \langle { v}_{\rm R}\rangle _{\rho } \right) \left( { v}_{\varphi } - \langle { v}_{\varphi }\rangle _{\rho } \right). \end{aligned} $$(18)

The gravitational stress comes from correlations between the R and φ gravitational accelerations and takes the form

G R φ = ( R Φ disk ) ( φ Φ disk ) 4 π R G · $$ \begin{aligned} {G}_{\rm R\varphi } = \frac{\left(\partial _{\rm R} \Phi _{\mathrm{disk} }\right) \left(\partial _{\varphi }\Phi _{\mathrm{disk} }\right)}{4\uppi R G}\cdot \end{aligned} $$(19)

We normalize the stress tensors by the vertically averaged gas pressure to obtain dimensionless viscosity coefficients:

α R R R φ P θ , φ ; α G G R φ P θ , φ ; α α R + α G . $$ \begin{aligned} \alpha _{\mathrm{R} } \equiv \frac{\mathrm{R} _{\rm R\varphi }}{\langle P\rangle _{\theta ,\varphi }}; \qquad \alpha _{\mathrm{G} } \equiv \frac{\mathrm{G} _{\rm R\varphi }}{\langle P\rangle _{\theta ,\varphi }}; \qquad \alpha \equiv \alpha _{\mathrm{R} } + \alpha _{\mathrm{G} }. \end{aligned} $$(20)

Assuming that the orbital energy extracted by the stress is locally dissipated into heat, the condition of thermal balance against β cooling (here denoted by LTE) uniquely determines the value of α (Gammie 2001). With the chosen normalization2, this value is:

α LTE = 1 ( d log Ω d log R ) ( γ 1 ) β · $$ \begin{aligned} \alpha _{\mathrm{LTE} } = \frac{1}{\left(-\frac{\mathrm{d}\!\log \Omega }{\mathrm{d}\!\log R}\right) \left(\gamma - 1\right) \beta }\cdot \end{aligned} $$(21)

We show in Fig. 7 the time-averaged radial profiles of stress in run M3B10. The total stress is positive and oscillates by less than thirty per cent of its mean value. The mean value α r , θ , φ ¯ = 0.10 ± 0.01 $ \overline{\langle\alpha\rangle_{r,\theta,\varphi}} = 0.10\pm 0.01 $ is compatible with the predicted αLTE = 0.1 assuming a Keplerian shear rate d​logΩ/d​log R = −3/2. In comparison, the hydrodynamic stress oscillates around zero with small amplitudes: αR ≲ 10−1αLTE. It is therefore the stress induced by fluctuations of the gravitational potential that predominantly transports angular momentum. This is consistent with previous global simulations (Lodato & Rice 2004; Boley et al. 2006; Forgan et al. 2011), although the exact contribution of αR is delicate to measure because of the large velocity deviations relative to their azimuthal average (Michael et al. 2012). We believe that our grid resolution is insufficient to capture small-scale parametric instabilities, which could invigorate the Reynolds stress up to αR ≈ αG/2 at most (Riols et al. 2017; Booth & Clarke 2019).

thumbnail Fig. 7.

Radial profiles of the normalized total stress (solid black) and gravitational stress alone (dashed red) after spatial and time averaging over 181tin in run M3B10. The dotted horizontal line marks the expected value αLTE assuming a Keplerian shear rate in Eq. (21).

In principle, angular momentum could also be transported vertically. The fluctuations of the gravitational potential have the same phase as the density fluctuations, so the vertical gradient ∂zΦdisk has a turning point above density extrema whereas the azimuthal gradient ∂φΦdisk changes sign. The component of the gravitational stress tensor therefore vanishes upon azimuthal averaging on the θ boundaries. In other words, the upper/lower boundaries exert no torque on the disk, so mass accretion is determined by the radial flux of angular momentum (α) alone.

Because α vanishes at the inner radius, the inner radial boundary effectively imposes a stress-free condition. The steep increase of α inside r/rin ≲ 1.4 drives a net inward mass flux in this region and it is responsible for most of the mass losses over time. In the rest of the domain, a constant α means that the angular momentum flux decreases radially as fast as the initially prescribed pressure profile P ∝ R−4, driving a net outward mass flux in agreement with the disk spreading shown in Fig. 1 (see Eq. (28) of Balbus & Papaloizou 1999). We stress that the redistribution of mass does not represent a steady decretion disk. Only local thermal balance is satisfied over the timescales of our simulation, while the global mass profile keeps evolving.

While the time and azimuthally averaged stress coefficient α matches the expected αLTE and can be related to the average mass accretion rate, the root-mean-square dispersion of the stress reaches several 102αLTE at any given radius and time. At any time in the equatorial plane, we find large areas spanning δφ ≳ π/4 and δr/h ≳ 6 exhibiting a negative G < 0 stress. It is remarkable that the stress fluctuations reduce precisely to αLTE upon averaging, but it is also expected since a viscous α disk theory should apply (only) in an average sense to quasi-steady disks (Balbus & Papaloizou 1999). In regards to the dynamics of objects embedded in the disk, they may actually feel a nonvanishing net torque upon orbit-averaging and migrate in a way that viscous disks models cannot account for (Baruteau et al. 2011; Michael et al. 2011; Stamatellos 2015).

3.7. Velocity dispersion

A common extrapolation of viscous disk theory is to connect the Mach number ℳ of turbulent velocity fluctuations with the stress parameter α ∼ ℳ2. This relation was used to extract a stress parameter α from turbulent line broadening (Flaherty et al. 2015, 2018) or from the vertical settling of dust grains against turbulent diffusion (Dubrulle et al. 1995; Pinte et al. 2016). It was also used to constrain the largest grain sizes surviving collisional fragmentation for a prescribed value of α (Birnstiel et al. 2009, 2010). However, the inequality α ≤ ℳ2 holds for a purely hydrodynamic stress R and does not necessarily apply if the stress is dominated by the gravitational contribution G. We verify the validity of this relation by examining the peak-to-peak range of all three velocity components in the disk midplane as a function of radius, as shown in Fig. 8.

thumbnail Fig. 8.

Radial profiles of the peak-to-peak range of the three velocity components in the midplane, relative to the azimuthally averaged midplane sound speed at time 500tin in run M3B10.

The velocity dispersion clearly satisfies the above inequality: the squared turbulent Mach number ℳ2 is two orders of magnitude larger than the average stress coefficient α, so in fact α ≪ ℳ2. This remains true by a factor ∼10 when taking standard deviations over φ of each velocity component instead of their peak-to-peak ranges. Considering that the hydrodynamic stress amounts to αR/α ≲ 10−1 of the total stress at most, we deduce that the vr and vφ velocity fluctuations are weakly correlated on average. The Reynolds stress tensor is essentially diagonal, so the transsonic velocity fluctuations mainly induce an effective pressure. The relation αR ≈ ℳ2 would thus under-estimate the turbulent velocity dispersion by more than two orders of magnitude for a given αR.

The extreme values reported in Fig. 8 most likely correspond to the largest structures developing in the disk, namely, the spiral wakes. From the equatorial velocity components alone we estimate the typical velocity jump Δ v r 2 + Δ v φ 2 4 c s φ $ \sqrt{\Delta \mathit{v}_r^2 + \Delta \mathit{v}_{\varphi}^2} \approx 4 \langle c_s\rangle_{\varphi} $ across the wakes, as previously estimated from Fig. 4. Meanwhile, the vertical velocity fluctuations also reach a significant fraction of the sound speed in the midplane. These vertical motions could promote vertical mixing (Riols & Latter 2018) and alter the morphology and propagation of the spiral wakes, which we examine in the following section.

3.8. Spiral wakes

In this section, we attempt to average the action of the small-scale turbulent fluctuations and wakes so as to examine their cumulative effect, and thus to better isolate the generic properties of large-scale spiral structures in our reference run M3B10.

3.8.1. Spiral pitch angle

We show in Fig. 9 the map of the relative density fluctuations ρ $ \tilde{\rho} $ after vertical averaging at time 500tin in run M3B10. The surface density fluctuates in the range ρ [ 0.9 , 4 ] $ \tilde{\rho} \in \left[-0.9,4\right] $ relative to its azimuthal average at every radius. The density maxima in Fig. 9 are located on the vr = 0 wake fronts visible in Fig. 4. The morphology of these wakes can be parametrized locally by Eq. (7), where the pitch angle i is allowed to vary with radius and from one wake to another. For a pitch angle independent of r, this equation describes logarithmic spirals log(r/r0) = −tan(i)(φφ0) and corresponds to straight oblique lines in Fig. 9.

thumbnail Fig. 9.

Vertically averaged density fluctuations ρ $ \tilde{\rho} $ at time 500tin in run M3B10. The slope of the dashed lines refers to the estimated mean pitch angle tan(i) ≈ 0.25, or i ≈ 14°.

As is clear from Fig. 9, the wake fronts appear to share similar pitch angles over the entire disk. As all the dimensionless numbers (Q, β and h/R) are roughly constant with radius, perhaps this was to be expected. In order to obtain a quantitative handle on these features we designed the following measure of the characteristic pitch angle. Instead of trying to fit lines over every visually identified wakes, as in the figure, we used the information contained in the entire plane by computing the 2D autocorrelation of the relative density fluctuations:

C [ ρ ] ( δ log r , δ φ ) = ρ ( log r , φ ) × ρ ( log r + δ log r , φ + δ φ ) d log r d φ , $$ \begin{aligned} C\left[\tilde{\rho }\right]\left(\delta \log r, \delta \varphi \right) = \int \!\!&\!\int \tilde{\rho } \left(\log r, \varphi \right) \\&\times \tilde{\rho }\left(\log r + \delta \log r, \varphi + \delta \varphi \right) \mathrm{d}\!\log r \mathrm{d} \varphi , \nonumber \end{aligned} $$(22)

which represents how values of ρ $ \tilde{\rho} $ between two points separated by (δ log r,δφ) are correlated on average. It is maximal when (δ log r,δφ) separates two points on density maxima (minima), either along one wake front or separating different fronts. It is minimal when (δ log r,δφ) separates a density maximum from a density minimum on average, for which the values of ρ $ \tilde{\rho} $ have opposite signs. The extrema of the autocorrelation are thus distributed on a series of oblique bands in the (δ log r,δφ) plane, whose slope is the desired pitch angle.

We fitted a line passing through the origin (δ log r,δφ) = (0,0) and maximized the integral of the autocorrelation:

S ( θ ) δ φ k = π π C [ ρ ] ( δ log r = tan ( θ ) δ φ k , δ φ k ) . $$ \begin{aligned} S\left(\theta \right) \equiv \sum _{\delta \varphi _k=-\uppi }^{\uppi } C\left[\tilde{\rho }\right]\left(\delta \log r = -\tan \left(\theta \right) \delta \varphi _k, \delta \varphi _k\right). \end{aligned} $$(23)

The curvature of this integral near the optimal pitch angle was used to estimate an uncertainty δ i = S / ( 2 d 2 S / d θ 2 ) $ \delta i = \sqrt{S / \left(2 \mathrm{d}^2 S / \mathrm{d} \theta^2\right)} $. We repeated this procedure over fifteen simulation snapshots spanning 7.5tin and assimilated the best fit values into a single estimate of the pitch angle and of its uncertainty. The dashed lines drawn in Fig. 9 illustrate the estimated mean pitch angle tan(i) ≈ 0.25. It fits the slope of the wake crests and troughs on average over the whole plane, supporting our assumption of self-similarity. It is over five times larger than the expected tan(i) = ⟨h/Rρ ≈ 0.05 for tidally induced spiral density waves in non-self-gravitating razor-thin disks (Ogilvie & Lubow 2002).

The pitch angle relates the wavenumbers (k,m) via Eq. (7). In run M3B10 we evaluate m ∈ {2,4} for the dominant spirals, with no clear radial trend over time. We use this information in Sect. 4.3 to connect back with the linear dispersion relation Eq. (8).

As all spirals share a similar pitch angle regardless of radius and time, they either propagate without being sheared out, or they emerge intermittently as locally unstable perturbations with the same narrow spectrum of wavenumbers. Should these be understood as standing waves or tightly wound traveling wave trains, their angular pattern speed would be independent of radius. To better understand the nature of these spiral structures we move on to examine their propagation in the disk plane.

3.8.2. Equatorial motions

As an overview, we show in Fig. 10 the relative density fluctuations in four consecutive snapshots spanning two inner Keplerian periods in run M3B10. At all radii the spiral wakes follow the green curves, which are corotating with the gas. Examining smaller radii r/rin ≲ 6, density peaks that are initially disjoint and sitting on the same vertical line eventually meet to form one larger canted wake, as highlighted by specimens in the colored boxes. Meanwhile, other wakes that are extended radially in the first snapshot are sheared apart into smaller disjoint fragments corotating with the flow. In other words, large-scale spiral arms only appear momentarily, when wakes at different radii and drifting at different velocities meet. It then follows that the logarithmic spiral fit illustrated in Fig. 9 describes these individual wake fragments, since large-scale spiral arms only exist in a statistical sense, as a product of many locally corotating wakes sharing the same wavenumbers.

thumbnail Fig. 10.

Relative density fluctuations at four consecutive times separated by 0.5tin in the equatorial plane of run M3B10. The gray scale covers ρ [ 0 , 1 ] $ \tilde{\rho} \in \left[0,1\right] $. The green curves are transported at the mean orbital velocity ⟨vφρ of the gas and serve as corotation references. The blue and red contours delimit initially disjoint wakes in the process of merging by the third snapshot.

To quantify how far the wakes propagate through the gas, we measure their angular pattern speed, assuming that it admits a representative value Ωs at every radius. Our first measurement method uses azimuthal Fourier transforms of the gas surface density. At a given radius we extract the azimuthal phase of a Fourier component, compute the phase speed from consecutive times, and divide it by the azimuthal order m to obtain a pattern speed. We computed the median pattern speed of the 2 ≤ m ≤ 6 Fourier components in 15 snapshots spanning 7.5tin. We denote the resulting pattern speed Ω2 ≤ m ≤ 6 and plot its radial profile in blue in Fig. 11. This curve fluctuates around the gas orbital frequency Ωρ ≡ ⟨vφ/Rρ (in red) with supersonic amplitudes on short length scales. It shows signs of coherent radial propagation only over narrow intervals where Ω2 ≤ m ≤ 6 roughly follows the dashed green lines of constant angular velocity (e.g., at r/rin ∈ [7,10]). Though somewhat noisy, this curve indicates that the disk splits into narrow strips of coherent wave activity, but that overall the pattern speeds track the orbital frequency.

thumbnail Fig. 11.

Radial profiles of angular velocities relative the Keplerian frequency Ω ∝ r−3/2 after averaging over 7.5tin in run M3B10. The red area encompasses the bulk orbital velocity plus or minus the bulk sound speed. The thin blue curve is the spiral pattern speed estimated from the 2 ≤ m ≤ 6 Fourier modes. The thick black curve is the pattern speed estimated by maximizing the time correlation Eq. (24). The dashed green lines are contours of constant angular velocity.

Concerned with the robustness of this estimator against spatial discontinuities and fluctuations in time, we designed an alternative method to measure the spiral pattern speed. At a given radius, we assume that the wake pattern is transported in the φ direction at a constant angular velocity independent of φ and time. The angular velocity Ωcorrelation is estimated by maximizing the correlation Z(ω) between N − 1 pairs of consecutive snapshots separated in time by δtn:

Z ( ω ) n = 1 N 1 0 2 π f ( φ , t ) f ( φ + ω δ t n , t + δ t n ) d φ , $$ \begin{aligned} Z \left( \omega \right) \equiv \sum _{n=1}^{N-1} \int _{0}^{2\uppi } f\left(\varphi , t\right) f\left(\varphi +\omega \,\delta t_n, t+\delta t_n\right)\mathrm{d} \varphi , \end{aligned} $$(24)

where we compressed the density fluctuations using

f tanh ( ρ std ( ρ ) ) $$ \begin{aligned} f \equiv \tanh \left(\frac{\tilde{\rho }}{\mathrm{std}(\tilde{\rho })}\right) \end{aligned} $$(25)

to smooth the extrema and focus the optimization problem on well defined wake patterns. While the Fourier method projects the superposition of nonlinear spirals onto low-order smooth functions of φ, the correlation method is sensitive to the wake pattern as a whole and yields its mean angular velocity. The resulting profile of Ωcorrelation is drawn in black in Fig. 11. The deviations from corotation are on the order of the bulk sound speed and biased toward negative values on average:

Ω correlation Ω ρ Ω ρ = ( 0.5 ± 0.7 ) h / R ρ 2.25 × 10 2 . $$ \begin{aligned} \frac{\Omega _{\mathrm{correlation} } - \Omega _{\rho }}{\Omega _{\rho }} = \left(-0.5 \pm 0.7\right) \langle h/R\rangle _{\rho } \approx 2.25 \times 10^{-2}. \end{aligned} $$(26)

We conclude that the spiral patterns move at the gas orbital velocity to better than ⟨h/Rρ ≈ 0.05 accuracy on average, meaning that they typically propagate at subsonic velocities. Intervals of roughly constant Ωcorrelation are present but they are still narrow and do not necessarily coincide with those of Ω2 ≤ m ≤ 6.

In contrast, Cossins et al. (2009) reported significantly larger values of ξ ≡ |(Ωs−Ω)/Ω|≳0.08 at lower disk masses M ≤ 0.125 (h/R ≲ 0.03), and found that ξ increases with M. At larger masses 0.5 ≤ M ≤ 1.5, Forgan et al. (2011) reported values as high as 0.5 ≤ ξ ≤ 1.5 throughout the disk. This is problematic because WKBJ waves far from corotation should be stable when Q ≈ 1, so the disk would nowhere be unstable. The tension presumably originates in their use of the dispersion relation Eq. (9) to deduce |Ωs − Ω| from measured wavenumbers (k,m), whereas we directly track the gas kinematics.

Following the arguments of Cossins et al. (2009), large values of ξ would imply that a substantial fraction of the energy extracted by the instability is carried away by waves before thermalization. This picture assumes a genuine propagation of the spirals through the gas, which would appear as oblique lines of constant Ωs in Fig. 11. Our profiles of Ωs are never constant over significant distances, so we are lead to believe that the spirals are neither global wave modes nor traveling WKBJ wave trains. The excitation of coherent waves might be facilitated if, unlike in our case, there existed a preferential forcing timescale (e.g., in a circumbinary disk, Desai et al. 2019) or corotation radius (Mejía et al. 2005; Boley et al. 2006), which gravito-turbulence does not seem to produce by itself. To complete our measurements of local corotation, we move on to examine the wakes in their proper corotating frame.

3.8.3. Wake proper motions

If the spiral wakes were strictly corotating with the gas at every radius then they would quickly be sheared out by the disk’s differential rotation. The most probable reason for all wakes to share similar pitch angles is that they form and vanish intermittently, and that we observe them only at their peak amplitude during a cycle of growth, swing and decay (Mamatsashvili & Chagelishvili 2007; Heinemann & Papaloizou 2009, 2012). To test this hypothesis, we select a radius in the disk and follow the flow at its local orbital velocity. Doing so, we can disentangle the proper evolution of a wake from its orbital motion.

Figure 12 shows the evolution of the mean radial velocity ⟨vrρ relative to the azimuthally averaged bulk sound speed ⟨csρ, φ ≡ ⟨⟨csρφ in a frame corotating with the gas at a radius r/rin = 5.2. By looking at a single colored curve, one can inspect the instantaneous azimuthal structure of the flow from left to right. Reciprocally, by picking a specific absissa, one can inspect the evolution of the flow in time (colors) as seen in a frame corotating with the gas. In this frame, the propagation of a traveling wave in the φ direction would appear as a steady horizontal translation of the pattern over time.

thumbnail Fig. 12.

Radial velocity fluctuations in a frame corotating with the gas in run M3B10 at r = 5.2rin. The horizontal axis φ − Ωρt corresponds to an azimuthal coordinate in the comoving frame. The vertical axis is the density-weighted mean radial velocity ⟨vrρ relative to the average sound speed ⟨csρ, φ. Different curves correspond to different times as color-coded in units of local orbital periods. The horizontal segment illustrates how far a signal would propagate in the φ direction at the sound speed ⟨csρ, φ over this time interval.

All times considered, rather than azimuthally translating, the azimuthal profile of ⟨vrρ oscillates around zero up to sonic amplitudes. The first and final times (dark and light, respectively) appear to be roughly anticorrelated along the horizontal axis: a maximum becomes a minimum over this time interval, and vice versa. This sign reversal occurs with moderate and irregular displacements in the horizontal direction, and suggests the presence of oscillation nodes reminiscent of a standing wave. The initial time displayed in Fig. 12 was chosen arbitrarily and we see similar reversals at different radii and in all simulations, pointing to a generic oscillatory phenomenon. We conclude that the spiral patterns are indeed corotating with the gas and do not exhibit the features of traveling wave trains. Since the compressive ⟨vrρ motions are responsible for concentrating the gas into density wakes (see Figs. 4 and 9), the reversal of ⟨vrρ implies their destruction over orbital timescales. The wakes are therefore short-lived intermittent structures, corroborating the statistical nature of the large-scale spiral arms in our simulations.

3.8.4. Instantaneous meridional structure

We have shown that the spiral wakes result from compressive motions in the disk midplane, we now examine the verticality of these motions. We cut a meridional slice through the disk at φ = π, that is, roughly transversal to the wakes. We focus on the front visible at R ≈ 7.5rin in Fig. 9, and represent the surrounding flow in Fig. 13.

thumbnail Fig. 13.

Instantaneous flow in a meridional slice at φ = π, centered on a spiral wake at R/rin ≈ 7.5 and time 500tin in run M3B10. The gas density is represented by filled contours spanning four orders of magnitude. The orientation of the meridional gas velocity is represented with colored arrows, blue (resp. green) being oriented away (resp. toward) the disk midplane. At the wake location, the midplane pressure scale height hmid/rin ≈ 0.3.

Away from the disk midplane (|z|≳0.6rin ≈ 2h) the flow is organized in large scale rolls, as apparent from the horizontally alternating sign (blue and green colors) of the vertical velocity component. These rolls stretch up to the upper/lower θ boundaries and must therefore be controlled to some extent by our boundary conditions. However, their mirror symmetry with respect to z = 0 points to a coherent behavior across the disk midplane, presumably correlated with the wake patterns.

Near the disk midplane (|z|≲0.3rin ≈ h) the flow is essentially horizontal and it converges toward the density maximum from both sides. The flow is smooth on scales ≲h, with no obvious discontinuities in the velocity and density fields, suggesting that shock heating is restricted to a fraction of the disk height. The flow in the vicinity of the wake is mainly oriented away from the midplane (blue color). This reorientation of the horizontal flow into vertical motions is smooth, with no indications of violent splashing nor shock bores (Boley & Durisen 2006). While several smaller rolls can be identified within 2h of the midplane, it is unclear whether the associated circulations are pinned to the density maximum, as the flow is rather disorganized. We stress that the disk is strongly subadiabatic and thus not susceptible to convection. Nonetheless, roll motions can be induced by the baroclinic production of vorticity (Riols & Latter 2018), and this is possibly what we are observing here. The reason that roll motions do not extend as far vertically as shown in Riols & Latter (2018) may be attributed to the sheer strength of the stable stratification for the numerical set-ups we consider.

3.8.5. Mean meridional structure

The superposition of turbulent motions makes it difficult to isolate the flow structure induced by the density wakes alone, so we designed the following filtering procedure. We normalized the density at every radius by its average value at that radius: ρ′ = ρ/⟨ρθ, φ. We also normalized the velocity fluctuations by the bulk orbital velocity: v′ = (v−⟨vρ)/⟨vφρ. Next, we defined the meridional cross-correlation of two quantities X and Y at a given azimuthal angle φ as:

C [ X , Y ] ( δ log r , δ θ ) = X ( log r , θ ) × Y ( log r + δ log r , θ + δ θ ) d log r d θ . $$ \begin{aligned} C\left[X,Y\right]\left(\delta \log r, \delta \theta \right)&= \int \!\!\!\int X\left(\log r, \theta \right)\nonumber \\&\quad \times Y\left(\log r + \delta \log r, \theta + \delta \theta \right) \mathrm{d}\!\log r \mathrm{d} \theta . \end{aligned} $$(27)

By choosing X = ρ′, only the regions of high density fluctuations such as spiral wakes contribute to the integral. We thus obtain the mean meridional distribution of Y as seen from the location a typical density maximum. We took Y = ρv′ as a measure of the momentum fluctuations, including a second ρ′ weight to filter out the large velocity fluctuations away from the disk midplane. Finally, we averaged this meridional cross-correlation over φ and over fifteen snapshots spanning 7.5tin.

The resulting cross-correlation map is shown in Fig. 14. The center (δ log r,δθ) = (0,0) represents the location of a typical wake at an arbitrary radius in the disk. Positive (resp. negative) horizontal coordinates δr/r are located outside (resp. inside) of the wake. The vertical axis −δθ measures an angle from the midplane and spans roughly ±2hmid, so it does not capture the largest rolls visible in Fig. 13 at high altitudes.

thumbnail Fig. 14.

Mean meridional cross-correlation between the density fluctuations ρ′ and the momentum fluctuations ρv′. The arrows give the orientation of the meridional components of the flow while the blue-red colors map its azimuthal component with an arbitrary scale. The right axis indicates the corresponding height measured from the midplane in units of mean midplane pressure scales.

Near the equatorial plane (δθ = 0), the flow converges horizontally toward the wake from both sides before being reoriented vertically away from the midplane. The vertical component C[ ρ , ρ v z ] $ C\left[\rho^\prime,\rho^\prime {\it v}_z^\prime\right] $ of the cross-correlation is in fact significant inside |δr/r|≲0.05, so the orientation of the flow changes mainly within one pressure scale of the wake. Looking at the azimuthal component C[ ρ , ρ v φ ] $ C\left[\rho^\prime,\rho^\prime {\it v}_{\varphi}^\prime\right] $, it is negative inside (δr < 0) of the wake and positive outside. This can be interpreted as a consequence of angular momentum conservation during the radial motions of the gas. The absence of apparent backflow returning to the midplane from high altitudes is caused by the density weights incorporated in Eq. (27), biasing the integral to keep the structures closest to the wake. In summary, Figs. 13 and 14 illustrate the genuinely 3D nature of the flow around spiral wakes and their ability to promote vertical mixing even in global simulations with a steep entropy stratification.

4. Varying the disk mass and cooling time

In this section we sample different dimensionless cooling times β and disk masses M relative to the central object. To control the disk mass, we rescale the gas surface density Σ while keeping the radial dependence Σ(R) ∝ R−2. We also keep the initial Toomre number Q ≈ 1 by varying the initial sound speed in the same proportion as Σ. As cooling brings the disk toward a gravitationally unstable state, we always find that the density-weighted mean ⟨Qρ ≈ 1 in the turbulent phase. All simulations are integrated over more than 300tin to characterize quasi-steady states. Since Σ evolves slowly over the duration of a simulation, the density-weighted mean sound speed remains close to its initial profile and different disk masses correspond to different aspect ratios ⟨h/Rρ ≃ πΣR2/m. For example, we verified that πΣR2/m = (1.09 ± 0.05)⟨h/rρ in run M3B10. The simulations are listed in Table 1 with their main parameters and diagnostics.

Table 1.

Simulation label, initial disk mass, cooling timescale, integration time, mean aspect ratio, mean dimensionless stress, mean spiral pitch angle and mean standard deviation of the gas surface density.

4.1. Low-mass disk

We start with run M10B10, which has the lowest mass considered, M = 1/10 and the same cooling timescale β = 10 as previously. For comparison, a global SPH simulation with the same parameters was presented by Cossins et al. (2009), albeit with a different slope for the surface density profile Σ ∝ R−3/2. Lodato & Rice (2004) also presented a simulation with M = 1/10 but with a surface density Σ ∝ R−1 and a shorter cooling timescale β = 7.5.

4.1.1. Overview

The disk of run M10B10 loses less than seven per cent of its mass over 500tin. After the initial transient, the Toomre number settles at ⟨Qρ ≈ 1 with fluctuations less than ten per cent across the disk. Given the smaller surface density Σ, marginal stability is attained for a smaller aspect ratio ⟨h/Rρ ≈ 1.4 × 10−2. With less than five grid cells per midplane scale height, the vertical structure of the disk is only marginally resolved in this case.

We show the mean vertical structure of this disk at time 500tin in Fig. 15. The gas density profile is strongly peaked in the midplane and nearly flat at high latitudes. The sound speed increases with height by up to a factor twenty relative to its midplane value. The temperature high above the disk is in fact comparable to that of run M3B10, confirming that this region is predominantly supported by pressure against the gravity of the central object both in the vertical and radial directions.

thumbnail Fig. 15.

Same as Fig. 2 at time 500tin in the least massive case M10B10.

4.1.2. Spiral patterns

We show in Fig. 16 the relative density fluctuations in the equatorial plane of run M10B10 at time 500tin. Spiral wakes fill the disk up to r/rin ≲ 22 while turbulence is still developing in the outermost regions. The peak-to-peak equatorial velocity dispersion Δ v r 2 + Δ v φ 2 5 c s φ $ \sqrt{\Delta \mathit{v}_r^2 + \Delta \mathit{v}_{\varphi}^2} \approx 5 \langle c_s\rangle_{\varphi} $ amounts to a similar turbulent Mach number as in run M3B10. Accordingly, the density fluctuations reach similar amplitudes as in run M3B10 with a standard deviation std(Σ)/Σ ≈ 0.63 (see Table 1).

thumbnail Fig. 16.

Same as Fig. 9 in the least massive case M10B10. The dashed lines refer to a pitch angle of tan(i) ≈ 0.23, or i ≈ 13°.

Compared to run M3B10, the wake fronts are narrower and closer to each other, with a typical azimuthal orders m ≈ 8. We employ the method described in Sect. 3.8.1 to measure the mean spiral pitch angle, using 21 snapshots spanning 10tin. We find that tan(i) ≈ 0.23, slightly smaller than in run M3B10. This mean pitch angle is illustrated by the dashed oblique lines in Fig. 16, which fit the mean slope of the density wakes at all radii.

Using the procedure described in Sect. 3.8.2 over the same 21 snapshots, we measure the mean angular velocity of the spiral density patterns. In this case, the deviations from corotation amount to

Ω correlation Ω ρ Ω ρ = ( 0.4 ± 1.0 ) h / R ρ 5.6 × 10 3 $$ \begin{aligned} \frac{\Omega _{\mathrm{correlation} } - \Omega _{\rho }}{\Omega _{\rho }} = \left(-0.4 \pm 1.0\right) \langle h/R\rangle _{\rho } \approx 5.6\times 10^{-3} \end{aligned} $$(28)

after radial averaging, smaller than in run M3B10.

Moving to a corotating frame as in Sect. 3.8.3, we retrieve standing oscillations of the radial velocity over orbital timescales, such as appear in Fig. 12, with little evidence for the propagation of the wake patterns. As in run M3B10, we conclude that the spiral patterns are nearly corotating with the gas and that they are destroyed and reborn before they can travel over significant radial distances.

4.2. High-mass disk

Our most massive disk M1B10 has an initial mass equal to the central mass (M = 1). Similarly massive disks were examined by Lodato & Rice (2005) and Forgan et al. (2011), again with different Σ profiles. After the Toomre number has saturated to ⟨Qρ ≈ 1, the average thickness of the disk ⟨h/Rρ ≈ 0.1. While the dominant spiral pattern has an azimuthal order m = 2, one-armed spirals also occasionally appear. Such asymmetries displace the barycenter of the disk and should in principle cause the central object to react, an effect we neglect here but briefly consider in Appendix B. Bearing this drawback in mind, we carry on to evaluate the pattern speed of the wakes in run M1B10.

We estimate the angular pattern speed Ωs of the wakes via the two methods explained in Sect. 3.8.2, using 41 simulation snapshots over 20tin. The corresponding pattern speed profiles are represented in Fig. 17 along with the bulk orbital velocity Ωρ ≡ ⟨vφ/Rρ of the gas. The two methods are in relatively good agreement: the wake pattern speed is significantly slower than the gas in the inner r/rin ∈ [2,8], corotating with the gas at r/rin ≈ 16 and faster than the gas further out. Moving to a frame corotating with the gas at r/rin ≈ 5 as in Sect. 3.8.3, we do see a translation of the wake pattern as expected for traveling waves. These could be indications of spiral waves being excited near r/rin ≈ 16 and propagating radially.

thumbnail Fig. 17.

Same as Fig. 11 after averaging over 41 snapshots spanning 20tin in the most massive case M1B10.

After saturation of the GI, we observe a redistribution of the gas vorticity ωz and surface density Σ such that their ratio exhibit maxima at r/rin ≈ 1.4 and 18, see Fig. 18. Such a vortensity profile only appears in run M1B10 and it can destabilize global waves even without invoking self-gravity (Papaloizou & Savonije 1991). Such instabilities should be attributed to the choice of initial and boundary conditions that lead to this vortensity profile, and could presumably be reproduced in lower-mass disks as well. In addition, we note that the Ωcorrelation and Ω2 ≤ m ≤ 6 curves only follow the lines of constant pattern speed near r/rin ≈ 16 but clearly not over r/rin ∈ [4,32]. The deviations from corotation are saturated at sonic amplitudes, suggesting that waves cannot freely propagate due to shock dissipation. In summary, if global waves are excited in run M1B10, they may retain a local character by dissipating energy close to where it is extracted – namely, near corotation.

thumbnail Fig. 18.

Radial profile of the vortensity in run M1B10. The maximum at r/rin ≈ 18 coincides with the corotation radius in Fig. 17.

4.3. Spiral spectrum

From now on we consider our simulation set as a whole. In every simulation m and tan(i) show no consistent variations with radius, so k ∝ R−1 on average and the spirals keep a logarithmic morphology. We list in Table 1 the pitch angles measured via the method described in Sect. 3.8.1. They differ from each other by twenty per cent overall, despite the cooling timescale β and the disk mass M varying by a factor ten. Considering only the simulations with a cooling time β = 10, we can fit a linear dependence of the pitch angle on the aspect ratio: tan(i) ≈ 0.222 + 0.602 × ⟨h/Rρ.

Using the measured tan(i) and the density-weighted average k 0  πG Σ/ c s,iso ρ 2 $ k_0 \simeq \uppi G \Sigma / \langle c_{s,\mathrm{iso}}\rangle_{\rho}^2 $, we estimate that 2D linear waves corotating with the gas at any location should possess m0 = k0Rtan(i) ≈ 6 in run M3B10, m0 ≈ 9 in run M5B10 and m0 ≈ 16 in run M10B10. Instead, we identify spirals of order m ∈ {2,4} in run M3B10, m ∈ {4,5} in run M5B10 and m ∈ {7,9} in run M10B10. Since the dominant azimuthal wavenumber is m ≈ m0/2, the radial wavenumber must be k ≈ k0/2 to account for the measured pitch angles. The discrepancy k ≠ k0 should come as no surprise given the several issues raised3 in Sect. 2.4.1. On the other hand, the consistent scaling of k with k0 across the various simulations supports the basic thrust of the quasi-linear theory.

The ratio of k/k0 obtained by Cossins et al. (2009) can be deduced from their Fig. 13. They show that ms−Ω)/kcs ≈ 1 regardless of disk mass and cooling time when injecting the measured mean wavenumbers (k,m) at each radius in the WKBJ dispersion relation Eq. (9). This implies that the first two terms on the right hand side of Eq. (9) cancel each other, or equivalently k = k0 when Q = 1. Their dominant radial wavenumber is thus twice larger than ours. It is possible that the discrepancy issues from the different vertical structures and discretizations of the disk. In particular, in their SPH simulations the disk scale height h ≲ 0.03R was resolved by roughly two smoothing lengths (see their Fig. B5). If the spirals observed by Cossins et al. (2009) were closely corotating with the gas – as we found in all cases – then the equality k = k0 could indicate that vertical structures and motions were not fully resolved, so that the flow was effectively 2D. It is also possible that their mean wavenumber k was biased to larger values by the skewed power spectrum which they used as an averaging weight (see their Fig. 10).

We may extend our results to parametrize the predominant wavenumbers (k,m) and describe the spiral structure analytically. For the azimuthal wavenumber:

m = k R tan ( i ) 0.12 ( h / R ) 1 , $$ \begin{aligned} m = k R \tan \left(i\right) \approx 0.12 \left( h / R \right)^{-1}, \end{aligned} $$(29)

with the limit m = 1 of one-armed spirals requiring a proper account of the reaction of the central object (Adams et al. 1989). The morphology of a spiral arm may then be described by

| d φ d R | = k m 1 2 κ 2 m π G Σ , $$ \begin{aligned} \left|\frac{\mathrm{d} \varphi }{\mathrm{d} R} \right|= \frac{k}{m} \approx \frac{1}{2} \frac{\kappa ^2}{m\uppi G \Sigma }, \end{aligned} $$(30)

which is meaningful only on short radial intervals (m must not change) and in a statistical sense given the local and transient nature of the spiral wakes. For our choice of Σ profile, the right hand side of Eq. (30) is ∝R−1 and, on integrating, we obtain the logarithmic spiral structure seen throughout all our simulations. As mentioned, the crude relations Eqs. (29) and (30) may be sensitive to the detailed stratification of the disk and we caution against using them to quantitatively identify self-gravitating spirals in observations until this issue is settled.

4.4. Density contrast

As is apparent in Figs. 9 and 16 the spiral wakes are narrow and separated by large areas of low density. At any given radius, the maximum surface density fluctuation is approximately five times larger than the standard deviation std(Σ) at that radius in run M3B10, six times in run M3B32 and seven times in run M3B100. In comparison, this ratio would be 2 1.4 $ \sqrt{2} \approx 1.4 $ for linear waves with a sinusoidal behavior. Despite nonlinearities, there exist a simple scaling between the density dispersion δΣ and the cooling timescale: δΣ/Σ ∝ β−1/2 (Cossins et al. 2009). This relation could be used to constrain the observability of spirals in self-gravitating disks (Hall et al. 2016; Cadman et al. 2020).

We find that the formula std ( Σ ) / Σ = 2 / β $ \mathrm{std}(\Sigma)/\Sigma = 2/\sqrt{\beta} $ is compatible with all the standard deviations gathered in Table 1. It is however twice larger than the density dispersion measured by Cossins et al. (2009) at similar aspect ratios h/R ≲ 0.03. In Appendix B we show that the density contrast does not change with the grid resolution, as reported by Cossins et al. (2009, see their Fig. A1). If the discrepancy is not a direct consequence of our different spatial discretizations, then it may result from the different thermal stratifications and vertical motions in the two codes.

4.5. Angular momentum transport

From the requirement of local thermal balance we expect the dimensionless stress α to take the constant value αLTE defined by Eq. (21) at all radii. The time and radially averaged values of α gathered in Table 1 are all compatible with αLTE to ten per cent accuracy. The largest deviation from αLTE is found in the least massive case M10B10. There, the gravitational stress αG ≈ 1.2 × 10−1 is larger than αLTE = 10−1, so the inequality α < αLTE comes from a negative contribution of the average Reynolds stress. The results presented in Sect. 4.1 support the idea of local thermal equilibrium in run M10B10, so this slight departure α ≠ αLTE likely comes from the marginal resolution of the disk in the vertical dimension (Michael et al. 2012). Overall, the agreement between the measured α and the theoretical αLTE supports the hypothesis of local energy conversion for all the disks considered.

5. Summary

We simulated the global 3D hydrodynamics of gravitational instability with a shock-capturing grid-based code in self-similar models of self-gravitating disks. We considered disk masses in the range M ∈ [1/10,1] relative to that of the central object, for which a Toomre number Q = 1 corresponds to aspect ratios h/R ∈ [0.01,0.1]. We prescribed a local cooling law with a characteristic timescale β ∈ [10,100] to force the disk into a gravitationally unstable state. We characterized the ensuing gravito-turbulence with a focus on the spiral patterns commonly observed in self-gravitating disks. Our main findings may be summarized as follows.

  1. Nonlocal energy transport is negligible in all our simulations: the conversion from orbital to thermal energy happens at the same radius. The angular momentum flux is dominated by the gravitational stress and it matches the value expected for viscous α disks in thermal equilibrium upon averaging.

  2. Gravito-turbulence generates spiral density wakes that intermittently form and vanish over orbital timescales. These wakes are corotating with the gas at every radius to a good approximation. Large-scale spiral arms only manifest transiently through the coalescence of several neighboring wakes which then are sheared apart.

  3. The spiral wakes result from nearly symmetric compressive motions of supersonic amplitude in the midplane, making them the main site of energy dissipation. Simultaneously, these wakes promote vertical mixing by lifting the gas off from the midplane over several disk scale heights.

  4. At every radius, the prominent spiral wakes are restricted to a narrow range of pitch angles which only weakly depend on the disk mass and cooling time. We can thus predict the morphology of large-scale apparent self-gravitating spirals for given disk rotation and mass surface density profiles.

The numerical experiments presented in this paper cover a limited parameter space with a prescribed disk structure and simplified thermodynamics. In particular, our self-similar disk model is designed to promote locality in GI turbulence, and to reduce global effects as might be generated by edges or radial substructures. Circumstellar disks are known to exhibit substructures (Huang et al. 2018) and to be gravitationally stable in their inner regions due to stellar irradiation (D’Alessio et al. 1998). Such radial structure would certainly affect the excitation and propagation of waves in massive disks (Lin & Shu 1964; Papaloizou & Savonije 1991; Lin & Papaloizou 2011). Global waves may also result from introducing a characteristic tidal forcing, as is the case in circumbinary disks for example (Marzari et al. 2009; Desai et al. 2019).

Concerning the gas thermodynamics, we excluded the β ≲ 3 regime of disk fragmentation to keep the disk relatively smooth and well-behaved. While disk fragmentation is an avenue of research on its own, it is also unclear how a more detailed radiative energy transport and the thermal stratification of the disk affect the spectrum of unstable perturbations and their nonlinear saturation in global geometry (Tsukamoto et al. 2015; Hirose & Shi 2017, 2019).

Finally, the vertical flows accompanying spiral wakes can entrain dust grains away from the disk midplane (Riols et al. 2020). GI turbulence may thus affect the vertical segregation of different dust sizes (Pinte et al. 2016; Villenave et al. 2020) or induce local enhancements of the dust disk thickness (Benisty et al. 2015; Dong et al. 2015). These 3D flows also open a new route of magnetic field amplification via a spiral dynamo in sufficiently ionized disks (Riols & Latter 2019; Deng et al. 2020).


1

The disk mass is only used as a scaling factor. With the imposed profile Σ ∝ R−2, the disk mass would not converge to a finite value when extending the computational domain radially to infinity.

2

The different choice of normalization used by Gammie (2001) introduces an additional γ(d​logΩ/d​log R) in the denominator.

3

From the vertical dilution of the potential alone, we would expect k/k0 ≈ 0.466 for the marginally stable modes of Eq. (9) at Q ≈ 0.647.

Acknowledgments

We thank Richard P. Nelson for his feedback on our preliminary results. We thank our referee, Doug Lin, for his questions and comments that helped clarify the final manuscript. William Béthune acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG) through Grant KL 650/31-1. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the DFG through grant no INST 37/935-1 FUGG.

References

  1. Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959 [NASA ADS] [CrossRef] [Google Scholar]
  2. ALMA Partnership, (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [Google Scholar]
  3. Baehr, H., Klahr, H., & Kratter, K. M. 2017, ApJ, 848, 40 [Google Scholar]
  4. Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. 1997, in Modern Software Tools in Scientific Computing, eds. E. Arge, A. M. Bruaset, & H. P. Langtangen (Basel: Birkhäuser Press), 163 [Google Scholar]
  5. Balay, S., Abhyankar, S., Adams, M. F., et al. 2019, PETSc Web page: https://www.mcs.anl.gov/petsc [Google Scholar]
  6. Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 [NASA ADS] [CrossRef] [Google Scholar]
  7. Baruteau, C., Meru, F., & Paardekooper, S.-J. 2011, MNRAS, 416, 1971 [Google Scholar]
  8. Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Bertin, G. 2014, Dynamics of Galaxies, 2nd edn. (Cambridge: Cambridge University Press) [Google Scholar]
  10. Béthune, W., & Rafikov, R. R. 2019, MNRAS, 487, 2319 [Google Scholar]
  11. Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn., Princeton Series in Astrophysics (Princeton: Princeton University Press) [Google Scholar]
  12. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Boley, A. C., & Durisen, R. H. 2006, ApJ, 641, 534 [Google Scholar]
  15. Boley, A. C., & Durisen, R. H. 2008, ApJ, 685, 1193 [Google Scholar]
  16. Boley, A. C., Mejía, A. C., Durisen, R. H., et al. 2006, ApJ, 651, 517 [Google Scholar]
  17. Booth, R. A., & Clarke, C. J. 2016, MNRAS, 458, 2676 [Google Scholar]
  18. Booth, R. A., & Clarke, C. J. 2019, MNRAS, 483, 3718 [Google Scholar]
  19. Boss, A. P. 1998, ApJ, 503, 923 [Google Scholar]
  20. Boss, A. P. 2000, ApJ, 536, L101 [CrossRef] [Google Scholar]
  21. Boss, A. P., & Durisen, R. H. 2005, ApJ, 621, L137 [Google Scholar]
  22. Cadman, J., Hall, C., Rice, K., Harries, T. J., & Klaassen, P. D. 2020, MNRAS, 498, 4256 [Google Scholar]
  23. Clarke, C. J., & Lodato, G. 2009, MNRAS, 398, L6 [Google Scholar]
  24. Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157 [Google Scholar]
  25. D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 [Google Scholar]
  26. Deng, H., Mayer, L., & Latter, H. 2020, ApJ, 891, 154 [Google Scholar]
  27. Desai, K. M., Steiman-Cameron, T. Y., Michael, S., Cai, K., & Durisen, R. H. 2019, MNRAS, 483, 2347 [Google Scholar]
  28. Dipierro, G., Pinilla, P., Lodato, G., & Testi, L. 2015, MNRAS, 451, 974 [Google Scholar]
  29. Dong, R., Hall, C., Rice, K., & Chiang, E. 2015, ApJ, 812, L32 [Google Scholar]
  30. Dong, R., Vorobyov, E., Pavlyuchenkov, Y., Chiang, E., & Liu, H. B. 2016, ApJ, 823, 141 [Google Scholar]
  31. Dong, R., Najita, J. R., & Brittain, S. 2018, ApJ, 862, 103 [NASA ADS] [CrossRef] [Google Scholar]
  32. Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] [Google Scholar]
  33. Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 607 [Google Scholar]
  34. Eisner, J. A., Hillenbrand, L. A., Carpenter, J. M., & Wolf, S. 2005, ApJ, 635, 396 [Google Scholar]
  35. Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99 [Google Scholar]
  36. Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117 [NASA ADS] [CrossRef] [Google Scholar]
  37. Forgan, D., Rice, K., Cossins, P., & Lodato, G. 2011, MNRAS, 410, 994 [Google Scholar]
  38. Forgan, D. H., Ilee, J. D., Cyganowski, C. J., Brogan, C. L., & Hunter, T. R. 2016, MNRAS, 463, 957 [Google Scholar]
  39. Forgan, D. H., Ilee, J. D., & Meru, F. 2018, ApJ, 860, L5 [Google Scholar]
  40. Fromang, S., & Papaloizou, J. 2007, A&A, 468, 1 [EDP Sciences] [Google Scholar]
  41. Gammie, C. F. 2001, ApJ, 553, 174 [Google Scholar]
  42. Goldreich, P., & Lynden-Bell, D. 1965a, MNRAS, 130, 97 [NASA ADS] [CrossRef] [Google Scholar]
  43. Goldreich, P., & Lynden-Bell, D. 1965b, MNRAS, 130, 125 [Google Scholar]
  44. Hall, C., Forgan, D., Rice, K., et al. 2016, MNRAS, 458, 306 [Google Scholar]
  45. Harten, A., Lax, P. D., & Leer, B. V. 1983, SIAM Rev., 25, 35 [Google Scholar]
  46. Heemskerk, M. H. M., Papaloizou, J. C., & Savonije, G. J. 1992, A&A, 260, 161 [NASA ADS] [Google Scholar]
  47. Heinemann, T., & Papaloizou, J. C. B. 2009, MNRAS, 397, 52 [NASA ADS] [CrossRef] [Google Scholar]
  48. Heinemann, T., & Papaloizou, J. C. B. 2012, MNRAS, 419, 1085 [Google Scholar]
  49. Hill, G. W. 1878, Am. J. Math., 1, 5 [Google Scholar]
  50. Hirose, S., & Shi, J.-M. 2017, MNRAS, 469, 561 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hirose, S., & Shi, J.-M. 2019, MNRAS, 485, 266 [Google Scholar]
  52. Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42 [Google Scholar]
  53. Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292 [Google Scholar]
  54. Kolykhalov, P. I., & Syunyaev, R. A. 1980, Sov. Astron. Lett., 6, 357 [Google Scholar]
  55. Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271 [Google Scholar]
  56. Laughlin, G., & Bodenheimer, P. 1994, ApJ, 436, 335 [Google Scholar]
  57. Laughlin, G., & Rozyczka, M. 1996, ApJ, 456, 279 [Google Scholar]
  58. Laughlin, G., Korchagin, V., & Adams, F. C. 1997, ApJ, 477, 410 [Google Scholar]
  59. Laughlin, G., Korchagin, V., & Adams, F. C. 1998, ApJ, 504, 945 [Google Scholar]
  60. Lin, M.-K., & Kratter, K. M. 2016, ApJ, 824, 91 [Google Scholar]
  61. Lin, M.-K., & Papaloizou, J. C. B. 2011, MNRAS, 415, 1445 [NASA ADS] [CrossRef] [Google Scholar]
  62. Lin, D. N. C., & Pringle, J. E. 1987, MNRAS, 225, 607 [Google Scholar]
  63. Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  64. Liu, H. B., Takami, M., Kudo, T., et al. 2016, Sci. Adv., 2, e1500875 [Google Scholar]
  65. Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630 [Google Scholar]
  66. Lodato, G., & Rice, W. K. M. 2005, MNRAS, 358, 1489 [Google Scholar]
  67. Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [Google Scholar]
  68. Lubow, S. H., & Ogilvie, G. I. 1998, ApJ, 504, 983 [Google Scholar]
  69. Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1 [Google Scholar]
  70. Mamatsashvili, G. R., & Chagelishvili, G. D. 2007, MNRAS, 381, 809 [Google Scholar]
  71. Mamatsashvili, G. R., & Rice, W. K. M. 2010, MNRAS, 406, 2050 [Google Scholar]
  72. Mann, R. K., Andrews, S. M., Eisner, J. A., et al. 2015, ApJ, 802, 77 [Google Scholar]
  73. Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Mejía, A. C., Durisen, R. H., Pickett, M. K., & Cai, K. 2005, ApJ, 619, 1098 [Google Scholar]
  75. Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJ, 839, L24 [Google Scholar]
  76. Michael, S., & Durisen, R. H. 2010, MNRAS, 406, 279 [Google Scholar]
  77. Michael, S., Durisen, R. H., & Boley, A. C. 2011, ApJ, 737, L42 [Google Scholar]
  78. Michael, S., Steiman-Cameron, T. Y., Durisen, R. H., & Boley, A. C. 2012, ApJ, 746, 98 [Google Scholar]
  79. Mignone, A. 2014, J. Comput. Phys., 270, 784 [Google Scholar]
  80. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
  81. Müller, E., & Steinmetz, M. 1995, Comput. Phys. Commun., 89, 45 [Google Scholar]
  82. Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [NASA ADS] [CrossRef] [Google Scholar]
  83. Noh, H., Vishniac, E. T., & Cochran, W. D. 1991, ApJ, 383, 372 [Google Scholar]
  84. Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950 [NASA ADS] [CrossRef] [Google Scholar]
  85. Paczynski, B. 1978, Acta Astron., 28, 91 [Google Scholar]
  86. Papaloizou, J. C., & Savonije, G. J. 1991, MNRAS, 248, 353 [Google Scholar]
  87. Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [Google Scholar]
  88. Podolak, M., Hubbard, W. B., & Pollack, J. B. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 1109 [Google Scholar]
  89. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
  90. Rafikov, R. R. 2015, ApJ, 804, 62 [Google Scholar]
  91. Rice, W. K. M., Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 1025 [Google Scholar]
  92. Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543 [NASA ADS] [CrossRef] [Google Scholar]
  93. Richardson, D. C. 1994, MNRAS, 269, 493 [Google Scholar]
  94. Riols, A., & Latter, H. 2018, MNRAS, 476, 5115 [Google Scholar]
  95. Riols, A., & Latter, H. 2019, MNRAS, 482, 3989 [Google Scholar]
  96. Riols, A., Latter, H., & Paardekooper, S. J. 2017, MNRAS, 471, 317 [Google Scholar]
  97. Riols, A., Roux, B., Latter, H., & Lesur, G. 2020, MNRAS, 493, 4631 [Google Scholar]
  98. Roe, P. L. 1986, Ann. Rev. Fluid Mech., 18, 337 [Google Scholar]
  99. Safronov, V. S. 1960, Ann. d’Astrophys., 23, 979 [Google Scholar]
  100. Salo, H. 1992, Nature, 359, 619 [Google Scholar]
  101. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
  102. Shi, J.-M., & Chiang, E. 2014, ApJ, 789, 34 [Google Scholar]
  103. Shlosman, I., & Begelman, M. C. 1987, Nature, 329, 810 [Google Scholar]
  104. Shu, F. H., Yuan, C., & Lissauer, J. J. 1985, ApJ, 291, 356 [Google Scholar]
  105. Shu, F. H., Tremaine, S., Adams, F. C., & Ruden, S. P. 1990, ApJ, 358, 495 [Google Scholar]
  106. Stamatellos, D. 2015, ApJ, 810, L11 [Google Scholar]
  107. Steiman-Cameron, T. Y., Durisen, R. H., Boley, A. C., Michael, S., & McConnell, C. R. 2013, ApJ, 768, 192 [Google Scholar]
  108. Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
  109. Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Tsukamoto, Y., Takahashi, S. Z., Machida, M. N., & Inutsuka, S. 2015, MNRAS, 446, 1175 [Google Scholar]
  111. Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399 [Google Scholar]
  112. van Leer, B. 1974, J. Comput. Phys., 14, 361 [Google Scholar]
  113. Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164 [EDP Sciences] [Google Scholar]
  114. Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137 [Google Scholar]
  115. Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115 [Google Scholar]
  116. Yang, L. T., & Brent, R. P. 2002, in Fifth International Conference on Algorithms and Architectures for Parallel Processing (IEEE), 324 [Google Scholar]

Appendix A: Self-gravity in PLUTO

A.1. Poisson solver

To determine the gravitational potential of the gas we solved Poisson’s Eq. (4) with a method similar to that presented in Appendix B of Béthune & Rafikov (2019). This equation is cast into a matrix inversion problem after spatial discretization. We used the IBiCGStab algorithm (Yang & Brent 2002) implemented in the PETSC library (Balay et al. 1997, 2019) to handle the parallel inversion of the Laplacian matrix. The convergence of the solver was characterized by a tolerance of 10−9 on the relative decrease of the residual norm between successive iterations.

The Laplacian operator was discretized by second-order finite differences, accounting for the nonuniform spacing and spherical geometry of the grid. Boundary conditions were applied in the rows corresponding to boundary cells in the matrix and right-hand-side vector. We used periodic boundary conditions in the azimuthal direction and Dirichlet (imposed Φdisk value, see below) on all the other boundaries. The Poisson problem was solved once per hydrodynamic timestep. Although the time-integration scheme formally becomes first-order accurate, Béthune & Rafikov (2019) showed that it is able to accurately reproduce the growth rates of Jeans instability.

A.2. Implementation of boundary conditions

In its differential form, the Poisson problem allows any choice of values for Φdisk on the (r,θ) boundaries. In principle, one could impose Φdisk so as to describe the gravitational influence of massive bodies external to the computational domain. However, there is a single boundary distribution of Φdisk corresponding to a disk isolated in vacuum. This distribution is straightforward to obtain from the integral from of Poisson’s equation:

Φ ( r ) = G ρ ( r ) | r r | d 3 r , $$ \begin{aligned} \Phi \left(\boldsymbol{r}\right) = - G \int \frac{\rho \left(\boldsymbol{r^\prime }\right)}{\vert \boldsymbol{r} - \boldsymbol{r^\prime } \vert } \mathrm{d}^3 \boldsymbol{r^\prime }, \end{aligned} $$(A.1)

where the integral spans the entire computational domain. While this integral can easily be implemented, it can also become prohibitively expensive to evaluate on all boundary cells.

We followed the method described by Müller & Steinmetz (1995) to evaluate the integral Eq. (A.1) at a reduced computational cost for a moderate accuracy loss. The Green’s function 1/|r − r| is expanded onto spherical harmonics Y l m $ Y_l^m $ and the series is truncated at a finite order l ≤ L much smaller than the number of cells in φ. After storing some integration weights, the algorithm reduces to two steps per Poisson problem. First, we compute the low-order moments of the density distribution by projecting it onto the Y l m $ Y_l^m $. Then, we evaluate Eq. (A.1) as a finite sum over these low-order moments.

We used the recursion relations highlighted by Müller & Steinmetz (1995) to further reduce the number of operations involved. For subdomains of size nr × nθ × nφ, parallel communications amount to two reductions (sums) of 2nr × (L+1)2 values in the (θ,φ) dimensions and two sums of 2np × (L+1)2 values across the np subdomains in the radial dimension. At the heart of this method, the spherical harmonics were evaluated via the stable recursion relation provided by Press et al. (2007).

Although this algorithm makes it computationally affordable to evaluate Eq. (A.1) on the domain boundaries, the number of operations involved with the management of large arrays restrict the maximum order L of the harmonic expansion. Luckily, this expansion quickly converges with order l for smooth density distributions. As a compromise between speed and accuracy, we found that a truncation order L = 4 was sufficient to capture the main features of the gravitational potential on the domain boundaries.

A.3. Tests

We designed two experiments to test the correct implementation of the boundary conditions using spherical harmonics. In the first experiment, we meshed the spherical domain (r,θ,φ) ∈ [1,32] × [0,π] × [0,2π] with 512 × 32 × 64 grid cells, with a logarithmic spacing in r and uniform spacings in θ and φ. We uniformly distributed a mass M = 1 inside the spherical shell 2 ≤ r ≤ 4 and kept the density ρ = 10−16 outside. The gravitational potential of the gas should then be constant inside r ≤ 2 and vary as −GM/r outside r > 4. The integral Eq. (A.1) was evaluated on both radial boundaries while the θ and φ boundary conditions respected the spherical topology of the domain. Even though the monopole L = 0 completely describes this density distribution, we expanded the series to L = 4 to confirm that multipolar orders do not manifest themselves in this situation.

We obtained a spherically symmetric potential Φdisk(r) as expected. Outside r > 4, the potential followed the expected −GM/r to better than 5 × 10−3 accuracy all the way to the outer boundary. Inside 1 ≤ r ≤ 2 the gravitational potential varied by less than 2 × 10−4 of its boundary value. At the inner radius rin = 1 the potential matched the expected value Φdisk(rin) ≈ −0.321GM/rin to 4 × 10−4 accuracy. These residual errors mainly come from the mass shell not being exactly bounded by 2 ≤ r ≤ 4 on the radial grid. We can thus confirm that the boundary potential satisfies Poisson’s equation in its integral form Eq. (A.1), with the correct prefactors and spatial scalings.

In the second experiment, we tested the convergence of the boundary conditions with the truncation order L in an asymmetric configuration. We used the same domain and boundary conditions as in the first test, but we increased the resolution to 256 × 192 × 384 so the grid cells are roughly cubic. We distributed a mass M = 1 uniformly within a sphere of radius 1 and center x = y = z = 16 / 3 $ x=\mathit{y}=z=16/\sqrt{3} $ while keeping ρ = 10−16 outside. With L = 8 the two arrays of integration weights already amount to

2 complex × 256 × 192 × 384 grid × ( 8 + 1 ) 2 moments × 8 bytes double 24 Gb $$ \begin{aligned} \underbrace{2}_\text{complex}\times \underbrace{256 \times 192 \times 384}_\text{grid} \times \underbrace{\left(8+1\right)^2}_\text{moments} \times \underbrace{8 \,\mathrm{bytes} }_\text{double} \approx 24\mathrm{Gb} \end{aligned} $$

each – in total, to be subdivided per process – when stored in double precision.

We show in Fig. A.1 a slice of the sphere intersecting both the origin and the massive ball in this test with L = 8. The key point is that the iso-potential contours remain circular when crossing the radial boundaries. This property is accurately satisfied in the vicinity of the mass: the gravitational potential behaves as if there were no boundaries at all, as intended. Deviations from circular contours only appear in the lower-left corner of Fig. A.1. These are expected because the density distribution is all but smooth in this test, leading to a slow convergence of the harmonic series with L. These two tests confirm that the boundary conditions are properly implemented and that subsequent errors would primarily come from the finite truncation order L.

thumbnail Fig. A.1.

Gravitational potential (color-filled contours) of a massive ball inside a spherical domain. The domain extends radially over r ∈ [1,32], the massive ball has a radius of 1 and is offset from the origin by a cartesian displacement x = y = z = 16 / 3 $ x=\mathit{y}=z=16/\sqrt{3} $. This slice intersects the origin and is normal to the vector (0,−1,1) so as to intersect the mass. The boundary values of the potential are evaluated using spherical harmonics of order up to L = 8. The white circles are centered on the mass and correspond to theoretical iso-potential contours. The white disk r < 1 at the center is excluded from the computational domain.

Appendix B: Reaction of the central object

In this paper, we considered disks orbiting a point-like massive object located at the origin of the frame of reference. By assuming that this frame was inertial, we neglected the acceleration of the central object due to the attraction of massive inhomogeneities in the disk. To test the validity of this approximation, we ran four additional simulations including noninertial effects. We kept the same simulation setup as before apart from dividing the grid resolution by two in every dimension for the sake of computational time. At every timestep we computed the acceleration vector

a = i , j , k [ G ρ δ V r ] i , j , k r i , j , k 3 $$ \begin{aligned} \boldsymbol{a}_{\star } = \sum _{i,j,k} \frac{\left[ G \,\rho \,\delta V \,\boldsymbol{r} \right]_{i,j,k}}{r_{i,j,k}^3} \end{aligned} $$(B.1)

of the central object by summing over every grid cell indices (i,j,k) in the computational domain. For the central object to remain at the origin of the frame, the disk must feel a uniform and opposite acceleration −a which we included as a source term in the momentum equation. We fixed the dimensionless cooling time β = 10 and sampled disk masses M ∈ {1/5,1/3,1/2,1}.

Because mass can flow out of the computational domain, it is not possible to track the exact position of the barycenter over time. Asymmetric accretion bursts can cause a momentary acceleration of the central object in one direction, with no simple means to guarantee global momentum conservation after the gas has left the computational domain. We found that such asymmetric accretion bursts mainly occur during the onset of gravito-turbulence. As shown in Fig. B.1, the total mass losses amount to almost 25% over 300tin in the most massive M = 1 case. While the net momentum lost tends to vanish upon time averaging, it is generally different from zero and introduces a bias in the acceleration Eq. (B.1) which we cannot simply correct for.

thumbnail Fig. B.1.

Total mass contained inside the computational domain relative to the initial mass as a function of time for different initial disk masses (see legend), with the central object reacting to the disk gravity (solid lines) and without (dashed lines).

Once in a steady gravito-turbulent phase, the disk remains Keplerian and develops no substantial eccentricity. Compared to the inertial simulations, the diagnostics gathered in Table B.1 are overall compatible with those of Table 1 for a given disk mass. The slight differences found for the least massive M = 1/5 case likely come from the reduced grid resolution. The accretion histories drawn in Fig. B.1 are also similar after the onset of turbulence t ≳ 100tin. We conclude that the reaction of the central object does not significantly alter our diagnostics regarding the disk dynamics, in agreement with Michael & Durisen (2010). This conclusion is also supported by the simulations of Forgan et al. (2011), who found that m = 1 perturbations remain subdominant even for disk masses M = 3/2 of the central mass. The stability analysis of Heemskerk et al. (1992) predicts an m = 1 instability only when M ≳ 1, and it does not account for the effective viscosity of the gravito-turbulent disk, which may explain why our results with and without the reaction of the central object are compatible.

Table B.1.

Simulations including the reaction of the central object: initial disk mass, mean dimensionless stress, mean spiral pitch angle and standard deviation of the gas surface density.

By integrating the individual components of the acceleration vector a over time, one can reconstruct the hypothetical trajectory of the central object around the barycenter. Since the exact position of the barycenter cannot be known, this trajectory only serves to estimate the typical displacements of the central object. Reversing the sign of the acceleration vector, we equivalently obtain the typical displacements of the disk. We begin the integration at 100tin to discard the initial transient phase and show the resulting trajectories in Fig. B.2. All the trajectories remain confined inside r/rin < 1, so the disk is never significantly displaced through the inner radial boundary. As expected, the displacements are increasingly large for larger disk masses, reaching up to r ≈ 0.6rin in the M = 1 case. The trajectories of the two most massive cases M ≥ 1/2 roughly describe orbits but with no clear cyclic pattern, as one would expect if the disk supported long-lived m = 1 density perturbations.

thumbnail Fig. B.2.

Trajectory of the central object reacting to the disk gravity from 100tin to 300tin in the equatorial plane for different disk masses (see legend). The gray arcs represent the inner boundary of the computational domain at a radius r = 1. The inset zooms on the innermost regions for the two least massive disks.

All Tables

Table 1.

Simulation label, initial disk mass, cooling timescale, integration time, mean aspect ratio, mean dimensionless stress, mean spiral pitch angle and mean standard deviation of the gas surface density.

Table B.1.

Simulations including the reaction of the central object: initial disk mass, mean dimensionless stress, mean spiral pitch angle and standard deviation of the gas surface density.

All Figures

thumbnail Fig. 1.

Evolution in time (horizontal axis) of the radial profile (vertical axis) of the gas surface density relative to its initial profile in run M3B10. The dashed curves on the left fit the initial transient fronts with r(t) ∝ t2/3, i.e., a radial velocity proportional to r−1/2.

In the text
thumbnail Fig. 2.

Radial and azimuthal averages of the density (blue, left axis) and sound speed (red, right axis) relative to their azimuthally averaged midplane value at time 500tin in run M3B10.

In the text
thumbnail Fig. 3.

Gas surface density Σ at time 500tin in run M3B10; the color scale spans three orders of magnitude from black to white.

In the text
thumbnail Fig. 4.

Radial velocity relative to the azimuthally averaged sound speed in the midplane of run M3B10 at time 500tin. The velocity fluctuations reach supersonic amplitudes both outward (red) and inward (blue).

In the text
thumbnail Fig. 5.

Radial profiles of the time-averaged Toomre number Q in run M3B10, using either the density-weighted sound speed ⟨cs, isoρ (solid black) or a simple vertical averaging (dashed red) in Eq. (6). The dotted horizontal line marks the threshold Q = 1 of marginal stability.

In the text
thumbnail Fig. 6.

Time-averaged radial profiles of three key quantities to assess the disk stability: the dimensionless shear-rate (upper panel), the function ℒ introduced by Lovelace et al. (1999) and defined by Eq. (17), and the vertical component of the vortensity.

In the text
thumbnail Fig. 7.

Radial profiles of the normalized total stress (solid black) and gravitational stress alone (dashed red) after spatial and time averaging over 181tin in run M3B10. The dotted horizontal line marks the expected value αLTE assuming a Keplerian shear rate in Eq. (21).

In the text
thumbnail Fig. 8.

Radial profiles of the peak-to-peak range of the three velocity components in the midplane, relative to the azimuthally averaged midplane sound speed at time 500tin in run M3B10.

In the text
thumbnail Fig. 9.

Vertically averaged density fluctuations ρ $ \tilde{\rho} $ at time 500tin in run M3B10. The slope of the dashed lines refers to the estimated mean pitch angle tan(i) ≈ 0.25, or i ≈ 14°.

In the text
thumbnail Fig. 10.

Relative density fluctuations at four consecutive times separated by 0.5tin in the equatorial plane of run M3B10. The gray scale covers ρ [ 0 , 1 ] $ \tilde{\rho} \in \left[0,1\right] $. The green curves are transported at the mean orbital velocity ⟨vφρ of the gas and serve as corotation references. The blue and red contours delimit initially disjoint wakes in the process of merging by the third snapshot.

In the text
thumbnail Fig. 11.

Radial profiles of angular velocities relative the Keplerian frequency Ω ∝ r−3/2 after averaging over 7.5tin in run M3B10. The red area encompasses the bulk orbital velocity plus or minus the bulk sound speed. The thin blue curve is the spiral pattern speed estimated from the 2 ≤ m ≤ 6 Fourier modes. The thick black curve is the pattern speed estimated by maximizing the time correlation Eq. (24). The dashed green lines are contours of constant angular velocity.

In the text
thumbnail Fig. 12.

Radial velocity fluctuations in a frame corotating with the gas in run M3B10 at r = 5.2rin. The horizontal axis φ − Ωρt corresponds to an azimuthal coordinate in the comoving frame. The vertical axis is the density-weighted mean radial velocity ⟨vrρ relative to the average sound speed ⟨csρ, φ. Different curves correspond to different times as color-coded in units of local orbital periods. The horizontal segment illustrates how far a signal would propagate in the φ direction at the sound speed ⟨csρ, φ over this time interval.

In the text
thumbnail Fig. 13.

Instantaneous flow in a meridional slice at φ = π, centered on a spiral wake at R/rin ≈ 7.5 and time 500tin in run M3B10. The gas density is represented by filled contours spanning four orders of magnitude. The orientation of the meridional gas velocity is represented with colored arrows, blue (resp. green) being oriented away (resp. toward) the disk midplane. At the wake location, the midplane pressure scale height hmid/rin ≈ 0.3.

In the text
thumbnail Fig. 14.

Mean meridional cross-correlation between the density fluctuations ρ′ and the momentum fluctuations ρv′. The arrows give the orientation of the meridional components of the flow while the blue-red colors map its azimuthal component with an arbitrary scale. The right axis indicates the corresponding height measured from the midplane in units of mean midplane pressure scales.

In the text
thumbnail Fig. 15.

Same as Fig. 2 at time 500tin in the least massive case M10B10.

In the text
thumbnail Fig. 16.

Same as Fig. 9 in the least massive case M10B10. The dashed lines refer to a pitch angle of tan(i) ≈ 0.23, or i ≈ 13°.

In the text
thumbnail Fig. 17.

Same as Fig. 11 after averaging over 41 snapshots spanning 20tin in the most massive case M1B10.

In the text
thumbnail Fig. 18.

Radial profile of the vortensity in run M1B10. The maximum at r/rin ≈ 18 coincides with the corotation radius in Fig. 17.

In the text
thumbnail Fig. A.1.

Gravitational potential (color-filled contours) of a massive ball inside a spherical domain. The domain extends radially over r ∈ [1,32], the massive ball has a radius of 1 and is offset from the origin by a cartesian displacement x = y = z = 16 / 3 $ x=\mathit{y}=z=16/\sqrt{3} $. This slice intersects the origin and is normal to the vector (0,−1,1) so as to intersect the mass. The boundary values of the potential are evaluated using spherical harmonics of order up to L = 8. The white circles are centered on the mass and correspond to theoretical iso-potential contours. The white disk r < 1 at the center is excluded from the computational domain.

In the text
thumbnail Fig. B.1.

Total mass contained inside the computational domain relative to the initial mass as a function of time for different initial disk masses (see legend), with the central object reacting to the disk gravity (solid lines) and without (dashed lines).

In the text
thumbnail Fig. B.2.

Trajectory of the central object reacting to the disk gravity from 100tin to 300tin in the equatorial plane for different disk masses (see legend). The gray arcs represent the inner boundary of the computational domain at a radius r = 1. The inset zooms on the innermost regions for the two least massive disks.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.