Dust dynamics in 2D gravito-turbulent discs

Shi, Ji-Ming; Zhu, Zhaohuan; Stone, James M.; Chiang, Eugene

doi:10.1093/mnras/stw692

Abstract

The dynamics of solid bodies in protoplanetary discs are subject to the properties of any underlying gas turbulence. Turbulence driven by disc self-gravity shows features distinct from those driven by the magnetorotational instability (MRI). We study the dynamics of solids in gravito-turbulent discs with two-dimensional (in the disc plane), hybrid (particle and gas) simulations. Gravito-turbulent discs can exhibit stronger gravitational stirring than MRI-active discs, resulting in greater radial diffusion and larger eccentricities and relative speeds for large particles (those with dimensionless stopping times t_stopΩ > 1, where Ω is the orbital frequency). The agglomeration of large particles into planetesimals by pairwise collisions is therefore disfavoured in gravito-turbulent discs. However, the relative speeds of intermediate-size particles (t_stopΩ ∼ 1) are significantly reduced as such particles are collected by gas drag and gas gravity into coherent filament-like structures with densities high enough to trigger gravitational collapse. First-generation planetesimals may form via gravitational instability of dust in marginally gravitationally unstable gas discs.

hydrodynamics, turbulence, methods: numerical, planets and satellites: formation, protoplanetary discs

1 INTRODUCTION

As protoplanetary discs are often thought to be turbulent (Armitage 2011; but see Flaherty et al. 2015), understanding how disc solids interact with turbulent gas is crucial to modelling the formation of planetesimals and planets (Weidenschilling & Cuzzi 1993; Chiang & Youdin 2010; Johansen et al. 2014; Testi et al. 2014) and to explaining observations of discs (Williams & Cieza 2011; Andrews 2015). Turbulence determines the spatial distribution of solid particles and their relative collision velocities (‘turbulent stirring’; Voelk et al. 1980; Cuzzi, Dobrovolskis & Champney 1993; Ormel & Cuzzi 2007; Youdin & Lithwick 2007). For example, dust particles can be concentrated in local pressure maxima at the interstices of turbulent eddies; pressure bumps can also be found in spiral density waves, anti-cyclonic vortices, or zonal flows associated with whatever mechanism drives disc transport (Maxey 1987; Klahr & Bodenheimer 2003; Rice et al. 2004; Barranco & Marcus 2005; Johansen, Youdin & Klahr 2009; Mamatsashvili & Rice 2009; Johansen, Klahr & Henning 2011; Pan et al. 2011). Turbulent density fluctuations can also exert stochastic gravitational torques on solid objects and alter their orbital dynamics (‘gravitational stirring’; Laughlin, Steinacker & Adams 2004; Nelson 2005; Ogihara, Ida & Morbidelli 2007; Ida, Guillot & Morbidelli 2008).

Disc self-gravity can drive turbulence, provided that discs are sufficiently massive and their cooling times are longer than their dynamical times (Paczynski 1978; Gammie 2001; Forgan, Armitage & Simon 2012; Shi & Chiang 2014). ‘Gravito-turbulence’ may characterize the early phase of star formation, when discs are still massive (Eisner et al. 2005; Rodríguez et al. 2005; Andrews & Williams 2007). Observations show signs of early grain growth in some very young stellar systems (Ricci et al. 2010; Tobin et al. 2013). Millimeter-sized chondrules in primitive meteorites indicate that they might once have been melted via strong shock waves in self-gravitating discs (Cuzzi & Alexander 2006; Alexander et al. 2008). Simulations suggest large dust concentrations via spiral waves and vortices present in gravito-turbulent discs, possibly leading to planetesimal formation via gravitational instability in the dust itself (Rice et al. 2004; Gibbons, Mamatsashvili & Rice 2015).

Many studies of dust dynamics to date focus on particles in discs made turbulent by the magneto-rotational instability (MRI; Balbus & Hawley 1991, 1998). Both analytical and numerical works have been carried out to obtain radial/vertical diffusivities and particle relative velocities due to turbulent stirring (Cuzzi et al. 1993; Ormel & Cuzzi 2007; Youdin & Lithwick 2007; Carballido, Cuzzi & Hogan 2010; Carballido, Bai & Cuzzi 2011; Zhu, Stone & Bai 2015). For large particles, gravitational forces by MRI-turbulent density fluctuations exceed aerodynamic drag forces by gas and generate relative particle velocities too high to be conducive to planetesimal formation (Nelson 2005; Johnson, Goodman & Menou 2006; Yang, Mac Low & Menou 2009, 2012; Nelson & Gressel 2010; Gressel, Nelson & Turner 2012). A few useful metrics common to many of these papers include as follows: (1) the diffusion coefficient, which characterizes how quickly solids random walk through the disc; (2) the particle eccentricity or velocity dispersion; and (3) the pairwise relative velocity, which is crucial for determining collision outcomes. Quantity (2) usually serves as a good proxy for quantity (3).¹

Relatively, fewer groups investigated the dynamics of dust in turbulence driven by disc self-gravity. Gibbons, Rice & Mamatsashvili (2012), Gibbons, Mamatsashvili & Rice (2014) and Gibbons et al. (2015) study particles with a range of sizes (with stopping times t_stop = 10⁻²–10²Ω⁻¹, where Ω is the local orbital frequency) accumulate in local, two-dimensional (in the disc plane) simulations. They find that intermediate-sized dust can concentrate by up to two orders of magnitude, and that the dispersion of particle velocities can approach the gas sound speed, consistent with the results of 2D global simulations (Rice et al. 2004, 2006). Britsch, Clarke & Lodato (2008) and Walmswell, Clarke & Cossins (2013) find strong eccentricity growth for large-sized planetesimals forced more by gravitational stirring than by gas drag. Boss (2015) investigates the radial diffusion process for particles 1 cm–10 m in size (t_stop ∼ 10⁻²–1Ω⁻¹ in their model), finding enhanced diffusion for m-sized or larger bodies.

However, no systematic study has yet been performed to directly measure the dynamical properties (diffusivities, eccentricities, and relative speeds as listed above) of solids in gravito-turbulent discs as has been done for MRI-active discs. Gravito-turbulence tends to produce relatively stronger density fluctuations (δρ/ρ ∼ 1 for a typical Shakura–Sunyaev turbulence parameter α ∼ 10⁻²; see Shi & Chiang 2014) than are seen in MRI turbulence (δρ/ρ ∼ 0.1 for α ∼ 10⁻²). The prominent spiral density features that characterize self-gravitating discs and that help trap dust particles are also absent in MRI-turbulent discs.

It is the goal of this paper to study the dynamics and spatial distribution of dust in gravito-turbulent discs in a systematic manner, placing our measurements into direct comparison with analogous measurements made for MRI-turbulent discs. We first describe our simulation setup in Section 2. Results are given in Section 3, where we describe the radial diffusion, eccentricity growth, and relative velocities of particles, and how these quantities are affected by gravitational stirring, particle stopping time, gas cooling rate, numerical resolution, and simulation domain size. In Section 4, we put our results into physical context, discuss their astrophysical implications, and make comparison with MRI-active discs. We conclude in Section 5.

2 METHODS

2.1 Equations solved and code description

We study the diffusion of solids in gravito-turbulent discs using hybrid (particle+fluid) hydro simulations in the disc plane.

For the gas, we solve the hydrodynamic equations governing 2D, self-gravitating accretion discs, including the effects of secular cooling. The disc is modelled in the local shearing sheet approximation assuming that the disc aspect ratio H/r ≪ 1. In a Cartesian reference frame corotating with the disc at fixed orbital frequency Ω, the equations solved are similar to those Shi & Chiang (2014), but restricted to be in the disc plane:

