Convective overstability in radially global protoplanetary disks I -- Pure gas dynamics

Protoplanetary disks are prone to several hydrodynamic instabilities. One candidate, Convective Overstability (COS), can drive radial semi-convection that may influence dust dynamics and planetesimal formation. However, the COS has primarily been studied in local models. This paper investigates the COS near the mid-plane of radially global disk models. We first conduct a global linear stability analysis, which shows that linear COS modes exist only radially inward of their Lindblad resonance (LR). The fastest-growing modes have LRs near the inner radial domain boundary with effective radial wavelengths that can be a substantial fraction of the disk radius. We then perform axisymmetric global simulations and find that the COS's nonlinear saturation is similar to previous incompressible shearing box simulations. In particular, we observe the onset of persistent zonal and elevator flows for sufficiently steep radial entropy gradients. In full 3D, non-axisymmetric global simulations, we find the COS produces large-scale, long-lived vortices, which induce outward radial transport of angular momentum via the excitation of spiral density waves. The corresponding $\alpha$-viscosity values of order $10^{-3}$ agree well with those found in previous 3D compressible shearing box simulations. However, in global disks, significant modifications to their radial structure are found, including the formation of pressure bumps. Interestingly, the COS typically generates an outward radial mass transport, i.e. decretion. We briefly discuss the possible implications of our results for planetesimal formation and for interpreting dust rings and asymmetries observed in protoplanetary disks.


INTRODUCTION
The proposed formation of km-sized planetesimals in protoplanetary disks (PPDs) via the core accretion model (Safronov 1972) faces several obstacles, such as the so-called bouncing and drift barriers, inhibiting efficient coagulation of dust grains (e.g.Blum 2018;Drazkowska et al. 2022).In order to overcome these, mechanisms to effectively concentrate dust grains in the mid-plane of PPDs are required.The most popular such mechanism is the streaming instability (SI: Youdin & Goodman 2005;Youdin & Johansen 2007;Johansen & Youdin 2007), which in its nonlinear saturation can produce dense clumps of dust particles, which can subsequently collapse under the action of selfgravity to form planetesimals (Johansen & Youdin 2007; Bai & Stone mlehmann@asiaa.sinica.edu.tw2010; Simon et al. 2016;Carrera et al. 2021).However, for small dust grains the clumping by the SI is strong enough only if the metallicity (or dust abundance) takes at least a few times the solar value (Li & Youdin 2021).
In the recent years, the possible occurrence of purely hydrodynamic instabilities in PPDs has received an increased amount of attention.The main reason is that the magneto-rotational instability (MRI: Balbus & Hawley 1991), which had long been considered to be the primary source for mass accretion in PPDs, is likely weak or absent in large parts of PPDs, due to low ionization rates (Gammie 1996;Turner & Drake 2009;Armitage 2011;Turner et al. 2014).While purely hydrodynamic instabilities, most prominently the vertical shear instability (VSI: Urpin & Brandenburg 1998;Urpin 2003;Nelson et al. 2013;Barker & Latter 2015;Lin & Youdin 2015), the convective overstability (COS: Klahr & Hubbard 2014;Lyra 2014;Latter 2016) and its nonlinear counterpart -baroclinic vortex amplification (Petersen et al. 2007a,b;Lesur & Papaloizou 2010;Lyra & Klahr 2011;Raettig et al. 2013) are likely too weak to explain global mass accretion rates in PPDs, they can strongly affect the concentration of dust grains, and can thus be active elements in the planetesimal formation process (Lesur et al. 2023) In this paper, we focus on the COS.Several studies have investigated this instability in terms of its linear growth phase and its nonlinear saturation, which includes vortex formation (Lyra 2014;Klahr & Hubbard 2014;Latter 2016;Teed & Latter 2021).In addition, the concentration of dust in COS-vortices was demonstrated in local shearing box simulations of gas including Lagrangian dust particles (Lyra et al. 2018;Raettig et al. 2021).The results of Raettig et al. (2021) suggest that vortices induced by the COS are suitable sites for planetesimal formation, as their simulations revealed dust densities exceeding the Roche density, depending on the disk radius under consideration.
However, the aforementioned studies are based on a fully local shearing box framework, where the radial buoyancy that powers the COS is treated as a constant forcing.This means that the possible effects of the ensuing instability on the disk structure are being ignored.The exception is the simulations of Klahr & Hubbard (2014), which were conducted in a global cylindrical setup, and more recently the simulations of Klahr et al. (2023), which are in principle fully global.However, the former simulations merely probed the linear growth of the instability, and, more importantly, they were 2Daxisymmetric, rather than 3D.The latter simulations were axisymmetric and adopted small radial and vertical extents in favor of resolution.
The goal of this paper is to study the COS in a radially global disk setup, by conducting both linear analyses and nonlinear hydrodynamic simulations (2D axisymmetric and 3D).In contrast to local shearing box simulations, global simulations self-consistently consider the disk's background structure, and possible changes thereof in response to the nonlinear evolution of instabilities and their entailing angular momentum transport.This is expected to be particularly relevant in 3D, where angular momentum can be transport radially due to non-axisymmetric structures.Also the concentration of dust grains should be affected by the global disk evolution.Furthermore, the most unstable linear COS modes found in fully local calculations actually possess a vanishing radial wavenumber (i.e. an infinite radial wavelength), which somewhat contradicts the local approximation in the first place.Radially global linear calculations presented below indeed confirm that the most unstable modes have a radially global extent, but these modes are not captured by a fully local model.
The paper is organised as follows.In Section 2 we introduce our hydrodynamic model of a PPD.In Section 3 we summarize the theoretical background of the COS, which will be useful in interpreting our results, and describe its potential relevance to disk observations.In Section 4, we perform a radially global linear analysis of the COS.In Section 5, we describe the setup and results of our hydrodynamic simulations.In Section 6, we summarise and discuss caveats of our work and prospects for future work.

Basic equations
We consider a global hydrodynamic model of a PPD consisting of a compressible, non-isothermal gas, governed by the set of fluid equations In these equations ρ, v and P are the gas volume mass density, the three-dimensional velocity, and pressure, respectively, at time t.We adopt a non-rotating frame with cylindrical coordinates (r, φ, z) and with origin on a central star of mass M * with gravitational potential Φ * .In what follows, we neglect the vertical variation of the potential and set Φ * = −GM * /r, where G is the gravitational constant, which implies that we consider a region close to the disk mid-plane such that |z| ≪ H, where H denotes the pressure scale height defined below.Furthermore, the indirect gravitational term, selfgravity, and magnetic fields are neglected.The remaining symbols in the above equations are explained below.
We assume an ideal gas equation of state with pressure where R is the gas constant and µ the mean molecular weight.Following Klahr & Hubbard (2014), we assume gas cooling to occur in the optically thin regime and define the cooling function where δT is the deviation of temperature from its equilibrium profile specified below.We subsume viscous heating into this cooling function.In what follows we work with the dimensionless constant cooling time where is the local Keplerian frequency.We fix the adiabatic index γ = 1.4.Finally, is the viscous stress tensor with the constant kinematic viscosity ν.The symbol † denotes the conjugate transpose and I stands for the unit tensor.In our nonlinear simulations, viscous terms are only included to ensure numerical stability, such that ν is chosen to be very small.In particular, ν is much smaller than the typical turbulent viscosities (defined in Section 5.1) measured in 3D simulations.Moreover, below, we perform linear calculations, including a kinematic viscosity, to describe the effect of external turbulence on the linear COS.

Disk equilibrium
We focus on the region close to the disk mid-plane z = 0 and neglect vertical gravity, so the equilibrium configuration does not depend on height z away from the mid-plane.In particular, vertical shear is absent, eliminating the possible occurrence of the VSI.We assume an axisymmetric disk equilibrium with radial power-law profiles for the gas volume mass density and temperature with constant powerlaw indices p and q, and reference density ρ 0 = ρ(r 0 ) and temperature T 0 = T (r 0 ) at the reference radius r 0 .
Although the disk model is unstratified, we can define the characteristic disk thickness or pressure scaleheight of a corresponding vertically stratified model via where is the sound speed.
The equilibrium velocity of the gas is given by with the orbital frequency where we defined the dimensionless radial pressure gradient (Youdin & Goodman 2005) with the disk aspect ratio h = H/r.We also define Ω 0 = Ω(r 0 ) and H 0 = H(r 0 ).Furthermore, the squared radial buoyancy frequency is given by (16) with the dimensionless entropy The second equality in Equations ( 15)-( 16) assumes the radial power-law profiles defined above.

