© The Authors 2026

1 Introduction

In the standard scenario for planet formation, the coagulation of micron-sized grains leads to the formation of millimeter-sized dust (Dullemond & Dominik 2005). However, any further growth of particles above the millimeter-size range is rendered difficult because of bouncing (Zsom et al. 2010) and fragmentation (Zsom et al. 2010) barriers. An emerging scenario that is capable of overcoming these barriers is based on the assumption that 100-km sized planetesimals form directly through the streaming instability (SI; Youdin & Goodman 2005; Youdin & Johansen 2007), which occurs as a result of angular momentum exchange between the gas and dust components, as well as rotation. It has been shown that during the nonlinear evolution of the SI, dust particles can eventually directly concentrate into clumps or filaments, which can subsequently become gravitationally unstable to form 100–1000 km-sized bodies (Johansen et al. 2009; Simon et al. 2016; Abod et al. 2019). The ability for solids to experience a regime of strong clumping depends mainly on the dimensionless stopping time St and local dust-to-gas ratio, ϵ, which can also be quantified through the dust-to-gas surface density ratio, Z = Σ_d/Σ_g, where Σ_d (resp. Σ_g) is the dust (resp. gas) surface density. For Z ~ 5 × 10⁻³, strong clumping requires St ≳ 0.1 in a disk where no external turbulence operates, whereas for St ≲ 0.01 strong clumping is triggered for a critical solid abundance of a few percent (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021). In any case, this suggests that a level of dust growth is needed for the SI to enter the strong clumping regime (Drążkowska & Dullemond 2014).

Nonetheless, taking into account the effect of dust growth on the SI in a self-consistent fashion is a difficult task. The process of dust growth is indeed expected to be a sensitive function of the relative collision velocity between grains, which depends on multiple processes taking place in the protoplanetary disk. As a consequence, the resulting dust size distribution is not known a priori and remains uncertain. The generalization of the classical, monodisperse SI to multiple dust sizes has been examined by a number of recent studies, but on the basis of prescribed dust size distributions. These have shown that the growth rate of the polydisperse, multispecies SI (PSI) tends to decrease as the number of species is increased (Krapp et al. 2019) and becomes very small in the limit of a continuous size distribution (Paardekooper et al. 2020; Zhu & Yang 2021). However, it was also found that a top-heavy size distribution, resulting for instance from a pure coagulation process, leads to PSI linear growth rates equivalent to those of the classical SI (McNally et al. 2021).

The influence of the process of coagulation itself on the nonlinear development of the monodisperse SI has been investigated by Ho et al. (2024) using stratified simulations. These authors found that coagulation tends to broaden the range of Stokes numbers and dust-to-gas ratios over which the SI can be triggered, due to ability for coagulation to promote dust growth. Conversely, it has even been proposed that dust growth can be significantly boosted in the strong clumping regime of the SI (Tominaga & Tanaka 2023, 2025). This arises because dust-loading increases the inertia of the dusty fluid, resulting in a reduction of the effective sound speed and lower collision velocities between dust grains. These results suggest that a feedback loop between coagulation and the development of the SI may exist and this was recently studied by Carrera et al. (2025).

In this paper, we examine the effect of the coagulation instability (CI) on the linear growth and nonlinear saturation phase of the SI (Tominaga et al. 2021). The CI is a consequence of the interplay between the dust coagulation process and radial drift of dust particles and primarily related to the dependence of mass growth rate on dust density. As described by Tominaga et al. (2021), if we assume a positive dust density perturbation, the resulting increased coagulation efficiency leads to a variation of the dust size. This subsequently causes an enhanced radial drift speed which tends to amplify the initial density perturbation, resulting in a feedback loop. It has been shown that the nonlinear evolution of the CI can enable the Stokes number of solids to rapidly reach unity, even for relatively modest dust-to-gas ratios ϵ ~ 10⁻³ (Tominaga et al. 2022a,b). Therefore, the CI might represent a promising mechanism to extend the parameter space (in terms of dust-to-gas ratios and Stokes numbers) over which the SI can operate efficiently.

To investigate the possible effect of the CI on the SI in more detail, we performed 2D axisymmetric, unstratified, shearingbox simulations in which coagulation was modeled using the single-size approximation. Our main aims are to (i) assess whether the CI can help the SI enter the regime of strong clumping and (ii) evaluate to what extent the turbulence driven by the CI subsequently impacts the nonlinear development of the SI.

The paper is organized as follows. In Sect. 2, we present the governing evolution equations of the gas and dust components. We then describe in Sect. 3 the numerical setup and the initial conditions that are used in the simulations, whose results are presented in Sect. 4. We discuss our results in Sect. 5 and present a summary of our findings in Sect. 6.

2 Physical model

2.1 Governing equations

We modeled a local patch in a protoplanetary disk within the shearing box framework. A Cartesian coordinate system with origin located at an arbitrary distance, R₀, from the central star corotates with angular velocity Ω₀ = Ω_k(R₀), with Ω_k representing the Keplerian angular velocity. The x- and y-axes are oriented radially outward and along the orbital direction respectively, while the z-axis is directed in the vertical direction. We assume that the domain is small compared to the orbital distance so that Keplerian rotation appears as a linear shear flow with velocity ${\mathbf{U}}_{\mathbf{k}}=-\frac{3}{2}x{\mathrm{\Omega }}_0\hat{\mathbf{y}}$. We assumed the system is axisymmetric and neglect the vertical component of stellar gravity. In these limits, the continuity and momentum equations for the gas and dust components are, respectively, given by $$\frac{\mathrm{\partial }{\rho }_g}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\rho }_g\mathbf{u}\right)=0,$$(1) $$\begin{array}{rl}\frac{\mathrm{\partial }\mathbf{u}}{\mathrm{\partial }t}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u}=& 2u_y\mathrm{\Omega }\hat{\mathbf{x}}-u_x\frac{\mathrm{\Omega }}{2}\hat{\mathbf{y}}-\frac{1}{{\rho }_g}\mathrm{\nabla }P+2\eta R{\mathrm{\Omega }}^2\hat{\mathbf{x}}+\frac{ϵ}{{\tau }_s}(\mathbf{v}-\mathbf{u})\\ & +\frac{1}{{\rho }_g}\mathrm{\nabla }\cdot \mathbf{T}\end{array}$$(2) $$\frac{\mathrm{\partial }{\rho }_d}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\rho }_d\mathbf{v}\right)=\mathrm{\nabla }\cdot \left(D{\rho }_d\mathrm{\nabla }ϵ\right),$$(3) $$\frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }\mathbf{v}=2v_y\mathrm{\Omega }\hat{\mathbf{x}}-v_x\frac{\mathrm{\Omega }}{2}\hat{\mathbf{y}}-\frac{1}{{\tau }_s}(\mathbf{v}-\mathbf{u}),$$(4)

where ρ_g and ρ_d are the gas and dust densities, respectively, and u, v are their velocities measured in relation to the background Keplerian shear. We adopted an isothermal equation of state, $P={\rho }_g{c}_s^2$, where c_s = H_gΩ₀ is the (constant) sound speed, with H_g as the gas pressure scale height. In Eq. (2), T is the viscous stress tensor, given by $$\mathbf{T}={\rho }_g\nu (\mathrm{\nabla }\mathbf{u}+(\mathrm{\nabla }\mathbf{u}{)}^{†}-\frac{2}{3}I\mathrm{\nabla }\cdot \mathbf{u}),$$(5)

where ν is the kinematic viscosity employed to model the effects of possible underlying turbulence. Regarding the resulting turbulence-induced particle stirring, it is modeled by the source term appearing in the dust continuity equation, where D is the dust diffusion coefficient that is set here to D = ν. It is important to note that viscosity and diffusion are generally not included in our simulations, except in those presented in Appendix B.

In Eq. (2), the term ∝ η corresponds to an outward force acting on the gas due to a background radial pressure gradient, and determined by the dimensionless parameter, $$\eta =-\frac{1}{2{\rho }_g{\mathrm{\Omega }}_0^2{R}_0}\frac{\mathrm{\partial }P}{\mathrm{\partial }r}.$$(6)