\begin{eqnarray} \frac{\mathrm{\partial} \Sigma _{\rm g}}{\mathrm{\partial} t} + \nabla \cdot (\Sigma _{\rm g}\boldsymbol {u}) = 0 \,, \end{eqnarray}

(1)

\begin{eqnarray} &&\frac{\mathrm{\partial} \Sigma _{\rm g}\boldsymbol {u}}{\mathrm{\partial} t} + \nabla \cdot \left(\rho \boldsymbol {uu} +P\boldsymbol {I} + \boldsymbol {T_{\rm g}} \right)\nonumber \\ &&\quad=\, 2q\Sigma _{\rm g}\Omega ^2 x \hat{\boldsymbol {x}} -2\Sigma _{\rm g}\Omega \hat{\boldsymbol {z}}\times \boldsymbol {u}\,, \end{eqnarray}

(2)

\begin{eqnarray} &&\frac{\mathrm{\partial} E}{\mathrm{\partial} t} + \nabla \cdot (E+P)\boldsymbol {u}\nonumber \\ &&\quad =\, -\Sigma _{\rm g}\boldsymbol {u}\cdot \nabla \Phi + \rho \Omega ^2\boldsymbol {u}\cdot \left(2 q x \hat{\boldsymbol {x}} - z\hat{\boldsymbol {z}}\right) -\Sigma _{\rm g}\mathcal {L}\,, \end{eqnarray}

(3)

\begin{eqnarray} \nabla ^2\Phi = 4\pi G (\Sigma _{\rm g}-\Sigma _0)\delta (z)\,, \end{eqnarray}

(4)

where |$\hat{\boldsymbol {x}}$| points in the radial direction, ρ is the gas mass density, |$\boldsymbol {u} = (u_{\rm x}, u_{\rm y})$| is the gas velocity relative to the background Keplerian flow |$\boldsymbol {u_0} = (0, -q\Omega x)$|⁠, P is the gas pressure, Φ is the self-gravitational potential of a razor-thin disc, q = 3/2 is the Keplerian shear parameter,

\begin{equation} E = {U} + {K} = \frac{P}{\Gamma -1} + \frac{1}{2}\Sigma _{\rm g}u^2 \end{equation}

(5)

is the sum of the internal energy density U and bulk kinetic energy density K for an ideal gas with 2D specific heat ratio Γ = 2, and

\begin{equation} \boldsymbol {T_{\rm g}} = \frac{1}{4\pi G}\left[\nabla \Phi \nabla \Phi - \frac{1}{2}\left(\nabla \Phi \right)\cdot \left(\nabla \Phi \right)\boldsymbol {I}\right] \end{equation}

(6)

is the gravitational stress tensor with identity tensor |$\boldsymbol {I}$|⁠.

We choose a very simple cooling function,

\begin{equation} \Sigma _{\rm g}\mathcal {L}= U/t_{\rm {cool}}= \Omega U / \beta \,\,\,\, {\rm (constant \,\, cooling \,\, time)} \end{equation}

(7)

with β ≡ Ωt_cool constant everywhere. The assumption of constant cooling time t_cool is adopted by many 2D (e.g. Gammie 2001; Johnson & Gammie 2003; Paardekooper 2012) and three-dimensional (3D; e.g. Rice et al. 2003; Lodato & Rice 2004, 2005; Mejía et al. 2005; Cossins, Lodato & Clarke 2009; Meru & Bate 2011) simulations of self-gravitating discs. This prescription enables direct experimental control over the rate of energy loss.

For the solids, we assume that the dust particles only passively respond to the GT turbulence via the aerodynamical drag and also the gravitational acceleration from the gas. No particle feedback is included in this study. In the 2D local approximation, the equation of particle motion reads

\begin{equation} \frac{\mathrm{d}\boldsymbol {v_i}}{\mathrm{d}t} = \frac{\boldsymbol {u}-\boldsymbol {v_i}}{t_{\rm {stop}}{_{,i}}} - \nabla \Phi -2\Omega \hat{\boldsymbol {z}}\times \boldsymbol {v_i} + 2q\Omega ^2 x\hat{\boldsymbol {x}} \,, \end{equation}

(8)

in which the first term is the drag force per unit mass, the second term −∇Φ is the gravitational pull from the self-gravitating gas. The particle velocity |$\boldsymbol {v_i} = (v_{\rm x},v_{\rm y})$| represents the ith particle specie relative to the background shear.

Our simulations are run with Athena (Stone et al. 2008) with the built-in particle module (Bai & Stone 2010). We adopt the van Leer integrator (van Leer 2006; Stone & Gardiner 2009), a piecewise linear spatial reconstruction in the primitive variables, and the HLLC (Harten–Lax–van Leer-Contact) Riemann solver. We solve Poisson's equation of a razor-thin disc using fast Fourier transforms (Gammie 2001; Paardekooper 2012)² Boundary conditions for our physical variables (ρ, |$\boldsymbol {v}$|⁠, U, and Φ) are shearing-periodic in radius (x) and periodic in azimuth (y). We also use orbital advection algorithms to shorten the timestep and improve conservation (Masset 2000; Johnson, Guan & Gammie 2008; Stone & Gardiner 2010).

2.2 Initial conditions and run setup

We start with pure gas simulations. At t = 0, we initialize a uniformly distributed gas disc and set Σ₀ = Ω = G = 1. The thermal energy is such thatQ = c_sΩ/(πGΣ₀) = 1.1 close to the critical Toomre Q-parameter for a razor-thin disc (≃1; Toomre 1964; Goldreich & Lynden-Bell 1965), where |$c_{\rm {s}}= \sqrt{\Gamma (\Gamma -1)U/\Sigma }$| is the sound speed. The velocity is |$\boldsymbol {u_0} = (\delta _x,-q\Omega x + \delta _y)$|⁠, where δ_x and δ_y are randomized perturbations at ∼1 per cent of the initial sound speed c_s₀ = 1.1π. Our simulation domain is a box which covers [−80, 80] GΣ₀/Ω² in both radial (x) and azimuthal (y) direction with 512² grid points. This amounts to a spatial span of 160H/(πQ) ≃ 51Q^{− 1}H and a resolution of ≃10Q/H, where H ≡ c_s/Ω = πQ(GΣ₀/Ω²) is the disc scale height.

We allow the disc to cool off immediately with a fixed cooling time t_coolΩ ∈ {5, 10, 20, 40}. After a short transient phase which normally takes about twice the cooling time, the disc settles to a quasi-steady gravito-turbulent state in which the heating from compression and shocks is balanced by the imposed cell-by-cell cooling. For example, we show the time evolution of Reynolds and gravitational stresses, surface density and velocity dispersions and Toomre's Q parameter from the t_coolΩ = 10 case in Fig. 1.

Figure 1.

Time history of the gravito-turbulent disc of pure gas case with t_coolΩ = 10. Top left: the gravitational (green) and Reynolds (red) stresses normalized with averaged pressure as a function of time. Top right: the density dispersion versus time. Bottom left: the velocity dispersion with (red) and without (black) density weight. Bottom right: the averaged Toomre Q-parameter using sound speed with (red dashed) and without (black solid)density weight. In our dust+gas hybrid simulations, we choose t = 50Ω⁻¹ (indicated by vertical dotted lines) as our initial state and distribute dust particles randomly in space.

Open in new tab Download slide