Theoretical aspects
The COS comprises an oscillatory destabilization of inertial waves of typical frequency ∼ Ω, induced by a convectively unstable radial entropy gradient, where epicyclic oscillations are amplified on account of a radial buoyancy force which results from the cooling of gas occurring on a typical time scale t c ∼ 1/Ω (see Latter 2016 for an illustrative explanation and Teed & Latter 2021 for more details).Under typical conditions where the radial pressure gradient ∂ r P < 0, the COS requires a negative radial entropy gradient, ∂ r S < 0, since then the squared radial buoyancy frequency (16), which is the central quantity for the linear COS, is negative.The linear instability mechanism behind the COS was first studied by Klahr & Hubbard (2014), Lyra (2014), andLatter (2016).
The nonlinear saturation of the COS was first explored in compressible 3D shearing box simulations by Lyra (2014).These simulations revealed the formation of large scale vortices and an entailing radial angular momentum transport characterized by α-viscous values on the order ∼ 10 −3 .The nonlinear saturation of the axisymmetric COS was studied in detail by Teed & Latter (2021) (TL21 henceforth) employing local incompressible shearing box simulations.These authors found that two dimensionless quantities largely determine the nonlinear outcome of the instability.These are the Reynolds number measuring the strength of viscous diffusion, and the (pseudo-)Richardson number measuring the strength of the destabilizing radial entropy stratification against the stabilizing effect of the radial angular momentum gradient.Here, L is a characteristic length scale of the flow and κ denotes the epicyclic frequency.Since our disk model is (as in TL21) vertically unstratified, this length scale should be the vertical extent of the simulation region.Our simulations below adopt a vertical domain L z = 0.5H (see §5), but for simplicity we set L ∼ H, and we use κ ≈ Ω K .Assuming this and considering the power-law disks adopted above gives the second equalities in Eqs. ( 18)-( 19).Note that in contrast to TL21, our model is radially global, such that Re and R are functions of disk radius.Specifically, R (Re) increases (decreases) with decreasing radius, implying a more vigorous COS at smaller radii.TL21 found that, depending on the values of R and Re, the COS saturates either into a weakly nonlinear state characterised by ordered nonlinear waves, a state characterised by wave turbulence, or a state supporting intermittent or persistent zonal flows and elevator flows (see §5.3).They found that R needs to exceed a critical value, roughly delineated by the relation1 in order for zonal flows to be formed in the saturated state.The point is that the COS is only expected to play a role in dynamical process in PPDs if zonal flows are formed, which for instance can concentrate dust grains.Furthermore, in 3D simulations zonal flows are expected to wrap up into vortices, which can induce radial angular momentum transport, and, which are likewise able to concentrate dust (Raettig et al. 2021).

Occurrence in PPDs
The requirement of N2 < 0 for the COS to operate in the disk mid-plane typically translates to a rather steep temperature gradient and a flat density profile (see Figure 1), which would be inconsistent with global disk pro-files inferred from early mm and sub-mm continuum observations of PPDs (e.g.Andrews et al. 2009;Isella et al. 2009).However, flat density profiles could be realized, for instance, near the edges of dead zones, planet-carved gaps, or pressure bumps in general.As such, COS-active regions would be expected to have typical radial widths ∆r ∼ H. On the other hand, more recent ALMA imaging of HD 163296 (Mathews et al. 2013) and several disks in Lupus (Tazzari et al. 2017) revealed very flat surface mass density profiles at disk radii around tens of AU, which should in principle increase the likelihood for the occurrence of the COS.
Whether the COS can operate also strongly depends on the gas' cooling timescale and is thus tied to the distribution of small dust grains tightly coupled to the gas (Malygin et al. 2017;Barranco et al. 2018;Pfeil & Klahr 2019;Fukuhara et al. 2021).The above studies attempted to map out the gas cooling time scale across the radial and vertical extent of a PPD.However, such calculations require multiple simplifications, such that the computed cooling times are subject to large uncertainties.Nevertheless, these estimates suggest that one or more purely hydrodynamic instabilities are likely active in planet-forming regions of the disk.In particular, the results of Pfeil & Klahr (2019) suggest that the COS may occur in the potential planet-forming region at r = 1 − 10AU if the structural requirement N 2 < 0 is fulfilled.

Linearised Equations
We now assume the disk equilibrium described in §2.2 and investigate the stability of axisymmetric radially global, vertically local perturbations of the form and similar definitions for the perturbations δv r , δv φ , δv z and δP .We define the complex mode frequency ω = ω R + iω I , such that ω R and ω I denote the realvalued mode's frequency and growth rate, respectively.Furthermore, k z is the real-valued vertical wavenumber and we take k z > 0 without loss of generality.We also define for convenience./ΩK(r0) 2 within a radial region r0 ± 2H0 in a power law disk with density and temperature slopes p and q, respectively.The diagonal lines delineate the regions where min N 2 takes negative or positive values, such that the COS is in principle expected to occur in the parameter-space labeled "COS".Open (solid) circles represent values of p and q used in our 2D-simulations (2D as well as 3D simulations) below.
The linearised cylindrical equations then read Note that the driving force of the COS is the radial buoyancy force, encapsulated in the term ∝ ∂ r P in the radial momentum equation.

WKB dispersion relation
We can easily derive the WKB dispersion relation corresponding to Eqs. ( 24)-( 28) by ignoring all radial variations of the equilibrium disk structure, as well as curvature terms.Furthermore, we introduce the complex radial WKB wavenumber k r (r), assuming that ∂ r → ik r (r).We then obtain the fifth order dispersion relation In the isothermal and adiabatic limits t c → 0 and t c → ∞ we find the fourth order dispersion relation and respectively.Note that (29) should not yield any growing solutions, which can easily be shown in the isothermal and adiabatic limits, since in those cases the dispersion relation yields a simple quadratic equation in ω 2 .Furthermore, in the limit of high mode frequencies |ω| ≫ Ω we isolate acoustic waves in the adiabatic limit, where k = k 2 r + k 2 z .The isothermal result can then be obtained by setting γ to unity.
On the other hand, in the incompressible limit c s → ∞ we find The solution of the latter is the cooling mode decaying on cooling time scale t c , as well as the pair of inertial waves Since the axisymmetric COS is captured within the incompressible Boussinesq approximation (Latter 2016;Latter & Papaloizou 2017;TL21;Lehmann & Lin 2023), we expect that short wavelength COS modes in our radially global model essentially are inertial waves, described locally by Eq. ( 35).This will indeed be confirmed below.The Boussinesq approximation can be obtained from Eqs. (1)-(3) by assuming that the length scale of the phenomena under investigation is much shorter than the radial and vertical disk stratification lengths, that the flow is strongly subsonic (|δv| ≪ c s0 ), and that the fractional perturbations of density δρ/ρ ≪ 1 and pressure δP/P ≪ 1, but with δρ/ρ ≫ δP/P .The latter assumption is necessary to retain the effect of the background entropy gradient in the model equations (see Latter & Papaloizou 2017 for details).

Numerical methods
When supplemented with two boundary conditions accounting for the two first order derivatives in r, such as Eqs. ( 24)-( 28) constitute a linear eigenvalue problem of the form Lb = iωb (37) with the linear differential operator L, eigenvector b = {Q, δv r , δv φ , δv z , W } T and complex eigenvalue iω.We use the spectral code Dedalus (Burns et al. 2020) to solve the linear eigenvalue problem (37) on a Chebyshev grid using a collocation point method (see, for instance, Lin 2021 for a similar treatment of a vertically stratified disk model).
Complementary to this, we verify the results obtained with the spectral code using an alternative algorithm based on the shooting method.To this end, we first recast the above equations as two first-order ordinary differential equations in r: which are then integrated using a fourth or fifth-order Runge-Kutta scheme with pre-specified outer (we shoot from the outer to the inner domain boundary) boundary values for δv r and W .The solution with the correct value for ω is then obtained via a Newton-Raphson iteration algorithm to match the corresponding inner two boundary values for one of the fields (we use δv z ), obtained first with the spectral method.We find that both methods generally are in excellent agreement.However, in certain cases modes possess very short radial wavelengths (k r H 0 ≳ 100), such that the spectral method becomes prohibitively slow to obtain fully resolved radial profiles of the eigenfunctions.Nevertheless, the method can still find the correct eigenvalue ω corresponding to such solutions, which we verified using example calculations with subsequently increased radial resolution.
For such cases, we present radial profiles obtained with the shooting method, which can be computed efficiently even for a large number of radial grid points n r ≳ 10 4 .