Finally, τ_s is the stopping time that we will characterize in the following in terms of the dimensionless Stokes number, St = τ_s Ω. Dust grains are assumed to be in the Epstein regime such that St is the related to the grain size, a, via $$\mathrm{S}\mathrm{t}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_{i}a}{{\rho }_g{c}_s}\mathrm{\Omega },$$(7)

where ρ_i is the internal density of dust grains, and m =4πρ_ia³/3 their mass. In this work, we allowed m (and consequently a) to increase as a result of dust coagulation, which we modeled using a single-size approximation (Sato et al. 2016; Tominaga et al. 2021). This consists of following the size evolution of grains whose mass, m_p, dominate the dust density distribution at each location. The evolution equation for the so-called peak mass, m_p, can be obtained by considering the first-moment of the Smoluchowski equation (e.g., Appendix A of Tominaga et al. 2021), expressed as $$\frac{\mathrm{\partial }{m}_p}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }\left(m_p\right)={\mathcal{S}}^{+}\left({\mathrm{S}\mathrm{t}}_p\right)-{\mathcal{S}}^{-},$$(8)

where ${\mathcal{S}}^{+}=4\pi a_p^2\mathrm{\Delta }{v}_{\mathrm{p}\mathrm{p}}{\rho }_d$, with a_p the grain size corresponding to m_p and Δv_pp the typical relative velocity between particles, is the mass growth rate resulting from coagulation. This value can generally be expressed as a function of the Stokes number St_p corresponding to m_p. Although dissipation effects were not included in our dynamical equations, we followed Tominaga et al. (2021) and assumed that relative velocities are dominated by the turbulence experienced by the disk. This allowed us to i) make the problem analytically tractable and ii) identify clearly the role played by coagulation. In that case and within the limit St << 1, we can further make use of a simplified expression for Δv_pp given by (Ormel & Cuzzi 2007): $$\mathrm{\Delta }{v}_{\mathrm{p}\mathrm{p}}=\sqrt{C{\mathrm{S}\mathrm{t}}_p{\alpha }_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}}{c}_s.$$(9)

Here, C is a numerical factor set to C = 2.3 and α is the dimensionless turbulent stress parameter (Shakura & Sunyaev 1973), which we fixed to α_coag = 10⁻⁴. In the following and for convenience, we simply drop the subscript p in the definition of St_p and a_p so that (unless otherwise stated) we have $\mathrm{S}\mathrm{t}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_{i}a}{{\rho }_g{c}_s}\mathrm{\Omega }$ with ρ_i = 3m_p/4πa³.

With respect to the original work of Tominaga et al. (2021), we included in the mass evolution equation (Eq. (8)) an additional, constant sink term, 𝒮⁻, whose value was set to the initial value of 𝒮⁺. The primary aim of adding this term is to ensure that the initial chosen particle mass corresponds to an equilibrium solution to the set of Eqs. (1)–(4), combined with Eq. (8). Given the radial and vertical dependences of equilibrium quantities are ignored in the box, the left-hand side of Eq. (8) becomes zero at steady-state, thereby leading to 𝒮⁻ = 𝒮⁺.

We notice that in principle, dust loading leads to a reduction of the sound speed by increasing the inertia of the dusty fluid. This can be expressed by defining an effective sound speed given by (Lin & Youdin 2017): $${\tilde{c}}_s=\frac{c_s}{\sqrt{1+ϵ}}.$$(10)

From Eq. (9), we see a reduction in the sound speed due to dust-loading results in a decrease in the collision velocity. Recent works have shown that this effect may cause rapid dust growth in the strong clumping regime of the SI, due to the significant increase in the local dust-to-gas ratio there (Tominaga & Tanaka 2023, 2025; Carrera et al. 2025). Here, we do not account for this potential important effect and simply employ the nominal sound speed when computing the collision velocity.

2.2 Equilibrium state

For constant values of ρ_g, ρ_d, and St, equilibrium solutions to Eqs. (1)–(4) can be obtained, leading to velocity deviations from Keplerian rotation given by $$u_x=\frac{2ϵ\mathrm{S}\mathrm{t}}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(11) $$u_y=-\frac{(1+ϵ+{\mathrm{S}\mathrm{t}}^2)}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(12) $$v_x=\frac{-2\mathrm{S}\mathrm{t}}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(13) $$v_y=-\frac{(1+ϵ)}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(14)

where Δ² = St² + (1 + ϵ)². For a fixed grain size, these solutions correspond to the traditional steady-state drift solutions for a dusty disk obtained by Nakagawa et al. (1986). As mentioned above, the inclusion of coagulation would not be expected to impact these equilibrium solutions given the presence of the sink term in the right-hand side of Eq. (8).

3 Methodology

3.1 Numerical method and model setup

We examined the nonlinear dynamical evolution of the system by performing 2.5D (axisymmetric), unstratified simulations in a shearing-box using the multifluid version of the publicly available code FARGO3D (Benítez-Llambay & Masset 2016; Benítez-Llambay et al. 2019). We considered single dust species runs in which the effect of coagulation was implemented in the aforementioned single-size approximation, albeit using a conservative version of Eq. (8), expressed as $$\frac{\mathrm{\partial }{m}_p{\rho }_d}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(m_p{\rho }_d\mathbf{v}\right)=p_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathcal{S}}^{+}-{\mathcal{S}}^{-}){\rho }_d.$$(15)

Using this form, the previous equation can be integrated using the operator splitting technique employed by FARGO3D, with a coagulation term solved as a new substep within the source step. Tests of the implementation are presented in Appendix A. With respect to Eq. (8), the previous expression includes an additional, sticking efficiency, term p_eff whose expression is similar to that used by Tominaga et al. (2022a): $$p_{\mathrm{e}\mathrm{f}\mathrm{f}}=min(1,-\frac{\ln \left(\mathrm{S}\mathrm{t}/{\mathrm{S}\mathrm{t}}_f\right)}{\ln 5}),$$(16)

which also prevents the value of the Stokes number from increasing above St_f. This enables us to (i) take into account the effect of the fragmentation process in the simulations and (ii) examine the evolution outcome as a function of parameter St_f.

In our shearing-box setup, the Cartesian coordinate system is oriented with x along the (outward) radial direction and z pointing along the (upward) vertical direction. The domain considered is such that x ∈ [−0.0125, 0.0125] and z ∈ [−0.0025, 0.0025] or, equivalently, x ∈ (−0.25, 0.25)H₀ and z ∈ (−0.05, 0.05)H₀. The appropriate size of the shearing-box was determined from convergence tests; however, typically, a large radial extent of the domain is needed to accommodate the most unstable wavelength expected from linear growth maps of the SI as St → St_f, whereas the small vertical extent of the domain is consistent with an unstratified setup, focusing on the disk midplane. The domain is covered by N_x × N_z = 2048 × 512 grid cells, with the radial resolution chosen specifically to capture the largest number of CI unstable modes, as these have growth rates increasing with the radial wavenumber in the absence of diffusive processes. As we do not consider the vertical gravity in this work, boundaries are taken to be periodic in the z-direction, while we use standard shearing periodic boundaries in x.

3.2 Parameters and initial conditions

We adopted a unique value for the pressure gradient parameter, η = 0.005, corresponding to a typical pressure gradient in a protoplanetary disk. Three values for the initial dust-to-gas mass ratio, ϵ, were considered, ϵ = 0.02, 0.2, 3, which are representative of low (ϵ = 0.02), moderate (ϵ = 0.2), and high dust abundances (ϵ = 3). For simulations where dust coagulation was discarded, the Stokes number was kept constant with values: St = 0.01, 0.1, 1; whereas in those including coagulation, the Stokes number was computed from Eq. (15), starting from an initial value of St_i = 10⁻³. The details of parameters for each simulation are presented in Table 1.