All quantities saturate after ≥20Ω⁻¹, and well-established turbulence sustains to the end of the simulations (hundreds of dynamical times). The time-averaged (t > 50Ω⁻¹) nominal α, i.e. stresses normalized with pressure, is α ≃ 0.020 for the sum of the Reynolds and gravitational stress. Toomre's Q is hovering ∼3 (Gammie 2001) which also sets the spatial averaged sound speed 〈c_s〉 = πQ(GΣ₀/Ω). The Q-parameter using density weighted sound speed |$\langle c_{\rm {s}}\rangle _{\scriptscriptstyle \Sigma }$| (red dashed curve in bottom-right panel) shows less fluctuation and slightly diminished (∼10 per cent) mean. We choose the density weighted measure hereafter and simply use 〈c_s〉 to represent the density weighted sound speed. But we do note that the spatial and temporal averaged velocity dispersion 〈u_x〉 (black) is about twice the density weighted value |$\langle u_{\rm {x}}\rangle _{\scriptscriptstyle \Sigma }$| (red) as shown in the bottom-left panel of Fig. 1. The cause and effects will be discussed in Section 3.2. We also emphasize that the density dispersion δΣ_g/Σ₀ ≃ 0.6 for α ∼ 0.02, much stronger than found in the MRI-driven turbulence case where |$\delta \rho /\rho \sim \sqrt{0.5\alpha }\sim 0.1$| for similar α values with or without net weak magnetic field (Nelson & Gressel 2010; Okuzumi & Hirose 2011; Shi, Stone & Huang 2016). We will discuss it further in Section 3.7.

We then randomly distribute dust particles in space at t = 50Ω⁻¹ (for t_coolΩ ≤ 20) or 100Ω⁻¹ (for t_coolΩ = 40), and evolve the particle+fluid system for another 200Ω⁻¹. We implemented seven types of particles with constant stopping time such that the Stokes number τ_s ≡ Ωt_stop ∈ [10⁻³, 10⁻², 0.1, 1, 10, 10², 10³], evenly spaced in logarithmic scale. For each type, we use 2¹⁹ (∼5 × 10⁵) particles, or ∼2 particles per grid cell on average. Their velocities follow the background shear initially. Since gas densities vary in time and space, it would be more physical to fix the size of each particle rather than its stopping time; nevertheless, our default simulations fix stopping times to more easily compare with previous simulations that do likewise. We also verify explicitly that our simulations with fixed stopping time agree well with a simulation that uses fixed particle sizes (see Section 3.6). For this latter run, we employ a stopping time based on the Epstein drag law:

\begin{equation} t_{\rm {stop}}{_{,i}} = \frac{f a^{\ast }_i}{\Sigma c_{\rm {s}}} \,, \end{equation}

(9)

where |$a^{\ast }_i = 10^{i-4}$| for i = 1–7 is the dimensionless size of the ith particle species. The converting factor |$f \simeq 3\pi (G\Sigma _0^2/\Omega ^2)$| is chosen such that |$f a^{\ast }_i/\langle \!\langle \Sigma c_{\rm {s}}\rangle \!\rangle _{\rm t}\sim a^{\ast }_i$|⁠, matching τ_s used in the fixed stopping time runs.

3 RESULTS

3.1 Standard run

We first present our standard run tc = 10, i.e. t_coolΩ = 10 run (see Table 1 for the gas properties). After distributing the particles randomly on the grid, the particles quickly adjust in response to the dynamical gas flows. After ∼20Ω⁻¹, the distribution of particles reaches a steady state. As an illustration, we have shown the gas and dust density in Fig. 2 at t = 180Ω⁻¹. Clearly, the small particles (τ_s ≤ 10⁻²) are nearly perfectly coupled to the gas and therefore share the same density structures of the gas. Particles with intermediate stopping time, for both τ_s = 0.1 and 1, appear to concentrate along the dense gas filaments and cause dust density enhancement of two orders of magnitude relative to the mean (note the colour bars for τ_s = 0.1 and 1 now extend to higher values). Transient vortices are also observable in the snapshots for gas and τ_s < 1 particles, but are probably underresolved, see Gibbons et al. (2015) for the effects of vortices on particle concentration. For large particles (τ_s ≥ 10), they are strongly disturbed by the gravitational stirring from the fluctuating gas (see further discussion of the effects of gravity in Section 3.5). After ∼20Ω⁻¹, they are completely redistributed and their end status recovers a random distribution similar to the initial.

Snapshots of gas and dust surface density distributions at the end of simulation. The domain size is Lx = Ly = 160GΣ0/Ω2 ≃ 17H, and cooling time is 10Ω−1. Values are normalized against initial values and colour coded in log scale. We see small tstop particles tracing the gas, particles with intermediate τs = 0.1–1 strongly clustering, and larger particles diffused across the domain.

Figure 2.

Snapshots of gas and dust surface density distributions at the end of simulation. The domain size is L_x = L_y = 160GΣ₀/Ω² ≃ 17H, and cooling time is 10Ω⁻¹. Values are normalized against initial values and colour coded in log scale. We see small t_stop particles tracing the gas, particles with intermediate τ_s = 0.1–1 strongly clustering, and larger particles diffused across the domain.

Open in new tab Download slide

In Fig. 3, we show the time-averaged fraction of cumulative mass for each type of particle, binned in local surface density (number of particles in each cell) normalized to the mean (2/cell). The running profiles of small-sized particles resemble that of the time-averaged gas (solid black curve). Those with large stopping times track the initial random distribution (dotted black). In contrast, the blue (τ_s = 1) and green (τ_s = 0.1) curves show that intermediate-sized particles exhibit relatively higher densities, Σ_p/〈Σ_p〉 ∼ 10–100.

$The time-averaged cumulative fraction of the total particle mass as a function of particle surface density. The surface density is normalized by the averaged value (equal to the initial value) and thus represents the concentration factor of the dust particles. We also show the time-averaged gas distribution (black solid) and the initial dust distribution (black dotted) for comparison.$

Figure 3.

The time-averaged cumulative fraction of the total particle mass as a function of particle surface density. The surface density is normalized by the averaged value (equal to the initial value) and thus represents the concentration factor of the dust particles. We also show the time-averaged gas distribution (black solid) and the initial dust distribution (black dotted) for comparison.

Open in new tab Download slide

After reaching quasi-steady state, we further evolve the dust+gas mixture for a total duration of 200Ω⁻¹. We then measure the radial diffusion coefficients, particle eccentricity and relative velocities averaged overall particles of the same type. The results are presented in the following subsections.

3.2 Radial diffusion

We utilize the following formula to derive the radial diffusion coefficients D_p,x of different types of particles in our simulations (Youdin & Lithwick 2007; Carballido et al. 2011):

\begin{equation} D_{\rm {p,x}}\equiv \frac{1}{2}\frac{\mathrm{d}\langle |x_p(t)-x_p(0)|^2\rangle }{\mathrm{d}t} \,, \end{equation}

(10)

where D_p,x is the diffusion coefficient for given particle stopping time, x_p(t) and x_p(0) are the radial position at time t and initial (taken to be 50Ω⁻¹ after injecting the particles to the gas disc). The radial coordinate is extended beyond the edges of the sheet so that particles moves on radially without periodic boundary conditions. The measurements are made at every time interval δt = 0.5Ω⁻¹, and are performed for a duration of 150Ω⁻¹ long. The average 〈〉 here is taken for all particles within the same type.

