Particle Dynamics in 3D Self-gravitating Disks I: Spirals

Protoplanetary disks are observed with a number of features and substructures, from the common rings of recent surveys (Andrews et al., 2018; Long et al., 2018) to the less common but equally intriguing spirals (Grady et al., 2013; Pérez et al., 2016; Huang et al., 2018a; Johnston et al., 2020), vortices (van der Marel et al., 2013, 2016) and potential warps (Benisty et al., 2015). Spirals in the gas can be caused by two distinctive mechanisms, either by the wakes produced by the superposition of higher order modes induced by a point source (i.e. a planet) (Zhu et al., 2015; Dong et al., 2015; Bae & Zhu, 2018a, b) or by non-axisymmetric modes of gravitational instabilities (Lin & Shu, 1964; Dong et al., 2016; Lee et al., 2019). Distinguishing between these two sometimes relies on the pitch angle of the spiral arms (Pohl et al., 2015; Hall et al., 2018), the number of arms (Bae & Zhu, 2018a, b) or their amplitude (Juhász et al., 2015). We focus on the pitch angle for the ease in comparing and constraining this parameter with observations. While the number of arms is easier to constrain observationally, most disks only show two which is in line with theoretical expectations for either cause (Bae & Zhu, 2018a). Due in part to observational constraints, determining which causes the spirals in images, however, has proven to be a challenge (Forgan et al., 2018b, a; Ren et al., 2018).

Observations come in three regimes, those which are observed in light scattered off the optically thick gas surface of the disk (Muto et al., 2012; Grady et al., 2013; Garufi et al., 2013; Stolker et al., 2016; Benisty et al., 2015; Dong et al., 2018b, etc.), those which are observed as spirals in the dust emission (Pérez et al., 2016; Huang et al., 2018a; Reynolds et al., 2020) with ALMA and molecular line emission (Boehler et al., 2018; Booth et al., 2019; Huang et al., 2020). By only observing the surface of the disk, spirals in scattered light images may not indicate the presence of a feature throughout down to the disk midplane. In TW Hydrae for example, there exists an apparent discrepancy between the spirals in gas observations (Teague et al., 2019) and the rings seen in the dust continuum (Andrews et al., 2016). On the other hand, although spirals in the dust may probe the midplane of the disk, they may not reflect the spirals in the gaseous disk (Pérez et al., 2016).

Gravitationally unstable disks are commonly characterized through the presence of spiral arms due to the sheared self-gravitating waves generated by the massive disk (Lin & Shu, 1964; Goldreich & Lynden-Bell, 1965). These are typically represented by the dominant $m=2$ spiral mode, but higher order spirals are common in simulations depending on disk mass and disk aspect ratio (Cossins et al., 2009; Kratter & Lodato, 2016). Comparatively, a spiral formed by a planet will feature a dominant $m=1$ exterior mode, with higher order interior modes, depending on the mass of the perturber (Bae & Zhu, 2018a). Thus a planet that induces multiple spirals are almost always considered to be located further out than the observed spirals.

However, recent images of potentially self-gravitating disks reveal either no signatures of spirals (Andrews et al., 2018) or only faint spirals (Pérez et al., 2016) and disks with spirals are not conclusively in the regime where self-gravitational effects should be relevant (Huang et al., 2018b). Similarly, most disks which show spiral features also do not have obvious planet candidates (Boehler et al., 2018; Dong et al., 2018a), although this may be due to detection limitations of low luminosity planets (Ren et al., 2018).

Thus, in the absence of more robust indicators, the characterization of the pitch angle perhaps makes it possible to disentangle the cause of spiral features. In the case of an embedded planet, the angle of the spiral is influenced by the mass of the planet as well as the viscosity and local temperature (Rafikov, 2006; Muto et al., 2012). Overall, the pitch angle is greater close-in to the planet and decreases further away from the planet. For studies which delve into the properties of the spirals generated by embedded planets we refer to Zhu et al. (2015); Meru et al. (2017); Hall et al. (2018); Dong et al. (2018b, b); Bae & Zhu (2018b).

The spiral structure of self-gravitating disks are often studied for their gas dynamics (Cossins et al., 2009; Michael et al., 2012), but the effect of self-gravity with respect to solid material is an often neglected part of mesh-based simulations, due in part to the computational expense and because of the minor effect in low-mass disks. Without considering the self-gravity effects from the gas onto the dust and between the dust particles themselves, larger dust species will remain decoupled from the gas and form axisymmetric structures, and this may be the expected outcome in many cases (Cadman et al., 2020). However, as we will show, when the effects of the gas self-gravity is included on the dust, even large species will form non-axisymmetric structure similar to the gas.