Regarding the initial conditions, radial and azimuthal components of dust and gas velocities were initialized according to Eqs. (11)–(14), while setting the vertical component to zero. In most runs, the initial dust velocities were then seeded with random noise of amplitude 10⁻²c_s. Because the effects of stratification were not considered, we used an initial gas density of ρ₀ = 1. We also used a system of units where c_s = Ω_k = H_g = 1.

3.3 Diagnostics

One aim of this work is to assess the level of turbulence generated by the coagulation and streaming instabilities. To quantify the flux of angular momentum carried by turbulence, we first define the Shakura–Sunyaev stress parameter (Shakura & Sunyaev 1973): $$<{\alpha }_{\mathrm{S}\mathrm{S}}>=\frac{<{\rho }_g{u}_x{u}_y>}{{\rho }_0{c}_s^2},$$(17)

where the symbol <> denotes averaging in the x–z plane. To obtain a reliable measure of α_SS, we further performed a time average of the latter quantity, expressed as $$\overline{<{\alpha }_{\mathrm{S}\mathrm{S}}>}=\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}<{\alpha }_{\mathrm{S}\mathrm{S}}>dt,$$(18)

which we typically compute over a time interval of Δt = t₂ − t₁ ~ 10Ω⁻¹ once saturation has been attained. We also measured the velocity dispersion of the gas δu whose ith component as $$\delta u_i=\sqrt{<u_i^2>-<u_i{>}^2},$$(19)

along with the autocorrelation time, t_c,i, of the velocity fluctuations, obtained by evaluating each component of the autocorrelation function given by (Yang et al. 2018): $${\mathcal{R}}_i(t)=\int [u_i(\tau )-\overline{u_i}][u_i(t+\tau )-\overline{u_i}]d\tau .$$(20)

In the previous expression, $\overline{u}$ is the mean velocity, whose components u_i are estimated by considering their time average in the interval [t₁, t₂], $$\overline{u_i}=\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}{u}_i(\tau )d\tau .$$(21)

Following Hsu & Lin (2022), we set t₂ − t₁ = 5P and recorded the data with a cadence of 0.01P (Yang et al. 2018) to compute Eq. (20). We then determined t_c,i by fitting an exponential function to the autocorrelation function (Fromang & Papaloizou 2006).

From Eqs. (19) and (20), we can subsequently compute the gas bulk diffusion coefficient, D_g,i, in the ith direction, expressed as $$D_{g,i}={\overline{\delta u_i}}^2{t}_{c,i}.$$(22)

In the following, we made use of the dimensionless bulk gas diffusion coefficient instead, defined by $${\alpha }_{g,i}=\frac{D_{g,i}}{c_s{H}_g}={\left(\frac{\overline{\delta u_i}}{c_s}\right)}^2{\tau }_{c,i},$$(23)

where τ_c,i is the dimensionless equivalent of the correlation time, $${\tau }_{c,i}=t_{c,i}{\mathrm{\Omega }}_0^{-1}.$$(24)

4 Results

In this section, we discuss the impact of dust coagulation on the nonlinear evolution of the SI by providing a systematic comparison between pure SI simulations and runs that include dust coagulation. Our comparisons rely on models in which the maximum Stokes number St_f coincides with the value for St set in standard SI simulations. For simplicity, we chose to simply denote St_f as St when presenting the results of the simulations.

In Fig. 1, we present an overview of the results, where we show the time evolution of the maximum dust-to-gas ratio, ϵ_max, and Stokes number, St, for each model we considered. From this figure, it is immediately obvious that in dust rich disks, the impact of coagulation is relatively weak as simulations that include this effect reach overdensities that are very similar to those obtained in standard SI simulations.

In dust-poor disks with ϵ < 0.2, however, including coagulation has a stronger effect, with overdensities that can be a factor of 5–30 higher than without coagulation for tightly coupled particles. In this regime, we find that the CI (Tominaga et al. 2021) plays an important role in the evolution, as it is found to dominate over other processes for ϵ = 0.02, or even interact with the SI for ϵ = 0.2. Below, we describe these two different modes of evolution found for Stokes numbers St < 0.1 in more details. We then discuss the case of marginally coupled particles with St ≈ 1.

Fig. 1 Time evolution of the maximum dust-to-gas ratio, ϵmax, and Stokes number for all simulations listed in Table 1. The left, middle, and right columns correspond to ϵ = 0.02, 0.2, 3, respectively, while the upper, middle, and lower panels correspond to St = 0.01, 0.1, 1 respectively. In each panel, we show the evolution of ϵmax for the case with (blue) and without (green) coagulation. For the evolution of St, the light colour traces the minimum and maximum values, whereas the dark colour traces the box averaged value.

4.1 Tightly coupled particles

4.1.1 ϵ = 0.02: Evolution mediated by Cl

The CI is by essence a dust-driven instability which involves two ingredients: (i) the dust-density density dependence of the collisional coagulation rate and (ii) the dependence of St with the particle size. When these are mixed, a dust surface density perturbation δΣ_d leads to a radial variation of St, resulting in a traffic jam which tends to enhance δΣ_d. Compared to other dust-driven instabilities, the CI can develop over a few orbital periods even for small dust-to-gas ratios ϵ ~ 10⁻³. More precisely, the growth rate of the CI is typically given by (Tominaga et al. 2021): $$\sigma =\sqrt{\frac{1}{2}\frac{ϵ}{3t_0}k\mid <v_x>}\mid {\mathrm{\Omega }}^{-1},$$(25)

where k is the radial wavenumber and v_x is given by Eq. (13). Then, $t_0=4/3\sqrt{2{\mathrm{S}\mathrm{t}}_i/\alpha C\pi }$ is the mass growth timescale of a dust particle and can be obtained from Eq. (8). For St_i = 10⁻³ and ϵ = 0.02, setting kH = 4 × 10³, which corresponds to the maximum wavenumber that is resolved by ≈6 grid cells, we get σ ≈ 2.10⁻³Ω⁻¹. The time evolution of the maximum dust-to-gas ratio and Stokes number is shown, for ϵ = 0.02 and St = 0.01 (Model E0v02S0v01c), in the upper-left panel of Fig. 1. Both quantities exhibit an initial rapid exponential growh phase and then reach nonlinear saturation at t ~ 1500Ω⁻¹. We note in passing that in the context of the CI, unstable modes are essentially incompressible such that the evolution of ϵ_max reflects the evolution of the maximum of the perturbed dust density, δρ_d. From this plot, we can estimate a linear growth rate of σ ~ 1.7 × 10⁻³Ω⁻¹ during the initial linear growth phase, which is in good agreement with our previous analytical estimation. This clearly indicates that the CI is responsible for the initial increase in δρ_d and St that are observed, although disentangling the effect of the CI from the evolution of the averaged Stokes number is more difficult. For such small values of ϵ and St, we also note that the effect of the SI is expected to be negligible. This can be confirmed by inspecting the evolution for the case where coagulation is switched off (green line) and for which ϵ_max remains close to its initial value.

Figure 2 shows snapshots of the dust density and Stokes number at saturation, for this model and in the case with ϵ = 0.02 and St = 0.1 (Model E0v02S0v01c) for which a similar mode of evolution is found. Since the CI is originally a one-dimensional instability with no dependence in the vertical direction, it is not surprising to observe the formation of vertically extended filaments, whose number tends to decrease as St is increased. For St = 0.1, we indeed ultimately obtain a single filament with ϵ_max ≈ 0.1; whereas for St = 0.01, five filaments remain in the system at the end of the simulation with consequently smaller ϵ_max ≈ 0.006. We find that these filaments are themselves formed through the merging of thinner ones. This is illustrated in Fig. 3, where we display the vertically averaged dust density as a function of time and radial position. From this figure, it is clear that filaments grow not only via merging events but also by simply collecting dust as they drift inward. Both effects then combine to make the amplitude of the perturbed dust density start level off and saturate at later times. Another mechanism that may contribute to saturation is the level of turbulence generated by the CI itself. This is suggested by looking back to Fig. 2, where it is evident that turbulence operates in the disk, especially for St = 0.1. In this case, the time evolution of the dimensionless bulk gas diffusion coefficients, α_g,i, is shown in the upper panel of Fig. 4. Not surprisingly, this reveals anisotropy of the CI with an enhanced dust diffusion coefficient α_g,z of 𝒪(10⁻⁶) in the vertical direction, which is larger by two orders of magnitude compared to α_g,y. A similar trend is also found for St = 0.01 (see Table 1), with α_g,z of 𝒪(10⁻⁷) and α_g,y of 𝒪(10⁻⁸) in that case. The bottom panel of Fig. 4 compares linear growth rates for α_g,i = 0 (left) and for α_g,i = 3 × 10⁻⁶ (right), and which have been obtained from the linearized equations presented in Appendix A using ν = D = α_g,ic_sH_g. We see that the turbulent diffusion induced by the CI itself can indeed significantly cut off the most unstable modes corresponding to α_g,i = 0. Although it is not shown here, we found a similar effect for St = 0.01, which confirms that the CI-induced dust diffusion can also contribute to the saturation of the instability.