Smooth disk region
Figure 2 shows frequencies ω R and growth rates ω I of linear COS modes for vertical wavenumbers 0.2 ≤ k z H 0 ≤ 1, 000, obtained with the spectral code.Here we consider a radial domain 5 ≤ r/H 0 ≤ 15.As outlined in §1, typical regions that are unstable to the COS are expected to be significantly smaller.However, in our simulations below we will adopt the same radial domain, spanning 10H 0 .The reason is that we wish to avoid spurious influences from the radial boundary conditions on the long term nonlinear evolution in our simulations.Moreover, the linear analysis of a smooth disk region does not depend qualitatively on the radial domain size (cf.§4.4.2).
For each value of k z we generally find a band of modes within a range of oscillation frequencies and growth rates, with the average frequencies and growth rates increasing with increasing k z .Moreover, larger growth rates generally correspond to larger oscillation frequencies.Overall, oscillation frequencies ranging from small values ω R ∼ 10 −2 Ω 0 to roughly the largest value of the Keplerian frequency within the considered domain ω R ≈ 2.7Ω 0 (i.e.typical frequencies of local inertial waves with k z ≲ k) can occur.Note that for smaller k z the lowest occurring frequency is set by the radial grid resolution of the spectral solver.Furthermore, modes with lower frequency generally extend further out through the domain, whereas with increasing frequency (and hence increasing growth rate) modes tend to be increasingly concentrated toward the inner domain boundary.This is illustrated in Figure 3, which shows Figure 3. Radial eigenfunctions corresponding to the fastest growing COS modes in the fiducial disk with p = 1.5 and q = 2 for increasing vertical wavenumbers kzH0 = 0.2, 2, 20, as obtained with the spectral method.Notably, the eigenfunctions become increasingly concentrated toward the inner domain boundary with increasing kz.The bottom panel shows the wavelet power corresponding to the radial profile of Re [δvz].The red dashed curves in the wavelet power spectra correspond to kr obtained from (35), with accordingly specified kz and using ωIW = ωR.For the modes with kzH0 = 2, 20 the Lindblad resonance defined via ( 41) is located within the computational domain, where the red dashed curves drop to zero.
the radial profiles of the various hydrodynamic quantities involved in the problem.
The behavior of the modes just described can be better understood if we consider the WKB dispersion relation in the isothermal or adiabatic limit (30)-( 31).The latter can be rearranged to yield This equation is analogous to Eq. ( 60) of Lubow & Pringle (1993) and Eq. ( 12) of Svanberg et al. (2022), the latter applying to vertically stratified, radially local disk models.According to Eq. ( 40) neutral 2 waves can only propagate in regions where the two conditions ω > (<)Ω and ω > (<) √ γc s k z are simultaneously fulfilled.The first (second) equality applies to acoustic (inertial) waves.
Since our linear analysis above reveals no modes with ω > Ω, we conclude that COS modes only exist radially 2 with real-valued kr and ω inwards from their Lindblad resonances (LRs), whose radii r LR are given by the condition This appears plausible, since in the Boussinesq approximation the linear COS is known to excite inertial waves (Klahr & Hubbard 2014;Lyra 2014;Latter 2016), and presumably there exist no linear (global) modes in our disk model which connect both inertial and acoustic waves.
Nevertheless, modes that can appear in our simulations are restricted to vertical wavenumbers k z H 0 ≳ 10, such that the modes in our simulations can only fulfill the condition ω < Ω and waves can only propagate inside of their Lindblad resonance radius.Note that generally, our vertically unstratified analysis is only meaningful for modes with k z H 0 > 1.
The above implies that lower frequency modes have their LR at larger radii, such that these modes possess smaller radial wavelengths within the computational domain.Consequently, their growth rates are smaller, as their coupling to the background entropy gradient becomes less efficient.
In addition, we compute a wavelet transform3 of the vertical velocity δv z , which enables us to estimate the radial wavenumber k r (r) of the eigenmodes.The bottom panel in Figure 3 shows the wavelet power spectrum as a function of radius and radial wavenumber k r H 0 .The mode with k z H 0 = 0.2 shows a radial change of its radial wavenumber that is characteristic of global waves, such as spiral density waves in planetary rings (Shu 1984).This radial change results from the radially varying Ω K .Note that the eigenfunctions for the cases k z H 0 = 2 and k z H 0 = 20 are actually not wavelike, i.e. they do not undergo at least one full oscillation in radial direction.Therefore, the wavelet transforms for these cases cannot be expected to yield a clear radial wavenumber.To improve the results we re-scaled the wavelet power such that its minimal value is set to zero (by manual subtraction) for these cases, as this procedure eliminates artificial residuals in the wavelet power spectrum.The red dashed curve corresponds to the radial wavelength of an inertial wave (35) with given vertical wavenumber k z and frequency ω R .
Finally, we note that we checked the robustness of the computed modes by using a number of different boundary conditions.Also we verified the existence of some (including the fastest growing modes) with the shooting method.

Connection to local COS modes
We can attempt to compare the linear growth rates of the global COS modes derived in the previous section with those of fully local modes.To this end, we consider the expression for the growth rates of local modes with wavenumber k = k r e r + k z e z (Klahr & Hubbard 2014;Lyra 2014;Lehmann & Lin 2023): where This expression has been derived using the incompressible Boussinesq approximation (Latter & Papaloizou 2017), which yields almost identical results to a local compressible calculation (see Appendix A).We expect that the global COS modes derived in the previous section increasingly resemble local COS modes as we decrease the radial domain size, since then the variation of background quantities becomes smaller.This is illustrated in Figure 4, where we present growth rates of global modes with k z H 0 = 10 for three different domain sizes ∆r = H 0 , ∆r = 5H 0 and ∆r = 10H 0 , as a function of radial wavenumber k r , represented by the solid lines.These are compared with local modes following from (42), represented by open circles.The values of k r for each mode are obtained using the WKB dispersion relation ( 29), where we assume values for c 2 s , t c and Ω corresponding to radial averages over the entire domain.Similarly, radial averages of N 2 are used in (42).Note that generally the WKB wavenumber k r is complex, and the values in Figure 4 are the corresponding real parts.
Thus, for domain sizes ∆r ≲ H 0 , short wavelength modes with k r H 0 ≫ 1 resemble local COS modes, as indicated by the good match in the corresponding growth rates.However, the fastest growing modes are still global, as they do not extent throughout the entire domain unlike local modes, but are concentrated toward the inner domain boundary.This is illustrated in the insert figure, which shows the corresponding vertical velocity field for ∆r = H 0 and ∆r = 5H 0 .The reason is that the fastest growing modes' LR lies within the computational domain, implying that their frequencies ω R ≳ Ω 0 .Again, we checked these modes' robustness by using several different boundary conditions.With increasing domain size the radial variation of the background quantities (particularly Ω K ) results in an increasing discrepancy between local and global growth rates, as k r exhibits an increasing variation throughout the domain.Obtaining local growth rates becomes less meaningful for larger domain sizes.

Dependence on disk parameters
For completeness, Figure 5 shows growth rates of the fastest-growing COS modes as a function of various disk parameters.The left panel compares different values of an included α-viscous stress 4 (see Appendix B).Here, we also indicate the region that is expected to be resolved in our hydrodynamic simulations below ( §5).The lower boundary is set by the vertical size of the simulation domain (∆z = 0.5H 0 ), whereas the upper boundary is set by the grid resolution, assuming that 10 grid points per wavelength are required to resolve a mode.Thus, based on the linear results, we expect the COS in our simulations to be extinguished for α ≳ 10 −4 , which we will find in reasonable agreement with the simulation results presented below.The behavior of the curves in the remaining panels largely agrees with the fully local calculations presented in previous works.That is, the growth rates in the three rightmost panels are well explained by Eq. (42) using Eq. ( 16) for the buoyancy frequency (its radial average).For reference, the open circles in the second panel represent fully local growth rates computed using the method explained in §4.4.2 for the fiducial case.It should be noted again that modes 4 The indicated values of α correspond to a disk radius r = r 0 .
k z H 0 ≲ 1 are likely to be affected by the neglect of vertical disk stratification.

Disk region containing a radial density structure
It is interesting to verify if the COS can operate in a more "realistic" disk model.By this we mean a model of a disk region with temperature and density slopes compatible with observational constraints, but containing a "special" structure that locally gives rise to the necessary temperature and density gradients to yield a negative value of N 2 across a certain radial domain.To this end we consider an initially convectively stable disk with p = 2 and q = 0.75 (cf. Figure 1) on which we superimpose either a density bump or a density step to render the region unstable to the COS.This is illustrated in Figure 6 (left panel) for an amplitude A = 0.4.The upper panel displays the effective density slope p eff = −∂ log ρ/∂ log r, whereas the lower panel shows N 2 for the three cases of a smooth disk, and a disk containing a density bump and a density step, respectively.A linear analysis adopting the latter two background disk structures yields growing modes for k z H 0 ≳ 5.For increasing k z bands of growing modes are found which become increasingly narrow and concentrated around the largest Keplerian frequencies within the unstable region.The fastest growing modes have their LR close to the outer boundary of the unstable region.However, note that even for large k z H 0 ≫ 1 the fastest growing modes possess growth rates which are substantially smaller than in the smooth disk case discussed above.It is interesting to note that in the present case most5 of the growing linear modes in this problem resemble to some extent, resonantly forced (spiral) density waves (e.g., Shu 1984) since they similarly appear to be excited within a narrow, unstable region and propagate away from it with radially decreasing wavelength.As in the case of the smooth disk, here we also find only modes radially inward of their Lindblad resonances.In Figure 7 we show the fastest growing mode with k z H 0 = 20.3) forward in time.As our disk is vertically unstratified, we apply periodic boundary conditions to all quantities in the azimuthal and vertical directions.Radial boundary conditions are such that equilibrium values of gas density and azimuthal velocity are extrapolated into the ghost zones, while gas radial and vertical velocities are set to zero at the boundaries.Simulations are carried out in cylindrical coordinates as defined in §2.Computational units are such that r 0 = M * = G = 1.We adopt h 0 = H 0 /r 0 = 0.1 in all runs.For our fiducial setup the numerical grid covers 0.5 ≤ r/r 0 ≤ 1.5 and −0.25H 0 ≤ z ≤ +0.25H 0 , and in 3D simulations 0 ≤ φ ≤ 2π.The grid resolution is n r × n z × n ϕ = 2000 × 100 × 628.Thus, the radial and vertical resolutions are 200/H 0 , whereas the azimuthal resolution ∼ 10/H 0 is significantly lower for reasons of computational costs 7 , and since we expect structures to form in our simulations to be predominantly axisymmetric or elongated in azimuthal direction.3D simulations are carried out on a GPU cluster, which is necessary due to the high spatial resolution.All 3D simulations cover 1,000 reference orbits.

Diagnostics
Following previous studies (e.g.Lehmann & Lin 2022), we describe the radial and vertical turbulent angular momentum transport in 3D simulations via the dimensionless quantities and respectively, where the brackets denote averaging over spatial dimensions as indicated in the subscript, and where δv φ is the azimuthal velocity deviation from its ground state value.Furthermore, we define the deviation ∆X of a quantity X(r, φ, z) from its vertical-azimuthal average across 7 The typical size of a single 3D simulation is ≈ 600Gb.(47) If not stated otherwise, our radial diagnostic domain is defined by r min /r 0 = 0.8 and r max /r 0 = 1.2.