In this paper, we investigate how the simulation parameters of domain size and cooling timescale affect the pitch angle of gas and dust spiral features in local self-gravitating disks. In Section 2 we outline the formation of spirals in self-gravitating disks and then describe the numerical setup in Section 3. We then present the results in Section 4 and discuss the implications and limitations in Section 5 before concluding in Section 6.

We use 3D hydrodynamic shearing box simulations to study the dynamics of self-gravitating disks with Lagrangian super-particles embedded in the Eulerian mesh using the Pencil code (Brandenburg, 2003). Pencil is a sixth-order in space and third-order in time finite difference code with MPI for high parallelization. The code also has robust routines for self-gravity and particle-gas interaction and has been thoroughly tested for a number of astrophysical contexts (Youdin & Johansen, 2007; Lyra et al., 2007; Yang & Johansen, 2016).

For our simulations we require thermal and hydrodynamic properties such that the disk is marginally Toomre stable, i.e. $Q>1$. That is to say, for the gravitational constant $G$ and $\Omega$ both equal to one, the sound speed $c_{s}$ and combined gas and dust surface density $\Sigma_{\mathrm{0}}=\Sigma_{\mathrm{g,0}}+\Sigma_{\mathrm{p,0}}$ are chosen such that the initial $Q_{\mathrm{0}}=1.02$ everywhere.

Local simulations allow the Toomre wavelength²²2Also known as the largest unstable wavelength $\sim 2\pi H$ to be well-resolved in the radial and azimuthal directions ($x$ and $y$ in the Cartesian coordinates, respectively) and keep the boundary conditions periodic. The Toomre wavelength needs to be resolved by at least four grid cells to avoid spurious errors growing into larger, unphysical overdensities (Truelove et al., 1997; Nelson, 2006). We easily satisfy this by resolving each pressure scale height by at least 10 grid cells and thus the Toomre wavelength by around 60 grid cells.

The shearing box simulations used here employ hydrodynamic equations for density, momentum and energy conservation which are linearized and transformed into co-rotating Cartesian coordinates, where $q=-\mathrm{ln}\Omega/\mathrm{ln}R=3/2$ is the shear parameter:

In equations (18) - (20), $\mathbf{u}=(v_{\mathrm{x}},v_{\mathrm{y}}+q\Omega x,v_{\mathrm{z}})^{T}$ is the gas flow plus shear velocity in the local box, $\rho_{\mathrm{g}}$ is the gas density, $\Phi$ is the gravitational potential of the gas and dust, and $s$ is the gas entropy. The gas and dust are vertically stratified with a sinusiodal gravity profile to avoid a discontinuity at the vertical boundary (cf. Baehr et al., 2017). Viscous heat is generated by the term $2\rho_{\mathrm{g}}\nu\mathbf{S}^{2}$, with rate-of-strain tensor $\mathbf{S}$. Hyperdissipation is applied with the terms $f_{D}(\rho_{\mathrm{g}})$, $f_{\nu}(\bm{u})$, $f_{\chi}(s)$ which for each has the form

with constant $\nu=2.5\,H^{6}\Omega$ (Yang & Krumholz, 2012). Heat is lost through $\beta$-cooling prescription

with $t_{\mathrm{c}}$ given by $t_{\mathrm{c}}=\beta\Omega^{-1}$ and background irradiation term $c_{\mathrm{s,irr}}^{2}$. This background term is set to the initial sound speed so that cooling does not bring the local temperature too close to zero. This cooling prescription has no dependence on the optical depth and thus all regions cool with the same efficiency. In reality, the opacity will be dominated by the small dust grains and an increase of the particle density will increase the local cooling timescale.

Self-gravity of the gas and dust is solved in Fourier space by transforming the density to find the potential at wavenumber $k$ and transforming the solution back into real space. The solution to the Poisson equation in Fourier space at wavenumber $\bm{k}=(k_{x},k_{y},k_{z})$ is

where $\Phi=\Phi_{\mathrm{g}}+\Phi_{\mathrm{d}}$ and $\rho=\rho_{\mathrm{g}}+\rho_{\mathrm{d}}$ are the potential and density of the gas plus dust particles combined. This is the default setup for the first six simulations in Table 1 and we distinguish these six from the next three by the inclusion of gas-particle and particle-particle gravitational interactions self-consistently. Among these three simulations, the simulation with ’low dust mass’ indicates that the particle mass is so low that it does not contribute to the gravitational potential $\Phi=\Phi_{\mathrm{g}}$ and thus every particle does not feel the gravity of other particles while still feeling the gravitational acceleration from the gas potential. For the simulation indicated as ’no particle self-gravity’, the particles neither contribute to the potential nor feel the gas potential and thus the particles only feel the aerodynamic drag force, the Shearing box inertial forces, and the gravitational force from the central star.