In summary, we find that for tightly coupled particles, the CI can operate in dust-poor disks and lead to the formation of dust rings with ϵ_max of 𝒪(0.1). However, this is not enough to subsequently trigger the SI within these rings because the dust backreaction onto the gas tends to significantly regulate the maximum dust concentration that can be achieved through the CI (Tominaga et al. 2022b). This is illustrated in Fig. 5, where we show the time evolution of the maximum dust-to-gas ratio and Stokes number for the same models that include dust coagulation, but in the case where the effect of the dust backreaction onto the gas is discarded. We find a significantly higher ϵ_max in that case, which is consistent with the results of Tominaga et al. (2022b, see their Fig. 3).

Fig. 2 Case of ϵ = 0.02, snapshots of the dust density (left) and Stokes number (right) for model E0v02S0v01c with St = 0.01 (left) and for model E0v02S0v1c with St = 0.1 (right). Both models include the effect of coagulation and result in the growth of the CI.

Fig. 3 Case of ϵ = 0.02, showing the time evolution of the vertically averaged dust density for model E0v02S0v01c with St = 0.01 (left), and model E0v02S0v1c with St = 0.1 (right).

Fig. 4 Case of ϵ = 0.2 and St = 0.01 (model E0v02S0v01c), showing the time evolution of the dimensionless diffusion coefficients, αg,i, in each direction (top). The lower panel shows the expected linear growth rates for the same model in the inviscid case (left) and for αg,i = 5 × 10−6 (right). These have been obtained from linear theory (see Appendix A) using D = ν = αg,icsHg. Here, ${\tilde{K}}_{x,z}=k_{x,z}\eta R_0$ are the dimensionless wavenumbers.

4.1.2 ϵ =0.2: Interplay between CI/SI and turbulence

In the case with ϵ = 0.2, we find that the difference in amplitude between the α_g,i coefficients is smaller. For simulations including coagulation, their time evolution is represented in Fig. 6. The underlying turbulence is more isotropic in that case, with α_g,y ~ α_g,z for both St = 0.01 (Model E0v2S0v01c) and St = 0.1 (Model E0v2S0v1c). This is also illustrated in Fig. 7 where we plot the kinetic energy spectrum at saturation for St = 0.1. The power spectrum is broadly consistent with a negative five third power law (dashed line), which suggests Kolmogorov-like turbulence. Isotropic turbulence is expected for classical SI-driven turbulence, although we found an early evolution driven by the CI. In Fig. 1, the top-center and central panels show the evolution for St = 0.01 and St = 0.1, respectively. For models with coagulation, the initial rise in ϵ_max (for t ≲ 100Ω⁻¹ typically) is caused by the CI, involving density growth of vertically extended filaments through dust capture and merging of these filaments into denser ones. Compared to the previous case with ϵ = 0.02, however, it appears that the value for ϵ_max is high enough to trigger the SI within the filaments, which leads to the subsequent increase in ϵ_max that is observed. This would be consistent with the isotropic diffusion mentioned above. Figure 8 shows the dust density at the onset of the secondary increase in ϵ_max. It is clear that nonlinear overdensities can grow in very localized regions within the filaments. This corresponds to unstable modes with vertical wavenumbers k_z ≠ 0, which favors an SI-like instability. To confirm that the SI is definitely at work inside the filaments, we conducted a test calculation by restarting the run for St = 0.1 at t = 200Ω⁻¹, but with the dust back-reaction onto the gas switched off. The results of this model indicate that the overdensities are quickly damped, which confirms the important role of the dust feedback for the formation of clumps in the filaments.

Interestingly, with respect to simulations without coagulation, we can notice that the overdensities that are produced by the SI in presence of coagulation are stronger. For St = 0.1 (middle-middle panel) and St = 0.01 (top-middle panel), models with coagulation yield ϵ_max ≈ 6 and ϵ_max ≈ 0.3–0.35, respectively, whereas the growth of the instability without coagulation is marginal. We argue that for this range of parameters, the SI may be assisted by dust coagulation. To support this argument, we present in Fig. 9 the growth rates of the SI (upper panel) and those obtained from our linear analysis with coagulation included (bottom panel), for St = 0.01. For the run with St = 0.1 (Model E0v2S0v1c), St = 0.01 corresponds roughly to the value for which overdensities start to grow inside the filaments (see middle-middle panel in Fig. 1). On this plot, the red-coloured rectangular region for ${\tilde{K}}_x\gtrsim 100$, with ${\tilde{K}}_x=k_x\eta R_0$ as the dimensionless wavenumber, corresponds to unstable CI modes as these modes have no vertical structure. The region enclosed between $10\lesssim {\tilde{K}}_x\lesssim 100$ and $10\lesssim {\tilde{K}}_z\lesssim 100$ could, rather, be identified as unstable SI modes. Based on a comparison with the upper panel, we see that these modes have larger growth rates when coagulation is taken into account, which suggests that coagulation can amplify unstable SI modes. From Fig. 8, we estimate the typical radial extent of a vertically extended filament to be 0.02H_g, which corresponds to ${\tilde{K}}_x\sim 15$. This is well within this region of coagulation-assisted SI, confirming that this mechanism can indeed be at work here.

For St = 0.1 (Model E0v2S0v1c), inspecting the two panels of Fig. 8 reveals that the vertical filaments broaden with time, which is simply a consequence of turbulent diffusion generated by the coagulation-assisted SI. As the filaments spead out, they can subsequently start to intersect each other, which leads to an isotropic turbulent configuration as suggested by the middle panel of Fig. 10. We note that a similar evolution is found for St = 0.01, for which we show in the upper panel of Fig. 10 a snapshot of the dust density when the instability has saturated.

Going back to the middle panel Fig. 10, it is interesting to note that the fully turbulent state consists of large voids and narrow particle streams, which is similar to the patterns developed by the SI for strongly coupled particles (Johansen & Youdin 2007). For standard SI, however, this regime is typically found for ϵ ≳ 1 whereas for ϵ = 0.2, the dust density is left almost unchanged (see bottom panel of Fig. 10). This implies that coagulation makes this mode of nonlinear saturation of the SI occur at much lower dust-to-gas ratios. This is not surprising as coagulation tends to increase radial drift speed by increasing the Stokes number, so that growing particles can fall more rapidly through the voids created by Poisson fluctuations (Johansen & Youdin 2007). We note, however, that the dust density enhancement is limited to ~30–40 with coagulation included; whereas for standard SI, this mode of nonlinear saturation generally corresponds to strong clumping with dust density enhancements ≳100.

Fig. 5 Case of ϵ = 0.2, showing the time evolution of the maximum dust-to-gas ratio, ϵmax, and Stokes number for St = 0.01 (top) and St = 0.1 (bottom), where the effect of the dust backreaction onto the gas is discarded.

Fig. 6 Case of ϵ = 0.2, showing the time evolution of the dimensionless diffusion coefficients, αg,i, in each direction for models E0v2S0v01c with St = 0.01 (top) and E0v2S0v1c with St = 0.1 (bottom).

