The role of the drag force in the gravitational stability of dusty planet-forming disc – II. Numerical simulations Free

ABSTRACT

Young protostellar discs are likely to be both self-gravitating, and to support grain growth to sizes where the particles decoupled from the gas. This combination could lead to short-wavelength fragmentation of the solid component in otherwise non-fragmenting gas discs, forming Earth-mass solid cores during the Class 0/I stages of young stellar object evolution. We use three-dimensional smoothed particle hydrodynamics simulations of two-fluid discs, in the regime where the Stokes number of the particles St > 1, to study how the formation of solid clumps depends on the disc-to-star mass ratio, the strength of gravitational instability, and the Stokes number. Gravitational instability of the simulated discs is sustained by local cooling. We find that the ability of the spiral structures to concentrate solids increases with the cooling time and decreases with the Stokes number, while the relative dynamical temperature between gas and dust of the particles decreases with the cooling time and the disc-to-star mass ratio and increases with the Stokes number. Dust collapse occurs in a subset of high disc mass simulations, yielding clumps whose mass is close to linear theory estimates, namely 1–10 M_⊕. Our results suggest that if planet formation occurs via this mechanism, the best conditions correspond to near the end of the self-gravitating phase, when the cooling time is long and the Stokes number close to unity.

1 INTRODUCTION

Direct imaging of protostellar discs, in scattered light and especially at sub-mm wavelengths, has opened debate as to how the time-scales for planet formation align with the established evolutionary sequence for young stellar objects (YSO; Adams, Lada & Shu 1987). Many discs show substructure (Andrews 2020) that is consistent with – and often interpreted as – the theoretically expected signature of planet–disc interaction (Dong, Zhu & Whitney 2015; Zhang et al. 2018; Lodato et al. 2019). A subset of these discs, including HL Tau (ALMA Partnership 2015), IRS 63 (Segura-Cox et al. 2020), and GY91 (Sheehan & Eisner 2018) are young, with ages estimated to be less than 1  Myr. Assuming a planet origin for the substructure, a substantial part of the planet formation process must overlap with the time when protostellar discs at 10–100 au scales are likely self-gravitating. How this occurs is unclear. Classical planetesimal-only variants of core accretion model (Safronov 1969; Pollack et al. 1996) are not viable; meeting the time-scale constraint at 10 au is already non-trivial, and at substantially larger radii sufficiently massive planet formation is not possible at all. Models with a dominant pebble component (e.g. Rosenthal & Murray-Clay 2018; Guilera et al. 2020; Morbidelli 2020; Chambers 2021; Cummins, Owen & Booth 2022; Jiang & Ormel 2023), or those with large-scale core migration (Levison, Thommes & Duncan 2010), are more promising, but the ability of these processes to form the inferred large-orbital-radius population of planets has not been fully established.

Alternatively, planets may form via gravitational instability (GI; Boss 1997) in a self-gravitating protostellar disc. Self-gravity is expected to be an important process during the Class 0/I phase of protostellar evolution, and is most vigorous at large orbital radii where the cooling time is short (Clarke 2009; Rafikov 2009) and infall continues. GI in a gaseous disc, however, is a poor candidate for forming the Atacama Large Millimetre Array (ALMA)-inferred planet population, because it substantially overshoots the inferred masses. The mass of a gaseous fragment created by GI is of the order of the Jeans mass of the spiral perturbation, and for typical protoplanetary discs its value is, assuming Toomre Q = 1 (Toomre 1964; Kratter & Lodato 2016)

typically between 1 and 10 M_J (where M_Jup is the mass of Jupiter) depending on the disc aspect ratio. However, this is the initial mass of the object: subsequently, the fragment would start accreting material belonging to the accretion disc, increasing its mass and likely becoming a brown dwarf (Kratter et al. 2010; Kratter & Lodato 2016).

A mechanism that may allow lower mass planets to form as a consequence of GI was proposed by Rice, Lodato & Armitage (2005). Concentration of solid particles in spiral arms, together with vertical settling, can lead to gravitational collapse in the solid component even in a non-fragmenting gravitationally unstable gaseous disc (Gibbons, Rice & Mamatsashvili 2012; Gibbons, Mamatsashvili & Rice 2014, 2015). It is known that gas spiral arms act as dust traps, since they are pressure maxima (Rice et al. 2005; Shi et al. 2016). Solid particles collect inside them (Dipierro et al. 2015a), reaching dust-to-gas ratio that can be of the order of unity. Conversely, the interaction between gas spiral arms and dust particles can excite them, imparting random motions that reduce the peak density and potential for collapse (Riols et al. 2020). Walmswell, Clarke & Cossins (2013) found that large dust particles experience gravitational scattering by the spiral arms, while Booth & Clarke (2016) related the level of excitation of solid particles to the aerodynamical coupling and the cooling factor. Marginally coupled solid particles are less excited by spiral arms, while, in rapidly cooled discs, the level of dust excitation is higher. Baehr & Zhu (2021) confirmed this trend through three-dimensional shearing box simulations.

In this work, we perform three-dimensional smoothed particle hydrodynamics (SPH) simulations using the phantom code (Price et al. 2018) of gas and dust in gravitationally unstable protostellar discs. The aim is to investigate the interplay between GI and the dynamics of particles that are aerodynamically coupled to the gas, and to understand under which conditions it is possible to form solid cores inside GI spirals. In Section 2, we recall the basics of the aerodynamical coupling between gas and dust, and the role of dust in GI. In Section 3, we briefly describe the SPH code phantom (Price et al. 2018) and the set of simulations we have performed for this work. In Section 4, we analyse the simulations, show our main results, and compare our findings with previous works. In Section 5, we discuss our findings in an evolutionary perspective. Finally, in Section 6, we present our conclusions.

2 DRAG FORCE AND GRAVITATIONAL INSTABILITY

2.1 Radial drift