Finally, we use an ideal equation of state, with internal energy $\varepsilon$, and specific heat ratio $\gamma$

with $\gamma=5/3$.

Particles are initially placed at rest with a vertical Gaussian distribution but are otherwise randomly distributed in $x$ and $y$ directions. The $1.5\times 10^{6}$ particles are evenly divided amongst six sizes $\mathrm{St}=[0.01,0.1,1,10,100,1000]$ in all simulations. The equations of motion for the particles are then

where $t_{s}$ is the particle stopping time, which we write normalized by the dynamical time $\Omega^{-1}$ as the Stokes number

In the Epstein regime, the stopping time is related to the particle size as in (Weidenschilling, 1977)

where $a$ is the particle diameter and $\rho_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}$ is the material density of a dust particle. A super-particle (aka swarm-particle) is a collection of identical dust particles that do not move with respect to the others in the swarm. Particle drag is calculated by interpolating particle positions to the grid using a second-order triangular-shaped cloud scheme (i.e. Yang et al., 2018). Larger particles are less coupled to small scale gas motions, retain their initial perturbations for longer and thus have higher Stokes numbers. On the other hand, smaller particles are well-coupled and will quickly match the velocity of the gas motions in the vicinity and thus have low Stokes numbers. Particle mass is calculated from the gas through the initial condition of the dust-to-gas ratio $\epsilon=\rho_{d}/\rho_{g}$.

Our full collection of simulations and relevant parameters is summarized in Table 1. We compare two different box sizes, $L_{x}=L_{y}=(80/\pi)H$ and $(160/\pi)H$, which we label medium and large respectively. Both sizes have the same vertical extent $L_{y}=(40/\pi)H$, but the effective resolution of the larger simulations is half that of the medium ones. Additional simulations at a higher grid resolution were run with $\beta=20$ and $\beta=30$ in order to explore trends to broader disk conditions.

Figure 1 shows the vertically integrated particle count for each of the particle species included in two of our simulations after gravitoturbulence has been established. All species trace the gas to some extent, matching the same features of the gas with particularly narrow structures in the case of $\mathrm{St}=1$. Since particles of this size should drift most efficiently towards pressure maxima, this is expected. This also compares well with similar simulations in Gibbons et al. (2012), particularly for particle sizes $\mathrm{St}=1$ and smaller (see their Figure 2). In Gibbons et al. (2012) however, particles of size $\mathrm{St}=10$ and $100$ show some axisymmetric structure, but do not correspond to the gas features, suggesting more inertial behavior. Discrepancies at these larger sizes can be attributed to two differences between how particles are treated. In our primary simulations, we include the effects of self-gravity of and on the particles, which facilitates dust concentration, and the 3D implementation of these simulations which considers sedimentation effects which proceed faster for moderate Stokes number particles, $\mathrm{St}=0.1$ and $\mathrm{St}=1$. The effects of self-gravity will be discussed later in this section, and the effects of settling will be discussed in Paper II.

It is from these particle distributions and the vertically integrated gas that we calculate the pitch angle $\theta$ in each simulation. To determine the average opening angle, we calculate the autocorrelation of the gas and particle surface densities separately (cf. Gammie, 2001; Michikoshi et al., 2015)

Thus each point on the autocorrelation map, $(\tilde{x},\tilde{y})$ represents a distance away from each physical point $(x,y)$ on the 2D density map, all superimposed. This allows for a characterization of the overall direction and shape of density structures at a given point in time. Calculating the autocorrelation on a $512^{2}$ grid was computationally expensive, so all density data at that resolution was interpolated to a $256^{2}$ grid for calculating the autocorrelation over the entire domain. To calculate the autocorrelation when the distance $(\tilde{x},\tilde{y})$ crosses a shear periodic boundary at $x=-L_{x}/2,L_{x}/2$ (cf. Hawley et al., 1995), the density is shifted in the azimuthal direction by $\Delta y=(3/2)\Omega x$ so that features are aligned correctly.

From the autocorrelation, the angle with respect to the azimuthal axis was measured via the ridge of maximum $\xi(\bar{x},\bar{y})$ tracing the diagonal feature, which we mark with a dotted line. The particle autocorrelation is typically less smooth, and thus a Gaussian smooth was applied to the data before determining the line which traces the diagonal feature. This dotted line is then fit with a straight line and the angle $\theta$ is defined as that opened counterclockwise with respect to the azimuthal ($y$) axis.