The squared displacements as a function of time on the right-hand side of equation (10) are shown in the left-hand panel of Fig. 4. Each curve represents the displacement of one distinctive type of particles. For particles of τ_s ≤ 0.1, the curves overlap. As they are well coupled with the turbulent gas flows, their displacements reflect the properties of gas diffusion. For larger particles, the gravitational stirring dominates the drag force which introduces some extra effective diffusion. As a result, the displacement curves of those particles have bigger slopes. For the largest two dust species we implemented, the curves show periodic oscillations which are due to the epicyclic motion of individual particles. The large amplitude of the epicyclic motion is a result of the strong gravitational forcing of the background gas.

Figure 4.

Left: the squared radial displacement versus time for all types of particles (τ_s ≡ t_stopΩ) in run tc = 10 (t_coolΩ = 10). The colour symbols (with solid curves) are measurement, the dashed curves are the linear fits to the data. Their slopes are then twice the diffusion coefficients based on equation (10). Right: symbols mark the radial diffusion coefficients for particles with different stopping time. At larger stopping time, the measured coefficients are largely different than models of either Youdin & Lithwick (2007) (dashed grey curve) or Cuzzi et al. (1993) (dotted) due to extra stirring from the gas self-gravity.

Open in new tab Download slide

Now with the help of equation (10) and Fig. 4, we can derive the radial diffusion coefficients based on the slope of the squared displacements using linear fitting. The results are shown in the right-hand panel of Fig. 4 and also recorded in Table 2. The radial diffusion coefficients D_{p, x} are normalized against the gas diffusion coefficient D_{g, x}; the latter is approximated using the D_{p, x} of particles with τ_s = 10⁻³. We find D_{p, x} ≃ D_{g, x} for τ_s ≤ 0.1. However, D_{p, x} exceeds D_{g, x} for particles with longer t_stop, and stay roughly constant, D_{p, x} ∼ 6–8D_{g, x}. For comparison, we also plot the diffusion coefficient predicted by models in Cuzzi et al. (1993, dotted grey curve) and Youdin & Lithwick (2007, dashed grey). Both predict small and decreasing values for larger particles based on homogeneous turbulence without extra forcing like self-gravity, in clear contrast with what we obtain in our gravito-turbulent disc. When gravitational stirring is artificially suppressed (see Section 3.5), we do recover similar relationship as they predicted.

We also check the validity of approximating D_{g, x} with D_{p, x} by measuring the auto-correlation function of the turbulent velocity field directly. In general,

\begin{equation} D_{\rm {g,x}}= \int _0^{\infty }R_{\rm xx}(\tau ) \rm{d}\tau \,, \end{equation}

(11)

where R_xx(τ) = 〈u_x(τ)u_x(0)〉 is the auto-correlation of the gas velocity at time τ. However, in gravito-turbulent disc, the spiral density shock waves cause low density (Σ/Σ₀ ∼ O(10⁻²)) valleys between high-density (Σ/Σ₀ ∼ O(1)) ridges. Most of the matter in the high-density regions has small velocity, but the gas in low-density region has very high velocity (∼c_s). The auto-correlation function defined above would strongly bias towards the low-density instead of the high-density region where most of the small dust particles reside. One way to remove the bias is to calculate the gas or dust³ density weighted auto-correlation function

\begin{equation} R_{\rm xx}(\tau ) = \frac{\int \! \Sigma (\tau )u_{\rm {x}}(\tau )\Sigma (0)u_{\rm {x}}(0) \rm{d}x\text{d}y}{\int \! \Sigma (\tau )\Sigma (0) \text{d}x\text{d}y} = \langle u_{\rm {x}}(\tau )u_{\rm {x}}(0)\rangle _{\scriptscriptstyle \Sigma }\,, \end{equation}

(12)

in which u_x(0) and Σ(0) are velocity and density at a reference time, u_x(τ) and Σ(τ) are measured at time τ from the reference point and are both sheared back to that point in order to calculate the correlation. Shown in Fig. 5 is the velocity auto-correlation calculated with (solid curve) and without (dashed curve) density weight. Integrating the density weighted R_xx, we get D_{g, x} ≃ 0.014〈〈c_s〉〉Ω⁻¹ close to the D_{p, x} ≃ 0.012〈〈c_s〉〉Ω⁻¹ measured with particles of τ_s = 10⁻³. However, using R_xx without density weight would overestimate the diffusion coefficient by a factor of 15.

Figure 5.

The radial velocity auto-correlation functions of the gravito-turbulent gas disc calculated with (solid) and without (dashed) density weight.

Open in new tab Download slide

3.3 Eccentricity growth

When particle eccentricity is small, we have the relation

\begin{equation} v_{\rm {x}}^2 \simeq 2 v_{\rm {y}}^2 \simeq {e}^2\Omega ^2 r^2 \,. \end{equation}

(13)

We can therefore measure the orbital eccentricity of each individual particle according to this relation, and obtain the evolution of the eccentricity as shown in Fig. 6. Although the initial e = 0, we find the mean eccentricity quickly saturates at e ≃ 0.2–0.3(H/R) for particles of τ_s ≤ 1. It then saturates slower and levels off at greater value for increasing τ_s > 1. The saturated values are e ≃ 0.7(H/R) and 1.3(H/R) for τ_s = 10 and 10² particles (see Fig. 7). The eccentricity of τ_s = 10³ particle keeps rising gradually through the end of the simulation, suggesting that a saturation level of ≳2(H/R) might be achieved for a longer simulation which is consistent with previous simulations of planetesimals in gravito-turbulent discs (Britsch et al. 2008; Walmswell et al. 2013). We also find that the particle eccentricity obtained in gravito-turbulent disc is in general much greater than that could be excited in the MRI-driven turbulent disc with similar turbulent stress-to-pressure ratio α. The latter usually gives e ∼ 10⁻²–10⁻¹(H/R)Q⁻¹ (Yang et al. 2009, 2012; Nelson & Gressel 2010), orders of magnitude smaller than what we have observed here in our simulations.

Figure 6.

Top: the time evolution of orbital eccentricities averaged overall particles of the same type. The dashed line is 0.15(Ωt)^1/2, which characterize the early growth of the eccentricity for large particles. Bottom: the relative number distribution of particle eccentricity at late time of simulation. We choose 100 bins evenly divided between log₁₀(eR/H) = −3 and 1 for the distribution plot. The dashed line shows a typical Rayleigh distribution for comparison.

Open in new tab Download slide

Figure 7.

The time-averaged (last 15 orbits) particle eccentricity (or equivalently the velocity dispersion δv/c_s) varies with stopping time. The standard run using t_cool = 10Ω⁻¹ are shown in asterisk. The label ‘high resolution’, ‘large domain’, and ‘no self-gravity’ represent results from tc = 10.hires using doubled resolution, tc = 10.dble with doubled sheet size, and tc = 10.wosg without gravitational forcing from the gas respectively, and are discussed in Sections 3.5 and 3.8.

Open in new tab Download slide

Since the particle eccentricity is excited by nearly random gravity field, the evolution will follow the general |$\propto\! \sqrt{t}$| law (Ogihara et al. 2007),

\begin{equation} e = C \left(\frac{H}{R}\right) (\Omega t)^{1/2} \,, \end{equation}

(14)

where the dimensionless coefficient C determines the growth rate and can be measured with our simulation. We fit the early growing phase of the both τ_s = 10² (cyan square) and 10³ (magenta cross) with the above relation and obtain C ≃ 0.15 as the best-fitting coefficient (see the black dashed curve in Fig. 6) The excitation time-scale of eccentricity could be estimated as t_exc ∼ e/(de/dt) ∼ 2e²(R/H)²C⁻²Ω⁻¹. We can therefore predict the saturated eccentricity by equating t_exc with the damping time-scale, in our case, is simply the stopping time t_stop. The eccentricity at saturation is therefore