Fig. 7 Spectrum of the vertically averaged gas kinetic energy density for ϵ = 0.2 and St = 0.1 (run E0v2S0v1c). Here, ${\tilde{K}}_x=k_x\eta R_0$ is the dimensionless radial wavenumber.

Fig. 8 For ϵ = 0.2 and St = 0.1 (run E0v2S0v1c, with coagulation), dust density snapshots at t = 120 Ω−1 (left) and t = 150 Ω−1 (right).

Fig. 9 Linear growth rates obtained from the linear analysis presented in Appendix A for ϵ = 0.2 and St = 0.01 for the case without coagulation (top) and with coagulation(bottom). Here, ${\tilde{K}}_{x,z}=k_{x,z}\eta R_0$ are the dimensionless wavenumbers.

Fig. 10 Dust density map during the saturation phase for model E0v2S0v01c (ϵ = 0.2, St = 0.01, with coagulation), model E0v2S0v1c (ϵ = 0.2, St = 0.1, with coagulation), model E0v2S0v1 (ϵ = 0.2, St = 0.1, without coagulation), shown from top to bottom.

4.2 Marginally coupled particles

For ϵ = 0.02, the evolution outcome of the model with final Stokes number St = 1 and including coagulation (run E0v02s1c, bottom left panel in Fig. 1) is found to be very similar to that of St = 0.1 particles (run E0v02s0v1c), involving the formation of a vertically extended filament inside which CI-induced turbulence operates.

For ϵ = 0.2, an interesting issue that emerges from the previous discussion (in Sect. 4.1.2) is whether the conditions needed to reach strong clumping are less stringent with coagulation. To investigate this further, we allowed, in the case where coagulation is included, the Stokes number to increase up to St = 1. For standard SI, a model with ϵ = 0.2 is expected to lie in the strong clumping regime for St ~ 1 and our main aim here is to examine if this threshold is diminished when coagulation is taken into account. In Fig. 1, models with ϵ = 0.2 and St = 1 correspond to the bottom middle panel where we compare the runs with and without coagulation. From this plot, it is immediately evident that: (i) the final outcome of the simulations does not depend on whether coagulation is included or not and (ii) including coagulation does not reduce the threshold in terms of Stokes number for the onset of strong clumping. Looking at the time evolution of ϵ_max for the model with coagulation, we see that the initial behaviour is consistent with the results presented in the previous section, involving increase in ϵ_max as coagulation-assisted SI grows and develops isotropic turbulence. At a later time, t ~ 200 Ω⁻¹, however, ϵ_max is observed to decrease before increasing again as St → 1. From the time evolution of St, we estimate the (mean) Stokes number associated with the observed drop of ϵ_max occuring at t ~ 200 Ω⁻¹ to be St ≈ 0.3. Interestingly, we find by running simulations of standard SI with different Stokes numbers that this value coincides approximately with the one above which significant growth of the SI occurs. To explore the origin of this decrease in ϵ_max in more details, we have conducted a number of test calculations, by restarting the simulation with coagulation (run E0v2S1c) at t ~ 200 Ω⁻¹, but with dust coagulation (feedback) switched on (off). The findings of this study can be summarised as follows:

• a restarted run with dust backreaction onto the gas switched off leads to a continuous increase in ϵ_max due to CI, which is the opposite behaviour to that found in model E0v2S1c. This demonstrates that dust feedback plays an important role in the decrease in ϵ_max seen in run E0v2S1c;

• a restarted run with dust coagulation switched off results in the development of the SI into stripes aligned in the vertical direction but tilted in the radial direction, similarly to the patterns formed by the SI in the case of large particles (Johansen & Youdin 2007). When applied to run E0v2S1c, this shows that from t ~ 200 Ω⁻¹ onward, the SI can indeed be triggered, consistenly with the results of our simulations of standard SI with various Stokes numbers. During the saturated state, we find ϵ_max ~ 0.5–1 for this test calculation, which is also in agreement with our standard SI calculation with St = 0.3;

• a restarted run with coagulation and dust backreaction switched off results in a maintained quasi-stationary turbulent state, with ϵ_max remaining close to its initial value.

Taken together, these results suggest that provided its growth rate is high enough, the SI can counteract the effect of the CI such that the final turbulent state is similar to that found in pure SI calculations. This is illustrated in Fig. 11, which shows the evolution of the dust density for the model with coagulation included (run E0v2S1c). The large stripes tilted toward the vertical direction that are expected to develop during nonlinear saturation of the SI for St = 1 particles can easily be recognized. For the particular model with ϵ = 0.2 and St = 1, examining Table 1 confirms that the values α_SS and α_g,i do not depend on whether or not coagulation is included.

To clearly demonstrate the finding that the SI can overcome the effects of the CI, we show in Fig. 12 linear growth rates of the instability for ϵ = 0.2 and two different values of St. Here, both the CI and SI operate, with the unstable CI (resp. SI) modes corresponding to the unstable modes that have little (resp. strong) dependence on ${\tilde{K}}_z$. We see that as coagulation proceeds and Stokes number transitions from St = 0.3 to St = 1, the maximum growth rate of the SI increases whereas unstable CI modes with radial wavenumbers in the range ${\tilde{K}}_x\in [5-100]$ are damped. This possibly arises because both CI and SI are driven by the radial drift of dust, which acts as a finite energy reservoir for the two instabilities such that any of these two instabilities can only grow at the expense of the other.

Fig. 11 Case of ϵ = 0.2 and St = 1 (model E0v2S1c, with coagulation), with dust density snapshots taken at different times.

Fig. 12 Case ϵ = 0.2, including the effects of coagulation, along with the linear growth rates obtained from the linear analysis presented in Appendix A for St = 0.3 (top) and St = 1 (bottom). Here, ${\tilde{K}}_{x,z}=k_{x,z}\eta R_0$ are the dimensionless wavenumbers.

5 Discussion

5.1 Impact of coagulation on the streaming instability

Although a model with ϵ = 0.2 and St = 1 leads to strong clumping, consistently with the results of Johansen & Youdin (2007), a necessary condition for particle clumping through the SI is generally approximated as ϵ ≳ 1 in the disk midplane, due to the strong increase in the SI linear growth rates at ϵ ≈ 1 for St ≲ 1. The right panels of Fig. 1 show the effect of coagulation on the SI for ϵ = 3. The SI linear growth rate for St = 0.01 is σ ~ Ω⁻¹ but this case only produced modest overdensities, whereas strong clumping occurred for St ≳ 0.1, as expected. Clearly, the effect of coagulation appears to be negligible for ϵ = 3, regardless the (final) value for St. Combined with our results for ϵ = 0.2 and St = 1, this suggests that for the optimally dust sizes and dust-to-gas ratios that maximize the SI growth rate, coagulation has an almost null effect. In that case, this means not only that i) coagulation has no impact on the SI linear growth rate, but also that ii) the dust concentration at saturation is not affected by coagulation. Physically, this arises simply when the coagulation timescale is longer that the SI growth timescale or, equivalently, when the particule coagulation growth rate is smaller than the SI growth rate. From Eq. (8), the inverse of the coagulation timescale is approximately given by $$t_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}^{-1}=\frac{1}{2}\sqrt{\frac{\pi \alpha C}{8\mathrm{S}\mathrm{t}}}ϵ.$$(26)

For St = 0.01 this leads to $t_{\text{coag\hspace{0.17em}}}^{-1}\sim 0.01{\mathrm{\Omega }}^{-1}$, which is indeed smaller that the aforementioned SI growth rate. Starting from an initial Stokes number St = 10⁻³, coagulation can of course promote SI by increasing the average dust size by one to two orders of magnitude, consistently with the results of Ho et al. (2024). However, we find that once the SI is triggered, its dynamics dominate over coagulation effects.

5.2 Effect of coagulation on turbulence and dust diffusion

