Mass and Angular Momentum Transport in a Gravitationally Unstable Protoplanetary Disk with Improved 3D Radiative Hydrodynamics

During early phases of a protoplanetary disk's life, gravitational instabilities (GIs) can produce significant mass transport, can dramatically alter disk structure, can mix and shock-process gas and solids, and may be instrumental in planet formation. We present a 3D grid-based radiative hydrodynamics study with varied resolutions of a 0.07 M ⊙ disk orbiting a 0.5 M ⊙ star as it settles over most of its radial extent into a quasi-steady asymptotic state that maintains approximate balance between heating produced by GIs and radiative cooling governed by realistic dust opacities. We assess disk stability criteria, thermodynamic properties, strengths of GIs, characteristics of density waves and torques produced by GIs, radial mass transport arising from these torques, and the level to which transport can be represented as local or nonlocal processes. Physical and thermal processes display distinct differences between inner optically thick and outer optically thin regions of the disk. In the inner region, gravitational torques are dominated by low-order Fourier components of the azimuthal mass distribution. These torques are strongly variable on the local dynamical time and are subject to rapid flaring presumably driven by recurrent swing amplification. In the outer region, m = 1 torques dominate. Ring-like structures exhibiting strong noncircular motions, and vortices develop near the inner edge between 8 and 14 au. We find that GI-induced spiral modes erupt in a chaotic manner over the whole low-Q part of the disk, with many spiral modes appearing and disappearing, producing gravitoturbulence, but dominated by fluctuating large-scale modes, very different from a simple α-disk.