\begin{equation} e \sim C (H / R) (\tau _s/2)^{1/2} \,. \end{equation}

(15)

For τ_s = 10² particles, this gives e ∼ 1.1(H/R) that matches what we observe in Fig. 6 very well. It also predicts a saturation level of e ∼ 4.7(H/R) if we extend the simulation for another ∼300Ω⁻¹.

Equation (14) also allows us to measure the dimensionless parameter γ which characterizes the amplitude of the fluctuating gravity field as defined in Ogihara et al. (2007). Following Ida et al. (2008) and Okuzumi & Ormel (2013), we write the eccentricity growth as

\begin{equation} e \simeq 1.6\gamma \left(\frac{M_{{\odot }}}{M_{\ast }}\right) \left(\frac{\Sigma }{10\,{\rm g}\,{\rm cm^{-3}}}\right) \left(\frac{R}{100\,{\text{au}}}\right)^2 \left(\Omega t\right)^{1/2} \,. \end{equation}

(16)

After comparing with equation (14), we find the dimensionless turbulent strength γ ≃ 0.01 in our gravito-turbulent disc with Ωt_cool = 10 or α ≃ 0.02, considerably larger than that in MRI discs with similar α (e.g. γ ∼ 10⁻⁴ in Yang et al. 2012).

We also calculate the saturated eccentricity distributions of particles. In Fig. 6, we find that particles of all types obey a Rayleigh-like distribution (Yang et al. 2009, 2012). The probability rises towards greater eccentricity, and then drops at roughly the mean eccentricity measured in the top panel of Fig. 6. Particles having τ_s < 1 share nearly the same properties as gas. As τ_s increases above unity, increasingly many particles obtain higher eccentricities owing to stochastic gravitational forcing by gas. This behaviour contrasts with that shown in fig. 6 of Gibbons et al. (2012), in which the particle velocity distribution narrows as τ_s exceeds unity; their simulations omit gravitational stirring.

3.4 Relative velocity

The high eccentricities obtained in these particles indicate large velocity dispersion, which leads to the question of what is the relative velocity between pair collision. We thus measure the relative velocity of the same type particle or the monodisperse case v_rel, and also the relative velocity with respect to the smallest grains (τ_s = 10⁻³) or a bidisperse velocity |$v_{\rm rel}^{\prime }$|⁠, in each grid cell assuming there is no sub-structure below this scale. The results are shown in Fig. 8.

$Average relative velocities versus particle stopping time (top row), and underlying probability density distributions of relative velocity (bottom row). Averages are taken over the final orbit of the simulations. The relative velocity vrel is evaluated between particles of the same type (stopping time), while $v^{\prime }_{\rm rel}$ is measured with respect to τs = 10−3 particles. Error bars show the standard deviation for run tc = 10 alone (asterisks). The probability density curve for vrel between τs = 0.1 particles (green squares in bottom-left panel) has a diminished amplitude because there is significant probability at vrel/〈cs〉 < 10−4, which is off the scale of the plot. For comparison, a Maxwellian distribution is shown in the bottom-right panel as a dashed curve.$

Figure 8.

Average relative velocities versus particle stopping time (top row), and underlying probability density distributions of relative velocity (bottom row). Averages are taken over the final orbit of the simulations. The relative velocity v_rel is evaluated between particles of the same type (stopping time), while |$v^{\prime }_{\rm rel}$| is measured with respect to τ_s = 10⁻³ particles. Error bars show the standard deviation for run tc = 10 alone (asterisks). The probability density curve for v_rel between τ_s = 0.1 particles (green squares in bottom-left panel) has a diminished amplitude because there is significant probability at v_rel/〈c_s〉 < 10⁻⁴, which is off the scale of the plot. For comparison, a Maxwellian distribution is shown in the bottom-right panel as a dashed curve.

Open in new tab Download slide

In general, the relative velocity (see the asterisk symbols in top two panels) v_rel and |$v_{\rm rel}^{\prime }$| increase from ∼10⁻²c_s to ≳c_s as the τ_s increases from 10⁻³ to 10³ which is indicated by the increasing velocity dispersion or eccentricity as discussed in previous section. Exceptions occur at τ_s = 0.1-1, where particles show the strongest clustering effects (coherent and therefore less relative motions in the filaments), the relative velocity of the same type particles is largely reduced. It even drops below 10⁻³c_s for τ_s = 0.1 case. The general approach of approximating particle relative velocity via the measurement of velocity dispersion fails here. Our v_rel and |$v_{\rm rel}^{\prime }$| for smaller particles (τ_s ≤ 0.1) seem similar to what are found in MRI-driven turbulent discs (Carballido et al. 2010). However, our v_rel and |$v^{\prime }_{\rm rel}$| increases with t_stop for larger particles that have τ_s > 0.1, in contrast to MRI-driven turbulence results (cf. figs 2 and 3 in Carballido et al. 2010) where turbulent torquing due to the fluctuating gravity field of the gas is not considered. The latter show v_rel falls continuously, and |$v^{\prime }_{\rm rel}$| stays roughly constant with t_stop, in agreement with the theoretical prediction in isotropic turbulence (Ormel & Cuzzi 2007). We will discuss the effects of gravitational stirring in Section 3.5.

We also show the probability density functions (PDFs) for v_rel and |$v_{\rm rel}^{\prime }$| in Fig. 8. The PDFs of larger particles (τ_s > 1) are similar to Maxwellian distribution (see the dashed curve in the bottom-right panel as an example of Maxwellian distribution; Windmark et al. 2012; Garaud et al. 2013). Particles smaller than τ_s ∼ 1, however, show broader distributions and smaller mean values than the bigger particles. They deviate significantly from a Maxwellian and resemble more of a log-normal distribution (Mitra, Wettlaufer & Brandenburg 2013; Pan, Padoan & Scalo 2014). Relative velocities with respect to the smallest particles, i.e. |$v^{\prime }_{\rm rel}$|⁠, have distributions that shift gradually to higher values as the size of the larger particle increases. Between τ_s = 0.01 and τ_s = 0.001 particles, the largest relative velocities |$v^{\prime }_{\rm rel} > 0.1 c_{\rm {s}}$|⁠, which might lead to collisional destruction, are mostly due to their different deceleration rates in post-shock regions (Nakamoto & Miura 2004; Jacquet & Thompson 2014). As plotted in Fig. 9, the locations where |$v^{\prime }_{\rm rel} >0.1 c_{\rm {s}}$| (white symbols) are found just behind shock fronts, and where gas radial velocities are discontinuous.

Snapshot of gas radial velocities (normalized by the volume-averaged sound speed) at t = 200Ω−1. Large relative velocities (>0.1cs) between τs = 10−2 and 10−3 particles are overlaid as white symbols. Evidently, large relative velocities between particles characterize regions just behind shock fronts; in post-shock regions, particles decelerate at different rates according to their different stopping times (Nakamoto & Miura 2004; Jacquet & Thompson 2014).

Figure 9.

Snapshot of gas radial velocities (normalized by the volume-averaged sound speed) at t = 200Ω⁻¹. Large relative velocities (>0.1c_s) between τ_s = 10⁻² and 10⁻³ particles are overlaid as white symbols. Evidently, large relative velocities between particles characterize regions just behind shock fronts; in post-shock regions, particles decelerate at different rates according to their different stopping times (Nakamoto & Miura 2004; Jacquet & Thompson 2014).

Open in new tab Download slide