Results of 2D axisymmetric simulations
In this section we present the results of 2D axisymmetric (radial-vertical) simulations of the nonlinear saturation of the COS.In §5.3.1 we consider simulations of a smooth power-law disk with N 2 < 0 for all radii.Next, in §5.3.2 we consider a disk containing a pressure bump which locally gives rise to the COS in an other-wise stable disk.The different values of p and q adopted throughout our simulations are shown in Figure 1 for illustration.

Nonlinear saturation of the COS in a smooth disk
Figure 8 illustrates the nonlinear saturation of the COS in our fiducial axisymmetric simulation.The upper left frame shows the evolution of initially small scale, quasi-local modes into radially global modes.The latter eventually develop persistent flow structures.The upper right frame shows the evolution of the power spectral density of the radial velocity field (contours), averaged over the radial diagnostic domain 0.8 ≤ r/r 0 ≤ 1.2, as a function of vertical wavelength λ z .Generally, we find that the prevalent vertical length scale of flow structures increases with time.Note that in the saturated state the prevalent vertical wavelength of the radial velocity, which should represent secondary parasitic modes (Latter 2016; TL21) is L z /2.This is consistent with our numerical convergence analysis, presented in Appendix C, which shows that smaller vertical domain sizes than our fiducial one result in a weakened nonlinear saturation.The over-plotted solid curves are the rms vertical (green) and radial (cyan) velocities.As explained in TL21, the radial (vertical) rms velocities represent the power in the primary COS (secondary parasitic) modes.The onset of parasitic modes eventually controls the nonlinear saturation amplitude of the COS.
Furthermore, the circles represent the turbulent radial angular momentum flux α r computed via (45), such that open (solid) circles denote negative (positive) values.The dashed line denotes its negative time-averaged (over the last 400 orbits) value ⟨α r ⟩ t ≈ −1.3 • 10 −6 , and is roughly in agreement with corresponding values reported by TL21 (their Figure 17) for R ≈ 0.05.Note that since our scaling length is H, rather than L z (as in TL21), our value is to be multiplied by a factor of 4 for comparison.
The bottom frames in Figure 8 show snapshots of the vertical, azimuthal and radial velocity field components, respectively, in the saturated state of the simulation.The left and middle panels show the presence of persistent elevator flows and zonal flows, respectively, which constitute radial patches of alternating deviations of the corresponding velocity components (see TL21 for details on these structures).Zonal flows are accompanied by significant undulations in the radial pressure gradient, the latter being required to establish geostrophic balance in these flow structures.This is verified by the over-plotted curves of 2Ωδv φ and 1/ρ ∂ r δP , which are practically identical, thereby demonstrating geostrophic balance between the Coriolis force and the radial pres- sure gradient.Note that here, pressure variations δP are caused mainly due to variations in the density, as temperature variations are found to be subdominant.
In Figure 9 we plot rms vertical velocities illustrating the dependence of the nonlinear saturation of the COS on various physical simulation parameters.These are the cooling time β (upper left panel), the input viscosity α defined via Eq.( 18) (upper right panel), the initial radial density slope p (lower left panel) and the initial radial temperature slope q (lower right panel).In all panels of Figure 9 the dashed curve represents our fiducial simulation.Note the different simulation times in the left and right panels .The fiducial values of β and p are chosen so as to maximize linear growth rates and the saturation level of the COS reached within 1, 000 orbits.The adopted viscosity α = 10 −7 is roughly the largest possible value which does not notably affect the evolution of the instability.All simulations are initialized with low amplitude white noise ∼ 10 −3 c s0 in all velocity components.For α = 10 −7 we have R crit = 0.0032 by Eq. ( 20).The simulations with varying density slope and temperate slope which develop COS8 yield values in the range 0.037 ≲ R ≲ 0.085 and 0.008 ≲ R ≲ 0.085 within the diagnostic domain, respectively, such that in all of these simulations persistent zonal flows are expected to develop, at least after sufficient simulation time (see §3).This is found to be largely the case, except for simulations with large cooling times β ≳ 100.We suspect that in the latter simulations initially another double diffusive instability develops, making use of viscous diffusion (see Appendix D for details).However, after roughly 2,000 orbits we find the emergence of elevator flows, indicating that the double diffusive instability eventually gets supplanted by the COS.In the simulations presented here the former instability is weak due to the small values of α adopted.We hypothesize that after sufficient integration time the COS eventually saturates to similar levels in all simulations shown in the upper left panel.
We also performed additional simulations using wave damping boundary conditions (not shown).In these simulations we find a very similar nonlinear saturation, albeit slightly delayed (by 100-200 orbits) compared to the simulations shown in Figure 9.However, we find that the COS is almost entirely suppressed (for at least 10,000 orbits) in simulations using damping boundary conditions with p = 4.125 and p = −1.125,and also in simulations with cooling times β ≳ 10 (with final rms vertical velocity ∼ 10 −6 which keeps decreasing).This is an interesting result, as it suggests that the COS developing in our simulations is indeed a radially global, rather than a local instability.

Nonlinear saturation of the COS in the vicinity of a pressure bump
In §4.4.4 we conducted a linear analysis of a disk region containing density structures which give rise to the COS via a local modification of N 2 (r) in an otherwise stable disk, corresponding to the left most circle in Figure 1. Figure 10 illustrates the nonlinear evolution of the COS in the vicinity of a pressure bump (the same as in Figure 6).In this simulation, we seeded the linear eigenfunction with k z H 0 = 20 (Figure 7) with a small initial amplitude.
A similar outcome of the nonlinear saturation is obtained by seeding low amplitude white noise as done before.However, seeding the fastest growing mode significantly speeds up the evolution.Moreover, we can test its growth rate against the linear analysis presented above.The vertical domain for this simulation is set to be ∆z = π/5H 0 , such that the seeded mode fits exactly two times in the vertical domain.The black curve in the upper left frame displays the rms vertical velocity, with a linear growth phase of about 1,000 orbits, and the subsequent nonlinear saturation.The dashed green line represents the linear growth rate as predicted by our linear analysis, and which is indicated also in Figure 7.The insert figure shows the initial and final radial profile of the vertical velocity (at z = 0).Thus, even though only the region 0.7 ≲ r/r 0 ≲ 1 is unstable to the COS, excited waves penetrate significantly into the stable region at smaller radii.This is further illustrated in the upper right frames, which show contour plots of the vertical and azimuthal velocity perturbations after 30,000 orbits.The over-plotted curves are the buoyancy and pressure profiles N 2 and P , respectively.

Results of 3D simulations
We now shift our focus to our 3D simulations.Given the significant computational resources needed, we present a concise qualitative summary instead of a comprehensive parameter exploration.In particular, we perform 3D simulations of smooth disks ( §5.4.2) where we consider three different values of the cooling time In the simulations adopting different p we slightly adjust the temperature slope q, such that the radial profile of N 2 is the same in all simulations.Thus, essentially q ≈ 2 in these simulations.In addition, we consider one semi-3D simulation in the vicinity of a pressure bump, which we discuss first in §5.4.1

Semi-3D simulation of the COS in the vicinity of a pressure bump
Before we describe our main results below, we briefly discuss one particular 3D simulation illustrated in the bottom panels of Figure 10.The latter shows the contours of the vorticity of a 3D simulation that adopted the axisymmetric nonlinear saturation shown in the top panels of the exact figure as the initial state.In particular, the zonal flow formed around r ≈ 0.9 becomes unstable to the Rossby wave instability (Lovelace et al. 1999) and spawns several small vortices.This is consistent with the occurrence of a maximum in the vortencity L (over-plotted dashed curve in all panels) corresponding to the zonal flow, where with the vertically integrate pressure Π.The vortices eventually merge and strengthen, and the vortensity maximum disappears.The vortices then excite (weak) spiral density waves (cf §5.4.4) and migrate radially inward.The corresponding angular momentum flux is minimal, though, taking values α r ∼ 10 −5 (not shown).
Note that this simulation's non-axisymmetric outcome is likely somewhat flawed due to the artificial 'switching on' of the azimuthal dimension only after 30,000 orbits.It is unclear if the zonal flow would have attained the same strength if the simulation had been 3D from the outset.

Nonlinear saturation of the COS in a smooth disk
Figure 11 illustrates the nonlinear saturation of the COS in our fiducial 3D simulation (p = 1.5, q = 2).Compared to Figure 8 we show here the evolution of the vertically averaged vertical velocity (left panels) in the rφ-plane.The vertical structure of the flow is similar to our 2D simulations (5.3), but no notable persistent zonal flows and elevator flows form, as these are rapidly turned into vortices.The right panel shows the normalized azimuthal power spectrum of the mid-plane radial velocity field where n is a positive integer denoting the azimuthal mode number.The curves describe the rms vertical (red) and radial (green) velocities, as well as α r (t) (blue).Most power is located in modes with azimuthal mode numbers n = 3 − 5, owing to the excitation of spiral waves by vortices, which are responsible for the radial angular momentum transport, measured through α r (see §5.4.4).It is noteworthy that the magnitude of the rms radial and vertical velocities is reversed as compared with our 2D simulations.This holds true for all simulations we conducted in this study.We speculate that this difference is related to the presence of largescale vortices in 3D simulations, in contrast to persistent elevator flows in 2D simulations. In