Nevertheless, we have also shown that under certain conditions, coagulation and SI can act in concert. This is particularly true in dust-poor disks where the SI can grow in the vertically extended filaments formed by the CI, leading to a regime of coagulation-assisted SI that ultimately results in a level of turbulence much higher compared to standard SI (see Sect. 4.1.2). For each model, the α_SS dimensionless parameter that quantities the radial flux of gas orbital momentum is listed in the first column of Table 1. Most runs have α_SS < 0, in agreement with previous estimations of the flux of angular momentum in SI-driven turbulent disks (Johansen & Youdin 2007; Hsu & Lin 2022). The fact that α_SS is more negative when coagulation is taken into account suggests that the radial transport of angular momentum induced by the CI is also inward. This is not surprising since, similarly to the SI (Johansen & Youdin 2007), CI is also powered by a global pressure gradient.

An interesting issue that needs to be investigated is how the collision velocity resulting from the turbulence compares with the fragmentation velocity, v_f. It can indeed not be excluded that the turbulence driven by the coagulation and streaming instabilities in turn hinders further coagulation by producing collision velocities that exceed the fragmentation threshold. To address this question, we can estimate the maximum Stokes number, St_f, achievable through coagulation. From Eq. (9), it is straightforward to show that St_f is given by $${\mathrm{S}\mathrm{t}}_f=\frac{v_f^2}{C{\alpha }_{\mathrm{S}\mathrm{S}}{c}_s^2}.$$(27)

The fragmentation velocity has been studied both numerically and experimentally (Blum & Wurm 2008; Zsom et al. 2010; Wada et al. 2013; Gundlach & Blum 2015; Ueda et al. 2019), and it has been estimated to be between 1 and 10 m s⁻¹ for silicate particles and between 10 and 80 m s⁻¹ for icy aggregates. Regardless the nature of dust aggregates, all but one simulation resulted in turbulent collision velocities smaller than v_f for typical disk parameters. Assuming that R₀ corresponds to 1 AU, Model E0v2S0v1c which enters the regime of coagulation-assisted SI is the only exception, as we find St_f ≈ 0.07 for v_f = 1 m s⁻¹. By preventing the formation of large dust grains, this implies that coagulation-driven turbulence may in that case inhibit the strong clumping phase of the SI.

A related effect of coagulation is the vertical diffusion of particles driven by the CI and/or SI. At a steady state (i.e., once turbulent diffusion counterbalances gravitational settling), it is expected that the particle profile is characterized by a dust scale height, H_d, given by (Dubrulle et al. 1995; Youdin & Lithwick 2007): $$\frac{H_d}{H_g}=\sqrt{\frac{{\alpha }_{g,z}}{{\alpha }_{g,z}+\mathrm{S}\mathrm{t}}}.$$(28)

Using the values for α_g,z reported in Table 1, we find that it is possible to differentiate between CI and SI in dust-poor disks with ϵ ≲ 0.2 and for tightly coupled particles with St ≲ 0.1, as the SI is ineffective at producing a significant level of turbulence under these conditions. In contrast, the CI on its own can generate a vertical bulk diffusion coefficient α_gz ~ 5 × 10⁻⁷ for ϵ = 0.02, which results in H_d ~ 0.006H_g. In the coagulation-assisted SI regime (run E0v2S0v1c), α_g,z can even reach 𝒪(10⁻⁴) for St = 0.1, corresponding to H_d ~ 0.01H_g.

5.3 Caveats

The results that we presented in this work may suffer from a number of assumptions and limitations, the main one being that we adopted a moment equation to describe dust coagulation, based on the vertically integrated Smoluchowski equation (Sato et al. 2016; Tominaga et al. 2021). Rather than focusing on the evolution of the full dust size distribution, this approach consists in following the evolution of the so-called peak mass, which represents the weighted average mass. Sato et al. (2016) have shown that provided the effect of finite size dispersion is taken into account when evaluating the collision velocities, good agreement is obtained between single-size and full size calculations of drifting and growing dust grains in protoplanetary disks. We caution, however, that the original model of Sato et al. (2016) does not take into account the effect of fragmentation and it is not clear whether or not the single-size approximation would still hold when fragmentation is considered, since this process would tend to broaden the initial dust size distribution. Relaxing the single-size approximation and considering a full, possibly wide, dust size distribution produced by fragmentation may have important consequences on our results. In the context of the SI, recent studies have shown that taking into account the effect of a dust size distribution can significantly modify the growth rate of (monodisperse) SI (Krapp et al. 2019; Paardekooper et al. 2020; Zhu & Yang 2021; McNally et al. 2021). Therefore, our model does not have the capacity to explore the effect of coagulation on this “polydisperse” SI and this should be examined in a future study.

We also considered inviscid disks, while various hydrodynamical instabilities can be triggered in the planet-forming regions of protoplanetary disks. This gives rise to turbulence, with a corresponding turbulent viscosity, typically leading to α_SS ~ 10⁻⁴. These levels of external turbulence may significantly impact the nonlinear evolution of both the SI (Chen & Lin 2020), and CI (Tominaga et al. 2021). Although models that include dissipation effects should be explored in more detail, we present in Appendix B the results of a few viscous simulations that also take dust diffusion into account. This suggests that the results presented in this paper could be affected by the external turbulence for α_SS ≳ 10⁻⁶.

Another limitation lies in the size of our vertical domain that is (based on the discussion in Sect. 5.2) systematically higher than the dust layer thickness. The consequence of this is that the effects of stratification might be important and these should be examined in more details in a future study. Another reason for considering vertically stratified disks is that the CI is characterised by vertically extended filaments and it is not clear whether or not these would persist in stratified models.

We note that these two aspects, namely, considering a full dust size distribution and vertical stratification, were recently studied by Ho et al. (2024). Although the potential important effect of the CI was not investigated in this work, they found that the dust size distribution can be significantly altered by coagulation within the dense dust clumps that are formed as the dust settles toward the midplane and the SI is triggered.

6 Summary

In this paper, we present the results of unstratified, axisymmetric, shearing-box simulations of SI including the effect of dust coagulation. We solved a simplified version of the coagulation equation by a moment method, which consists of following the size evolution of solids that dominate the total dust surface density. Our primary aim was to examine to what extent dust growth impacts the linear growth of the SI and to get insight into its impact on the nonlinear saturation phase of the SI. We are particularly interested in the role played by the CI (Tominaga et al. 2021), which results from the dust-density dependence of coagulation efficiency and the grain size dependence of radial drift speed. Our motivation stems from the fact that the CI can be triggered for dust-to-gas ratios as small as ϵ ~ 10⁻³ and make the Stokes number of solids rapidly reach unity, leading thereby to optimal conditions for efficient growth of the SI. We studied how the CI and the SI can interact with each other, as well as how the interplay between these two instabilities depends on the initial dust-to-gas ratio and on the maximum Stokes number St allowed before fragmentation sets in. Our key findings are summarized as follows:

1. In dust-poor disks with ϵ ~ 10⁻² and for tightly coupled particles with St ≲ 10⁻¹, we find that evolution is mainly driven by the CI. Since the latter instability is by essence an axisymmetric instability, early evolution involves the formation of vertically extended filaments that tend to merge themselves more rapidly for larger Stokes numbers. Within the filaments, the nonlinear saturation of the CI gives rise to anisotropic turbulence with an enhanced dust diffusion coefficient in the vertical direction;

2. In dust-rich disks with ϵ > 1, evolution is mainly driven by the SI simply because in that case, the growth timescale of the SI is smaller than the coagulation timescale. In this regime, coagulation is observed to have only little effect, not only on the initial linear growth rate of the instability, but also on the final dust concentration values reached once the SI has saturated;

3. For intermediate dust-to-gas ratios of ϵ ~ 0.1 and Stokes numbers St ≲ 0.1, we observe the onset of the SI within the filaments formed originally by the CI. An examination of linear growth rate maps indicates that for this range of parameters, the linear growth rate of the SI is significantly increased by coagulation effects, which we refer to as a regime of coagulation-assisted SI. In this regime, nonlinear evolution leads to isotropic turbulence with dust diffusion coefficients having similar value in each direction. The generated turbulence causes the vertically extended filaments to spread out until they start to intersect each other, ultimately leading to a turbulent state characterised by large voids and narrow particle streams. The patterns that develop look similar to those formed by the SI in the strong clumping regime, albeit with a smaller dust concentration, on the order of 𝒪(30–40), compared to that typical of the strong clumping regime;