Relative velocities v_rel between particles of the same type can be separated into two groups. One group corresponds to the small-size particles (τ_s ≤ 1), for which relative velocities peak at v_rel ≃ 0.01c_s. The other group corresponds to large-sized particles (τ_s > 1), for which peak velocities are sonic. We note that τ_s = 0.1 (green squares) particles show diminished amplitude and a significant shift towards smaller relative velocity due to the clustering effect: ∼80 per cent of contribution from v_rel/c_s ≪ 10⁻⁴ is not shown in this plot. While for τ_s = 1 particles, there is also a second, near sonic, contribution, which shows the transitional behaviour from small size (friction dominated) to large size (gravity dominated) particles.

3.5 Effects of gravitational stirring

The gravitational stirring effect has been studied in MRI-driven turbulent discs and shown to induce high-velocity dispersion for solid bodies bigger than ∼ 100 m (or τ_s ≳ 10³ at radius R ∼ au) (Nelson & Gressel 2010; Yang et al. 2009, 2012; Okuzumi & Ormel 2013; Ormel & Okuzumi 2013). But we will show the dominating gravitational forcing kicks in for even smaller particles owning to its much stronger density fluctuation, |$\delta \Sigma /\Sigma \sim 5\sqrt{\alpha }$| in gravito-turbulent discs (see Section 3.7) rather than |${\sim } \sqrt{0.5\alpha }$| in MRI-driven turbulent discs (Nelson & Gressel 2010; Yang et al. 2012).

As the stopping time increases, the first term in equation (8), the drag force, becomes less important compared to the second term, the gravitational acceleration. Assuming the characteristic length scale of the spiral density features is l, then |∇Φ| ∼ Φ/l ∼ GδΣ, where δΣ ∼ Σ in our gravito-turbulent disc is the density fluctuation. The length scale can be approximated with the most unstable wavelength for axisymmetric disturbances, or simply the disc scaleheight, |$l \sim c_{\rm {s}}^2/G\Sigma \sim Q^2 (G\Sigma /\Omega ^2) \sim H$|⁠. This is also supported by our simulations; the density waves have typical radial wavelength ∼H in the top-left panel of Fig. 2. The drag force |u_x − v_x|/t_stop ≲ Ωl/t_stop ∼ c_s/t_stop. The ratio of these two thus gives

\begin{equation} \frac{|\nabla \Phi |}{|u_{\rm {x}}-v_{\rm {x}}|/t_{\rm {stop}}} \sim \frac{1}{Q} \frac{\delta \Sigma }{\Sigma } \tau _s \,, \end{equation}

(17)

the gravitational stirring becomes more important than aerodynamical drag when τ_s > Q ∼ 1 in gravito-turbulent discs. However, for non-self-gravitating turbulent discs such as MRI-driven turbulent discs, δΣ/Σ is one order of magnitude smaller, and Q ≫ 1. As a results, the gravitational stirring only affect particles with τ_s > 10Q ≫ 1. In top panel of Fig. 10, we show the time variation of radial positions of two representative particles in the standard simulation. The particle with larger τ_s = 10³ oscillates with a peak-to-peak amplitude of ∼40GΣ₀/Ω² at a frequency of the orbital frequency Ω after it is released in the turbulent disc. However, the small particle, τ_s = 10⁻³, scatters randomly over time and does not show strong periodic variations.

Examples of radial displacement for two particles with fixed stopping time τs = 10−3 (black) and 103 (magenta) in run with (top panel) and without (bottom panel) gravity pull from the gas in the calculation. The gravitational stirring of the τs = 103 particle causes large amplitude oscillations which are absent in the bottom panel when gravity from the gas is removed.

Figure 10.

Examples of radial displacement for two particles with fixed stopping time τ_s = 10⁻³ (black) and 10³ (magenta) in run with (top panel) and without (bottom panel) gravity pull from the gas in the calculation. The gravitational stirring of the τ_s = 10³ particle causes large amplitude oscillations which are absent in the bottom panel when gravity from the gas is removed.

Open in new tab Download slide

To reveal the self-gravity effect to the dust dynamics, we perform a rerun (tc = 10.wosg) of the standard simulation but suppress the gravitational forcing manually, i.e. removing −∇Φ term in equation (8). The radially projected trajectories of the same two representative particles are plotted in the bottom panel of Fig. 10 for comparison. Clearly, the small particle is not affected by the missing of gravitational acceleration and appears to behave similarly as in the standard run. However, the larger particle in tc = 10.wosg barely moves radially when the only acceleration/deceleration is friction.

Without the gravitational stirring, the diffusion of particle also look significantly different than in Fig. 4. In Fig. 11, we present the diffusion coefficients measured in run tc = 10.wosg using the same method as in Fig. 4. Again the small t_stop particles are unaffected and therefore give nearly the same D_{p, x} as in the standard run. D_{p, x} of particles with τ_s > 1 drops rapidly as D_{p, x} ∝ (τ_s)⁻² which manifests the lack of gravitational forcing. We note that the exact relation between diffusion coefficient and stopping time slightly deviates from models of Cuzzi et al. (1993), Youdin & Lithwick (2007). We speculate that is a result of different power spectra in gravito-turbulent discs than normal homogeneous turbulence model adopted in these models.

$Similar to Figs 3 and 4 , we show radial diffusion coefficients (top panel) and cumulative fraction of dust mass (bottom panel) for run tc = 10.wosg, in which the gravitational forcing from the gas is artificially removed. At larger stopping time, the measured Dp,x without gravitational stirring are now closer to models of Youdin & Lithwick (2007) (dashed grey curve) or Cuzzi et al. (1993, dotted). The clustering effect for particles of τs = 0.1 observed in Fig. 3 is absent when gravitational stirring is removed.$

Figure 11.

Similar to Figs 3 and 4 , we show radial diffusion coefficients (top panel) and cumulative fraction of dust mass (bottom panel) for run tc = 10.wosg, in which the gravitational forcing from the gas is artificially removed. At larger stopping time, the measured D_p,x without gravitational stirring are now closer to models of Youdin & Lithwick (2007) (dashed grey curve) or Cuzzi et al. (1993, dotted). The clustering effect for particles of τ_s = 0.1 observed in Fig. 3 is absent when gravitational stirring is removed.

Open in new tab Download slide

In Fig. 11, we also compute the cumulative mass fraction for run tc = 10.wosg. As expected, the Ωt_stop = 1 particles show the strongest concentration and more contribution towards the higher concentration (≳100) than the standard tc = 10 run. Both Ωt_stop = 0.1 and Ωt_stop = 10.0 show similar secondary clustering but are much weaker than the concentration of Ωt_stop = 0.1 in the case where gravity is present.

Without the stochastic gravitational stirring, the averaged particle eccentricity for all types of particles never exceed e ∼ 0.3(H/R) as shown in Fig. 7 (green squares). The small particles (τ_s ≤ 1) have similar eccentricities as in the standard run; however, those larger particles are dynamically unimportant, and e approaches nearly zero (initial value) as τ_s increases. The gravitational stirring is also responsible for the increasing relative velocities for particles of τ_s ≥ 1. As shown in top row of Fig. 8 in green squares, instead of rising up, v_rel decreases while |$v^{\prime }_{\rm rel}$| stays constant with increasing τ_s when gravitational forcing is not included, which is consistent with the results from previous study (Ormel & Cuzzi 2007; Carballido et al. 2010).

3.6 Fixed particle size

In case that the fluctuations of the density and sound speed are large (≲O(1)), as is the case in gravito-turbulent discs, the constant stopping time assumption for individual particles might not be valid as the particles would have varying t_stop as they travel around different regions of the disc. We therefore, in this section, test if our results still hold statistically when we fix particle size instead of the stopping time.