Mass and angular momentum transport
In this section, we study the radial transport of mass and angular momentum resulting from the COS in 3D simulations in more detail.Note that the negative (timeaveraged) angular momentum flux induced by the nonlinear COS in our axisymmetric simulations is minimal (cf.§5.3.1),such that the disk's radial density structure is not altered to any significant extent (not shown).This is consistent with the results of TL21.Despite its small magnitude though, the axisymmetric radial angular momentum flux is the driving element for the formation of zonal flows, and hence vortices (see §6 of TL21 for details).As mentioned above, in our 3D simulations we generally find a significant outward radial transport of angular momentum.Interestingly, in most cases this goes along with an outward radial transport of mass.For the case p = 1.5 this is true at least for cooling times β = 0.1 and β = 1.Indeed, using Eq. ( 50) we find timeaveraged accretion rates ⟨ Ṁ ⟩ t = −8.3• 10 −10 , 1.6 • 10 −9 , 2.3 • 10 −7 , 1.5 • 10 −8 , for simulations adopting the initial slopes p = −0.25,0.5, 1.5, 2.5, in units M * Ω K0 .We also find ⟨ Ṁ ⟩ t = 1.6 • 10 −8 , −8.9 • 10 −8 , for simulations adopting the cooling times β = 0.1, 10 (and p = 1.5).The finding of a positive accretion rate ( 50) is contrary to the case of a laminar viscous accretion disk with a constant kinematic viscosity.Nevertheless, our results are consistent with axisymmetric viscous accretion disk theory (Lynden-Bell & Pringle 1974) as we show in the following.
Using (50) we can define the effective accretion velocity, as obtained from our simulations, through  (53) where in the current situation (i.e. for an α-viscous prescription of a turbulent, non-axisymmetric accretion disk) we define the axisymmetric kinematic shear viscosity via In particular, α r (r, t) and thus ν may vary with radius (and time), which explains our results, as shown below.Note that the angular frequency Ω ≈ Ω K and the sound speed c 2 s given through ( 12) and ( 10) do not change significantly during our simulations.Using ( 51) and ( 53) we find (for an axisymmetric laminar, viscous Keplerian accretion disk) where we assumed for simplicity.Note that a(t) and p(t), and hence also m(t) vary with time in our simulations.In contrast, as mentioned above, q remains essentially unchanged.In actuality, the density profile substantially deviates from a simple power law at later stages of our simulations (cf.Figures 14 and 16 below), but for simplicity we ignore this here.In the current approximation a steady state density profile (such that the accretion rate is radially constant) strictly requires m(t) = 0. Furthermore, the pre-factor [m(t) − 1 2 ] determines the sign and overall magnitude of the accretion rate.
Figure 14 shows time-azimuth-and height-averaged radial profiles of α r (top panel) and ρ (bottom panel) for simulations adopting different initial density slopes p.The dashed curves are power law profiles that have been fitted to the solid curves.The corresponding averaged power law slopes are a(t) = 0.74, 1.40, 3.68, 0.69 and p(t) = −0.68,−0.23, −0.32, 0.68 for simulations adopting the initial slopes p = −0.25,0.5, 1.5, 2.5.Using these values we find m(t) = 0.41, 1.29, 3.36, 1.46.Thus, once the COS saturates, none of our simulated disks is predicted to remain in steady state, which is indeed what we find in our simulations.Note that a larger absolute value of m corresponds to an increasingly radially varying accretion rate, which should result in an increasingly large modification of the radial density profile.Indeed, the above m-values are qualitatively compatible with the amount of radial mass-redistribution in our simulations, which is shown in Figure 16 (right frames).The left frames show simulations with different cooling times.As seen in the panels for N 2 , large portions of the simulation region are stabilized with respect to the COS at later stages in some of the simulations.Moreover, almost all simulations form at least one pressure bump, such that ∂ r ⟨P ⟩ φz = 0.
The aforementioned values of m(t) are also in qualitative agreement with the values of the time-averaged accretion rate ⟨ Ṁ ⟩ t given above.Indeed, using the values of m(t) just derived, only the simulation with p = −0.25 is predicted to have a negative accretion rate, whereas the other simulations are predicted to yield a positive accretion rate.Moreover, the largest value is predicted for the case p = 1.5, followed by the case p = 2.5, etc. , which is in agreement with the values obtained from our simulations.Similarly, for the simulation with p = 1.5 and the longer cooling time β = 10 we find a(t) = −2.03,p(t) = 0.93 and m(t) = −1.10,which is also in qualitative agreement with the measured accretion rate, which is negative.
Finally, Figure 15 shows the time evolution of the radial accretion velocity, obtained from Eqs. ( 51) and ( 53).Based on the degree of agreement between the two quantities, we again conclude that the mass transport can qualitatively be described by an axisymmetric laminar viscous accretion disk model.

Vortex migration and spiral density waves
Angular momentum transport in our 3D simulation is presumably realized through the propagation and (shock-)damping of spiral density waves, excited by sufficiently strong vortices (see Paardekooper et al. 2010 for an explanation for the excitation mechanism).This is illustrated in Figure 17, where two strong vortices are present in the computational domain, both of which excite spiral waves.The left panel shows the vertically We compare the two different expressions for the radial accretion velocity ( 51) and ( 53), labeled "sim" and "axi", respectively.The former is rigorously derived from our simulations, while the latter assumes a laminar axisymmetric viscous disk with kinematic viscosity (54).
averaged vertical component of the vorticity where δv is the velocity on top of the equilibrium velocity ( 13).Furthermore, we consider the quantity (cf.§5.2) representing the Reynolds stress induced by nonaxisymmetric structures.This quantity is displayed in the middle panel of Figure 17, showing that radial angular momentum transport is entirely due to vortices.Moreover, the plot in the right panel shows a vertical cut through the upper left vortex, illustrating that strong COS-vortices can possess a significant vertical flow structure.
The positive radial flux of angular momentum induced by the vortices leads to radial inward migration of the latter (Paardekooper et al. 2010).This is illustrated in Figure 18.Here, we follow Manger et al. (2020) and display the quantity min[⟨ω z ⟩ z − ⟨ω z ⟩ zφ ], where we additionally apply a smoothing filter to this quantity to suppress the visibility of sharp features related to spiral waves, and to increase the visibility of larger vortices.The dark diagonal streaks in the figure represent vortices, which undergo rapid inward migration with typical speeds of 0.01H per orbit.The faint horizontal bands visible in some panels are spiral waves that survived our smoothing procedure.It appears that vortices are more numerous in our simulations with larger p ≥ 1.5.Also, in general they more rapidly grow strong in the same simulations, which is reflected in the earlier inward mi-gration after their formation.This is consistent with the values of α r and Ṁ in figure 13.Furthermore, merging of vortices can be seen in many cases.In addition, we find that the COS can, in principle, generate large-scale, long-lived vortices that may be observable in real PPDs based on their size.An example vortex with a radial width ≈ 3H is presented in Figure 19.
Finally, we note that the somewhat unexpected behavior of the simulation with the longer cooling time β = 10 (for which we expect the COS to develop significantly later than in the case with β = 1) could be related to the fact that in the nonlinear state, sub-critical baroclinic vortex amplification (Petersen et al. 2007a,b;Lesur & Papaloizou 2010) should be active within the vortices in our simulations, and the fact that the latter typically rotate at a much slower rate than the local disk rotation rate, with a turnover frequency ∼ Ω K /χ, with the vortex aspect ratio χ, denoting the ratio of the vortex semi-major and semi-minor axes.Therefore, the larger the vortex aspect ratio, the longer the rotation time, and as such the longer the optimal cooling time for baroclinic amplification (Raettig et al. 2013).Indeed, we find the occurrence of strong and large vortices in the early stages of this simulation before 400 orbits (not shown), which produce spiral density waves which undergo strong shocking.The gradual decay of α r , Ṁ and α z for the simulation with β = 10 after 400 orbits is likely related to the strong modification of the radial disk structure.Figure 16 shows that N 2 > 0 almost everywhere within the diagnostic domain, such that the COS, and thus the production of new vortices, should be strongly suppressed.Nevertheless, a detailed understanding of the effect of the cooling time scale on COS vortices is beyond the scope of this paper.

Elliptic instability of vortices
Our simulation results suggest that vortices with a sufficiently small aspect ratio χ are subjected to the elliptic instability (Lesur & Papaloizou 2009), which should reduce their lifetimes.To illustrate this, we compute the vertically averaged z-component of the vorticity (59) and relate it to the vortex aspect ratio via (Lesur & Papaloizou 2010) which is valid for a Kida vortex (Kida 1981) if the constant dimensionless parameter C = 1.Note that, strictly, this simplified correspondence assumes that the minimum vorticity in our simulations occurs within the core of a dominant (or representative) vortex at each time during the simulation.This is a strongly idealized picture for two reasons.For one, several vortices  are present in the simulation domain during the nonlinear COS stage.Furthermore, these vortices in general excite spiral waves.The latter form shocks in many cases, which can dominate the minimum and maximum values of the vorticity.To determine the proportionality factor C, we measure the aspect ratio of a strong and resolved vortex in each of our simulations, using at least three snapshots, and relate it to min [⟨ω z ⟩ z ] via Eq.( 61).To suppress the influence of spiral waves, we apply the same smoothing procedure as in Figure 18 above.We find C = 3.0, 2.8, 2.4, 1.8 for density slopes p = −0.25,0.5, 1.5, 2.5, respectively.For illustration Figure 20 shows the results for the cases p = 0.5 and p = 2.5.In this figure, the solid curves represent χ, whereas the dashed curves show the quantity max [v z ] − min [v z ], which should be an indicator for the onset of the elliptic instability, which is expected to attack vortices with χ ≲ 4 (Lesur & Papaloizou 2010).Indeed, in many cases where the computed values of the effective vortex aspect ratio drop below this value in our simulations, we do find notable increases (decreases) in the maximum (minimum vertical velocity).Note that the correspondence mediated by Eq. ( 61) is in some cases compromised by the presence of spiral wave shocks that survive the above-mentioned smoothing procedure and which artificially reduce our computed χ-values.The insert figure in each panel shows an example of a vortex with χ < 4, which contains fine-scale structures within its core, accompanied by strong vertical velocity fluctuations, indicative of the elliptic instability.