To zeroth order, gas dynamics is dictated by the hydrodynamical equations. Neglecting the role of the disc self-gravity, in centrifugal equilibrium, the azimuthal component of the gas velocity is

where |$v_\mathrm{ k}^2 = GM_\star /R$| is the Keplerian speed and the second term is the negative contribution of the pressure gradient. It is possible to explicitly compute the pressure gradient knowing the disc temperature structure: the gas azimuthal speed is slightly sub-Keplerian, and the correction term is of the order of (H/R)² ∼ 0.01, where H is the disc thickness. In addition, the gas has also a negative radial velocity, given by viscosity.

Conversely, solid particle dynamics is dictated by both gravitational and drag force. Since dust is not supported by pressure, the basic velocity of gas and dust is different, but they interact through a drag force, that modifies the relative velocity. Writing Δu = v_d − v_g, the drag force F_d depends upon the relative velocity as

where C_D is a coefficient that depends on the drag regime, s is the dust particle size, and ρ is the total density. The different aerodynamical coupling regimes depend upon the size of dust particles relative to the gas mean free path. It is possible to define a dimensionless parameter, called the Knudsen number (Paardekooper & Mellema 2006), that measures this property

where s is the dust particle size and λ_g is the gas particles’ mean free path, given by

where σ_coll = 2.367 × 10⁻¹⁵cm² is the cross-section of gas particles, μ = 2.1 is the mean molecular weight, ρ_g is the gas density, and m_p is the mass of the proton. Typically, in protoplanetary discs, Kn > 1, meaning that the dust particles’ size is smaller than the gas mean free path: this is called the Epstein regime. However, when the disc is very massive, there is a transition between Epstein and Stokes regime (Kn < 1). The main impact of this transition is that in the inner denser region of the disc, where particles’ size is comparable with the gas mean free path, the Stokes number is higher compared to Epstein regime. The C_D coefficient is given by (Fassio & Probstein 1970)

where c_g is the gas sound speed and the Reynolds number based upon the difference of velocity is given by

and the last equivalence is true for collisional viscosity. We also define the stopping time, i.e. the time needed to modify the relative velocity between gas and dust. The longer it is, the less the particles are coupled. For spherical grains, it is given by

where m_d is the dust grain mass and ρ₀ is the dust grain’s intrinsic density. In order to measure the strength of the aerodynamical coupling, we define a dimensionless Stokes number as the ratio between the stopping time and the dynamical one: in general, its expression is

and, in the Epstein regime, its value is

The smaller the Stokes number is, the more tightly the particles are coupled: so, for St → 0, dust dynamics follows the gas dynamics, while for St → ∞, the two fluids do not influence each other aerodynamically. In a smooth axisymmetric disc, the drag force has dramatic physical consequences. Because of the azimuthal difference of speed between gas and dust, solid particles experience a headwind that slows them down, giving them a negative radial velocity. This effect is called radial drift, and its time-scale for St ∼ 1 is of the order of ∼100 yr, several orders of magnitude shorter than the disc lifetime. This is the ‘metre-sized barrier’, so called because St ∼ 1 corresponds to roughly this physical size in the inner disc (Weidenschilling 1977). Radial drift can be stopped in discs where the presence of a gas pressure maximum creates a ‘dust trap’. Indeed, in a gas pressure maximum, the pressure gradient is zero, and hence the velocity difference between gas and dust vanishes and no radial drift occurs. Several mechanisms have been proposed to form dust traps, such as gaps made by planets (Paardekooper & Mellema 2004, 2006; Ayliffe et al. 2012; Pinilla et al. 2012), zonal flows (Johansen, Youdin & Klahr 2009; Simon & Armitage 2014), or spirals induced by gravitational instabilities (Dipierro et al. 2015a). Concentration of solid material by aerodynamic effects must be considered in assessing models for planet formation in early protoplanetary stages.

2.2 Gravitational instability

In an axisymmetric thin disc composed of a single fluid, the dispersion relation for linear, tightly wound density waves is

where ω is the perturbation frequency, Ω(R) is the angular velocity, m/R is the azimuthal wavenumber of the perturbation, c(R) is the sound speed of the fluid, k is the perturbation radial wavenumber, κ(R) is the epicyclic frequency, and Σ(R) is the disc-surface density. The well-known instability criterion for axisymmetric disturbances is (Toomre 1964)

which identifies the parameter space for which ω² < 0. It is possible to generalize the instability taking into account a second fluid component, i.e. particles in the protostellar case. In the context of galactic dynamics, the gravitational role of a second component has been investigated by Jog & Solomon (1984) and Bertin & Romeo (1988). Longarini et al. (2023) applied this method to protoplanetary discs, taking into account also the role of the drag force between gas and dust. When the second component is considered, the outcome of GI can significantly change. In the protostellar scenario, the three fundamental parameters are the relative concentration of the dust component to that of the gas

the dust relative temperature

where c_d is the dust dispersion velocity, and the Stokes number

Two regimes of instability can be identified, depending on whether the instability is triggered by the gas or the dust. When the instability is driven by the gas component, the most unstable wavelength is close to the gas-only one, and the role of dust is negligible. Conversely, when the instability is controlled by the dust, the perturbation wavelength is much smaller. This different kind of instability happens when the dust is sufficiently cold, decoupled, and abundant: the threshold between these two regimes is given by (Bertin & Romeo 1988)

or, taking into account the drag interaction (Longarini et al. 2023)

Dust-driven instability could have important consequences for planet formation because when the most unstable wavelength (i.e. the Jeans wavelength) is smaller, the Jeans mass will be small as well. In the gas-only model, the Jeans length is

where we used |$\rho _\mathrm{ g} = \Sigma _\mathrm{ g}/\sqrt{2\pi }H$|⁠, and the Jeans mass is