We set up seven types of particles with fixed size, dimensionless size |$a_i^{\ast } = 10^{i-4}$| for i = 1–7, as described in equation (9) and Section 2.2. Everything else is kept the same as the standard run tc = 10. For each type of particles, we define an effective stopping time |$\langle \langle t_{\rm {stop}}{_{,i}}\rangle _{\scriptscriptstyle \Sigma }\rangle _t$| by first measuring the stopping time of individual particle according to equation (9), and then averaging over all particles of the same type and time. We find the averaged stopping time 〈〈τ_s_{, i}〉〉_tΩ ≃ [1.04 × 10⁻³, 1.01 × 10⁻², 0.086, 0.827, 13.5, 2.10 × 10², 2.31 × 10³] for |$a_i^{\ast } = 10^{-3}$|–10³. For small particles (⁠|$a_i^{\ast } < 0.1$|⁠), |$\langle \langle \tau _s{_{,i}}\rangle \rangle _t \simeq a_i^{\ast } \Omega ^{-1}$| as designed because we choose the converting factorf in equation (9) to match a* with τ_s of the fixed stopping time case. For intermediate size, the averaged stopping time is only ≲20 per cent smaller than the designed values. For larger particle with |$a_i^{\ast } > 1$|⁠, the effective stopping time is roughly |${\sim } 2 a_i^{\ast } \Omega ^{-1}$|⁠.

After converting a* to the effective stopping time, we find that the results stay the same as the simulations using particles of constant stopping times. In Fig. 12, we show the diffusion coefficient, particle eccentricity, and relative velocities in large coloured symbols. They almost follow the results using fixed t_stop (grey squares). Therefore, statistically, these results validate the usage of fixed t_stop particles even in the highly fluctuated density field of a gravito-turbulent disc. In general, the results of fixed t_stop can approximate the results of fixed size particles. But this approach (constant t_stop) might underestimate the clustering effect for the intermediate size particle τ_s = 1 as shown in Fig. 13. This also leads to an overestimate of the relative speed v_rel for intermediate size particles of τ_s ∼ 1. The actual relative speed for intermediate size particles would become even smaller than what we have obtained.

$Diffusion coefficients (top left), particle eccentricity (top right), relative velocities vrel (bottom left), and $v^{\prime }_{\rm rel}$ (bottom right) versus the effective stopping time for fixed size simulation tc = 10.size. The effective stopping times are calculated by time and particle averaging of tstop of individual particle. Results of tc = 10 from Fig. 4, which are shown with smaller grey squares, almost match fixed size results and therefore validate the usage of fixed stopping time particles in this work. We note that the constant tstop approach could overestimate the relative speed for intermediate size particles with a* = 1 or effective τs ∼ 1.$

Figure 12.

Diffusion coefficients (top left), particle eccentricity (top right), relative velocities v_rel (bottom left), and |$v^{\prime }_{\rm rel}$| (bottom right) versus the effective stopping time for fixed size simulation tc = 10.size. The effective stopping times are calculated by time and particle averaging of t_stop of individual particle. Results of tc = 10 from Fig. 4, which are shown with smaller grey squares, almost match fixed size results and therefore validate the usage of fixed stopping time particles in this work. We note that the constant t_stop approach could overestimate the relative speed for intermediate size particles with a* = 1 or effective τ_s ∼ 1.

Open in new tab Download slide

Figure 13.

How particles concentrate in a run with fixed stopping time τ_s = 1 (from run tc = 10) versus a run with fixed particle size a* = 1 (from run tc = 10.size). The latter case, having an effective τ_s ≃ 0.83, shows stronger concentration.

Open in new tab Download slide

3.7 Effects of cooling time

The dust dynamics would also be affected by the strength of the turbulence. In a gravito-turblent disc, the strength of the turbulence, e.g. in term of α, is inversely proportional to the cooling time (Gammie 2001; Shi & Chiang 2014). We can therefore study how the dust dynamics changes with the turbulent strength by varying t_cool.

We perform simulations with t_coolΩ = 5 (run tc=5), 20 (run tc=20), and 40 (run tc=40) using the same particle distribution as in the t_coolΩ = 10 standard run. We then measure the diffusion coefficient in each t_cool case and present the results in Fig. 14. We find D_{p, x} of various cooling time share rather similar profile. It is flat at τ_s < 0.1 end when aerodynamical coupling is strong. D_{p, x} traces D_{g, x} in this regime, and the Schmidt number, which measures the ratio of angular momentum transport and mass diffusion, |${\rm Sc} \equiv \alpha c_{\rm {s}}^2\Omega ^{-1}/D_{\rm {g,x}}\sim 1.4$|–2.4, stays roughly constant. The diffusion coefficient rises up for 0.1 ≤ τ_s ≤ 1, indicating the increasing effect of the gravitational stirring. In the very long stopping time regime, τ_s ≥ 10, the profile levels off again as the dominant gravitational acceleration does not depend on t_stop. The magnitude of D_{p, x} is usually 5–9 times of the diffusion of small particle or gas. For given particle size, we find D_{p, x} scales roughly as ∝ 1/t_cool which indicates stronger diffusion in stronger turbulence.

$The dependence on cooling time (different symbols) of the radial diffusion coefficient (top left), saturated particle eccentricity (top right), relative velocity between particles of the same type (bottom left), and relative velocity with respect to τs = 10−3 particles (bottom right). The measured Dp,x is nearly inversely proportional to tcool, while particle eccentricity and relative velocities roughly follow $t_{\rm {cool}}^{-1/2}$.$

Figure 14.

The dependence on cooling time (different symbols) of the radial diffusion coefficient (top left), saturated particle eccentricity (top right), relative velocity between particles of the same type (bottom left), and relative velocity with respect to τ_s = 10⁻³ particles (bottom right). The measured D_p,x is nearly inversely proportional to t_cool, while particle eccentricity and relative velocities roughly follow |$t_{\rm {cool}}^{-1/2}$|⁠.

Open in new tab Download slide

Varying the cooling time also affects the particle concentration. In Fig. 15, we show the density snapshots of simulations using different cooling times at the top row. As t_cool gets longer, we find the turbulent amplitude diminishes, so does the tilt angle of the shearing wave features. This leads to less frequent interactions between separated shearing waves and longer lifetime for existing density features. In the bottom row of Fig. 15, we find strong interactions of shearing waves in Ωt_cool = 5 run leads to plume-like structure of τ_s = 1 particle with typical concentration of ∼10–100. The concentration is greater, ≳100 − 10³ for Ωt_cool = 40 case, and particles are more spatially confined in very narrow streams. In Fig. 16, we quantitatively confirm this result by showing the cumulative mass fraction of the τ_s = 1 particles for various cooling times. As t_cool rises, more and more mass actually concentrate towards larger dust density, although the density dispersion of gas decreases.

Figure 15.

The logarithmic values of the gas (top) and Ωt_stop = 1 particle density (bottom) at the end of simulations with various cooling time as shown in the top panels of each column. Colour bars are on the left of each row. As t_cool increases, the fluctuation of gas density becomes weaker, the tilt angle of shearing waves also gets smaller, both indicates relatively weaker angular momentum transport as the cooling time-scale gets longer. However, stronger clustering effects are found for longer cooling cases as the velocity dispersion from the gas gets weaker.

Open in new tab Download slide

$Similar to Fig. 3, the cumulative mass fraction of the particles with Ωtstop = 1 for various cooling times. Longer tcool results in stronger concentration and clustering.$

Figure 16.

Similar to Fig. 3, the cumulative mass fraction of the particles with Ωt_stop = 1 for various cooling times. Longer t_cool results in stronger concentration and clustering.