Summary
We studied the convective overstability near the midplane of radially global protoplanetary discs.We first conducted a linear stability analysis.Generally, we find that linear COS modes can only exist radially inward of their corresponding Lindblad resonance.In particular, the fastest growing modes are those whose Lindblad resonances lie within the computational domain, such that their growth rate increases with decreasing resonance radius.Since the modes are not resonantly excited, Lindblad resonances merely serve as turning points in the present context.
The radially global, nonlinear saturation state of the COS in axisymmetric simulations is similar to that in local incompressible shearing box simulations (TL21).We find that the flow is dominated by intermittent or persistent zonal flows and elevator flows, for sufficiently large values of pseudo-Richardson number R. However, differences occur.In §5.3.1 we found that for a viscosity α ≳ 10 −5 the COS is entirely suppressed in our simulations.For Re = 10 5 we find R crit ≈ 0.1 via Eq.( 20), which is only marginally largr than the largest values of R (i.e.close to he inner radial domain boundary) for our fiducial parameters.However, based on the results of TL21 we would still expect the COS to generate some hydrodynamic activity in form of a (disordered) nonlinear wave state in this simulation.The reason for this discrepancy is unclear, but could be related to the (inevitably) different boundary conditions applied in our simulations, as compared to TL21.Nevertheless, the close resemblance of the results for sufficiently large values of R indicates that the nonlinear saturation process of the axisymmetric COS in this parameter regime, which is the most relevant one to planetesimal formation in PPDs, is well captured by a radially local, incompressible model.
On the other hand, the radially global nature of the COS is not captured in local shearing box models.This is expected to be relevant in more realistic, non-smooth disc setups containing density features such as pressure bumps or gaps.Here, we find, in agreement with our linear calculations, that the COS may be locally excited but that its influence penetrates stable regions radially inwards.Moreover, in simulations applying radial damping zones, we find that the COS can be suppressed entirely (for at least 10,000 orbital periods), which is not expected if it were a radially local instability.
The nonlinear evolution revealed by our 3D simulations appears to be rather complex, and is -in contrast to axisymmetric simulations -no longer governed only by the two quantities Re and R ( §3).If the COS is sufficiently strong, vortex formation ensues, which can result in significant radial mass and angular momentum trans-port via the excitation of spiral density waves.This transport modifies the radial disc structure, including forming pressure bumps.Moreover, our results show that the (initial) radial disc structure (described by the two slopes p and q) affects the intensity of vortices and spiral waves, and therefore also the magnitude of the ensuing radial mass and angular momentum transport.This is another aspect that local shearing box simulations cannot capture.Notably, the COS can induce outward radial mass transport, i.e., decretion, which is found to be the case in most of our 3D simulations.In contrast, in simulations of VSI-active discs, one typically finds an inward mass flux (Manger & Klahr 2018;Lehmann & Lin 2022).Note, though, that such simulations are vertically stratified, unlike those considered here.
Our 3D simulations reveal the emergence of largescale, long-lived vortices.These generally migrate radially inward with rates ∼ 0.01H per orbit, similar to vortices generated by the VSI (Manger et al. 2020;Lehmann & Lin 2022).We also find that strong vortices with small aspect ratios are subject to elliptic instability, which affects their lifetimes.Moreover, intense vortices, such as the one in Figure 17, can contain significant vertical structure in their flow field, whereas sufficiently weak vortices, such as the one in Figure 19, lack any notable In both cases, the displayed radial region is 0.7 < r/r0 < 1.1 and the azimuthal extent ∆φ = 2.Note that vortices around 300 orbits in the simulation with p = 0.5 are tiny, such that the elliptic instability likely cannot be resolved, explaining the absence of notable increases of vz.vertical structure.These differences are expected to affect the concentration of dust by such vortices.We also find that the COS can, in principle, spawn long-lived vortices (Figure 19) that should be large enough to be observable by ALMA and could, therefore, be a possible explanation for observed asymmetries in PPDs.

Implications for dust substructures
Small dust grains tightly coupled to the gas are expected to possess a velocity where τ s denotes the particle stopping time and ρ tot is total density of the dust-plus-gas mixture (Youdin & Goodman 2005).Most of our 3D simulations yield at least one radial location where the radial pressure gradient vanishes (Figure 16).For sufficiently small |v|, dust grains respond to the local pressure gradient and are thus expected to drift radially into these regions and, if the pressure bump remains axisymmetric, produce dust rings.Dust rings may also result from 'traffic jams' produced from differential dust drift owing to a non-uniform radial pressure gradient, which is almost always the outcome of the COS.
On the other hand, our 3D simulations also reveal that the nonlinear COS gives rise to an averaged (turbulent) radial gas velocity (cf. Figure 15), which, as explained above, in most cases results in mass decretion.Whether this velocity is large enough to overcome pressure gradients and affect the dust concentration in the pressure bumps described remains to be seen in future work.
Finally, the COS readily produces vortices, which can trap dust radially and azimuthally (Lyra & Lin 2013;Raettig et al. 2015Raettig et al. , 2021)).The COS may, therefore, explain dust rings and asymmetries observed in PPDs (Andrews et al. 2018;Long et al. 2018;van der Marel et al. 2021).However, dust trapping cannot be too significant; otherwise, planetesimal formation would ensue, which is not observable.

Non-local particle stirring
In the more probable scenario where N 2 < 0 is only realized near particular radii, our results indicate that the COS can still drive activity interior to such regions due to its radially global character.This provides a background turbulence to stir up dust grains.From the vertical Reynolds stress rms(α z ) (Figures 12-13), we can estimate the expected height of a passive dust layer H d to be (e.g.Lin 2021) with the particle Stokes number St = τ s Ω K , where it is assumed that St ≫ α z .Considering a standard Minimum Mass Solar Nebula (Chiang & Youdin 2010), pebbles with an internal density of 1 gcm −3 and a size of 3mm would have St ∼ 10 −3 at 1AU, which are then predicted to possess a scale height of H d ≈ 0.03H.If the verticallyintegrated metallicity Z = 0.01, then the resulting midplane dust-to-gas ratio ϵ ∼ 0.3 would be insufficient for the streaming instability to drive strong clumping and hence planetesimal formation should not proceed (Johansen & Youdin 2007).Furthermore, turbulent stirring in of itself stabilizes the streaming instability (Umurhan et al. 2020;Chen & Lin 2020).We can, therefore, expect a finite-sized pebble layer to be maintained radially inward from where the COS is locally excited.