4. However, as the Stokes number increases and approaches unity, SI can overcome the effect of the CI. Thus, the final turbulent state ends up being very similar to the results of pure SI calculations without the effect of coagulation included.

Taken together, these findings suggest that the formation of dense clumps through the SI may be rendered easier by coagulation, simply because the process of dust growth enables the SI to reside in a regime where it can operate efficiently. The CI on its own only has a slight impact on the onset of the SI, but it can nonetheless be a non negligible source of turbulence in protoplanetary disks, especially when it acts in concert with the SI.

Acknowledgements

Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Universite de Bordeaux and by HPC resources of Cines under the allocation A0170406957 made by GENCI (Grand Equipement National de Calcul Intensif). MKL is supported by the National Science and Technology Council (grants 113-2112-M-001-036-, 114-2112-M-001-018-, 113-2124-M-002-003-, 114-2124-M-002-003-) and an Academia Sinica Career Development Award (AS-CDA-110-M06). This work was partly supported by the “Action Thématique de Physique Stellaire” (ATPS) of CNRS/INSU PN Astro co-funded by CEA and CNES, and by ORCHID program (project number: 49523 PG, co-PI. Dutrey-Tang), funded by the French Ministry for Europe and Foreign Affairs, the French Ministry for Higher Education and Research and the National Science and Technology Council (NSTC).

References

Appendix A Linearized equations

To linearize Eqs. 1–4, we write X as a quantity of the system as the sum between the equilibrium state X₀ (Eqs. 11–14), and a perturbation δX << X₀ in Fourier modes with k_x and k_z the spatial modes and σ the growth rate of the instability. $$X(x,z,t)=X_0(x,z)+\delta X\times e^{i(k_{x}x+k_{z}z)+\sigma t}.$$(A.1)

By injecting this into the mass continuity equations, Eqs. 1 and 3, we obtain $$\sigma \delta {\rho }_g=-{\rho }_{g,0}(i{k}_x\delta u_x+i{k}_z\delta u_z)-i{k}_x{u}_{x,0}\delta {\rho }_g,$$(A.2a) $$\begin{array}{rl}\sigma \delta {\rho }_d=& -{\rho }_{d,0}(i{k}_x\delta v_x+i{k}_z\delta v_z)-i{k}_x{v}_{x,0}\delta {\rho }_d\\ & -{\rho }_{d,0}D(k_x^2+k_z^2)(\delta {\rho }_d/{\rho }_{d,0}-\delta {\rho }_g/{\rho }_{g,0}).\end{array}$$(A.2b)

Similarly, with the equations of conservation of momentum for the gas Eq. 2 $$\begin{array}{rl}\sigma \delta u_x=& -u_{x,0}i{k}_x\delta u_x+2\mathrm{\Omega }\delta u_y-c_s^{2}i{k}_x\frac{\delta {\rho }_g}{{\rho }_{g,0}}+\delta F_x^{visc}\\ & -C_{DF}{m}_p^{-1/3}\{[\delta {\rho }_d-\frac{1}{3m_p}{\rho }_d\delta m_p](u_{x,0}-v_{x,0})\\ & +{\rho }_d(\delta u_x-\delta v_x)\},\end{array}$$(A.3a) $$\begin{array}{rl}\sigma \delta u_y=& -u_{x,0}i{k}_x\delta u_y-\frac{1}{2}\mathrm{\Omega }\delta u_x+\delta F_y^{visc}\\ & -C_{DF}{m}_p^{-1/3}\{[\delta {\rho }_d-\frac{1}{3m_p}{\rho }_d\delta m_p](u_{y,0}-v_{y,0})\\ & +{\rho }_d(\delta u_y-\delta v_y)\},\end{array}$$(A.3b) $$\begin{array}{rl}\sigma \delta u_z=& -u_{x,0}i{k}_x\delta u_z-c_s^{2}i{k}_z\frac{\delta {\rho }_g}{{\rho }_{g,0}}+\delta F_z^{visc}\\ & -C_{DF}{m}_p^{-1/3}{\rho }_d(\delta u_z-\delta v_z),\end{array}$$(A.3c)

with $C_{DF}=\mathrm{\Omega }\sqrt{\frac{8}{\pi }}\frac{H_g}{{\rho }_i}{\left(\frac{4\pi {\rho }_i}{3}\right)}^{1/3}$ a coefficient related to drag force and $\delta F_{(x,y,z)}^{visc}$ the linearized viscous forces from Chen & Lin (2020): $$\delta F_x^{visc}=-\nu (k_z^2+\frac{4}{3}{k}_x^2)\delta v_x-\frac{1}{3}v{k}_x{k}_z\delta v_z,$$(A.4a) $$\delta F_y^{visc}=-\nu (k_z^2+k_x^2)\delta v_y,$$(A.4b) $$\delta F_z^{visc}=-\nu (\frac{4}{3}{k}_z^2+k_x^2)\delta v_z-\frac{1}{3}\nu k_x{k}_z\delta v_x,$$(A.4c)

Fig. A.1 Perturbation $\frac{max(|{\rho }_d-{\rho }_{d,0}|)}{{\rho }_{d,0}}$ and $\frac{max(|m_p-m_{p,0}|)}{m_{p,0}}$ obtained through hydrodynamic simulations with coagulation, in blue and orange. In green, the theoretical growth rate obtained analytically.

and for the equations of conservation of momentum for the dust Eq. 4, $$\begin{array}{rl}\sigma \delta v_x=& -v_{x,0}i{k}_x\delta v_x+2\mathrm{\Omega }\delta v_y\\ & -C_{DF}{m}_p^{-1/3}\{[\delta {\rho }_g-\frac{1}{3m_p}{\rho }_g\delta m_p](v_{x,0}-u_{x,0})\\ & +{\rho }_g(\delta v_x-\delta u_x)\},\end{array}$$(A.5a) $$\begin{array}{rl}\sigma \delta v_y=& -v_{x,0}i{k}_x\delta v_y-\frac{1}{2}\mathrm{\Omega }\delta v_x\\ & -C_{DF}{m}_p^{-1/3}\{[\delta {\rho }_g-\frac{1}{3m_p}{\rho }_g\delta m_p](v_{y,0}-u_{y,0})\\ & +{\rho }_g(\delta v_y-\delta u_y)\},\end{array}$$(A.5b) $$\sigma \delta v_z=-v_{x,0}i{k}_x\delta v_z-C_{DF}{m}_p^{-1/3}{\rho }_g(\delta v_z-\delta u_z).$$(A.5c)

And finally, the perturbed equation for the coagulation equation, expressed in Eq. 8, $$\sigma \delta m_p=-v_{x,0}i{k}_x\delta m_p+D_{coag}{m}_p^{5/6}{\rho }_g^{-1/2}{\rho }_d[\frac{5}{6}\frac{\delta m_p}{m_p}-\frac{1}{2}\frac{\delta {\rho }_g}{{\rho }_g}+\frac{\delta {\rho }_d}{{\rho }_d}],$$(A.6)

with $D_{coag}={\left(\frac{\pi }{2}\right)}^{5/12}\times 3^{5/6}{\rho }_i^{-1/3}\sqrt{2.3\alpha }{H}_g^{-1/2}{c}_s$.

This system of nine equations is rewritten in its matrix form as sV = AV, using the dimensionless quantities $\hat{v}=\frac{v}{\eta R\mathrm{\Omega }},K=k\eta R$ and $s=\frac{\sigma }{\mathrm{\Omega }}$, with the perturbation vector $\mathbf{V}={[\frac{\delta {\rho }_g}{{\rho }_g},\frac{\delta {\rho }_d}{{\rho }_d},\delta {\hat{u}}_x,\delta {\hat{u}}_y,\delta {\hat{u}}_z,\delta {\hat{v}}_x,\delta {\hat{v}}_y,\delta {\hat{v}}_z,\frac{\delta m_p}{m_p}]}^{\mathrm{\top }}$ and the matrix of the linearized system A.