Open in new tab Download slide

The particle relative velocity and eccentricity drops as the cooling time gets longer (see Fig. 14). In general, from Ωt_cool = 10 to 40, v_rel, |$v^{\prime }_{\rm rel}$|⁠, and e diminish by less than a factor of ∼2–3, roughly scales as ∝ α^1/2, similar to the way velocity/density dispersion scales. Exceptions occur at the intermediate size particles (as found in Section 3.4), τ_s = 0.1 and 1, where we find the relative velocity of the same kind v_rel drops rapidly with increasing cooling time. It diminishes ∼30 times from Ωt_cool = 5 to 40 for τ_s = 0.1 particles and ≲10³ times for particles with τ_s = 1. As the intermediate-sized particles show the strongest clustering, the rather slow relative velocity would favour gravitational collapsing.

Based on our results in Table 1, the density dispersion usually follows δΣ/Σ ≃ 2(Ωt_cool)^−1/2 ≃ 5α^1/2 in gravito-turbulent discs. As already mentioned in Section 3.5, this is much stronger fluctuation than found in MRI discs (∼≲α^1/2). Similar as in Section 3.3, we also measure the dimensionless turbulent strength γ based on the eccentricity growth for different t_cool runs and list them in Table 1. The strength parameter γ ≃ 0.07α^1/2 seems about one order of magnitude greater than those found in MRI-driven turbulent discs (Nelson & Gressel 2010; Yang et al. 2012) which also explains why we observe much larger eccentricity and relative speed for large particles (τ_s > 1) in our gravito-turbulent discs.

Table 1.

Open in new tab

Gas properties in gravito-turbulent discs.

Name	Ωt_cool	Q	δΣ_g/Σ_g^(a)	γ^(b)	δu_x/c_s^(c)	α_tot^(d)	\|$D_{\text{g,x}}\,(c_{\rm {s}}^2\Omega ^{-1}){}^{(e)}$\|	Sc^(f)
tc = 5	5	3.35	0.82	0.013	0.38 (0.91)	0.035	0.022	1.6 (1.7)
tc = 10	10	3.19	0.63	9.4 × 10⁻³	0.31 (0.61)	0.020	0.014	1.4 (1.7)
tc = 10.dble^(g)	10	3.15	0.61	9.0 × 10⁻³	0.34 (0.64)	0.022	0.016	1.4 (1.5)
tc = 10.hires^(h)	10	3.09	0.59	9.5 × 10⁻³	0.33 (0.59)	0.024	0.015	1.6 (1.7)
tc = 20	20	3.16	0.52	8.3 × 10⁻³	0.26 (0.38)	0.011	0.0061	1.8 (2.0)
tc = 40	40	2.99	0.32	5.2 × 10⁻³	0.20 (0.23)	0.006	0.0025	2.4 (2.9)

Name	Ωt_cool	Q	δΣ_g/Σ_g^(a)	γ^(b)	δu_x/c_s^(c)	α_tot^(d)	\|$D_{\text{g,x}}\,(c_{\rm {s}}^2\Omega ^{-1}){}^{(e)}$\|	Sc^(f)
tc = 5	5	3.35	0.82	0.013	0.38 (0.91)	0.035	0.022	1.6 (1.7)
tc = 10	10	3.19	0.63	9.4 × 10⁻³	0.31 (0.61)	0.020	0.014	1.4 (1.7)
tc = 10.dble^(g)	10	3.15	0.61	9.0 × 10⁻³	0.34 (0.64)	0.022	0.016	1.4 (1.5)
tc = 10.hires^(h)	10	3.09	0.59	9.5 × 10⁻³	0.33 (0.59)	0.024	0.015	1.6 (1.7)
tc = 20	20	3.16	0.52	8.3 × 10⁻³	0.26 (0.38)	0.011	0.0061	1.8 (2.0)
tc = 40	40	2.99	0.32	5.2 × 10⁻³	0.20 (0.23)	0.006	0.0025	2.4 (2.9)

Notes.^(a)Time- and spatial averaged density dispersion, i.e. |$\langle \!\langle \delta \Sigma _g^2\rangle \!\rangle _{t}^{1/2}/\Sigma _g$|⁠.

^(b)Dimensionless parameter used in Ida et al. (2008) which characterize the amplitude of the random gravity field.

^(c)Velocity dispersion with density weight, i.e. |$\langle \!\langle u_{\rm {x}}^2\rangle _{\scriptscriptstyle \Sigma }\rangle _{t}^{1/2}/c_{\rm {s}}$|⁠, where c_s is the density weighted sound speed with Γ = 2. The dispersions calculated without density weight are in the parentheses.

^(d)α_tot = α_R + α_g, the total internal stress normalized with averaged pressure.

^(e)Gas diffusion coefficient calculated using the auto-correlation function of equation (11).

^(f)Schmidt number |${\rm Sc}=\alpha c_{\rm {s}}^2\Omega ^{-1}/D_{\rm {g,x}}$|⁠. The values in parentheses are calculated assuming D_g,x ≃ D_p,x for τ_s = 10⁻³ particles (see Table 2).

^(g)Similar to run tc = 10 but using doubled domain sized and fixed grid resolution and particle number.

^(h)Similar to run tc = 10 but using doubled grid resolution and particle number.

Table 1.

Open in new tab

Gas properties in gravito-turbulent discs.

Name	Ωt_cool	Q	δΣ_g/Σ_g^(a)	γ^(b)	δu_x/c_s^(c)	α_tot^(d)	\|$D_{\text{g,x}}\,(c_{\rm {s}}^2\Omega ^{-1}){}^{(e)}$\|	Sc^(f)
tc = 5	5	3.35	0.82	0.013	0.38 (0.91)	0.035	0.022	1.6 (1.7)
tc = 10	10	3.19	0.63	9.4 × 10⁻³	0.31 (0.61)	0.020	0.014	1.4 (1.7)
tc = 10.dble^(g)	10	3.15	0.61	9.0 × 10⁻³	0.34 (0.64)	0.022	0.016	1.4 (1.5)
tc = 10.hires^(h)	10	3.09	0.59	9.5 × 10⁻³	0.33 (0.59)	0.024	0.015	1.6 (1.7)
tc = 20	20	3.16	0.52	8.3 × 10⁻³	0.26 (0.38)	0.011	0.0061	1.8 (2.0)
tc = 40	40	2.99	0.32	5.2 × 10⁻³	0.20 (0.23)	0.006	0.0025	2.4 (2.9)

Name	Ωt_cool	Q	δΣ_g/Σ_g^(a)	γ^(b)	δu_x/c_s^(c)	α_tot^(d)	\|$D_{\text{g,x}}\,(c_{\rm {s}}^2\Omega ^{-1}){}^{(e)}$\|	Sc^(f)
tc = 5	5	3.35	0.82	0.013	0.38 (0.91)	0.035	0.022	1.6 (1.7)
tc = 10	10	3.19	0.63	9.4 × 10⁻³	0.31 (0.61)	0.020	0.014	1.4 (1.7)
tc = 10.dble^(g)	10	3.15	0.61	9.0 × 10⁻³	0.34 (0.64)	0.022	0.016	1.4 (1.5)
tc = 10.hires^(h)	10	3.09	0.59	9.5 × 10⁻³	0.33 (0.59)	0.024	0.015	1.6 (1.7)
tc = 20	20	3.16	0.52	8.3 × 10⁻³	0.26 (0.38)	0.011	0.0061	1.8 (2.0)
tc = 40	40	2.99	0.32	5.2 × 10⁻³	0.20 (0.23)	0.006	0.0025	2.4 (2.9)