When the instability becomes dust driven, the value of the most unstable wavelength changes dramatically. It is possible to compute it through the dispersion relation of the two-fluid-component model with drag force (Longarini et al. 2023), and its value depends on the Stokes number, the relative temperature, and the dust-to-gas ratio. In particular, the gas-only limit is obtained for St → 0, ϵ → 0, and ξ → 1. For convenience, we define the ratio between the Jeans wavelength of one-fluid-component model and of the two-fluid one with drag force as Λ(ϵ, ξ, St): the value of Λ is between 0 and 1, where Λ = 1 is the one fluid limit. Indeed, by computing the most unstable wavelength, it is possible to extract Λ. Hence, the Jeans mass of the two-component-fluid model is

In Appendix  A, we show plots of Λ³ as a function of (ϵ, ξ, St).

3 NUMERICAL SIMULATIONS

In this work, we perform numerical SPH simulations of gas and dust protostellar discs using the code phantom (Price et al. 2018). This code is widely used in the astrophysical community to study gas and dust dynamics in accretion discs (Ragusa et al. 2017; Ceppi et al. 2022), both in a single fluid mixture (Veronesi et al. 2020) or dust-as-particles approach (Aly et al. 2021). In this work, we use the dust-as-particles formulation.

3.1 Two-fluid gas and dust mixtures

The dust formulation is based on the continuum fluid equations in the form

where the subscripts g and d refer to gas and dust properties. We define the stopping time, that is given by

where the drag coefficient K depends on the physical drag regime the system is in. In general, it is given by

where C_D is defined as before.

In the phantom implementation, the two phases are modelled as two distinct sets of particles: hereafter, we adopt the convention from Monaghan & Kocharyan (1995), and refer to gas particles with the subscripts a, b, c and to dust particles with i, j, k. Gas and dust densities are computed by weighted summation over the particles of the same type according to

where the kernel W is the same for gas and dust, h is the smoothing length, and h_fact = 1/2. The drag terms of the equations of motion are discretized by using a ‘double hump’ kernel (Fulk & Quinn 1996) that, instead of the bell-shaped kernel, goes to zero at r = 0 and has a peak at r/h ≲ 1.¹ In these simulations, we are using the velocity reconstruction procedure presented in Price & Laibe (2020).

3.2 Heating and cooling

In order to take into account the effect of cooling or heating phenomena, we write the complete equation for the evolution of gas internal energy e

where the first term on the RHS is the PdV work, the second is a heating term due to the shock viscosity, the third is the cooling of the disc, and the last term is the drag heating term. In this work, we assume an adiabatic equation of state. For an ideal gas, it is possible to link pressure and density as follows:

where γ = 5/3 and c_g is the gas sound speed, that is initialized as a power law c_g ∝ R^−0.25.

The shock viscosity term can be written as

where c₁, c₂ are two numerical factors, whose dimension is a specific energy per unit time, α^AV and β^AV are respectively the linear and the quadratic viscosity coefficients, and h_g/H_g is related to the numerical resolution.² The viscosity term is dissipative, so it heats the disc. In the phantom simulations, we are not using the disc-viscosity flag, meaning that the shock capturing viscosity is not described by an α_SS (Shakura & Sunyaev 1973) prescription. We did so since in these systems the main driver of angular momentum transport is GI.

For the cooling, we use the prescription from Gammie (2001) and Rice et al. (2004), in which the cooling time t_cool is proportional to the dynamical time, with a factor of proportionality β_cool

Under the assumption that the transfer of angular momentum driven by gravito-turbulence occurs locally (Lodato & Rice 2004; Béthune, Latter & Kley 2021), we can relate the cooling parameter to an effective α-viscosity parameter

We need to introduce a cooling prescription for the system in order to trigger GI. Finally, the drag heating term is

Currently, phantom neglects any thermal coupling between the dust and the gas, aside from the drag heating.

3.3 Numerical set-up

We perform simulations with three different disc-to-star mass ratios M_d/M_⋆ = {0.05, 0.1, 0.2}, three different cooling times β_cool = {8, 10, 15} and two different dust particle sizes, for a total of 18 simulations. See below for further details. We first initialize a gas-only disc around a solar mass star, with R_in = 0.25 au, R_out = 25 au, and Σ_g ∝ R⁻¹. We set the aspect ratio so that Q_ext = 2 initially and because of the cooling, it decreases, eventually reaching Q = 1. The shock viscosity coefficients α^AV = 0.1, β^AV = 0.2 as in Rice et al. (2005). As a test, we performed a simulation with β^AV = 2, and we did not find any differences in terms of the relevant quantities in this work. The self-gravity of both gas and dust is taken into account, as well as the dust back-reaction. We performed simulations at two different numerical resolutions: the standard runs are performed with N_g = 10⁶ and the high-resolution ones with N_g = 2 × 10⁶. In both cases, N_d = N_g/5, where N_g and N_d are the number of gas and dust particles, respectively. We verify that the results are consistent with the different resolutions.

We let the system evolve for an outer thermal time (β_coolΩ⁻¹ at the outer radius), and then we add dust particles and evolve for a further five outer dynamical times (i.e. 5 × 10³ inner dynamical times). The dust particles are added proportional to the gas distribution.³ Distributing dust particles proportional to the gas is a valid assumption only for St < 10. Shi et al. (2016) and Baehr & Zhu (2021) showed that for uncoupled particles gravitational interactions and stirring become quite relevant, and hence dust distribution is closer to a uniform one. This is particularly clear for St ∼ 100, where, in our simulations of large dust grains, although the initial distribution is proportional to the gas one, we observe that the spiral is less prominent compared to smaller grains.