Caveats and outlook
In this study, we applied a simple optically thin gas cooling law with a single parameter β.However, recent studies indicate that cooling times favorable for the COS most likely imply that gas cooling is in the optically thick regime close to the disk mid-plane.Future simulations should, therefore, consider more realistic gas cooling prescriptions, incorporating thermal diffusion.Furthermore, as our simulations are vertically unstratified, they omit the possible presence of the VSI, which may occur for the parameter regimes considered here, if vertical stratification and more realistic cooling were accounted for.Fully global models are required to study the possible interplay of the COS and the VSI in a realistic manner.
Our example simulation in which the COS results from the presence of a pressure bump shows the formation of weak non-axisymmetries and vortices only after ten thousands of orbits.It is, therefore, speculative whether the required time scales are relevant to planetesimal formation in PPDs.It further remains speculative if the resulting vortices can efficiently concentrate dust.Also, more detailed observational constraints on the radial density and temperature structure of PPDs are required to determine if the COS is potentially relevant to planetesimal formation and whether or not it is only expected to occur near certain density structures such as pressure bumps.
Due to limitations in computational resources, our results on the radial mass transport in 3D simulations are not converged with respect to spatial resolution (see Appendix C), and, perhaps more importantly, vertical domain size (not shown).The latter implies that a quantitative study of the radial mass transport due to the COS formally requires a vertically stratified disk model.Such simulations will also need to adopt even higher grid resolutions than in this work, which poses a great computational challenge.
We note that we do not fully understand the reasons behind the exceptionally large radial mass flux occurring in our fiducial 3D simulation with p = 1.5, as compared to the other simulations.We find (not shown) that in the former simulation the χ-values computed via (61) seem to hover on the threshold for marginal stability (i.e.χ ≈ 4), unlike the other simulations , where χ drops well below the threshold.Perhaps the vortices in the simulation with p = 1.5 are marginally unstable to the elliptic instability, and might therefore survive longer as those in the other simulations.Moreover, we find a relatively large pileup of mass in the diagnostic domain in this simulation, compared to the other simulations (see Figure 14).This explains to some extent the large mass flux.However, it remains unclear how this behaviour is connected to the equilibrium density slope p.We speculate that the radial vortencity profile ρ/ω z , which yields a constant in the special case p = 1.5 plays a role in the evolution of vortices in our simulations.We therefore suggest that future work should study the evolution of vortices in radially global disk models in a more systematic manner.
Note also that our discussion on the elliptic instability in §5.4.5 is not rigorous.As mentioned earlier, the correspondence mediated by Eq. ( 61) is in some cases compromised by the presence of spiral wave shocks that survive our smoothing procedure and which artificially reduce our computed χ-values.Moreover, the measurements of vortex aspect ratios are done by hand and are, therefore, subject to significant potential inaccuracy.Furthermore, vortices in our simulations are not necessarily Kida vortices, which are defined by a constant patch of relative vorticity.Nevertheless, our results strongly suggest that the elliptic instability occurs in the cores of COS-vortices of sufficiently small aspect ratio.
In a follow-up study, we will investigate the impact of the nonlinear COS on dust concentration in PPDs.In particular, we will clarify if COS vortices in radially global PPD models can result in dust-to-gas density ratios large enough to trigger vigorous streaming instability.
We thank an anonymous reviewer for useful comments.This work is supported by the National Science and Technology Council (grants 112-2112-M-001-064-, 112-2124-M-002-003-, 113-2124-M-002-003-) and an Academia Sinica Career Development Award (AS-CDA-110-M06).Simulations were performed on the Kawas cluster at ASIAA and the Taiwania-2 cluster at the National Center for High-performance Computing (NCHC).We thank NCHC for providing computational and storage resources.Figure 21 shows results of 2D simulations with different values of numerical parameters.These are the spatial resolution (upper panel), the radial domain size (middle panel) and the vertical domain size (bottom panel).The plots in the upper panel demonstrate that our 2D simulations, using a resolution of 200/H are converged.The plots in the middle panel reveal that the saturation level of the COS generally reduces with decreasing radial domain size.Based on the results of our linear analysis in §4.4.1 this is expected, and is due to the fact that the maximum value of R defined in ( 19) is reduced (the value at the inner domain boundary), as well as our finding that for vertical wavenumbers k z H 0 ≳ 10 the fastest growing linear modes (i.e. the modes that can "access" the largest R-values) are strongly concentrated towards the inner disk boundary.Thus, with decreasing radial domain size these modes are gradually removed.
We also find a significant drop in the saturation level when decreasing the vertical domain size as ∆z = 0.5H 0 → 0.25H 0 , whereas the corresponding drop is only modest for the reduction ∆z = H 0 → 0.5H 0 .These reductions imply that we increase the minimum vertical wavenumbers k z,min = 4π → 8π and k z,min = 2π → 4π, respectively.Considering Figure 2 (left panel) these findings are plausible, since the largest growing modes in our simulations possess 10 ≲ k z H 0 ≲ 100.

C.2. 3D simulations
As pointed out in §6.1, the radial angular momentum flux in our 3D simulations is not converged with respect to resolution in all spatial dimensions, which is illustrated in Figure 22 (lower panel).This can likely be traced back to the formation of spiral shocks (which generate very sharp gradients in the hydrodynamic quantities), which are ultimately responsible for the deposition of angular momentum in the disk.On the other hand, the plot in the upper panel shows that rms radial (dashed curves) and vertical (solid curves) velocities resulting from the COS are actually converged.Our chosen resolution of 200/H in both radial and vertical directions proves computationally demanding in 3D global simulations, making it currently impractical to pursue an even higher resolution.Nevertheless, as already noted in §6.1, a detailed study of the radial mass flux requires a vertically stratified disk model.
Figure 23 illustrates the robustness of our (3D) simulation results with respect to the radial boundary conditions.This is not a straight forward issue in the present study, as the fastest growing COS modes are concentrated towards the inner domain boundary.The plots show the time evolution of the averaged density slope ∂ ln ρ/∂ ln r (upper panel) and the α-viscosity (lower panel).Compared are fiducial simulations with and without wave damping radial boundary conditions.We find that the application of a damping zone near the boundaries merely results in a delayed nonlinear saturation of the COS, of equal magnitude as in the simulation with default boundaries.The insert plots in the upper panel show the initial and final radial profiles of the squared buoyancy frequency N 2 and pressure P , in good agreement between the two simulations.In §5.3.1 we found that rms vertical velocities in 2D simulations do not show the expected behavior of the COS with increasing cooling time.That is, for large cooling times β ≳ 100 they exhibit an unexpected rise.We suspect this to be caused by the presence of another instability.In the following we briefly show that a "double-diffusive" instability (not the COS) can operate in the adiabatic limit and in the presence of an unstable entropy gradient, as well as a considerable shear viscosity.It is double-diffusive in that is relies on a large (viscous) diffusion of angular momentum and small thermal diffusion.
We start by presenting in Figure 24 rms vertical velocities (right panel) of adiabatic simulations (β = 10 4 ) The simulations in the top panel were conducted using a vertical domain size ∆z = H.Note that all of these simulations adopted wave damping boundaries.and constant viscosity taking values α = 10 −5 − 10 −2 .In these simulations we use q = 2 and p = 0.5, which results in N 2 < 0, and a negligible equilibrium viscous accretion flow according to Eq.(B5).We see that for α ≳ 10 −4 the system settles on appreciable rms velocity values.Furthermore, the left panels show from top to bottom radial profiles of the vertically averaged radial velocity, squared buoyancy frequency and entropy, respectively.As seen, the instability gives rise to a ra- dially varying radial outward velocity, such that (due to the lack of cooling) entropy can be transported outwards, which results in a flattening of the entropy profile.The buoyancy frequency develops strong, stochastic deviations from the ground state profile.However, the averaged value of N 2 (over height and radius) saturates to a finite negative value (not shown), explaining the persistence of the instability in the nonlinear regime.Note that the instability similarly develops for different values of p.The difference is that the radial viscous accretion velocity (B5) adds up to the radial velocity field induced by the instability.Depending on the sign of the former (which depends on p) this somewhat speeds up or slows down the flattening of the entropy profile.

D.2. Local incompressible description
The existence of such an instability was already predicted by Latter et al. (2010), who studied a "resistive  39)-( 41)), for simplicity presented in the limit10 k x → 0 and β → ∞ (adiabatic limit), as well as vanishing ground state velocity of the gas: D20) where the radial entropy stratification is described by the quantity 1 and N 2 is given by ( 16).Furthermore, α is a dimensionless viscosity parameter.The above equations are defined in a Cartesian frame (x, y, z) within small rectangular patch of the disk, centered on a fiducial point (r 0 , ϕ 0 − Ω 0 t, z = 0), that rotates around the star at approximately the local Keplerian frequency.The unit vectors ⃗ e x , ⃗ e y , ⃗ e z in the box corresponding to the radial, azimuthal, and vertical directions in the global disk.
In the inviscid case α → 0, by taking the time derivative of Eq. (D21) we find which yields stable, buoyancy-adjusted epicycles.However, if we assume that the viscosity is non-zero and consider an initial radial velocity perturbation δv x > 0.
Then, for sufficiently large k z the viscous term in the azimuthal momentum equation prevents the azimuthal velocity δv y from becoming large enough to provide a restoring motion via the Coriolis term in the radial momentum equation.At the same time, the buoyancy term is found to be positive via (D20) and needs to exceed the viscous damping term in (D21) in order for δv x to grow.The ensuing instability is thus expected to be non-oscillatory, which can easily be confirmed as follows.The dispersion relation for local modes ∝ exp (ik x x + ik z z − iωt) with k x = 0 following from above equations reads ω 3 +ω 2 2ik 2 z αh 2 0 −ω 1 + N 2 + k 4 z α 2 h 4 0 −ik 2 z αh 2 0 N 2 = 0. (D25) If we assume purely growing modes with ω ≡ iω I we find that all coefficients of the resulting dispersion relation are positive, such that the condition is a sufficient condition for instability, expressing the necessity of a non-vanishing viscosity and an unstable entropy gradient with N 2 < 0. Nevertheless, even though the linear instability can be understood within the local Boussinesq framework, its nonlinear saturation involves considerable changes of the background disk structure, as shown above.Furthermore, the relevance of this instability to PPDs is questionable, since the considered values of α would require a significant level of turbulence generated by some other source in the first place.

Figure 1 .
Figure1.Contour plot showing the minimal value of N 2 /ΩK(r0) 2 within a radial region r0 ± 2H0 in a power law disk with density and temperature slopes p and q, respectively.The diagonal lines delineate the regions where min N 2 takes negative or positive values, such that the COS is in principle expected to occur in the parameter-space labeled "COS".Open (solid) circles represent values of p and q used in our 2D-simulations (2D as well as 3D simulations) below.

Figure 4 .
Figure 4. Linear growth rates of radially global COS modes (solid curves) as function of radial wavenumber kr for fixed kzH0 = 10 and for different radial sizes of the computational domain ∆r are compared with fully local growth rates (circles), computed via Eq.(42).The radial wavenumbers kr are obtained from (29) and correspond to radial averages.See the text for details.The insert figures illustrate contours of the vertical velocity perturbation of the fastest growing modes in the cases ∆r = 5H0 and ∆r = H0.

Figure 5 .
Figure 5. Maximal linear growth rates of radially global COS modes in a smooth disk.Compared is the effect of a varying α-viscosity (left panel), cooling time (second panel), temperature slope (third panel) and density slope (rightmost panel).If not indicated otherwise parameters correspond to the fiducial disk (q = 2, p = 1.5, β = 1, α = 0).Open circles in the second panel from the left correspond to growth rates obtained from a fully local model.