For a given spatial mode (K_x, K_z), the largest eigenvalue of A corresponds to the growth rate σ of the most unstable mode, while the associated eigenvector can be used to study the linear case of this mode numerically. This enables a verification of the coagulation implementation, as well as the analytical development growth rates derived from the matrix’s eigenvalues via a comparison of the two, as done in Fig. A.1 in the linA case of St = 0.1, ϵ = 3 and (K_x, K_z) = (30, 30) with coagulation. The analytical solution gives s = 0.3602, identical to the numerical growth rate derived from the particle density and particle mass until the perturbation is no longer negligible (δX ~ 0.1).

Appendix B Effect of gas viscosity and dust diffusion

In the inviscid limit that we considered in this work, the numerical dissipation may impact our results, such that these could be dependent on the employed resolution. This is because the growth rate of CI increases with radial wavenumber, such that the maximum wavenumber that is resolved by the grid is expected to be unstable. Here, we present the results of additional simulations with nonzero gas viscosity, parameterized as ν = αc_SH_g, and with the dust diffusion coefficient set to the same value, namely, D = ν. As in inviscid simulations presented above, the α_coag parameter introduced into Eq. 9 is set to α_coag = 10⁻⁴. The aim here is to examine: i) how the results are changed when dissipation processes are controlled and ii) how these would be modified in presence of external turbulence, namely when the turbulence originates from other sources of instabilities than the CI and/or SI individually. For different values of the standard α parameter, the time evolution of ϵ_max and st for these runs is presented in Fig. B.1. Not surprisingly, the linear growth rate of the instabilities is reduced as α is increased. This is a consequence of the stabilization of the CI by dust diffusion (Tominaga et al. 2021). Unstable high wavenumbers modes are also suppressed, as illustrated by Fig. B.2, which shows the dust density at the end of the linear growth phase. From Fig. B.1, we estimate the value for the numerical viscosity to correspond to α ~ 10⁻⁹–10⁻⁷. This is in line with the expectation that unstable modes with wavenumbers such that $kH\sim 1/\sqrt{\alpha }$ are damped by turbulence (Dubrulle et al. 1995). Using α = 10⁻⁷, this gives kH ≈ 3 × 10³, which is close to the value mentioned in Sect. 4.1.1. Going back to Fig. B.1, we can observe growth of the SI for each value of α we considered, except for α = 10⁻⁶. In that case, however, St ≲ 0.5 at the end of the simulation due to slower coagulation, leading to a reduced SI growth rate. As St → 1, the SI should nevertheless enter a strong clumping regime, as it is expected in that case the SI growth rate to be only impacted by diffusion for α = 10⁻⁶ (e.g., see the upper-left panel of Fig. 3. in Chen & Lin 2020).

Fig. B.1 For ϵ = 0.2, St = 1 and with coagulation included, time evolution of the maximum dust-to-gas ratio ϵmax and Stokes number St in viscous simulations where the α viscous stress parameter is taken to be equal to the dimensionless dust diffusion coefficient.

Fig. B.2 For ϵ = 0.2, St = 1, dust density snapshots for the reference inviscid run (top-left) and for viscous simulations with α = 10−9 (top-right), α = 10−7 (bottom-left), and α = 10−6 (bottom-right).

All Tables

All Figures

Fig. 1 Time evolution of the maximum dust-to-gas ratio, ϵmax, and Stokes number for all simulations listed in Table 1. The left, middle, and right columns correspond to ϵ = 0.02, 0.2, 3, respectively, while the upper, middle, and lower panels correspond to St = 0.01, 0.1, 1 respectively. In each panel, we show the evolution of ϵmax for the case with (blue) and without (green) coagulation. For the evolution of St, the light colour traces the minimum and maximum values, whereas the dark colour traces the box averaged value.

In the text

	Fig. 2 Case of ϵ = 0.02, snapshots of the dust density (left) and Stokes number (right) for model E0v02S0v01c with St = 0.01 (left) and for model E0v02S0v1c with St = 0.1 (right). Both models include the effect of coagulation and result in the growth of the CI.
In the text

	Fig. 3 Case of ϵ = 0.02, showing the time evolution of the vertically averaged dust density for model E0v02S0v01c with St = 0.01 (left), and model E0v02S0v1c with St = 0.1 (right).
In the text

Fig. 4 Case of ϵ = 0.2 and St = 0.01 (model E0v02S0v01c), showing the time evolution of the dimensionless diffusion coefficients, αg,i, in each direction (top). The lower panel shows the expected linear growth rates for the same model in the inviscid case (left) and for αg,i = 5 × 10−6 (right). These have been obtained from linear theory (see Appendix A) using D = ν = αg,icsHg. Here, ${\tilde{K}}_{x,z}=k_{x,z}\eta R_0$ are the dimensionless wavenumbers.

In the text

	Fig. 5 Case of ϵ = 0.2, showing the time evolution of the maximum dust-to-gas ratio, ϵmax, and Stokes number for St = 0.01 (top) and St = 0.1 (bottom), where the effect of the dust backreaction onto the gas is discarded.
In the text

	Fig. 6 Case of ϵ = 0.2, showing the time evolution of the dimensionless diffusion coefficients, αg,i, in each direction for models E0v2S0v01c with St = 0.01 (top) and E0v2S0v1c with St = 0.1 (bottom).
In the text

	Fig. 7 Spectrum of the vertically averaged gas kinetic energy density for ϵ = 0.2 and St = 0.1 (run E0v2S0v1c). Here, ${\tilde{K}}_x=k_x\eta R_0$ is the dimensionless radial wavenumber.
In the text

	Fig. 8 For ϵ = 0.2 and St = 0.1 (run E0v2S0v1c, with coagulation), dust density snapshots at t = 120 Ω−1 (left) and t = 150 Ω−1 (right).
In the text

	Fig. 9 Linear growth rates obtained from the linear analysis presented in Appendix A for ϵ = 0.2 and St = 0.01 for the case without coagulation (top) and with coagulation(bottom). Here, ${\tilde{K}}_{x,z}=k_{x,z}\eta R_0$ are the dimensionless wavenumbers.
In the text

	Fig. 10 Dust density map during the saturation phase for model E0v2S0v01c (ϵ = 0.2, St = 0.01, with coagulation), model E0v2S0v1c (ϵ = 0.2, St = 0.1, with coagulation), model E0v2S0v1 (ϵ = 0.2, St = 0.1, without coagulation), shown from top to bottom.
In the text

	Fig. 11 Case of ϵ = 0.2 and St = 1 (model E0v2S1c, with coagulation), with dust density snapshots taken at different times.
In the text

	Fig. 12 Case ϵ = 0.2, including the effects of coagulation, along with the linear growth rates obtained from the linear analysis presented in Appendix A for St = 0.3 (top) and St = 1 (bottom). Here, ${\tilde{K}}_{x,z}=k_{x,z}\eta R_0$ are the dimensionless wavenumbers.
In the text

	Fig. A.1 Perturbation $\frac{max(|{\rho }_d-{\rho }_{d,0}|)}{{\rho }_{d,0}}$ and $\frac{max(|m_p-m_{p,0}|)}{m_{p,0}}$ obtained through hydrodynamic simulations with coagulation, in blue and orange. In green, the theoretical growth rate obtained analytically.
In the text

	Fig. B.1 For ϵ = 0.2, St = 1 and with coagulation included, time evolution of the maximum dust-to-gas ratio ϵmax and Stokes number St in viscous simulations where the α viscous stress parameter is taken to be equal to the dimensionless dust diffusion coefficient.
In the text

	Fig. B.2 For ϵ = 0.2, St = 1, dust density snapshots for the reference inviscid run (top-left) and for viscous simulations with α = 10−9 (top-right), α = 10−7 (bottom-left), and α = 10−6 (bottom-right).
In the text