Dust back-reaction and dust self-gravity are always taken into account. Since the aim of the work is to study the effect of the drag force in GI, we use two different dust sizes: a larger one, to reproduce weakly coupled solid particles, and a smaller one, to study marginally coupled particles. Since the disc mass is different across our set of simulation, we chose to adapt the particles’ size in order to obtain the same Stokes number distribution. To do so, we computed the radially averaged Stokes number as a function of the particle size for different disc-to-star mass ratio, starting from the initial conditions of the gas disc, taking into account the transition between Epstein and Stokes regime. We decided to choose the small particles’ size so that the radially averaged Stokes number 〈St〉 = 8, and for the large ones 〈St〉 = 40: in this way, we effectively cover a Stokes number range from 1 to 10 with smaller grains, and from 10 to 100 for larger ones. One exception is that the small dust particles for the highest disc-to-star mass ratio simulations are chosen to have 〈St〉 = 16 for computational reasons. In addition, as we will show in the next section, a higher disc-to-star mass ratio makes the dust unstable, and collapse can happen: for that reason, dust in simulations S16, S17, and S18 is evolved only for an outer dynamical time. In every simulation, dust intrinsic density is fixed ρ₀ = 5 g cm⁻³.

Self-gravitating discs may be more radially extended than our models, which can be rescaled. If we rescale the outer radius of a factor λ, how does the dust particles’ size need to be rescaled? Since |$\text{St}\propto {s}/{\Sigma }\propto sR_\text{out}^2$|⁠, if we change the outer radius according to |$R_\text{out}^{\prime } = \lambda R_\text{out}$|⁠, in order to maintain the same Stokes number, the corresponding rescaling for dust particles size should be s′ = λ⁻² s. Hence, if we consider a larger disc with |$R_\text{out}^{\prime } = 10R_\text{out}=250$| au, the corresponding particle sizes should be rescaled as s′ = s/100. Dust properties are summarized in Table 1, where we also include the rescaled dust particles’ size. Snapshots of the hydrodynamical simulations are shown in Figs 1–3.

4 ANALYSIS AND RESULTS

In this section, we present the analysis of the numerical simulations. Since the simulations simulate gas and dust as two different sets of particles (Laibe & Price 2012a, b), gas and dust particles are treated as two different sets of particles; thus, they occupy different locations and carry their own physical information. In order to obtain properties that depend on both gas and dust, such as dust-to-gas ratio, or Stokes number, we interpolate gas properties to the location of dust particles. In addition, since dust is modelled as a pressure-less fluid, it has no internal energy and no thermal sound speed. However, stirring phenomena induce a velocity dispersion onto dust particles (Youdin & Lithwick 2007): to obtain this quantity, we compute it with an SPH interpolation over neighbouring dust particles, via

For our analysis, we will mainly focus our attention on the dust-to-gas ratio ϵ, the relative temperature between gas and dust ξ = (c_d/c_g)², the Stokes number St, the cooling factor β_cool, and the disc-to-star mass ratio M_d/M_⋆. These parameters are related to different physical phenomena: the dust-to-gas ratio and the relative temperature trace dust trapping and dust excitation, respectively, the cooling factor is linked to the gas spiral amplitude (and hence to the strength of GI), the Stokes number determines the power of the aerodynamical coupling, and the disc-to-star mass ratio is connected to the spiral morphology.

4.1 Dust trapping: the ϵ parameter

GI spiral arms trap dust particles (Dipierro et al. 2015a), since they are both pressure maxima and gravitational potential minima. We use the simulations to quantify how this phenomenon depends on the model parameters. Fig. 4 shows the distribution of the dust-to-gas ratio for different β_cool and M_d/M_⋆, for a set of simulations with large dust particles. The initial value of the dust-to-gas ratio is 10⁻². In Fig. 4, the higher tail of the distributions reaches values of ≳10⁻¹, implying that dust concentration by up to approximately an order of magnitude is happening. Fig. 5 shows a comparison between the dust-to-gas ratio distributions of large and small particles. As expected, dust concentration in spiral arms is stronger for smaller particles, and it can approach values of the order of unity, since their aerodynamical coupling with gas is stronger. The strength of dust trapping is determined by both the aerodynamic coupling between gas and dust and the gravitational potential of gas spiral arms. The combined effect of gravitational and drag interaction is maximized when St ≃ Q ≃ 1 (Baehr & Zhu 2021); thus, in our simulations, smaller particles reach higher values of the dust-to-gas ratio. No particular correlations are found between ϵ and the disc-to-star mass ratio, while there is a slight dependence on the cooling factor. In order to understand this relationship, Fig. 6 shows the quantity δ Σ/Σ₀ for gas (orange dots) and dust (large particles – blue dots, small particles – green dots) as a function of the cooling factor β_cool. The quantity Σ₀ is the azimuthally averaged surface density at a fiducial radius of 10 au, and the quantity δΣ is its standard deviation. For the gas, it is known that |$\delta \Sigma _\mathrm{ g}/\Sigma _{\mathrm{ g}0}\propto \beta _\text{cool}^{-1/2}$| (Cossins, Lodato & Clarke 2009), and we recover this behaviour in our simulations. For the dust, the situation is different. We do not find any evident correlation between the density contrast and the cooling factor. So, if we assume that

the ratio between these two quantities is

where Σ_d0/Σ_g0 = 1/100, meaning that |$\epsilon \propto \beta _\text{cool}^{1/2}$|⁠. This is the positive correlation we found before. Why do we not find any evident correlation between the dust density contrast and the cooling factor? Physically, the dust experiences the effect of the gas cooling through gravitational and drag forces. When St ≪ 1, dust and gas particles are indistinguishable, and so |$\delta \Sigma _\mathrm{ g}/\Sigma _{\mathrm{ g}0} = \delta \Sigma _\mathrm{ d}/\Sigma _{\mathrm{ d}0} \propto \beta _\text{cool}^{-1/2}$|⁠. For higher Stokes number, the drag force is weaker, and the dust is less influenced by the gas cooling. In this case, we expect the relationship between δΣ_d/Σ_d0 and β_cool to be flatter than |$\beta _\text{cool}^{-1/2}$|⁠: if this condition is respected, the correlation between ϵ and β_cool will be positive. In general, δΣ_d/Σ_d0 is a function of both the cooling factor and the Stokes number.