Figure 6 .
Figure 6.Illustration of the linear COS induced by a density structure in an otherwise linearly stable disk.The Left panel shows radial profiles of density ρ (and the effective density slope p eff ) and squared buoyancy frequency N 2 of the smooth disk, and disks containing either a density bump or a density step with amplitude A = 0.4.The right panel shows linear growth rates of COS modes for the latter two cases for different vertical wavenumbers kz.

Figure 7 .
Figure7.Hydrodynamic quantities (their real part) associated with the fastest growing COS mode with kzH0 = 20 in a disk containing a pressure bump (43) with A = 0.4, as obtained with the shooting method.The mode is excited in the region 0.7 ≲ r/r0 ≲ 1 where N 2 < 0 (see Figure6) but penetrates into the stable region at smaller radii.The bottom panel shows the wavelet power spectrum of δvz as a function of radial wavenumber kr.The red dashed curve is the radial wavenumber resulting from (35).The Lindblad resonance location where ωR = Ω is the location where the modes' radial wavenumber drops to zero.

Figure 8
Figure 8. a): The left panels show snapshots of the vertical velocity in the course of the axisymmetric nonlinear saturation of the COS, where brighter (darker) colors correspond to larger (smaller) values of δvz.Note that the color scaling of the bottom right frame at 1000 orbits is the same as in the left panel in the bottom frame b).The right frame shows the power spectral density of the radial velocity in vertical wavelength space, and illustrates the increased vertical elongation of structures, eventually saturating into elevator flows and zonal flows with almost vanishing vertical structure.The curves represent rms radial and vertical velocities as indicated.The circles and the dashed line represent αr(t) and its time-average, respectively.See the text for details b): Saturated state of the COS after 1,000 orbital periods, exhibiting persistent elevator flows (left panel) and zonal flows (middle panel).The dashed curves in the middle panel illustrate the geostrophic balance of the displayed zonal flows as explained in the text.

Figure 9 .
Figure 9. rms vertical velocities in 2D axisymmetric hydrodynamic simulations.Compared are different cooling times β (upper left), initial temperature slopes q (lower left), initial density slopes p (lower right) and values of the dimensionless viscosity α defined via Eqs.(18) and (8) at r = r0.The dashed line near 100 orbits in the upper left panel represents the maximum growth rate predicted by the linear radially global analysis above.It matches well with the largest linear growth rate seen in the simulation for the same disk parameters.Note that the COS in simulations with α ≳ 10 −4 is fully damped for at least 10,000 orbits (not shown).

Figure 10
Figure 10.a): Nonlinear saturation of the COS in the vicinity of a pressure bump.The upper left frame shows the evolution of the rms vertical velocity.The dashed line indicates the linear growth predicted by our linear analysis.The insert figure displays initial and final vertical velocities (at z = 0).The upper right frame shows contours of the vertical and azimuthal velocity perturbations.The over-plotted curves are the initial radial buoyancy and pressure profiles.b): Contours of the vorticity (59) defined below of a 3D simulation for which we used the axisymmeteric saturation depicted in the top panels as initial state.The dashed curves represent the vortencity L, defined in Eq. (48).

Figure 11 .
Figure 11.Nonlinear evolution of the COS in the fiducial 3D simulation.The left panel shows snapshots of the vertically averaged vertical velocity.In the right panel the contour plot represents the time evolution of the azimuthal power spectral density of the radial velocity (planar cut at z = 0).The curve labeled "⟨P (n)⟩t" is the time-averaged normalized function (49), corresponding to the left plot-axis.The remaining curves are the time evolution of the rms vertical velocity, rms radial velocity and αr, corresponding to the right plot-axis.β = 0.1, β = 1 and β = 10, as well as four different values of the initial density slope p = −0.25,p = 0.5, p = 1.5 and p = 2.5.The latter values are also shown in Figure 1.In the simulations adopting different p we slightly adjust the temperature slope q, such that the radial profile of N 2 is the same in all simulations.Thus, essentially q ≈ 2 in these simulations.In addition, we consider one semi-3D simulation in the vicinity of a pressure bump, which we discuss first in §5.4.1

Figure 12 .
Figure 12.Radial transport of angular momentum (top left), rms vertical velocities (bottom left), radial mass accretion rate (top right), and rms vertical transport of angular momentum (bottom right) in 3D simulations with varying cooling time β.The quantity M d denotes the initial mass contained within the diagnostic domain in each simulation.
Figures 12 and 13  we show the time evolution of α r (upper left panels), the rms vertical velocities (lower left panels), the radial mass accretion rate panels), and rms(α z ) (lower left panels), for different cooling times β, and density slopes p, respectively.Our measured α r -values are in the range 10 −4 − 10 −3 , roughly consistent with values obtained from local compressible shearing box simulations(Lyra 2014).

Figure 13 .
Figure 13.Same as Figure 12 but now for simulations with varying initial density slope p.

Figure 14 .
Figure 14.The top and bottom panels show radial profiles (solid curves) of the time-averaged αr(r, t) as defined in (55) and mass density ρ (averaged over time, azimuth and height), respectively, for simulations adopting different initial density slopes p.The dashed lines are corresponding linear fits.Note that the plots are double logarithmic.

Figure 15 .
Figure15.Testing the applicability of axisymmetric viscous accretion disk theory to our 3D simulations of the COS.We compare the two different expressions for the radial accretion velocity (51) and (53), labeled "sim" and "axi", respectively.The former is rigorously derived from our simulations, while the latter assumes a laminar axisymmetric viscous disk with kinematic viscosity (54).

Figure 16 .
Figure16.Effect of the nonlinear COS on the radial disk structure.From top to bottom we show initial (dashed curves) and final (after 1,000 orbits, solid curves) profiles of the radial pressure gradient (the vertically and azimuthally averaged pressure), the density and the squared radial buoyancy frequency.The left (right) panels compare different cooling times β (initial density slopes p).

Figure 17 .
Figure 17.Illustration of the excitation of spiral density waves by strong vortices induced by the COS in the simulation with p = 2.5.a): Vertically averaged vertical component of vorticity defined in (59), as well as the quantity ∆Rrφ defined in (60).b): Vertical structure of the upper left vortex in the left frames.The arrows indicate velocities in the rz-plane and the dashed curve denotes the pressure deviation ∆P , showing the pressure maximum associated with the vortex.

Figure 18 .
Figure 18.Illustration of the formation and subsequent migration of large scale vortices in 3D simulations with different initial density slopes p.The plots show contours of the quantity min[⟨ωz⟩z − ⟨ωz⟩zφ] (we plot its square root for improved visibility).The dark small-scale features, notably seen in the lower left panel, result from spiral wave shocks which survived our applied smoothing algorithm.

Figure 19 .
Figure 19.Appearance of a large scale (radial width ≈ 3H) long-lived vortex with aspect ratio χ ≈ 7 in the 3D simulation adopting a short cooling time β = 0.1 (after ≈ 1, 300 orbits).This vortex exhibits practically no vertical structure and vertical flow field (not shown), in contrast to the vortex shown in Figure 17.The left panel shows the vertically averaged z-component of vorticity (59).The arrows illustrate the retrograde flow field associated with the vortex.The over-plotted dashed curve is the scaled pressure deviation ∆P (at z = 0), clearly showing a pressure maximum associated with the vortex core.The right hand panel shows a polar plot of the scaled density deviation ∆ρ .Spiral waves excited by the vortex are clearly visible in the right panel.

Figure 20 .
Figure 20.Illustration of the elliptic instability in vortices resulting from the COS.The solid curves show the effective aspect ratio χ of vortices (61) occurring in 3D simulations with different initial density slopes p.The dashed curves represent the maximum vertical velocity.The insert figures show examples of vortices where the elliptic instability is suspected to be active.In both cases, the displayed radial region is 0.7 < r/r0 < 1.1 and the azimuthal extent ∆φ = 2.Note that vortices around 300 orbits in the simulation with p = 0.5 are tiny, such that the elliptic instability likely cannot be resolved, explaining the absence of notable increases of vz.
D. DOUBLE-DIFFUSIVE INSTABILITY IN SLOWLY COOLED VISCOUS DISKS D.1.Occurrence in hydrodynamic simulations

Figure 21 .
Figure 21.Testing convergence of 2D axisymmetric simulations with respect to resolution (top panel), radial domain size (middle panel) and vertical domain size (bottom panel).The simulations in the top panel were conducted using a vertical domain size ∆z = H.Note that all of these simulations adopted wave damping boundaries.

Figure 22 .
Figure 22.Testing convergence with respect to spatial resolution in 3D simulations.We consider convergence simultaneously in radial and vertical direction (labeled "XZ") as well as azimuthal direction (labeled "Y").The upper panel shows rms vertical (solid curves) and radial (dashed curves) velocities.The lower panel shows the time evolution of αr.The curves in the lower panel have been smoothed for clarity.The dashed curve in the lower panel describes the fiducial simulation.

Figure 23 .
Figure 23.The plots illustrate the robustness of our simulation results with respect to the radial boundary conditions.

Figure 24 .
Figure 24.Illustration of the nonlinear saturation of the double diffusive instability powered by viscous momentum diffusion and a negative radial entropy gradient.The right panel shows rms vertical velocities in simulations with increasing kinematic shear viscosity α, defined via (18).The left panels show the evolution of the vertically averaged radial velocity, squared buoyancy frequency and entropy, in the simulation with α = 10 −3 .