We focus our efforts at a time where gravitoturbulence has been established in the simulation, such that non-axisymmetric (non-linear) features have been allowed to evolve over at least a few dynamical timescales. In Figures 2 and 3, we show the gas and total dust autocorrelation calculations at $t=60\,\Omega^{-1}$. Overall, gas pitch angles $\theta_{g}$ are within a few degrees of $\theta=12^{\circ}$, with an apparent decreasing trend with increasing $\beta$ for the larger simulations.

The particle autocorrelation is partially complicated by the dust clumping due to particle drift and particle-particle self-gravity at high particle concentrations. Numerous bright spots make it hard to find the overall dust pitch angle, particularly for $\mathrm{St=1}$. However, there is less clumping at other sizes such that the overall dust picture is not disrupted significantly.

We further break down the autocorrelation analysis by particle size. Figure 4 shows two medium-sized simulations at $t=60\,\Omega^{-1}$. Particles with size $\mathrm{St}=1$ are especially prone to clumping, such that densities along filaments are higher and the autocorrelation shows a starker contrast. This makes determining the pitch angle at this species more difficult than the rest and this is reflected in the larger variation of pitch angles at this size. Overall, the spirals shown by different sized particles have surprisingly consistent pitch angle. The pitch angle is almost unchanged even if the $\mathrm{St}$ of the dust changes by 5 orders of magnitude.

In Figure 5, we plot the radial-azimuthal distribution of particles from two simulations identical to our initial six simulations, aside from how particle self-gravity is treated. The left panel includes particle density as part of the potential, but the particles only have $10^{-4}$ times the mass. Thus, they contribute little to the overall potential, even at high local concentrations. However, the particles still feel the acceleration of the gas potential and thus are still affected by dense gas structures. As the autocorrelations show in the left side of Figure 6, particles still retain effectively the same non-axisymmetric structure at all sizes, although it does not perfectly correlate to the gas.

In the other case, shown in the right panel of Figure 5, particle density is not included in calculating the self-gravity potential and the acceleration of the potential is not included in Equation (25). In this situation, particle-gas interaction is through the aerodynamic drag only. Thus, larger dust particles become increasingly axisymmetric with increasing dust size ($\theta$ approaching $0^{\circ}$), until they are randomly distributed at $\mathrm{St}=1000$. These axisymmetric features are similar to the results shown in Gibbons et al. (2012). On the other hand, the left panels in Figure 5 are consistent with previous 2D simulations by Shi et al. (2016) who have shown spiral features for large particles. By comparing the left and right panels of Figure 5 and 6, we identify that non-axisymmetric structure can occur for $\mathrm{St}>1$ since the self-gravity of the gas is acting on the dust.

We investigated the spiral arms in self-gravitating disks and analyzed a number of factors which could affect what form they take and how they are observed. Using local 3D shearing box simulations of gravitoturbulent disks we focused on the influence of the gas cooling rate and simulation domain and how different sized particles presented in each case. Self-gravitating disks are a possible cause of spiral features which have been observed in both gas and dust disks.

In our simulations, We characterize their features in not only gas properties, but in the dust as well, which allows for a better comparison with observations. We summarize our results in the following points:

• •

  The opening angles of gas spirals in GI disks are universal ($\sim 10^{o}$), which are not significantly affected by the size of the computational domain. In our simulations with the biggest domain size, the gas spirals become slightly more open with increased cooling efficiency.

• •

  Particles of different sizes have similar responses to the gas structure as measured by the pitch angle when gas self-gravity acts on the dust. While it is expected that smaller particles are more coupled and have pitch angles that are more similar to the gas, we find that this is also true for the larger particles which are less affected by aerodynamic drag but more affected by the gas gravity.

• •

  The above point suggests drag and self-gravity may be comparable influences on the shape spirals take. When particles do not feel the potential of the gas at all, the ability of large particles to form features corresponding to the gas is drastically diminished. This suggests that self-gravity plays a role in shaping the structure spirals in gravitoturbulent disks, particularly for larger particles. Even very low mass particles behave very similar to more massive particles, but will not collapse into small dense clouds unless particle-particle self-gravity is included.

• •

  While large particle species show very broad non-axisymmetric features, they become more refined with decreased cooling efficiency (increasing $\beta$). This occurs even for a vertically unsettled distribution of particles which are still undergoing overdamped oscillations about the midplane.

Overall, both gas and dust at different sizes show spirals with $\sim 10^{\circ}$ pitch angle (almost all from 6^∘ to 14^∘ ), independent on the simulation domain size. This pitch angle is roughly consistent with current protoplanetary disk observations. We observed some weak trends regarding particle sizes and cooling time. On the other hand, considering the global nature of GI, it is desired to investigate both gas and dust spirals in global simulations in future.