4.2 Dust excitation: the ξ parameter

To investigate dust excitation by spiral arms, we study how the relative temperature ξ = (c_d/c_g)² varies as a function of the simulation parameters. Fig. 7 shows the distribution of the dust relative temperature ξ for different values of β_cool and M_d/M_⋆ for a set of simulations with large dust particles. We observe that for the simulations in which dust collapse is not happening, the dust dispersion velocity reaches very quickly (<1 outer orbital time) a steady value. The relative temperature shows a negative correlation with both the disc-to-star mass ratio and the cooling factor. For the disc-to-star mass ratio, this can be understood by considering the relationship with the spiral morphology: M_d/M_⋆ is inversely proportional to the azimuthal wavenumber m (Cossins et al. 2009); hence, massive discs have fewer spiral arms. Since dust is excited because of ‘spiral kicks’ (Walmswell et al. 2013), the lower the number of spiral arms, the less the dust is excited. For the cooling factor, the negative correlation can be understood in two ways. First, the cooling rate β_cool is linked to the amplitude of the spiral perturbation according to equation (36). Since gas spiral arms excite dust particles by kicking them every passage, the higher is the perturbation, the stronger is the kick, and so the excitation. Secondly, in gravito-turbulent regime, transport of angular momentum is driven by the spiral perturbation, and it is possible to define an α-viscosity coefficient related to the cooling rate (equation 33). The height of a dust layer is determined by the interaction with the gas and by the vertical diffusion. In the hypothesis that the vertical diffusion coefficient is equal to the azimuthal one, we can obtain the height of the dust layer as

If we assume that α = α_GI, the dust layer height, and thus the dust dispersion velocity, is inversely proportional to β_cool.

Fig. 8 compares the relative temperature distributions of large and small dust particles. Small particles are colder than the large ones, and their relative temperature is almost completely enclosed in the interval ξ_small ∈ [10⁻², 1], meaning that their random motions are subsonic. On average, the distribution of ξ_small is shifted by one order of magnitude compared to ξ_large. This behaviour is in agreement with what Booth & Clarke (2016) found: larger particles tend to be dynamically hotter, because the kicks of the gas spiral are more effective, while if the coupling with the gas is stronger, the kicks are damped because of the drag force.

In principle, the observed trend would imply that increasing β_cool would lead to an arbitrarily thin dust layer, eventually causing gravitational collapse. This would be a direct analogue of the classical Goldreich & Ward (1973) mechanism for planetesimal formation in a (weakly) self-gravitating context, and in a more realistic model it is likely to be limited by the excitation of shear turbulence (Cuzzi, Dobrovolskis & Champney 1993). However, there is an upper limit for β_cool set by the value above which the transfer of angular momentum would be driven by some process other than GI. To compute an estimate of the maximum cooling time, we require α_GI to be larger than 10⁻³. Observations of protoplanetary discs that are not expected to be self-gravitating suggest that this is a reasonable upper limit to the strength of turbulence (Flaherty et al. 2017). Thus, we obtain

which corresponds to a minimum density perturbation

Table 2 summarizes the relationships we have discussed in these paragraphs. A broader collection of histograms can be found in Appendix  C.

Booth & Clarke (2016) studied the relationship between the dust excitation and both the cooling and the Stokes number. They found that c_d ∝ β^−1/2St^1/2v_k, where v_k is the Keplerian speed. To compare with Booth & Clarke (2016), we use our data to reproduce figs 7 and 13 of their paper, which show a relationship between the dust velocity dispersion and β_cool and St. To do so, we divided the particles into equally spaced intervals of Stokes number and, for each particle, we computed the mean value of |$c_\mathrm{ d}/c_\mathrm{ g} = \sqrt{\xi }$|⁠. The comparison is shown in Fig. 9, where we show the results of simulations with M_d/M_⋆ = 0.05, for both standard- and high-resolution cases. The previously derived relationships with the Stokes number (the left-hand panel) and with the cooling factor (the right-hand panel) are well recovered. Using our two-fluid algorithm, it is too computationally expensive to analyse properly the case of St < 1. However, in this case we expect that as the aerodynamical coupling with the gas is stronger, c_d/c_g should increase, eventually reaching c_d = c_g for St → 0. This growth for St < 1 has been shown by Booth & Clarke (2016).

4.3 Two-fluid instability

In this section, we apply the two-fluid instability theory presented in Section 2. The gas-only and gas-and-dust models for GI have both been developed within a linear framework; hence, in principle, the quantities ϵ, ξ, and St should be evaluated in the unperturbed state. However, in this work, we are evaluating them in the perturbed one. Although not completely self-consistent, it gives us an idea of the most unstable regions of the disc. Fig. 10 shows the distribution of large (blue) and small (orange) dust particles in the (ξ, ϵ) diagram: the black lines correspond to the dust-driven GI threshold for St = ∞ (equation 16, solid line) and for St = 0.5 (equation 17, dashed line). We choose St = 0.5 as a minimum value since the number of particles with Stokes number lower than this is negligible. The particles above the region are in a dust-driven GI regime. We find that the number of small particles for which the instability is dust driven is greater compared to large ones: this is because small particles have both larger dust-to-gas ratio and lower dispersion velocity. In addition, the number of dust-driven particles increases with the cooling factor and with the disc-to-star mass ratio, as already discussed in previous sections. To understand the spatial location of these particles in the disc, Fig. 11 shows the particles that satisfy condition 16 superimposed on the total density map. As expected, the most unstable regions of the disc are not randomly distributed, but correspond with the spiral arms.