INTRODUCTION
Recent capabilities provided by the Atacama Large Millimeter/submillimeter Array (ALMA) and other instruments have revealed beautifully detailed structures, including rings, spiral arms, and forming gas giants, in protoplanetary disks (PPDs) (ALMA Partnership et al. 2015;Huang et al. 2017Huang et al. , 2018aHuang et al. ,b, 2020;;Paneque-Carreño et al. 2021;Currie et al. 2022).During the early phases (Class 0 and 1) of a PPD's life, gravitational instabilities (GIs) can produce significant mass transport, dramatically alter disk structure, mix and shock-process gas and solids, and be instrumental in planet formation.
The existence, persistence, and characteristics of GIs in PPDs are driven primarily by thermal and accretion processes.When cooling is sufficiently rapid, fragmentation can occur (Boss 1997;Gammie 2001;Mayer et al. 2002).Such disk fragmentation can further lead to gravitationally bound clumps, which themselves offer a wide range of evolutionary possibilities.For example, the clumps can undergo nontrivial disk migration (Vorobyov & Basu 2010;Baruteau et al. 2011;Michael 2011;Michael et al. 2011), grow to become brown dwarfs or stars (Kratter et al. 2010a,b), or remain as giant planets (Boley et al. 2010).Moreover, clumps that are initially bound are not guaranteed to remain so, and through disk migration, scattering, or clump-clump interactions, they may become tidally disrupted (Boley et al. 2010;Nayakshin 2010).This disruption could in turn serve as a mechanism for processing solids.For lower cooling rates, GI-active disks can instead develop quasi-steady structures or bursts of activity that dramatically transform the disk without forming fragments.For reviews, see Durisen (2011), Armitage (2011), Kratter & Lodato (2016), Rice (2016), and Armitage (2019).
Detailed three-dimensional (3D) hydrodynamic modeling of GI-active PPDs, as reported herein, can provide considerable insight into the physical and thermal states of these disks and their shortand long-term evolution.An important consequence of GIs in nonfragmenting PPDs is the ability of GI-generated spiral arms to drive angular momentum transport, a fact initially recognized by Lynden-Bell & Kalnajs (1972) in the context of galactic dynamics.A disk's susceptibility to GIs can be parameterized by the Toomre Q parameter, where c s is the sound speed, κ is the epicyclic frequency, and Σ is the disk surface mass density (Toomre 1981).In a disk subject to GIs, small density perturbations grow exponentially on a timescale comparable to the local dynamical time when Q < ∼ 1.7 in the linear regime (Durisen et al. 2007).In the nonlinear regime, perturbations can grow for even larger Q.These perturbations manifest themselves as predominantly trailing multiarm spiral density waves over a broad range of radii; see Papaloizou & Savonije (1991), Laughlin et al. (1998), Nelson et al. (1998), Nelson et al. (2000a), Pickett et al. (1998), Pickett et al. (2000), Pickett et al. (2003), and Michael et al. (2012).
Gravitational torques arising from these spiral structures enable the disk to tap the free energy associated with the rotational shear.Some of this energy is then returned as heat when waves steepen into shocks.This heating, along with net inward transport of mass, pushes the disk back toward stability.At the same time, radiative cooling acts in the opposite sense, pushing the disk toward continued instability.
Spiral structure is not the only type of morphology expected in disks.Indeed, while it is well established that rings can arise from disk interactions with embedded objects (Goldreich & Tremaine Transport in Gravitationally Unstable Disks 3 1980; Lin & Papaloizou 1986;Paardekooper & Mellema 2006;Zhu et al. 2011;Zhang et al. 2018), they can also emerge in a disk without an embedded perturber (e.g., Lubow 1991;Ogilvie 2001;Takahashi & Inutsuka 2014;Tominaga et al. 2019;Lee et al. 2019a,b;Riols & Lesur 2019;Riols et al. 2020;Li et al. 2021).In particular, eccentric modes, corresponding to perturbations with azimuthal wavenumber m = 1, have received particular interest in the context of PPDs because of their global nature.A large corpus of work has examined the development and sustenance of these modes in fluid disks (Adams et al. 1989;Shu et al. 1990;Hirose & Osaki 1990;Lubow 1991;Heemskerk et al. 1992;Laughlin & Korchagin 1996;Ogilvie 2001;Tremaine 2001;Papaloizou 2002;Tominaga et al. 2019;Lee et al. 2019a,b;Li et al. 2021;Béthune et al. 2021).It has been shown that almost any disk with a realistic density profile can sustain long-lived eccentric modes (Lee et al. 2019b).Moreover, these modes can grow in amplitude via the sling mechanism that amplifies an eccentric perturbation through the wobble of the central star (Adams et al. 1989;Shu et al. 1990;Lin 2015).Ring formation can follow via angular momentum exchange with the unstable eccentric mode (Lubow 1991;Ogilvie 2001;Lee et al. 2019a,b).A recent 2D study of an eccentric spiral instability in a self-gravitating disk with prescribed cooling by Li et al. (2021) found that a trapped one-arm instability forms early in the simulation and evolves into a set of axisymmetric rings.
For these reasons, we expect that rings formed by a process not involving embedded objects may be a common product of early PPD evolution.Durisen et al. (2005), in an analysis of ring-like structures that developed in their PPD hydro simulation, proposed a hybrid avenue for planet formation where even if instabilities due to disk self-gravity do not produce gaseous protoplanets directly, they may create persistent dense rings that are conducive to accelerated growth of gas giants through core accretion.They suggested that even if GIs do not lead to permanent clump formation, they may significantly accelerate core accretion by creating persistent dense gas rings near boundaries between GI-active and GI-inactive regions (see also Haghighipour & Boss 2003a,b).
Ring features have been seen in previous 3D hydro simulations of self-gravitating PPDs carried out by our group and collaborators (Pickett et al. 1996(Pickett et al. , 2003;;Mejía 2004;Mejía et al. 2005;Durisen et al. 2005;Cai 2006; Boley et al. 2006Boley et al. , 2007a;;Cai et al. 2006Cai et al. , 2008;;Michael et al. 2012;Steiman-Cameron et al. 2013;Desai et al. 2019).These studies found that rings form early in the simulation, well before disks settle into an asymptotic state where heating and cooling balance.The simulations described in the current work also lead to multiple rings, which we explore in §3.8.
Another important question arises if a GI-active disk settles into a quasi-steady saturated state.Specifically, in this state does thermal balance of heating and cooling apply locally or only in a global, long-term, time-average sense?Torques due to spiral waves involve long-range interaction for low-order (few-armed) spirals, and the wave nature of GI activity opens the possibility for wave transport of energy (Laughlin & Rozyczka 1996;Balbus & Papaloizou 1999).On the other hand, it has been proposed by several authors that GI transport and evolution can be captured by a turbulent α-disk formulation by applying a Shakura & Syunyaev (1973) α that can be obtained in the case that the energy balance is precisely local; see Paczynski (1978), Pringle (1981), Gammie (2001), Vorobyov (2010), andZhu et al. (2010).In his razor-thin shearing box simulations, Gammie (2001) found good agreement between such a local derivation and the effective α eff computed from the observed stresses in his simulations.However, these calculations were local by their very nature.Full 3D hydrodynamics calculations have given somewhat mixed results on this issue depending on the disk mass, numerical resolution, numerical techniques, and the nature of the cooling (e.g., Lodato & Rice 2004, 2005;Boley et al. 2006;Cossins et al. 2009;Forgan et al. 2011;Michael & Durisen 2010;Michael et al. 2012;Steiman-Cameron et al. 2013;Evans et al. 2015a;Béthune et al. 2021).
To address further some of the issues outlined above, we report here results of a grid-based, finitedifference, 3D radiative hydrodynamics convergence study of a gravitationally unstable PPD where cooling is allowed to adjust naturally by radiative transport using realistic dust opacities and where star-disk interactions are explicitly modeled.This study builds on the earlier works of Michael et al. (2012) (hereafter Paper I) and Steiman-Cameron et al. (2013) (hereafter Paper II) using code improvements in radiative transfer and the incorporation of star-disk interactions.The number of azimuthal elements used in calculations is especially important because it is nonaxisymmetric structures, i.e., spiral waves, that produce the gravitational torques leading to mass and angular momentum transport.We follow the time evolution of four simulations of a PPD which differ only in the azimuthal resolution of their computational grid, allowing us to ensure that results have converged to a solution that is not affected by the size of the azimuthal mesh.
Simulations are run for a time period sufficient for the disk to settle into a long-lived, statistically quasi-steady, asymptotic state, allowing GIs and their consequences to be characterized and quantized."Quasi-steady, asymptotic state" refers here to the evolutionary phase of a gravitationally unstable disk that has settled into a long-lived, quasi-steady balance between radiative cooling and disk heating at a relatively constant, but unstable, value of Q.
Paper II and this work examined the same disk with the same initial conditions and hydrodynamics code, but here with important improvements to the code (Section 2), including the implementation of a subcyling approach to better control heating and cooling and the inclusion of an indirect potential approach to self-consistently account for star -disk interactions.In contrast to Paper II, great convergence is found here, demonstrating the importance of doing the radiative physics well.
The balance of this paper is organized as follows.Section 2 provides the details of the numerical approach and defines the models.Results of the simulations are reported in Section 3, with physical and thermodynamic properties of the converged asymptotic disk described in Sections 3.1, 3.2, and 3.3; the time dependence of nonaxisymmetric modal properties presented in Section 3.4; computation of the gravitational torques arising from these structures are examined in Section 3.4; and angular momentum transport, mass flux, and time variability are described in Sections 3.5 and 3.6.The locality/nonlocality of mass transport due to GIs and the applicability of an effective α-based viscosity is discussed in Section 3.7.The development of ring-like structures and their impact on disk evolution is found in Section 3.8.Convergence is discussed in Section 3.9, followed by a discussion section in Section 4. A summary and conclusions are found in Section 5.

Hydrodynamics
We seek to understand the physical and dynamical characteristics of a resolution-independent PPD simulation at a time when the disk has relaxed into a quasi-steady asymptotic state characterized by a statistically constant unstable Q(r).To this end, we conduct multiple hydrodynamic simulations of a PPD using CHYMERA, an explicit, second-order, 3D Eulerian code that self-consistently solves the hydrodynamic equations of motion, Poisson's equation, and the energy equation, on a uniform cylindrical grid (Boley 2007;Boley et al. 2007a).The number of grid elements in the r-, z-, and ϕ-directions are specified by j max , k max , and l max , respectively, and mirror symmetry is assumed about the equatorial plane.
The equation of state used in this work takes into account contributions of the translational, rotational, and vibrational states of H 2 and assumes a fixed H 2 ortho-to-para ratio of 3:1 (Boley et al. 2007b).For the temperature range in our simulations (Section 3.2), the gas is well approximated by an adiabatic index γ = 5/3.

Radiative Cooling
The disk is embedded in a 3 K background, a boundary condition in the radiative scheme for the z-direction of the hydro code for the simulations run here.Disk cooling is implemented using the radiative energy transfer scheme of Boley et al. (2007a) that combines flux-limited diffusion for optically thick regions in the r-and ϕ-directions and a single-ray discrete-ordinate radiative transfer solver in the z-direction that treats both optically thick and thin regions.This scheme produces smooth temperature profiles across the disk photosphere.The opacity tables of D' Alessio et al. (2001) are used, assuming minimum and maximum grain sizes of 0.005 and 1.0 µm, respectively, and a power-law size distribution with index -3.5 within this range (for opacity details, see Appendix A of Boley et al. 2006).
Because the code explicitly solves for radiative transport, it can become unstable if the radiative time step becomes smaller than the hydrodynamic time step producing excessive heating or cooling of the gas that lead to unphysical results and numerical instability.To avoid this situation, in Paper II limiters were placed on the local cooling and heating rates to prevent the computational time step ∆t from becoming shorter than ∼ 0.03 outer rotation periods (ORPs), where ORP is defined as the initial orbital period at radial grid element j = 200 (r ≈ 33 au, 1 ORP ≈ 255 yr; see Section 2.3).While this eliminated the numerical instability, it artificially set a computational time step that might be unrealistically large.In the current work, a subcycling approach is used to control heating and cooling.At each hydrodynamic time step, the radiative routine calculates a separate radiative time step size from the ratio of the internal energy density to the divergence of the radiative flux for each cell and compares this with the hydrodynamic time step.If the radiative time step size is smaller than the hydrodynamic time step, radiative transfer is subcycled, i.e., multiple calls are made to the radiative transfer routine for that hydrodynamic time step.Details are described in Shabram & Boley (2013) and Evans et al. (2015a).Boley et al. (2007a) demonstrated that, for numerical stability and accuracy, the optically thick portion of the disk must be resolved by a minimum of about five to seven vertical cells.When the condition is not satisfied, vertical oscillations that are purely numerical can occur that could, in turn, lead to "artificial heating" by an uncontrolled numerical effect.The vertical resolution of our simulations is too small to satisfy this requirement interior to ∼ 8 au.Thus, we lack confidence in the simulations in this region.

Star-Disk Interactions
In Papers I and II, the star was represented by a point-source gravitational field held fixed at the center of the computational grid and star-disk interactions were not modeled.In fact, star-disk interactions will inevitably displace the star from the system center of mass (Rice et al. 2003).Here we account for the star's acceleration using the indirect potential method (Adams et al. 1989;Nelson et al. 2000b), as discussed in Michael & Durisen (2010).In this approach, the star remains fixed at the grid center for computational convenience while the reference frame of the star plus grid is accelerated through inclusion of fictitious forces that self-consistently account for the gravitational interactions between the star and disk.

Initial Conditions
To assess mesh convergence, four simulations were run following the evolution of a 0.07 M ⊙ PPD surrounding a 0.5 M ⊙ central star with a background temperature of 3 K.These simulations have identical initial conditions and differ only in the azimuthal resolution of the computational grid.The numbers of grid elements in the r-and z-directions are given by j max = 512 and k max = 64, respectively, for all simulations with each increment in j and k corresponding to 0.167 au.This j max provides a sufficiently large outer computational grid radius to keep all material on the grid when the disk expands during the violent onset of nonlinear GIs and throughout the entire disk evolution.The number of grid elements in the ϕ-direction is given by l max = 64, 128, 256, or 512.Simulations will hereafter be referenced by their l max .
Outflow conditions are used for the vertical and radial boundaries.These are chosen, as opposed to allowing inflow and outflow, to ensure that artificial mass streams do not form at the boundaries during fluxing.Mass that flows out of the grid is removed from the simulation.In the vertical direction, only the top boundary uses the outflow condition, while the midplane boundary assumes mirror symmetry.
Like the outer radius, the inner radius also assumes an "outflow" condition.However, for the inner boundary, instead of removing mass from the grid entirely, mass that passes through the inner boundary is added to the star's mass, thus assuming that accretion has taken place.
Time is expressed in units of ORP, defined by the orbital period in the initial (t = 0) disk at radial grid element j = 200 (r ≈ 33 au); 1 ORP ≈ 255 years.Simulations were followed through ∼ 20 ORPs, a time when the disk has settled into a quasi-steady thermodynamic state where cooling and heating are in balance and GI-induced structural and thermal properties are approximately constant.
The initial disk configuration that serves as the basis for these simulations was developed by Mejía et al. (2005); (see also Pickett et al. 1996Pickett et al. , 2003;;Mejía 2004;Mejía et al. 2005;Cai 2006).At t = 0 the initial state of the Mejía disk is esentropic with inner and outer radii of 2.3 and 40 au, respectively, and a surface density profile Σ ∝ r −1/2 within this radial range.The initial thermodynamic state of this disk was set using an equilibrium star plus disk model generated by a modified Hachisu (1986) self-consistent field relaxation method, where random ∼ 10 −4 cell-to-cell density perturbations were introduced to seed the growth of GIs.The Mejía t = 0 disk has defined the initial state of disks in a number of hydrodynamical studies (e.g., Mejía et al. 2005;Boley et al. 2006Boley et al. , 2007a;;Cai et al. 2006Cai et al. , 2008;;Michael et al. 2012;Steiman-Cameron et al. 2013;Desai et al. 2019).In these studies, the unstable disk passes through several phases of evolution during which the disk's mass distribution and thermal state are significantly modified.Ultimately, the disk settles into a quasi-steady, long-lived asymptotic state of sustained GI activity over a large part of the disk, with an approximate overall balance between heating and cooling.Details of the asymptotic phase and resultant asymptotic disk structures in these works are strongly dependent on the thermodynamical properties, i.e., heating and cooling, of the disk.
The Paper II simulations used the state of the Mejía disk as its initial condition, defined at t = 0 ORPs.The calculations reported here begin with initial conditions defined by the Paper II disk state at 7 ORPs, at the time when the disk is still in transition toward its asymptotic phase.The star's motion was not accounted for up until this point, but is in the calculations that follow.This is part of the reason that we allow a lot of time for transients to decay.Specifically, we follow the evolution of these disks through t = 20 ORPs for all l max but limit most of our analyses to t > 16 ORPs to allow transients to fully decay and the disk to fully settle into an asymptotic phase.As seen in the following sections, all four simulations settle into an asymptotic phase by ∼ 16 ORPs.

The Asymptotic Converged Disk
Our goal is to understand the "asymptotic converged disk configuration" of a gravitationally unstable PPD and how this drives disk evolution."Asymptotic" refers here to the evolutionary phase of a gravitationally unstable disk that has settled into a long-lived, quasi-steady balance between radiative cooling and disk heating at a relatively constant, but unstable, value of Q (sometimes referred to as "saturated GIs").Four simulations, each with differing azimuthal resolution, are followed to their own asymptotic states.The detailed configurations of these four asymptotic states are then used to establish grid convergence of the asymptotic disks and define the asymptotic converged disk configuration (hereafter ACDC).
Figure 1 shows mass densities in the disk midplane (top panel) and above the midplane along an azimuthal cut through the midplane (bottom panel) at t = 20 ORPs for all four simulations.As described below, by this time each simulated disk has settled into a quasi-steady state.Figure 2(a) shows azimuthally averaged radial surface density profiles, Σ(r), for each l max simulation at t = 20 ORPs.These profiles are very similar between ∼ 7.5 and 52 au and well described between 8 and 40 au by the exponential fit to the l max = 512 profile shown in the figure by the dashed line.Superposed on this general trend are persistent local maxima at ∼ 8, 11, and 14 au containing "excess" masses of ∼ 6M J , 18M J , and 10M J , respectively.As will be seen, these radii correspond with notable physical, kinematic, and thermodynamic ring-like features in the disk.Figure 2(b) displays the mass enclosed on cylindrical shells with thickness of one radial cell length (∼ 0.1667 au) as a function of radius at this same time for the 512 simulation.Broad bumps are visible in the mass distribution centered around ∼31-32 au and 48 au.These features fluctuate with time but are persistent.The location of the 31-32 au bump corresponds with a strong local maximum in time-averaged gravitational torques (Section 3.5).The 48 au bump arises from a strong one-arm spiral in the outer disk.Masses interior to 10, 20, 30, 40, and 50 au are approximately 5%, 20%, 40%, 62%, and 86% of the disk mass, respectively.
Detailed nonaxisymmetric density structures can readily be seen by subtracting the exponential fit of Figure 2(a) (dashed line) from the full 2D surface density distribution.The resultant enhanced spiral density structures are shown in Figure 3 at t = 20 ORPs.The structure is similar at all resolutions but, as expected, not precisely identical.With higher resolution, the fine structure of the ring region is more pronounced and spiral structures tend to be sharper and more pronounced as l max increases.New structure emerges as the resolution increases, but by 256 the global structure is roughly consistent.
A direct means of ascertaining convergence is provided by the characteristics of Q(r) with time and azimuthal resolution.A GI-active disk in its asymptotic state should be characterized by a quasi-steady Q(r). Figure 4(a) shows azimuthally averaged Q(r) for the l max = 512 simulation at t =12, 16, and 20 ORPs.Details describing how Q was evaluated are described in detail in Section 3.1 of Paper II.
The inner disk relaxes faster than the outer disk.At 12 ORPs, the 512 disk is still relaxing toward lower Q over most of its radial extent.By 16 ORPs, the disk has settled into an asymptotic state to ∼ 40 au.We find that the entire 512 disk has reached this state by ∼ 17 ORPs (not shown in the figure) as represented by Q(r) at 20 ORPs.At a given radius, Q displays ∼5-25% variability on the local dynamical timescale, with the amplitude of variability larger at larger radii.This was previously noted in Paper II (see Figure 4 of that paper).
The middle panel of Figure 4 shows Q(r) at 20 ORPs.By this time, all four simulations have relaxed into similar, but not identical, asymptotic Q(r) over the entire disk with clear convergence toward the l max = 512 curve.In what follows, we consider the state of the 512 simulation at t = 20 ORPs as representing the quasi-steady asymptotic state for the system studied here, and we will refer to it as the ACDC, as defined at the beginning of Section 3. As seen in the top two panels of Figure 4, two radial regions exist in the ACDC within which Q is essentially constant.Region 1, where Q ≈ 1.4, lies between ∼11 and 32 au while Region 2, characterized by Q ≈ 2.1, lies between ∼ 40 and 50 au.Between these regions, Q increases in a roughly linear fashion with radius.A local enhancement in Q between 11 and 14 au in Region 1 will be discussed in detail later.Q = 1.4 at the inner and outer edges of this feature, consistent with Q in Region 1. Between the two regions, Q(r) increases in a roughly linear fashion from ∼ 1.4 to ∼ 2.1.Beyond 50 au, Q rises sharply near the disk's outer edge as Σ(r) goes to zero.
Figure 5 shows face-on color-coded optical depths of the 512 disk at t = 20 ORPs with the inner and outer Q-defined Regions 1 and 2 depicted.The optical depth distribution is also shown in Figure 6, which presents azimuthally averaged optical depths (solid black line) for the ACDC disk (panel (a)) and the analogous disk of Paper II (panel (b)).Pixel representations in Figure 6 of individual azimuthal values of optical depth that are included in the average are displayed as blue pixels (less than average) and red pixels (larger than average).Inspections of these two figures show that Region 1 in the ACDC disk is essentially optically thick while Region 2 is optically thin.Between these regions the optical depth is a roughly equal mix of optically thick and thin cells.
The bottom panel of A.U.
A.U. Implementing a realistic, self-consistent, cooling approach in Paper II led Q(r) to increase by ∼ 0.5-1 relative to the values in the constant cooling simulations of Paper I. However, as described below, this was not due to the cooling approach itself, but rather to heating and cooling limiters implemented to maintain numerical stability.
In Paper II, only in regions interior to ∼ 18 au did the four simulations with differing l max settle to approximately the same Q(r) by 20 ORPs and the outer disk of the 512 model may not have fully relaxed by 20 ORPs.We also note that Figure 9  on the hydro calculation of that work.However, as discussed in the next section, the use of limiters on cooling rates appears to have had adverse consequences on the outcomes of Paper II.

Cooling Times and Temperatures
Thermal radiative cooling times t cool are calculated on each cylindrical shell using an averaging-like scheme.Specifically, the cooling time for the jth shell is given by which is the total internal energy ϵ in a cylindrical shell divided by the radiative energy-loss rate in that shell.The ∇ • F terms represent 3D radiative transport, as calculated in the code, but due to the cylindrical averaging, each t cool is based on vertical and radiative energy transport only.The cooling time is useful for characterizing the disk, but it is not used directly by the code in any of the radiative gravito-hydrodynamic calculations.Rather, we use t cool in this work primarily to compare the disk evolution here with the prescription used by Gammie (2001) discussed in Section 3.7, which uses a simple cooling time to parameterize energy loss.
The resulting normalized cooling times for the ACDC disk are shown in Figure 7, i.e., the cooling times given above, normalized to the local orbital period, P orb .In this case, the normalized cooling times are also time-averaged over 17-20 ORPs.For comparison, we also show the prescribed cooling times used in Paper I (constant cooling t cool = 2 ORPs, no star-disk interaction) and the calculated cooling times found in Paper II (limiters on cooling rates, no star-disk interaction).In the ACDC disk, the normalized cooling times rise rapidly in the inner disk, peaking at 33 around 11 au.Due to known uncertainties in the calculations interior to 8 au (see Section 2.2), the cooling times inside this region should be interpreted with caution.At radii larger than 11 au, the normalized cooling times fall steadily, with an approximately linear decrease between about 15 and 40 au.These radii include the optically thick Region 1 and most of the transition to the optically thin Region 2.
One might worry that the cooling times in the outer disk are problematically long.However, the divergence of the radiative flux includes the 3 K heating from the background.Normally, this is not relevant for the t cool calculations because the effective temperatures are usually much higher than this.But in Region 2, the average effective temperature of the disk drops below 5 K and approaches 3 K.Note that the "effective temperature" used here is really a brightness temperature determined from the outgoing radiative intensity in the vertical direction.It is only a measure of the radiation field, and not directly used in the simulation.It is nonetheless useful for understanding the behavior of t cool .In short, the temperature of the background radiation field contributes nontrivially to t cool in the outer disk.With this, a long t cool does not necessarily mean that the disk will require an equally long time to reach that state.
We can develop a more complete picture of the heating and cooling in the outer disk by instead looking at temperature and shock structures.We do not understand the "glitch" in the bottom panel that shows temperatures from Paper II, but note that this occurs at the same radius where the spread of optical depth about the mean at that radius undergoes an abrupt transition (see Figure 6).17 and 20 ORPs.There are noticeable spikes in temperature at 20 ORPs that are not seen in the average, because such temperature spikes cool relatively quickly.The figure also shows that the effective temperature exhibits the same behavior.Any given snapshot can have sudden variations due to the spiral arms, but those temperature variations do not persist away from the shocks.For further context, the 3 K background is met at 58.5 au.
For comparison, the bottom panel of Figure 8 shows the results from Paper II, which highlights that the ACDC disk is cooler over most of its radial extent.Because the surface densities are very similar between the ACDC and Paper II disks, the hotter temperatures in the Paper II disk explain the higher Q values shown in Figure 4. We suspect this difference is primarily from the use of limiters, as already discussed in section 2.2.
One final point regarding the cooling times is that approximately 20% of the computation cells at any moment have a negative ∇ • F , which represents heating for our sign convention.An example of the locations of these cells is shown in Figure 9.They are strongly correlated with the strong spiral waves seen in Figure 3 and occur more commonly in optically thicker portions of the disk.Their sharpness further suggests that these regions are associated with shocks, which are expected to have nontrivial radiative transport in all directions.
Transport in Gravitationally Unstable Disks 17

Nonaxisymmetric Structures
As demonstrated in Figure 4(b), the ACDC is subject to GIs between 10 and 50 au, a region encompassing ∼76% of the disk's total mass, although, based on their asymptotic Q values, we expect that these instabilities will be manifested differently in Regions 1 and 2.Although the region interior to 10 au is not unstable to GIs, it is still affected by nonaxisymmetric density structures generated by GIs in Region I.
Figure 10 shows the ACDC with spiral structures accentuated by plotting the difference between Σ(r, ϕ) and the exponential fit to azimuthally averaged Σ(r) shown in Figure 2 Visual inspection of Figure 10 shows more complex spiral structures in the inner optically thick Region 1 than those seen in the optically thin Region 2. Indeed, structures with up to sixfold symmetry can be seen in Region 1 of the figure, while only a one-arm spiral is readily apparent in Region 2. Steiman-Cameron et al.
To quantify these nonaxisymmetric structures, we examine the Fourier amplitudes The limiting azimuthal resolution for each l max is given by m = l max /2.GI-active disks display power at all resolvable m values in their asymptotic phase (e.g., Lodato & Rice 2004;Mejía et al. 2005;Boley et al. 2006;Cossins et al. 2009;Michael et al. 2012;Kratter & Lodato 2016).
For a disk-to-star mass ratio of the simulations presented here (m d /M * = 0.14), dispersion studies predict that the power spectrum for a disk subject to quasi-steady GIs will be dominated by relatively loosely wound, low-m waves (Lin & Shu 1964;Vandervoort 1970;Lau & Bertin 1978;Bertin 2000;Kratter & Lodato 2016;Béthune et al. 2021).
Figure 11(a) shows global time-averaged Fourier amplitudes ⟨A m ⟩ for each simulation integrated over 10 -50 au and averaged from 16 to 20 ORPs.As expected, the global power spectra are dominated by low-order Fourier components, and the spectrum continues to grow in amplitude at high m as l max increases.In this sense, while the spectra converge well by l max = 512, we cannot say that the spectra are converged.There is a residual uncertainty in our work in that we cannot be certain that, if the resolution were increased dramatically, we might actually see contributions from higher m to the resultant transport.
The inset in Figure 11(a) compares the lower-order m-terms ⟨A m ⟩ of the converged disk of this paper (solid line) with the corresponding converged disk of Paper II (dashed line).Because star-disk interactions were not included in Paper II, the m = 1 component of the azimuthal mass distribution in that work was not accurately treated, leading to an artificially large ⟨A 1 ⟩.Indeed, the global power spectra of Paper II are dominated by the m = 1 and 2 Fourier components, much different than seen here.This led to unrealistically large m = 1 torques in the Paper II converged disk.
Figure 11(b) shows power spectra of the converged disk plotted separately for Regions 1 and 2 and the global results of panel (a).As expected from the visual appearance of Figure 10, Region 2 is dominated by m = 1 with moderately strong m = 2. Region 1 has strong contributions from m = 1 to 5. For small asymptotic Q (Region 1), many Fourier components, i.e., azimuthal symmetries, are important.For larger asymptotic Q, only the lowest-order components appear dynamically important.
From the perspective of numerical convergence, discussed in more detail in Section 3.9, the spectrum of amplitudes for high m values in Figure 11(a) is not converging.However, the spectrum is converging for the lower-order m values.As shown in Section 3.5, these dominate in the production of gravitational torques.
We stress here that we have been examining time-averaged properties.As will be demonstrated below, time averaging hides dynamically important time variabilities that affect the disk.

Time Variability
While the time-averaged Fourier analysis above gives insight into time-averaged, nonaxisymmetric disk structures, the radial and temporal stabilities of the Fourier components are important in determining their dynamical effects on the disk.For example, radial incoherence in a Fourier component will diminish its gravitational effects and lead the component to have more importance on a local scale than a global scale.
Power in a specific A m does not imply the existence of an (eigen)mode for that m, nor does it necessarily represent the strength of an eigenmode that truly exists (e.g., Michael et al. 2012;Steiman-Cameron et al. 2013).For example, a disk with a single m = 2 eigenmode growing to nonlinear amplitudes will exhibit power at all even values of m, while a disk with two nonlinear m = 2 and 3 eigenmodes will exhibit power at all m values.We will avoid referring to Fourier components as modes except for those cases where we have determined that a mode exists.To this end, Figure 12 displays periodograms (Scargle 1982;Horne & Baliunas 1986;Mejía et al. 2005;Boley et al. 2006) of the converged disk for the m = 1-6 Fourier components during the same 16-20 ORP time frame as Figure 11.Locations of the corotation, inner Lindblad, and outer Lindblad resonances (CR, ILR, OLR) are displayed in each panel of Figure 12.Constructions of these periodograms use only the phase information ϕ m from the Fourier decompositions of ρ(r, ϕ) in the midplane.Power spectra of cos(ϕ m ) are generated for each r-value using a large number of time steps over the time range 16-20 ORPs.If cos(ϕ m ) is strictly periodic, i.e. ϕ m is linear in t, then there will be a strong spike at the corresponding pattern frequency.These periodograms for all radii are combined into one plot in which isocontours of spectral power are traced.If a pattern with a well-defined pattern frequency is present over a range of radii, it will produce a vertical stripe in the contour diagram.Periodograms only measure the coherence, not the amplitude, of patterns present.Strong phase coherence combined with significant amplitude at the same m value over the same radial range demonstrates that a dynamically significant m-armed wave is present.The m=2 mode with ILR, CR at (22,38) would have an OLR around 48 au but does not have power in the periodogram extending to the OLR.There are several other instances at all m ̸ = 1 with strong power between an ILR and CR radius that do not extend to the OLR.In these cases, the OLR would fall in Region 2.
There are a number of well-defined, densely packed, high-m stripes for CRs outside about 35 au and pattern frequencies ∼(1/ORP) or less.These are ignored here owing to confusion and the fact that many have pattern periods roughly comparable with the 4 ORP time frame used in constructing the periodograms.
The strong ring-like structure at 8 au lies at the ILR of an m = 2 mode, the 11 au ring lies at the ILR of an m = 2 mode and the innermost m = 4 mode, and the 14 au ring is at the ILR of both m = 3 and m = 4 modes.Durisen et al. (2005) previously noted similar ILR overlaps with positions of rings in their work on a hybrid theory of gas giant giant formation.
In addition to these well-defined modes, there are indications of transitory structures between the Lindblad resonances for a number of pattern frequencies and components.For example, several of these can be seen as faint vertical strips in the m = 4 panel of Figure 12 at pattern frequencies between 1.4 and 2.6 ORP −1 and spanning 10-30 au.Similar features are visible for several mvalues.These features represent density structures that have phase coherence for some period but less than the full time window of the periodogram.One might think of them as modes that pop into and out of existence.As will be seen, these ephemeral modes are dynamically important.
The discussion of gravitational torques and mass transport in what follows (Sections 3.5 and 3.6) will make it clear that there are strong m = 2 and 3 effects for periods close to where the red blobs are in Figure 12.Hence, even if there are no pure sustained modes, the disk seems susceptible to bursts of global m = 2 and 3 waves, probably swing amplified, and this is happening near pattern periods of about an ORP with CRs near 25-30 au.
The ephemeral modes are consistent with the description of Béthune et al. ( 2021) that gravitoturbulence generates spiral wakes that intermittently form and vanish over orbital timescales while, at the same time, large-scale spiral arms only manifest transiently through the coalescence of several neighboring wakes that are then are sheared apart.However, as will be shown below, in addition to the ephemeral modes, we find recurrent (also on an orbital timescale) swing-amplified bursts that correspond with the modes listed above.These bursts extend radially over the full ILR to OLR resonances.In short, these are coherent modal structures.

Mass and Angular Momentum Transport
Mass motions in accretion disks arise from stresses embodied in the stress tensor where T Rey , T grav , and T mag represent the stresses arising from hydrodynamic (Reynolds), gravitational, and magnetic interactions.Magnetic stresses fall outside the purview of this work and thus will not be considered further (see Deng et al. 2020;Béthune & Latter 2022).The Reynolds stress tensor is defined by where ρ is the mass density and u ′ r and u ′ ρ are the fluctuations in the radial and azimuthal field components, respectively.For the i-component of the velocity field, these fluctuations are defined by u ′ i = u i − u i , where u i is the instantaneous velocity and u i represents the "mean" (bulk) flow.Unfortunately, it is not easy to properly determine Reynolds stresses in 3D nonlinear global simulations because of difficulties inherent in evaluating the local mean fluid flow.In particular, the precise methods used to determine bulk flows are problematic.In Paper I, we attempted to measure Steiman-Cameron et al.
Reynolds stresses in similar 3D simulations using several different averaging schemes to evaluate the mean flow.We found that different approaches yielded dramatically different results and there were no obvious criteria for selecting one approach over another.Several previous global 3D studies have reported that Reynolds stresses are small relative to gravitational stresses (e.g, Lodato & Rice 2004;Boley et al. 2006;Michael et al. 2012;Steiman-Cameron et al. 2013;Bae et al. 2016;Béthune & Latter 2022).For these reasons, we omit Reynolds stresses in the calculations of angular momentum and mass transport.
The global torque, C, acting on a cylindrical section of the disk at radius r can be calculated by integrating the stress tensor T over the surface of the cylinder (Lynden-Bell & Kalnajs 1972), i.e., C = r × T dS. (7) Since the stress tensor includes only gravitational stresses, the surface integral in Equation 7 can be replaced with the volume integral where Φ is the gravitational potential of the disk.Here we are interested only in the z-component of torque, as only this component drives mass and angular momentum transport.The torque can be deconvolved into contributions from each Fourier term by replacing ρ in Equation ( 9) with the density distribution reconstructed from that Fourier component, i.e., ρ m = a ϕm cos(mϕ) + b ϕm sin(mϕ), (10) where a ϕm = 1 π ρ cos(mϕ)dϕ, b ϕm = 1 π ρ sin(mϕ)dϕ, and only the gravitational potential produced by the mass distribution given by ρ m is included in Φ.The total torque is then the sum of these torque components.
Figure 13 displays the total time-averaged gravitational torque, ⟨C Z (tot)⟩, and time-averaged torques summed over a number of low-order Fourier components, m 1 ⟨C Z(m) ⟩, m = 1 -4, for each of the four simulations.All torques are time-averaged over 240 equally spaced times between 17 and 20 ORPs to suppress short-timescale fluctuations.
With the exception of the l max = 64 simulation, time-averaged torques are dominated by several low-order (m ∼ 2 -6) components throughout the optically thick Region 1 and much of the transition region, with m = 2 -6 providing > 95% of the total torque and no one component dominating the time average.In the optically thin Region 2, m = 1 and 2 dominate.Some m = 1 strength may be due to beating of m = 2 and m = 3 but we suspect that most of its strength arises from sling amplification, a type of eccentric GI in nearly Keplerian disks (Shu et al. 1990;Ostriker et al. 1992;Kratter & Lodato 2016).
These results, where several low-order Fourier components dominate the optically thick regions and low-order terms dominate optically thin regions, are consistent with several other studies (see Section 3.3).The l max = 64 disk is distinctly different in that it is dominated by m = 3 torques for r ≤ 20 au, m = 2 torques in the r = 25 -35 au region, and m = 1 torques outside 40 au.

17-20 ORP Average
Gravitational Torque (x 10 39 ergs) Clear convergence toward the 512 disk is visible.Torque profiles of the 64 and 128 simulations show very clear differences with the 256 and 512 disks while the 256 and 512 torque profiles are essentially identical.
The results found here for the converged disk are similar to those found in the constant cooling study of Paper I but considerably different from those of Paper II.In Paper II, a very strong m = 2 mode exists and is the dominant torque over most of the disk with significant, but not dominant, contributions from m = 1 in selected radial regions centered around 28 and 38 au.We attribute most of the problematic results of Paper II to the use of cooling limiters in that study.The results in this paper using radiative subcycling supersede the results of Paper II.
Radial mass fluxes arising from the gravitational torques of Equation 9 can be written as (see also Balbus & Papaloizou 1999, Eq.  ulations.These radial mass flows correspond to evolutionary times ∼ 10 6 yr.It is important to note that there are several bands with different directions of radial flow at different ranges of r, very different from a simple α-disk. For comparison, the panel shows "measured" fluxes obtained directly from the output products by measuring how the mass within a cylindrical shell one radial grid element thick changes with time.Instantaneous Ṁ (r) are determined at the same 240 times used in determining the gravitational torques in panel (a) and then time-averaged.
Predicted and measured mass fluxes agree well for radii larger than ∼ 22 au but poorly at smaller radii.The disagreement interior to ∼ 22 au suggests that some stresses or transport terms are not properly accounted for, possibly Reynolds stresses.
Figure 15.Time variability of the mass in shells at six radii.At each radius (labeled in the upper left corner of the panels), the solid line shows the change in the cylindrical mass, ∆M cyl (t), measured relative to the mass at 14 ORPs.Cylindrical masses vary by ∼ 2-10% in a roughly cyclical manner on the local dynamical time.GI-active disks slosh around considerably more mass on short timescales than the net transport, which is indicated in the diagrams by the dashed linear fits.This sloshing can be seen in the animation associated with Figure 16.
We discuss below in Section 3.6, but note here, that the region in the inner disk where the red curve shows poor agreement is a region of strong torque outbursts and a region of strong shocks, which suggests strong systematic (not turbulent) flows in the spiral arms, which are Reynolds stresses when rotation is used as the "mean flow."We also note that masses on shells significantly change, in both the positive and negative sense, on local dynamical time scales or less, as shown in Figure 15.These short-timescale changes, in turn, follow longer-term trends.These can make measured dM/dt difficult to accurately assess.These variations are readily apparent in the animation of Figure 16, which shows how the radial mass distribution profile M cyl (r) changes with time.On timescales of the local dynamical time, the mass distribution in the disk displays radial oscillations giving the appearance that the disk "sloshes."This animation can also be seen at https://www.dropbox.com/s/q4p1a1lrfc0kp44/Fig14-Mass_on_Cylinders_Animation.mp4?raw=1 3.6.Short-timescale variability GI-induced mass transport arises from nonaxisymmetries in the mass distribution and hence nonaxisymmetries in the gravitational potential of the disk plus star.These nonaxisymmetries change continuously with time, as seen by changes in the geometry and strengths of spiral structures in the animation associated with Figure 10.Fluctuations in strength and coherence manifest themselves as the ephemeral modes seen in the periodogram of Figure 12 leading to the expectation that gravitational torques and mass flows also display considerable variations on short timescales.
Figure 17 shows instantaneous gravitational torques, including low-order Fourier components and total gravitational torques, at nine points in time between 18 and 20 ORPs, at intervals of 0.25 ORPs

Radius (AU)
Torque (10 39 ergs) Figure 17.Instantaneous gravitational torques.This figure is available as an animation.Panels are labeled with the times, in ORPs, for which the torques are plotted.The elapsed time between panels (0.25 ORPs) is approximately 64 yr.Total torques (the sum of all m-terms) and summations of low-m torque components (m = 1, 1+2, 1+2+3, and 1+2+3+4) are depicted.The key is the same as that used in Figure 14(a).As is clear from the figure, torque components and total torques are highly variable on short timescales.The animation runs from 17.00 to 19.93 ORPs in 0.012 ORP intervals.The real-time duration of the animation is 40 s.
(∼ 64 yr) with low-order torque components labeled the same as in Figure 14(a).An animation of Figure 17 can be accessed at: https://www.dropbox.com/scl/fi/uj8l9nimwes3lvxbul55r/Fig17-Low_m_torques_NEW.mp4?dl=0 Instantaneous torques vary dramatically on timescales less than the local dynamic time, often by 50% or more, and sometimes exhibit brief changes in the sign of local torque.Large local increases in the total torque are typically produced by sudden strengthening in an individual particular Fourier component, strongly suggesting recurrent bursts of swing amplification.
This can also be seen in Figure 18 which shows the instantaneous radial mass distribution, total z-torques, and torques from low-order Fourier components at nine points in time.Individual time frames are selected to show the wide range of burst activity seen at different times.Three vertical Steiman-Cameron et al.  frames are associated with each point in time.These show: (1) the mass enclosed within a cylindrical shell one radial cell width wide (∼ 0.1667 au), (2) azimuthally averaged total torques, and (3) torque contributions from m = 1, 2, ..., 6 Fourier components.The bottom right panel (t =19.32 ORPs) displays a time with relatively small torques over the full disk interior to 40 au with no Fourier component dominating the total torque.All other frames show strong bursts arising from one or more components, with individual component torques often changing sign.
The animation associated with Figure 18 shows that, with the exception of m = 1, bursts typically have durations on the order of the local rotation period as expected for swing amplification (m = 1 is not subject to swing amplification).These recurrent amplifications are readily apparent in Figure 19, which shows the time-progression of torques associated with the m = 1 -5 components at six radii Transport in Gravitationally Unstable Disks 31 between 10 and 45 au.An animation showing the time progression of torques and torque components between 17 -20 ORPs in the same format as Figure 18 is available and can also be found at https://www.dropbox.com/scl/fi/vs5i1odswfxrga7b1oqxq/Fig18-Animation.mp4?dl=0 Figure 20 overplots the total torque and torque contributions from m = 1 -5 components at 41 equally spaced intervals of ∼ 25.5 yr between 17 and 20 ORPs.Readily noticeable in Figures 17 -20 and associated animations is that total torques and torque components vary greatly over short timescales.Total torques are usually positive but occasionally reverse sign in a limited radial range as a result of large variability in low-order m contributions to the total torque.
Significant bursts in low-order (particularly m = 2, 3) gravitational torque components sometimes range from 11 to 30 au or more (see Figure 16 and the animation linked to that figure).Sudden increases in m = 1 and, to some extent, m = 2 and m = 3 torques sometimes range over a significant fraction of the disk's radial extent.Torque outbursts arising from higher-order symmetries, like m = 4 and 5, are more local, but each still shows radial extents of 6-30 au with the radial range of efficacy and extent changing with time, even on short timescales.An example of this is seen in the the upper rightmost set of panels in Figure 18 (17.85ORPs), where the blue curve for m = 4 shows an amplitude burst extending over more than 20 au.
Although instantaneous torques vary considerably from their time-averaged values, the overall characteristics of broad humps in the surface density near 9, 11, and 14 are preserved.

Effective α
It has long been recognized that molecular viscosity alone is not sufficient to account for angular transport in accretion disks (e.g., Durisen 2011).Some additional process must augment or dominate molecular viscosity.Shakura & Syunyaev (1973) proposed turbulence in the gas as the source of this viscosity.Assuming subsonic turbulence and an upper limit on the size of eddies, they proposed the ansatz ν = αc s H, where ν is the coefficient of viscosity responsible for carrying angular momentum outward, c s is the sound speed, and H is the disk scale height.The free parameter α provides the link between known quantities and turbulent viscosity.Given this equation, α and thus ν are locally defined.Following the development of Gammie (2001) (see also Lodato & Rice 2004;Boley et al. 2006;Michael et al. 2012), an effective α eff arising from gravitational stresses alone can be written as where Ω is the azimuthally averaged rotation speed, T grav rϕ is the gravitational stress tensor, and ρ is the volume density; angle brackets indicate azimuthal averages of vertically integrated stresses.The integrated divisor is also azimuthally averaged.
The gravitational stresses of Equation ( 13) can be evaluated in a straightforward manner using the gravitational torque of Equation ( 9), i.e., Using a 2D shearing box approach, which formally included both gravitational and hydrodynamic stresses and enforced balance between heating and cooling, Gammie formulated an effective α (see also Pringle 1981) expressed in 3D by Eq. ( 21) of Béthune et al. (2021) as where β = Ωt cool is the normalized local cooling time and γ is the adiabatic index.In this formulation, for fixed γ in a Keplerian disk, α eff is a function only of β and thus the effective viscosity depends only on the cooling time.As described in Section 2.1, for the temperature range in our simulations the gas is well approximated by an adiabatic index γ = 5/3.Figure ( 21) shows time-averaged α eff (r) for the converged disk, determined using the gravitational stresses of Equations ( 13) and ( 14), averaged at 240 equally spaced times between 17 and 20 ORPs along with a subset of 25 instantaneous measures of α eff (r) equally spaced in time.Because the m = 1 Fourier component of the mass distribution does not directly arise from GIs ( §3.5), we show time-averaged α eff (r) inclusive of all Fourier components and minus the m = 1 Fourier component.For comparison, local predictions of Equation ( 15) averaged over the same limiting times and using instantaneous β from the ACDC are shown.
Based on evolutionary lifetimes of protostellar disks, Hartmann et al. (1998) estimated α ≈ 0.01 at disk radii between 10 and 100 au.Depending on disk parameters, saturated GIs provide transport at rates such that 10 −2 < α eff < 1 (Kratter & Lodato 2016).Our time-averaged curve (solid blue line) falls between ∼ 10 −2.2 and 10 −1.6 , consistent with this result.In contrast, instantaneous α eff vary by greater than an order of magnitude over much of the disk and are, at times, considerably lower than the time-averaged values.This reflects the great deal of variability in the torques discussed in the previous section.The deep minima of instantaneous α eff seen in Figure ( 21) occur because the α eff are subject to sign reversals due to sign reversals in the torque.When the torques pass through zero, one sees deep minima in log(|α eff |).Hence, reader needs to realize that, although the averaged α(r) values look relatively smooth in Figure 21, there are several sign reversals leading to several radial ranges of negative α.This does not much resemble the picture of a simple α-disk.
ACDC torque-based α eff and the local dissipation predictions of Equation ( 15) track each other over the full radial range of the disk and are intertwined such that they precisely agree at six different radii and differ by up to a factor of two between these radii of agreement.Better agreement is seen in the outer disk when m = 1 torques are not included.
Between 16 and 25 au, the local dissipation curve is above the torque-based curve, suggesting that a stress or energy transport term is missing.Where does this energy go?It seems to heat the disk inside 16 au and the region outside 25 au.How important is it?It is about 50% of the heating by gravitational stress.
We suspect that this energy flow is due to Reynolds stresses caused by the correlated deviations from mean circular motion in global spiral modes.In other words, the low-order global spiral waves, especially m = 2-4, are transporting energy nonlocally.It is also possible that there is radiative radial transport caused by the large temperature gradients associated with shocks.These sources of energy flow are difficult to tease out of our current simulations.The spiral modes erupt in a chaotic manner over the whole low-Q part of the disk.The deviations of the red and blue curves are the combined effect over time of many spiral modes appearing and disappearing.This is truly gravitoturbulence, but dominated by fluctuating large-scale modes.It would be useful to study the Reynolds stresses and radial radiative transport due to erupting global spirals in more detail.Such an analysis goes beyond the scope of the current paper, but Figure 21 suggests that such a study might be interesting and fruitful.
We further note that 16, 25, and 42 au, radii where both local and torque-predicted α eff curves agree in Figure 21, are the same radii where torque-predicted and measured dM/dt = 0 in Figure 14.
Total torques at 30 au, where time-averaged torque-based α eff peak, never drop to zero and always are fairly large.This can be seen in the instantaneous curves of Figure (21).As seen in Figure 19, there are frequent m = 2, 3, and 4 outbursts centered around 30 au.In addition, as shown in Figure 16, a persistent "bump" is seen in the radial mass distribution at this same radius.

Ring-like Structures
Prominent and persistent ring-like structures are present at 8, 11, and 14 au in the converged disk (locations shown by dashed lines in the Figure 22).These rings contain "excess" masses of ∼ 6M J , 18M J , and 10M J , respectively, and correlate strongly with physical characteristics of the A. U.
Right edge = 10.18"Bottom edge = 3.48" 1.0" = 6.7 AU Ring AU Ring dia (inches) 3.8 1.13" 5.3 1.58" 7.9 2.36" 10.9 3.25" 13.9 4.15" Ring locations correspond with several low-order orbital resonances.The 8 au ring serves as the ILR of an m = 2 mode and also the inner Q-barrier.Torques arising from m > 1 structures are negligible interior to the 8 au ring (Figure 20), but this is not the case for m = 1.In the parlance of Durisen et al. (2005), the 8 au ring corresponds with the "active boundary ring" (ABR), because it occurs at the boundary between GI-active and GI-inactive regions and because it displays active nonaxisymmetric dynamics.The 11 au ring serves as the ILR of m = 2, 4, and 6 modes, while the 14 au ring serves as the ILR of m = 3, 4, and 6 modes, as well as the CR radius of m = 4 and m = 6 modes.
Similar ring features have been seen in previous 3D hydro simulations of self-gravitating PPDs carried out by our group and collaborators (Pickett et al. 1996(Pickett et al. , 2003;;Mejía 2004;Mejía et al. 2005;Durisen et al. 2005;Cai 2006;Boley et al. 2006Boley et al. , 2007a;;Cai et al. 2006Cai et al. , 2008;;Michael et al. 2012;Steiman-Cameron et al. 2013;Desai et al. 2019).These studies find that rings form early in the simulation, well before disks settle into their asymptotic state, and persist throughout the settled asymptotic state.More recently, a 2D study of an eccentric spiral instability in a self-gravitating disk with cooling by Li et al. (2021) found that a trapped one-arm instability forms early in the  2), Q (Figure 4), time-averaged gravitational torques (Figure 14), cooling times (Figure 7), and temperatures (Figure 8).Curves are rescaled and vertically offset to aid readability.
simulation and evolves into a set of axisymmetric rings.Their disk was not subject to GIs, but, on a deeper level, the dynamics involved may be similar.
The rings found in this study correspond with local time-averaged torque maxima.Because both Σ and temperatures are local maxima at ring locations, rings are also pressure maxima (see also Carrera et al. 2021).Assuming that these pressure maxima rings are, roughly speaking, in radial force balance, the rotation curve Ω(r) will be supra-Keplerian on the inner slope and sub-Keplerian on the outer slope of rings.We have confirmed that this is the case in the simulation.
Measured fluid flows in the 8 -14 au region that encapsulates the rings are complicated and nonsteady, with distinct radial flows and anticyclonic vortices around local clumps, as shown in Figure 24.Velocity vectors are shown with small red circles at their heads, and the general bulk rotation is counterclockwise.Large red circles represent the centers of two vortices.These vortices are almost certainly caused by the Rossby wave instability (Li et al. 2000(Li et al. , 2001)).Disk vortices have long been heralded as promising routes for planet formation owing to their ability to trap significant Figure 24.Velocity field at t = 20 ORPs between 9.8 and 11.8 au relative to the azimuthally averaged velocity at 10.8 au, the center of the 11 au ring.The radial appearance of the annular region has been stretched to facilitate visibility.Velocity vectors are shown with small red circles at their heads.The sense of bulk rotation for the ring is counterclockwise.Strong noncircular motions and vorticities, marked by large red circles, are prominent.

Grid Convergence and l max
Convergence in this study is quite good.The results cited above permit an assessment of appropriate azimuthal grids for the problem under consideration and similar problems.
The l max = 512 disk ("converged disk") simulation has clearly achieved grid convergence in the ϕ-direction and settled into an asymptotic phase before 17 ORPs for regions interior to ∼50-52 au.Most of the analyses presented in this paper have been performed in the 17-20 ORP time window.All four simulations settled into similar surface densities Σ(r) (Figure 2), with differences of 25% or less between 10 and 52 au in their thermodynamic states, as represented by their asymptotic Q(r) (Figure 4), where l max = 64 differs from the 512 disk by 15-30% and the 128 disk by roughly half that amount.In contrast, Q for the 256 disk varies from the 512 disk by an amount comparable to, or less than, the size of cyclical variations in Q(r) Gravitational torques for the 256 and 512 disks are dominated by m = 1 -6 Fourier components.Power spectra of these components (Figure 11) and the resultant time-averaged torques display very good convergence for this range of Fourier components, with minimal differences between the 256 and 512 simulations.Thus, for a quick examination of radially averaged Q(r), surface densities, and torques, a modestly low resolution of l max = 128 may suffice in disks dominated by only the lowest-order spiral structure.For most purposes, a simulation with l max = 512 appears to offer no material advantage over an l max = 256 calculation.
It has been argued that prompt fragmentation is more likely to occur as the resolution of 3D simulations increases.For the range of resolutions presented here, there is no evidence that this is the case.It must be noted, however, that the full time duration of the simulations reported here is only ∼ 5000 yr, with the converged asymptotic state followed for ∼ 800 yr.Therefore, depending upon one's definition of "prompt", these simulations cannot preclude or confirm the possibility that limited resolution (however that is defined) might prevent fragmentation in some instances.
Paper II and this work both examined the same disk using identical initial conditions and the same hydrodynamics code, albeit with two important modifications.Specifically, this paper includes the implementation of a subcyling approach to better control heating and cooling (Section 2.2) and the inclusion of an indirect potential approach to self-consistently account for star -disk interactions (Section 2.3).In contrast to Paper II, where convergence was lacking, great convergence is found here.The state to which a nonfragmenting disk settles is fairly sensitive to the accurate treatment of radiative losses, leading to higher (and more nearly constant) α values, constant Q over the optically thick disk (much in line with approximations used by others), higher Q' values in the optically thin regions (at least for our disk).The results of this work demonstrate the importance of doing the radiative physics well.In addition, the greater prominence of m = 1 here, compared with our earlier papers, indicates that one needs to include star-disk interaction (see also Section 6.5 of Elbakyan et al. 2023).

DISCUSSION
The existence of pressure maxima at ring centers may have important implications for planetesimal/planet formation.Weidenschilling (1977) was the first to point out that solid particles orbiting in a gaseous disk drift radially in the direction of a radial pressure gradient (see also Cuzzi et al. 1993).Haghighipour & Boss (2003a,b, hereafter HBa, HBb) subsequently studied the motions of small solids in the vicinity of local pressure enhancements of a gaseous nebula.They showed that the combined effects of gas drag and pressure gradients lead solids to accumulate at the locations where the pressure of the gas maximizes.
We can envision a path where solid particles accumulate in the rings owing to radial drift toward pressure maxima in rings.These particles can then be trapped by vortices within the rings which could, in turn, accelerate the growth of protoplanets.Durisen et al. (2005) used the results of HBa and HBb to estimate drift times in rings that appeared in an earlier disk simulation, leading them to suggest that even if instabilities due to disk self-gravity do not produce gaseous protoplanets directly, Steiman-Cameron et al.
they may create persistent dense rings that are conducive to accelerated growth of gas giants through core accretion.
The coincidence of strong resonances with rings, as described in Section 3.4, is an argument for the probable role of resonances with GI waves in ring formation (Durisen et al. 2005).Eccentric modes, corresponding to perturbations with azimuthal wavenumber m = 1, have also received particular interest in the context of PPDs because of their global nature.A large corpus of work has examined the development and sustenance of these modes in fluid disks.These have shown that almost any disk with a realistic density profile can sustain long-lived eccentric modes (Lee et al. 2019b) and that, once initiated, a global eccentric mode can grow its amplitude via the sling mechanism that amplifies an eccentric perturbation through the wobble of the central star and instantaneous cooling (Adams et al. 1989;Shu et al. 1990;Lin 2015).Ring formation then follows via angular momentum exchange with an unstable eccentric mode (Lubow 1991;Ogilvie 2001;Lee et al. 2019a,b;Li et al. 2021).
For all these reasons, we expect ring formation to be a common product of PPD evolution and it may play an important role in giant planet formation.

Rings and Closely Spaced Giant Planets
Given the discussions above and in Section 3.8, we speculate on a possible mechanism for forming closely spaced giant planets.The Nice model for the early dynamical evolution of the solar system (Tsiganis et al. 2005;Morbidelli et al. 2005;Gomes et al. 2005;Levison et al. 2008;Nesvorný & Morbidelli 2012;Brasser & Morbidelli 2013;Nesvorný et al. 2013) proposed the migration of the giant planets from an initial compact configuration into their present positions, long after the dissipation of the solar PPD.Among other things, it successfully explains the late heavy bombardment of the inner solar system, the formation of the Oort Cloud, and the existence of populations of small solar system bodies.A critical aspect of the Nice model is that the four giant planets were originally in much more closely spaced orbits than today.
At the earliest phases, PPDs have spiral waves.If the growth of ring-like enhancements in Σ(r) is a natural process and persists long enough, then the disk will have true rings left after GIs largely shut down.Extending this line of thought, if the azimuthal mass concentrations (clumps) along with associated sustained vortices in these rings persist, this could lead to the growth of multiple Jovian and/or ice giant planets within a relatively small radial extent.Interactions between these newly formed massive planets would then redistribute them into a more stable orientation and move icy planetesimals into the inner disk.
We note that issues may exist with this scenario, in particular, the delayed action between the ring and planet formation in closely spaced rings.The spacing of our rings may give highly unstable planets.Of course, with different disk parameters and following the whole rapid infall phase, any possible set of rings might look quite different.In addition, the structure of the inner disk could differ a lot if other methods for generation of turbulence were included in regions too hot to produce GIs.These are serious concerns and should be mentioned as caveats to this idea.
The Durisen et al. (2005) paper about a hybrid theory of planet formation attributes the rings in part to ILRs of high-m modes, and the authors proposed that heating of the inner disk might be due to dissipation in the innermost disk of GI waves generated in the central low-Q region of the disk.Alternately, the rings may be due instead to Rossby wave Instabilities, which in the nonlinear regime produce vortices, as shown by Li et al. (2021).In reality, perhaps both mechanisms operate.It is also interesting that broad surface density bumps appear and persists near 30 and 48 au in our simulation.Hence, a GI-active disk can generate long-lived radial structures in which solids may accumulate, not just near edges.

SUMMARY AND CONCLUSIONS
We present results of a 3D grid-based radiative hydrodynamics convergence study of a 0.07 M ⊙ protoplanetary disk subject to GIs surrounding a 0.5 M ⊙ star.The disk is evolved with a significantly improved radiative transport scheme using realistic dust opacities.This work supersedes the work of Steiman-Cameron et al. (2013).Both works examined the same disk, but here with important improvements to the hydrodynamics code, including the implementation of a subcyling approach to better control heating and cooling and the inclusion of an indirect potential approach to selfconsistently account for star -disk interactions that inevitably displace the star from the center of mass.
Our goals include determining cooling times experienced by the disk, characterizing the spiral density perturbations produced by the GIs and the gravitational torques produced by these structures, understanding disk processes on both time-averaged and instantaneous timescales (time averaging can hide dynamically important time variabilities that affect both the short-and long-term properties of the disk), evaluating the level to which transport can be represented as a local or nonlocal process, and understanding ring-like structures in the inner disk and their possible role in planet formation.
Four simulations were conducted and followed through to the time when the disks have settled into an asymptotic state where heating and cooling are roughly in balance.These simulations were identical except for the number of azimuthal computational grid points, thus allowing for calculations to establish mesh convergence in the ϕ-direction.
1.The primary messages of this work are that GI-active disks are awash in vigorous dynamic behavior, sloshing around considerably more mass on short timescales than the net transport, exhibiting rapidly changing gravitational torques with recurrent amplitude bursts.Edges of various sorts, e.g., optically thick to optically thin, variable cooling times, real physical edges, radially varying Q at those edges, rings, etc., are important.Edges affect many physical processes that vary in time and space and thus underlie the global nature of GI-active PPDs.
2. Accurate treatment of radiative cooling is critically important.With the heating and cooling limiters used in Paper II, the disk remained significantly too hot, leading to higher Q values, higher optical depths, and mostly incorrect cooling times.Moreover, there was no convergence in Paper II even with l max = 512.Results in this paper, based on subcycling of the energy equation for radiative cooling, are thermally smoother and converge quickly as the azimuthal resolution is increased.
3. All simulations settled into similar radial surface density Σ(r) and Toomre Q(r) profiles.Σ(r) is well fit by an exponential profile over most of the disk.Superposed on this general trend are persistent local ring-like maxima at ∼ 8, 11, and 14 au that contain "excess" masses of ∼ 6M J , 18M J , and 10M J , respectively.In addition, broad fluctuating but persistent bumps in the radial mass distribution are seen around 30 and 48 au.The former location corresponds with a maxima of the time-averaged gravitational torques.These torques serve to maintain the density maximum.The latter is associated with a strong one-arm spiral feature in the outer disk.
4. The ring-like features correspond with local torque and pressure maxima.Fluid flows in the region encapsulating the rings are complicated, with distinct radial flows and anticyclonic vortices around local clumps within in the rings.
5. Two distinct radial regions were found with essentially constant Q, with each displaying different convergence characteristics.Region 1 lies between ∼ 11 and 32 au and is defined by Q ≈ 1.4.Region 2 is bounded by 40-50 au and is characterized by Q ≈ 2.1.Between these regions, Q increases in a roughly linear fashion with radius.These two regions were initially identified and defined by the Q behavior only.Further analysis determined that Region 1 was optically thick, Region 2 was optically thin, and the transition between them contained the vertical τ = 1 transition region.
6. Region 1 spans the fully optically thick portion of the disk.Cooling times peak at the inner edge of the region with t cool /P orb ∼ 33 and decrease to ∼ 14 orbital periods at the outer limit of the region.Optically thick spiral waves, arising from GIs, are embedded in the optically thick background of Region 1. Shock-heated cells are common in locations that strongly correlated with the spiral waves.Torques in Region 1 are dominated by the m = 2 -5 Fourier components of the azimuthal mass distribution.
7. Region 2 spans the fully optically thin portion of the disk.Here only low-order Fourier components appear dynamically important.Optically thick spiral waves dominated by m = 1 with moderately strong m = 2 are embedded in the optically thin background of Region 2. Local cooling times are smaller than in Region 1, dropping from t cool /P rot ∼ 10 at the inner edge to 5 at the outer edge.
8. We find several low-m (2-6) modes in Region 1 and the transition region, with ILR ranging from 8 to 32 au and OLR ranging from 18 to 40 au.Gravitational torques arising from these modes vary in strength by up to a factor of 10 in a roughly cyclical manner.We attribute this to recurrent swing-amplified bursts that cycle on approximately the local dynamical time.This activity is particularly strong over the radial range between 8 and 18 au, a region dominated by ring structures between 8 and 14 au.In addition to these persistent well-defined modes, numerous transitory ephemeral modes are found.While these have well-defined orbital resonances (ILR,CR,OLR), they are not persistent throughout the simulation, but rather come and go.
9. The agreement between α values due to gravitational torques alone and expectations from measured cooling times is fairly good and gives values ∼ 10 −2 and hence mass transport rates ∼ ±10 −7 M ⊙ yr −1 and evolutionary time scales ∼ 10 6 yr.Deviations of a factor of two about this agreement suggest that some stresses or transport mechanism are not accounted for.What this and our earlier papers have shown, however, is that the disk behavior is not at all well characterized by a simple α-disk.This GI-active disk produce bands of inflow and outflow of mass, associated with significant persistent and/or episodically eruptive spiral structures.The disk is highly dynamic on a range of time scales in a nonlocal manner.

Figure 1 .
Figure 1.Volume mass densities in the midplane (top of each panel) and out of the plane along an azimuthal cut (bottom of each panel) at t = 20 ORPs for each of the four simulations.Panels are labeled with l max , the number of azimuthal grid elements used in the simulation.The color scale is logarithmic in code units, and axes units are given in au.

Figure 2 .
Figure 2. (a) Azimuthally averaged surface densities at t = 20 ORPs, measured in g cm −2 , as a function of radius for the l max = 512, 256, 128, and 64 simulations.The dashed line depicts the best-fit exponential to the 512 Σ(r) between 8 and 40 au.(b) Masses on cylindrical shells one radial cell width wide (0.1667 au) at t = 20 ORPs for l max = 512.Note ring-like features at 8, 11, and 14 au and broad bumps around 31-32 and 48 au.

Figure 3 .
Figure 3. Enhanced density structures at t ≈ 20 ORPs for the four simulations.Panels are 120 au on a side and labeled with l max , the number of azimuthal grid elements used in the simulation.Differences between the local surface density, Σ(r, ϕ), and the exponential fit to the azimuthally averaged surface density of Figure 2(a), Σ f it (r), are represented by color contours.The color scale is logarithmic in code units, and axis units denote au.

Figure 4 .
Figure 4. (a) Azimuthally averaged Q(r) for the l max = 512 simulation, as a function of radius at 12, 16, and 20 ORPs.By 16 ORPs, the disk has achieved a quasi-steady Q interior to ∼ 40 au and is near that state out to the disk's outer edge.By 20 ORPs, the disk has settled over its entire radial extent.(b) Azimuthally averaged Q(r) for the four simulations, at t = 20 ORPs.By this time, all simulations have converged toward the settled 512 Q(r) profile.(c) Azimuthally averaged Q(r) for the converged disks of Paper I (constant cooling), Paper II (no star-disk interaction, limiters on heating/cooling times), and the work reported here.
Figure 4 compares Q(r) for the ACDC with the corresponding states of Papers I and II.Interior to 30 au (Region 1), the ACDC Q(r) looks very similar to the constant cooling

Figure 5 .
Figure 5. Optical depths normal to the disk plane of the 512 disk at t = 20 ORPs.Colors show log(τ ).Dashed circles at 32 and 40 au delineate the outer and inner radii of Q-defined Regions 1 and 2, respectively.

Figure 6 .
Figure 6.Optical depths normal to the disk plane at t = 20 ORPs for all radial and azimuthal grid centers in the disk midplane.Solid black curves trace the azimuthally averaged optical depth as a function of radius, red pixels correspond with grid cells with larger than the average at that radius, and blue pixels have smaller than the average at that radius.Panel (a) shows the 512 disk of this work, and panel (b) shows the 512 disk of Paper II.Note that the optical depths shown in Figure 9 of Paper II are in error.Panel (b) has the corrected information and replaces that figure.

Figure 7 .
Figure 7. Azimuthally averaged cooling times normalized to the local orbital time for the ACDC disk of this work, Paper I, and Paper II. Cooling times are averaged between 17 and 20 ORPS.

Figure 8 .
Figure 8. Azimuthally averaged midplane and effective temperatures for the 512 disk of this work timeaveraged over 17-20 ORPs (top), the 512 disk of this work at t = 20 ORPs (middle), and the 512 disk of Paper II at t = 20 ORPs (bottom).We do not understand the "glitch" in the bottom panel that shows temperatures from Paper II, but note that this occurs at the same radius where the spread of optical depth about the mean at that radius undergoes an abrupt transition (see Figure6).

Figure 9 .
Figure 9. Locations of cells with negative divergent fluxes ∇ • F at 20 ORPs.These are strongly correlated with the locations of strong spiral waves seen in Figure 3.
(a).Outer and inner boundaries of Regions 1 and 2, respectively, are delineated in the figure.An animation showing the time-evolution of Figure10can be accessed at https://www.dropbox.com/s/7htomn38tb5vw0b/Fig9_Enhanced_Spirals.mp4?raw=1

Figure 10 .
Figure 10.Enhanced density structures in the converged disk at 20 ORPs.This figure is available as an animation showing the evolution of these structures between 17 and 20 ORPs.Dashed circles are added at 32 and 40 au to delineate the outer and inner radii of Q-defined Regions 1 and 2, respectively.Color contours show fractional differences between the local surface density, Σ(r, ϕ), and the linear fit to the azimuthally averaged surface density at that radius, Σ f it (r), shown by the dashed line in Figure 2(a).The color scale is logarithmic in code units, and axis units denote au.

Figure 11 .
Figure 11.(a) Strengths of the time-averaged global Fourier coefficients ⟨A m ⟩ integrated over 10-50 au for each azimuthal resolution.Limiting resolutions are given by m = l max /2.The inset compares ⟨A m ⟩ for the converged disk of this work (solid line with triangles) with the converged disk of Paper II (dashed with crosses).(b) ⟨A m ⟩ calculated separately for Regions 1 and 2 in the converged disk along with the 512 ⟨A m ⟩ from panel (a).

Figure 12 .
Figure 12.Periodograms for the m = 1-6 Fourier components of the converged disk during the 16-20 ORP time frame.Red lines denote the corotation, inner Lindblad, and outer Lindblad resonances.Zigzags in the ILR curves for m = 2 and 3 are due to the non-Keplerian rotational frequencies, Ω(r), near the inner rings.

Figure 13 .
Figure 13.Time-averaged gravitational torques as a function of radius for each of the four simulations.The total gravitational torque (m = All curve) is shown along with contributions to the total arising from loworder (m = 1 -4) Fourier components of the azimuthal mass distribution.These are shown as summations of the individual components.Torques are averages of instantaneous torques at 240 equally time-spaced times between 17 and 20 ORPs.

Figure 14 MFigure 14 .
Figure14(b)  shows mass fluxes predicted by applying this equation to the time-averaged instantaneous torques used in generating Figure14(a).These gravitational torque-predicted rates range between ∼ ±10 −7 ∼ M ⊙ /yr −1 , rates comparable to accretion rates reported in other global 3D sim-

Figure 16 .
Figure 16.Mass on cylinders.This figure, available as an animation, shows masses in 0.167 au wide cylindrical shells, M cyl .Each frame is labeled with the time in ORPs.The animation runs from 17.00 to 19.96 ORPs in 0.03 intervals.The real-time duration of the animation is 8 s.

Figure 18 .
Figure 18.Total torques, low-m torques, and masses on cylinders as a function of radius at six points in time.This figure is available as an animation.Note recurrent m = 2, 3, and 4 amplitude bursts, which we attribute to swing amplification, and bursts of m = 1 amplitude at larger radii, which we attribute to sling amplification.Individual frames in the still figure are selected to display a wide range of different torque activity more readily visible in the animation.The animation runs from 17.00 to 19.98 ORPs in 0.01 ORP intervals.Time is shown at the top of each panel.The real-time duration of the animation is 48 s.

Figure 19 .
Figure 19.Strengths of low-order Fourier component torques at selected radii as a function of time.Note recurrent amplifications of m = 2-5 torques at radii interior to ∼18 au.

Figure 20 .
Figure20.Overplotted instantaneous gravitational torques arising from m = 1-5 Fourier components and the total torque at 41 equally spaced times between 17 and 20 ORPs.This demonstrates how variable torques are on short timescales.It also illustrates how there are strong bursts in the torques, especially for m = 2 and 3 in Region 1 and m = 1 and 2 in Region 2.

Figure 21 .
Figure 21.Effective α for the converged disk.Light gray lines show instantaneous α eff at a subset of 25 equally spaced time steps between 17 and 20 ORPs calculated using Equations (13) and (14).The time average of instantaneous α eff is displayed by the solid blue line, while the dashed blue line displays the time average minus contributions from the m = 1 Fourier component.Predictions for α ef f from Equation (15), are shown by the thick red line.

Figure 22 .
Figure 22.The central regions of Figure 10 with locations of the 8, 11, and 14 au "ring" structures shown.Careful inspection reveals tightly wound spiral arms originating in the 11 au feature.This radius corresponds with the ILR of m = 2, 4, and 6 modes.