Fig. 12 shows the value of the Jeans mass for M_d/M_⋆ = 0.05 as a function of the Stokes number. The curve can be divided into two parts: for St ∈ (0, 5), the Jeans mass is decreasing with the Stokes number, reaching its minimum at about St ∼ 3. This happens because for St → 0, the particles are indistinguishable, and the instability is gas driven. It is important to notice that the number of particles with small Stokes number is low, hence the binning in Stokes number presents a considerable scatter. By increasing the Stokes number, ξ decreases and the two fluids behave more and more differently. When the Stokes number is approximately 1, the relative temperature is a minimum and the dust-to-gas ratio is high, so the instability becomes dust driven. Otherwise, for St > 5, the Jeans mass increases with the Stokes number. Indeed, the relative temperature increases, and the dust-to-gas ratio decreases. Hence, the system transitions from dust-into-gas driven instability, again, eventually reaching the gas-only component model value.

5 DISCUSSION

5.1 Dust collapse

The simulations with M_d/M_⋆ = 0.2 and small dust particles do not reach five outer orbits, since the simulation stops due to the onset of dust collapse. This happens because the stopping time of collapsing dust particles becomes smaller than the time-step of the code. Indeed, the stopping time is inversely proportional to the total density ρ_tot = ρ_g + ρ_d, and since the dust density is increasing, because of the collapse, the stopping time tends to zero. The top panel of Fig. 13 shows the maximum dust density as a function of time for small and large dust particles, in the run with M_d/M_⋆ = 0.2 and β = 15. While large dust particles reach a quasi-steady state, the small particle density exponentially increases in the first orbit. This is the signature of dust collapse. This phenomenon is visible in simulations S16, S17, and S18, and happens only in the dust component. The bottom panel of Fig. 13 shows a comparison between gas and small dust averaged density as a function of time, for M_d/M_⋆ = 0.2 and β = 15. At t = 0, 〈ρ_d〉 = 10⁻²〈ρ_g〉. Whereas the gas average density is constant with time, the dust density increases because of dust trapping, eventually exceeding that of the gas. This means that any clumps forming from this mechanism would be substantially made up of solids, and would likely be identified with the rocky core of a giant planet. However, in this work, we do not want to characterize the outcome of this collapse, which is a complex topic. Indeed, simulations of planetesimal collapse (Nesvorný et al. 2021) show that a rotating self-gravitating cloud of dust does not monolithically collapse, meaning that it is not possible to directly equate the cloud mass with the planetary core one.

To identify and analyse dust clumps in more detail, we define the numerical conditions that should be respected for a clump to be physical and not affected by resolution. For a clump radius r_clump, the smoothing length of the dust particle h_i should be less than a fraction of the clump radius, in order to be resolved. This condition translates into h_i < ηr_clump, where η is less than unity. We take η = 1/2. Physically, a gravitationally bound clump is a collection of particles whose thermal support does not balance the gravitational one. We define the thermal and gravitational energy of particles inside r_clump as follows:

where φ is the gravitational softening kernel and r_ij = |r_i − r_j|. We used a cubic spline softening kernel, and the detailed expression can be found in the appendix of Price & Monaghan (2007). According to the virial theorem, if the force between any two particles can be described in terms of a potential energy Φ ∝ rⁿ, where r is the distance between two particles, the equilibrium state respects the following condition:

where T is the kinetic energy. In the case of gravitational interaction, the virial theorem reads 〈T〉/〈Φ〉 = −1/2. We define a clump as a region of the space where the dimensionless quantity α_J = −e_th/e_gr < 1/2. Then, in order to be sure that the collapse is physical, and not artificial, we verify that h_g < h_d in the region where there is the dust clump. The last condition requires that the resolution of the gas in the region where there is a possible clump should be higher compared to the dust resolution.

To summarize, the conditions under which we define a dust clump are the following:

• h_{g, i} ≤ h_{d, i},

• h_{d, i} < ηr_clump,

• α_J < 1/2.

In simulations S16, S17, and S18, there are particles that respect the previous conditions, implying that dust collapse has happened. Table 3 shows the mass of the clumps, obtained by summing the mass of each particle that is gravitationally bound and the one obtained from analytical theory from equation (20). The masses are broadly in agreement with analytical expectations. For S16, for example, the derived Jeans mass is of the order of an Earth mass. In general, the mass of the clump computed from the simulation is smaller compared to the one expected from the analytical theory: this is not surprising, since with the simulation we are only able to appreciate the initial phase of the collapse. Indeed, as soon as the stopping time is smaller than the time step, the simulation stops. To avoid this problem, one could decrease the time step of the code, but it is computationally expensive. Otherwise, one could not consider the dust density when computing the stopping time: this procedure has been applied in previous works to increase the velocity of the simulations (Poblete, Cuello & Cuadra 2019; Longarini et al. 2021) but this approximation is valid only for small dust-to-gas ratios. While there this was justified, here this is not possible, since dust is collapsing and the dust-to-gas ratio becomes higher than unity. Finally, to study the early evolution of clumps with an SPH code, it would be possible to simulate them as sink particles: so far, the creation of dust sink particles is not possible in phantom.

5.1.1 Gas–dust coupling during the collapse

The bottom panel of Fig. 13 shows that the collapse happens only in the dust component, and the gas is not influenced. This could sound surprising: indeed, in high-density regions, we expect the stopping time (equation 25) to be small. So, why is the dust collapse not influencing the gas? The degree of coupling is measured with the Stokes number, which compares the strength of the drag force with the ones that are acting on the particle. Usually, in a flat protoplanetary disc, the Stokes number is computed as the ratio of the stopping time and the dynamical time, which is the typical time-scale of a particle orbiting around a central object at a distance R. However, in this situation, the dust clump is driving the dynamics of the surrounding particles, and hence we should compare the stopping time with the free fall time in order to understand the degree of coupling of particles. The typical time-scale of the infall of a spherically symmetric distribution of mass is

Comparing the stopping time and the free fall time in the simulations where dust collapse is happening (Table 4), we obtain that the particles in the collapsing region are uncoupled, since the ratio between the two time-scales is higher than one.

It is possible to quantify the critical density a clump should reach so that t_ff < t_dyn, that is

when a clump reaches this density, the evolution of the surrounding particles is determined by the clump, and not by the star anymore. From that point, the aerodynamical coupling should be quantified by taking the ratio between the stopping time and the free fall time. Hence, for ρ < ρ_crit, St ∝ ρ⁻¹, since the dynamical time does not depend on the density. For ρ > ρ_crit, the scaling changes since the free fall time depends on the density St ∝ ρ^−1/2, and so does the degree of coupling.

5.2 Gravitational instability in the context of protoplanetary disc evolution

In this work, we focused on cooling-driven GI, but it is also possible to trigger it through infall. Kratter et al. (2010) found that in the infall-driven case the strength of the spiral perturbation is controlled by two dimensionless parameters, a thermal one, which relates the infall mass accretion rate |$\dot{M}_\text{inf}$| to the characteristic sound speed of the disc, and a rotational one, which compares the relative strength of rotation and gravity in the core. Obviously, the higher the accretion rate, the stronger is the spiral perturbation in the disc. We can naively associate the dust evolution we have modelled in the limit of fast (slow) cooling with high (low) infall rate. This association is expected to be qualitatively correct; however, a detailed comparison between these two regimes is given by Kratter & Lodato (2016).

By studying the non-linear evolution of gas and dust in protoplanetary discs, we found that the instability conditions for the two components are different. It is well established that spiral fragmentation occurs in fast cooling gas discs, and it is possible to define a critical value of β_cool below which fragmentation occurs. Simulations of cooling-driven fragmentation (Gammie 2001; Rice et al. 2005; Lodato & Clarke 2011; Meru & Bate 2012) currently suggest that β_min ≃ 3 (Deng, Mayer & Meru 2017). For the dust, Booth & Clarke (2016) found that dust becomes more unstable for higher β_cool, and we confirm this trend with 3D simulations. The differing behaviour of gas and dust suggests an interesting evolution in the outcome of GI within protoplanetary discs. During a first stage, at the beginning of the disc lifetime, we expect a very massive disc system characterized by strong GI, caused by the high infall rate from the molecular cloud. If conditions allow gas fragmentation, because of the high Jeans mass, low mass stellar companions can be formed (Kratter et al. 2010). A second stage is characterized by a less-massive disc, with lower infall rate, or equivalently longer cooling time. If conditions during this epoch trigger GI, it will lead to dust-driven fragmentation that could be responsible for the formation of rocky cores of giant planets. Then, in the third stage, there is a protostar surrounded by a planet-hosting disc, characterized by substructures such as gaps, rings, or planetary spirals. In this stage, GI is not effective anymore because the disc mass is small and the transport of angular momentum is controlled by disc winds (Tabone et al. 2022) or other, non-self-gravitating, sources of turbulence (Lesur et al. 2022). A schematic view of these stages is given by Fig. 14.

5.3 Application to an actual case: HL Tau

HL Tau is a young (<1 Myr) protostellar system that shows axisymmetric structures (gaps and rings) in dust continuum emission (ALMA Partnership 2015). The origin of rings and gaps is usually attributed to planet–disc interaction (Lin & Papaloizou 1986), with Dipierro et al. (2015b) finding that three protoplanets with masses M_p1 = 61 M_⊕, M_p2 = 83 M_⊕, and M_p3 = 170 M_⊕ at R_p1 = 13.2 au, R_p2 = 32.3 au, and R_p3 = 68.8 au best match the observations. The formation of super-Earth mass planets at large radii, in such a young system, is a challenge for core accretion theory. If planets can form via the dust-induced collapse mechanism we have discussed in this work, HL Tau is a plausible example of what the resulting planetary system might look like. We note that the inferred trend of increasing mass with radius can be interpreted within our model. In a simplistic way, we assume that St ∝ Σ⁻¹ ∝ R, and that |$\xi \propto \beta _\text{cool}^{-1}\text{St}\propto R^{4.5}R=R^{5.5}$|⁠, where we have supposed a realistic cooling law, according to which β_cool ∝ R^−4.5 (Rafikov 2005; Clarke & Lodato 2009). Supposing that the dust-to-gas ratio is constant, the Jeans mass of the gas-and-dust fluid model will be an increasing function of the radius, since the relative temperature is smaller for particles closer to the central star. So, within this hypothesis, we expect that the mass of the dust-driven GI protoplanets will be an increasing function of the radius. We want to point out that this is just a qualitative argument. Indeed, to thoroughly investigate whether HL Tau planets can be formed through dust-driven GI, one should properly model the system. In addition, planet migration and planet accretion should also be considered.

6 CONCLUSIONS

Self-gravitating gas discs may be ubiquitous during the Class 0/I phases of YSO evolution (e.g. Lin & Pringle 1990; Xu & Kunz 2021). Observations suggest that the self-gravity can in some systems be strong enough as to trigger fragmentation (Tobin et al. 2016), while in other cases, such as the massive disc of the more evolved IM Lup system (Lodato et al. 2023) the instability is expected to be more gentle. If that fragmentation occurs in the gas, as in the L1448 IRS3B system, the outcome is typically a star or brown dwarf formation. Lower mass objects can be formed if the fragmenting fluid is instead the solid component of a two-fluid self-gravitating disc, given disc conditions that allow the collisional growth of dust to small macroscopic dimensions while the disc remains self-gravitating. Analytical estimates suggest that planetary cores of ∼1−10 M_⊕ could form from this mechanism at large orbital radius, with properties that could be identified with the population of ALMA-inferred disc-embedded planets (Andrews 2020).

In this paper, we presented the results of SPH simulations of gas and dust in protoplanetary discs, studying the role of aerodynamic coupling in the context of GI. We analysed our results in the framework of two-fluid GI, and compared our findings with previous numerical works, finding generally good agreement.

Our main results can be summarized as follows:

• We studied the relationship between the dust-to-gas ratio ϵ, the relative temperature between gas and dust ξ, and the cooling factor β_cool, the disc-to-star mass ratio, and the Stokes number. We found that the dust-to-gas ratio increases with the cooling factor and decreases with the Stokes number, and that the relative temperature increases with the Stokes number and decreases with the cooling factor and the disc-to-star mass ratio. It is possible to explain these relationships by considering the interaction between dust particles and gas spiral arms. We compared our findings with Booth & Clarke (2016) and found good agreement.

• We investigated the role of dust in GI, and found that the most unstable regions of the disc are the spiral arms, where the instability tends to be dust driven. In addition, we studied the relationship between the theoretical Jeans mass and the Stokes number. We found that the Jeans mass – when the instability is dust driven – can reach values of the order of the Earth mass.

• We observe three cases of dust collapse in our set of simulation, which occur (as expected) in the high disc mass models. The values of the clump masses obtained numerically are close to those predicted by linear theory.

In applying our results to the possibility of early planet formation in Class 0/I discs, the main prerequisite is the requirement that dust is able to grow via coagulation to a large enough Stokes number, with a top-heavy particle mass function, in a short enough time. Our simulations that explicitly exhibited dust collapse had solid particles with an average Stokes number 〈St〉 = 16, which would correspond (rescaling our simulations to a disc size of R_out = 250 au) to particles with sizes between a few cm and a few metres. Fragmentation and radial drift pose barriers to growth to the required sizes (Birnstiel, Dullemond & Brauer 2010), and further work is needed to assess whether there are circumstances where the required Stokes numbers can be reached. Simulations with a constant Stokes number, along with runs with St ∼ 1 (a regime which is numerically difficult to access using our code), would also help to better define the regime where dust can fragment in a gas disc that is itself stable against fragmentation.

In an evolutionary context, our results imply that planet formation – if it is able to occur via the mechanism of dust-dominated gravitational disc instability – is likely to occur towards the end of the self-gravitating phase. This is when the competing effects of particle trapping and particle excitation are jointly most favourable for collapse. The masses and orbital radii of the planets formed via dust collapse are qualitatively in agreement with those inferred for the HL Tau system, and we speculate that they may form from the collapse of solids in the spiral arms of a formerly self-gravitating protostellar disc.

ACKNOWLEDGEMENTS

The authors thank the referee for useful comments and suggestions that significantly improved the quality of the work. In this work, we used splash to create figures of the hydrodynamical simulations (Price 2007). This project and the authors have received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no.823823 (DUSTBUSTERS RISE project). CL acknowledges support from Fulbright Commission through VRS scholarship. PJA acknowledges support from NASA TCAN award 80NSSC19K0639, and from award 644616 from the Simons Foundation. DJP acknowledges Australian Research Council funding via DP220103767 and DP180104235. The authors thank Sahl Rowther, Benedetta Veronesi, Cathie Clarke, and Richard Booth for useful discussions. We thank Stony Brook Research Computing and Cyberinfrastructure, and the Institute for Advanced Computational Science at Stony Brook University for access to the SeaWulf computing system, supported by National Science Foundation grant # 1531492.

DATA AVAILABILITY

The data presented in this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: VALUE OF Λ

In this appendix, we discuss the value of Λ used to compute the Jeans mass for the two-fluid-components model. The value of Λ can be obtained through the dispersion relation of the two-fluid-component model with drag force (Longarini et al. 2023), by computing the most unstable wavelength having fixed (ϵ, ξ, St). Fig. A1 shows the value of Λ³ as a function of (ϵ, ξ, St). The curves show a sharp fall: this is the point where instability becomes dust driven. When the instability is gas driven, the value of Λ³ decreases with the dust-to-gas ratio. Then, it starts increasing since the mass of the dust becomes higher and higher, significantly contributing to the total density. Conversely, the value of Λ³ increases with the relative temperature both in gas and dust-driven regimes. Finally, Λ³ decreases with the Stokes number.

APPENDIX B: RESOLUTION REQUIREMENTS

In this appendix, we compute the minimum number of dust particles to resolve the expected Jeans mass: to study this problem, we refer to Lodato & Clarke (2011). In order to resolve the Jeans length, we require the smoothing length to be smaller than the disc height: so, the condition for dust particles is

Dust smoothing length and height can be written as a function of gas properties as follows:

where we assumed M_dust/M_gas = 1/100. By considering the gas disc to be marginally unstable (Qg = 1), condition B1 becomes

where ϵ is the dust to gas ratio and m(r) = Σ_gr²/M_⋆. If we write the disc mass as |$M_\mathrm{ d} = \pi r_\text{out}^2 \Sigma _\mathrm{ g}$|⁠, m(r) ≃ 1/π(M_d/M_⋆), thus the resolution requirement is

Fig. B1 shows the minimum number of dust particles required to resolve the Jeans length as a function of the relative temperature ξ for M_d/M_⋆ = {0.05, 0.1, 0.2}, for a fixed dust-to-gas ratio ϵ = 0.1: even for extreme values of ξ, the standard resolution we have chosen (N_d = 2 × 10⁵) is sufficient to resolve the Jeans length. One should note that we never reach ξ ∼ 10⁻², and in addition, we have underestimated the dust-to-gas ratio: indeed, inside spiral arms it could reach values >0.1, making our N_d choice even safer.

APPENDIX C: LARGE AND SMALL DUST PARTICLES COMPARISON

Fig. C2 shows the complete set of histograms of the dust-to-gas ratio and the relative temperature, for small and large dust particles.