© The Authors 2024

1 Introduction

The early stages of planet formation are not well understood. Observations of protoplanetary and debris discs, along with the discovery of numerous exoplanets, have shed light on the efficiency of mechanisms that allow tiny dust particles to aggregate into objects several thousand kilometres in size (Dominik & Tielens 1997; Dullemond & Dominik 2005). The process of dust particle growth involves several stages. Initially, monomers grow thanks to Brownian motions, sticking together via Van der Waals forces (Cuzzi et al. 1993; Stepinski & Valageas 1997). However, after growing by just one order of magnitude, this mechanism becomes ineffective, giving way to turbulence within the gas, which promotes collisions between aggregates (Weidenschilling & Cuzzi 1993). Thanks to this turbulence, grains can grow to sizes ranging from millimetres to metres, depending on their porosity (Okuzumi et al. 2012; Garcia & Gonzalez 2020). At this point, they face various barriers to further growth, including radial drift and fragmentation.

Radial drift is caused by the friction between gas and dust (Whipple 1972; Weidenschilling 1977) resulting from the difference in orbital velocities due to gas pressure support. This friction causes dust particles to lose angular momentum and drift inward towards the star (Nakagawa et al. 1986). This phenomenon is especially rapid for intermediate-sized grains, posing a challenge to their growth (Weidenschilling 1977). In addition to radial drift, other barriers to dust growth are associated with the physics of grain interactions (Weidenschilling & Cuzzi 1993; Stepinski & Valageas 1997). These interactions depend on the relative velocities of grains during collisions, which can lead to coagulation, bouncing (elastic or plastic deformation), or fragmentation (Blum & Wurm 2000). The bouncing or fragmentation thresholds refer to the relative velocity at which grains begin to bounce or fragment, respectively. In any case, these thresholds, ranging from a few metres per second to a few tens of metres per second, are quickly reached when sizes are on the order of millimetres to metres. These barriers are called the bouncing barrier (Zsom et al. 2010; Windmark et al. 2012) and the fragmentation barrier (Weidenschilling & Cuzzi 1993; Dominik & Tielens 1997; Blum & Wurm 2008). Blum & Wurm (2000, 2008) and Güttler et al. (2010) showed in their experiments that adhesion and fragmentation are linked to the composition of the grains and the volatile materials covering them. Recent studies (Yamamoto et al. 2014; Kimura et al. 2020; San Sebastián et al. 2020) have shown that silicates are more resistant than suspected.

More recently, additional barriers to dust growth have been identified, including collisional erosion (Schräpler & Blum 2011) and aeolian erosion (Paraskov et al. 2006; Rozner et al. 2020; Grishin et al. 2020; Michoulier et al. 2024). Collisional erosion occurs when larger aggregates eject grains from their surface through successive collisions with smaller grains. Aeolian erosion is caused by gas friction, leading to the detachment of loosely bound grains from larger aggregates in the innermost region of the disc. Another mechanism that can destroy grains is rotation disruption (Tatsuuma & Kataoka 2021), where the gas flow torques the grains, causing them to rotate. When the centrifugal force exceeds the maximum tensile strength of the grains, they shatter into fragments. Tatsuuma & Kataoka (2021) and Michoulier & Gonzalez (2022a) studied this process for porous aggregates in protoplanetary discs.

All of these barriers theoretically hinder dust growth. In order to understand how the barriers are bypassed in nature to form planetesimals from micron-sized grains, various solutions have been proposed. Some of them involve creating dust traps with local pressure maxima, such as vortices (Barge & Sommeria 1995; Meheut et al. 2012); baroclinic instability (Klahr & Bodenheimer 2003; Loren-Aguilar & Bate 2015); planet gaps (Paardekooper & Mellema 2004; Fouchet et al. 2007, 2010; Gonzalez et al. 2012; Zhu et al. 2014); or snow lines where different chemical species sublimate (Kretke & Lin 2007; Brauer et al. 2008; Drążkowska et al. 2014; Drążkowska & Alibert 2017; Hyodo et al. 2019; Vericel & Gonzalez 2020). Self-induced dust traps, driven by dust feedback, have also been proposed by Gonzalez et al. (2017). Other mechanisms enable direct formation of planetesimals, such as the streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007; Youdin & Johansen 2007; Jacquet et al. 2011; Carrera et al. 2015; Yang et al. 2017; Schäfer et al. 2017; Auffinger & Laibe 2018; Li et al. 2019; Zhu & Yang 2021; Schaffer et al. 2021). The streaming instability arises from the weak coupling of large dust grains and gas in dust enriched regions, and it concentrates dust in filamentary structures, which can ultimately collapse gravitationally.

One solution that has recently gained attention is grain porosity (Ormel et al. 2007; Suyama et al. 2008; Okuzumi et al. 2009, 2012; Suyama et al. 2012; Kataoka et al. 2013; Garcia 2018; Garcia & Gonzalez 2020). Porosity, often overlooked for simplicity in dust modelling (Weidenschilling 1977; Nakagawa et al. 1986; Barrière-Fouchet et al. 2005; Laibe et al. 2008; Drążkowska et al. 2014; Gonzalez et al. 2015; Vericel et al. 2021), changes gas-dust coupling and grain evolution. Porous grains, due to their larger collision cross-sections for a given mass, grow faster and decouple from the gas at larger sizes, enhancing their survival in the disc. This has been observed in studies using local or 1D disc models to track the evolution of mass and filling factor for single grains (Ormel et al. 2007; Suyama et al. 2008, 2012; Okuzumi et al. 2009, 2012; Kataoka et al. 2013) and more recently in a 3D disc model evolving an entire population of porous grains (Garcia 2018; Garcia & Gonzalez 2020). They are also less susceptible to fragmentation, which contributes to planetesimal formation through coagulation (Garcia 2018). Laboratory experiments studying the growth of highly porous aggregates are restricted to low velocity (Blum & Schräpler 2004), thus making numerical experiments necessary (Blum & Wurm 2000; Dominik & Tielens 1997; Wada et al. 2007, 2009; Suyama et al. 2008, 2012; Seizinger et al. 2012).

By looking at dust in comets, Güttler et al. (2019) classified aggregates into three classes: fractal, porous, and compact. They also showed that porous grains are common, supporting the need to account for porosity in models. Observations of discs have shown that porosity must be considered in order to explain spectral energy distributions (SEDs, Kataoka et al. 2016; Zhang et al. 2023), polarization (Kataoka et al. 2015, 2016, 2019; Tazaki et al. 2019), and the amount of dust settling (Pinte et al. 2019; Verrios et al. 2022). The monomer size is key to explaining polarized light (Tazaki & Dominik 2022), and Tazaki et al. (2023) identify smaller grains on disc surfaces as fractal or porous aggregates.

In this paper, we focus on the impact of porosity on dust evolution by modelling mechanisms impairing the growth of grains in the frame of the mono-disperse approximation. Sections 2 and 3 present our model and numerical simulations, respectively. In Sect. 4, we show the impact of porosity on dust evolution and study the role of compaction during fragmentation using 1D and 3D simulations. We discuss our results and implications as well as the limitations of our models in Sect. 5, and we conclude in Sect. 6.

2 Dust evolution model

2.1 Grain growth

Dust coagulation results from collisions between grains with a certain relative velocity. Different velocities or mass ratios can result in sticking, mass transfer or penetration (Güttler et al. 2010). To treat dust growth, Stepinski & Valageas (1997) assumed that locally, the mass distribution of grains is peaked around a single value, or ‘mono-disperse’, where collisions occur between grains of identical mass. In the following, we adopted their formalism, following the implementations by Laibe et al. (2008) or, for the PHANTOM code, Vericel et al. (2021). The reader is referred to these two papers for discussion of its implications. In each collision, the mass m of a dust grain of size s doubles over a characteristic time τ_coll

$${\left(\frac{\text{d}m}{\text{d}t}\right)}_{\text{grow }}\approx \frac{m}{{\tau }_{\text{coll }}}.$$(1)

The characteristic time can be expressed in terms of the dust number density n_d, collision cross-sectional area (σ_d = 4π s²), and relative velocity υ_rel during collision as τ_coll = (n_dσυ_rel)⁻¹. The mass change rate can then be computed as

$${\left(\frac{\text{d}m}{\text{d}t}\right)}_{\text{grow }}=4\pi {\rho }_{\text{d}}{S}^2{\upsilon }_{\text{rel }}\text{, }$$(2)

where ρ_d is the dust volume density, while the relative velocity υ_rel can be expressed as

$${\upsilon }_{\text{rel }}=\sqrt{2}{\upsilon }_{\text{t}}\frac{\sqrt{\text{Sc}-1}}{\text{Sc}}.$$(3)

The turbulent velocity υ_t is given by

$${\upsilon }_{\text{t}}=\sqrt{2^{1/2}Ro\alpha }{c}_{\text{g}}\text{, }$$(4)

where Ro is the Rossby number, which is a constant equal to 3.0 (Stepinski & Valageas 1997), α the turbulent viscosity parameter (Shakura & Sunyaev 1973) and c_g the sound speed. The Schmidt number Sc of a dust grain measures the grain’s coupling to the vortex

$$\text{Sc}=(1+\text{St})\sqrt{1+\frac{\text{Δ}{\upsilon }^2}{{\upsilon }_{\text{t}}^2}},$$(5)

where ∆υ is the difference in velocity between grains due to friction with the gas. The Stokes number St measures the coupling between gas and dust, and it is expressed as

$$\text{St}=t_{\text{S}}{\Omega }_{\text{K}}\text{, }$$(6)

where t_s is the stopping time, that is, the time needed for a grain to reach the gas velocity, and Ω_K the Keplerian frequency. More details about the implementation and the growth model are provided in Vericel et al. (2021).

2.2 Porosity evolution during growth

In this section, we describe the modelling of porosity evolution during the growth of an aggregate. Suyama et al. (2008, 2012), and Okuzumi et al. (2012) have developed a porosity evolution model. Suyama et al. (2008) and Okuzumi et al. (2012) used the Smoluchowski equation to evolve an aggregate, which makes the porosity model discrete – the filling factor computed at iteration n directly depends on the quantities from iteration n − 1. While this allows for a detailed modelling of porosity evolution in an N-body simulation of grain collisions, the model cannot be used directly in global 3D simulations involving dust growth that rely on other methods. The continuous model developed by Garcia (2018) enabled its implementation in the LYONSPH code (Barrière-Fouchet et al. 2005; Laibe et al. 2008; Garcia & Gonzalez 2020).

An aggregate is a collection of several monomers, which we considered as compact spheres with mass m₀, size a₀, and intrinsic density ρ_s. Experiments (Blum & Wurm 2000; Shimaki & Arakawa 2012a), numerical simulations (Seizinger et al. 2012; Ringl et al. 2012; Wada et al. 2013; Planes et al. 2021), and observations (Güttler et al. 2019) show that dust grains are not perfect spheres but rather have ovoid or irregular shapes. For simplicity, we considered grains to be spherical aggregates with volume V. Not considering the grains as spherical would require knowing the collision history of the grain, which would greatly increase the complexity of the modelling.

We defined the filling factor ϕ related to the porosity 𝔭 as the ratio between the volume occupied by the monomers V_mat and the volume of the aggregate V:

$$\varphi =\frac{V_{\text{mat}}}{V}\text{\hspace{1em} and \hspace{1em}}\mathfrak{p}+\varphi =1.$$(7)

The mass of an aggregate with mass m and size s can thus be simply computed as

$$m={\rho }_{\text{s}}\varphi \frac{4\pi }{3}{s}^3.$$(8)

Two energies can be associated with an aggregate. The first is the kinetic energy upon impact when two identical grains with mass m collide with a relative velocity υ_rel. In the referential of the centre of mass, the kinetic energy is expressed as

$$E_{\text{kin}}=\frac{m}{4}{\upsilon }_{\text{rel}}^2,$$(9)

where the factor of 1/4 arises from the reduced mass for two grains of identical mass.

The second energy is the rolling energy E_roll, which corresponds to the energy required to rotate a monomer by 90° around a connection point. The ability of monomers to reorganize leads to internal rearrangement of the aggregate (Dominik & Tielens 1997)

$$E_{\text{roll }}=6{\pi }^2{\gamma }_{\text{s}}{a}_0{\xi }_{\text{crit }},$$(10)

where γ_s is the surface energy (energy per unit area) of a monomer, ξ_crit is the critical rolling displacement, and a₀ is the monomer size. The value of ξ_crit is still poorly constrained. According to Dominik & Tielens (1997), the critical separation δ_c between two monomers before they separate,

$${\delta }_{\text{c}}={\left(\frac{27{\pi }^2}{2}\frac{{\gamma }_{\text{s}}^2{a}_0}{{ℰ}^2}\right)}^{\frac{1}{3}},$$(11)

is of the same order of magnitude as ξ_crit (Chokshi et al. 1993). Thus, we can rewrite Eq. (10) as

$$E_{\text{roll }}={\left(2916{\pi }^8\frac{{\gamma }_{\text{s}}^5{a}_0^4}{{ℰ}^2}\right)}^{\frac{1}{3}}\approx 302{\left(\frac{{\gamma }_{\text{s}}^5{a}_0^4}{{ℰ}^2}\right)}^{\frac{1}{3}},$$(12)

with ɛ as Young’s modulus. Depending on the ratio between these two energies, two different growth regimes can be distinguished with distinct porosity evolution depending on the grain mass m (Suyama et al. 2012; Okuzumi et al. 2012).

• In the ‘hit & stick’ regime, the grains are small and thus coupled to the gas. Collisions occur at very low relative velocities, with a condition on kinetic energy (E_kin < 3b E_roll, Suyama et al. 2008), where b is a numerical factor equal to 0.15 (Okuzumi et al. 2012).

• As the grains grow, the kinetic energy during impacts increases, becoming much larger than the rolling energy. As a result, energy is dissipated by internal restructuring of the grain structure, leading to compaction. This is referred to as the collisional compression regime.

We recall in Appendix A the equations describing the porosity of dust grains in the different expansion and compression regimes presented in Garcia (2018) and Garcia & Gonzalez (2020). Figure 1 summarises the evolution of the filling factor ϕ of a grain, starting as a compact (ϕ = 1) monomer (s = a₀ = 0.2 µm) as its size s increases, for different distances to the star.

Fig. 1 Filling factor ϕ as a function of size s for different distances from the star R0 for silicate grains composed of 0.2 µm monomers in the case of pure growth in our standard disc model (see Sect. 3.3).

2.3 Bouncing

In addition to growth, grains can also bounce. We have improved upon the model initially developed by Garcia (2018). Wada et al. (2011), Shimaki & Arakawa (2012a), and Arakawa et al. (2023) have numerically and experimentally shown that grains are no longer capable of bouncing when the filling factor drops below ϕ = 0.3. This limit corresponds either to small grains composed of a few monomers or to large aggregates spanning several kilometres, when considering only growth – see Fig. 1. Once below this filling factor, grains can only grow or fragment when colliding.

During bounce, some of the energy is used to bind the grains together, and another part is dissipated in a deformation wave (Thornton & Ning 1998). Several velocities can be distinguished. The first is the sticking velocity υ_stick, corresponding to the velocity at which grains no longer stick together systematically. Johnson et al. (1971) and Thornton & Ning (1998) provided the sticking velocity υ_stick for the collision of two identical grains

$${\upsilon }_{\text{stick }}=4.23{\left(\frac{{\gamma }_{\text{s}}^5{s}^4}{m^3{ℰ}^2}\right)}^{1/6}\text{. }$$(13)

Dominik & Tielens (1997) provided another formula for υ_stick based on the mass and size of monomers. Their formula differs by a factor of 2. This factor arises from the medium in which the wave propagates, which is the monomer in the case of Dominik & Tielens (1997), and a collection of monomers for Eq. (13). However, Shimaki & Arakawa (2012b, 2021) have shown from experiments that the Young’s modulus ɛ of water ice varies with porosity by about an order of magnitude. Its value as a function of ϕ is not known for silicates. Since υ_stick ∝ ɛ^−1/3, by using the value of ɛ for monomers and retaining the factor of two, a reasonable compromise is achieved. Therefore, we kept

$$\begin{array}{l}{\upsilon }_{\text{stick }}=8.46{\left(\frac{{\gamma }_{\text{s}}^5{s}^4}{m^3{ℰ}^2}\right)}^{1/6}\hfill \\ =2.82{\text{mms}}^{-1}{\left(\frac{{\rho }_{\text{s}}}{1000{\text{kgm}}^{-3}}\right)}^{-1/2}{\left(\frac{{\gamma }_{\text{s}}}{0.1{\text{Jm}}^{-2}}\right)}^{5/6}\hfill \\ \times {\left(\frac{ℰ}{10\text{GPa}}\right)}^{-1/3}{\left(\frac{s}{1\text{mm}}\right)}^{-5/6}{\varphi }^{-1/2},\hfill \end{array}$$(14)

where we have used Eq. (8) to obtain a dependence on size and filling factor. With the material properties we adopted for silicate and water ice grains (see Sect. 3.3), we obtained υ_stick = 2.89 mm s⁻¹ and 5.26 mm s⁻¹, respectively, for millimeter-sized grains with ϕ = 0.3. Weidling et al. (2012) argued that instead of γ_s, one should use an effective surface energy, which accounts for porosity and is smaller, together with a smaller value of the Young’s modulus for porous aggregates. Both compensate in our expression for υ_stick, leading to values that are within one order of magnitude of their experimental data.

However, in their study of collisions of less massive silicate aggregates, Kothe et al. (2013) found values of υ_stick of several to tens of cm s⁻¹. Additionally, other studies suggested that, instead of the Young’s modulus of the monomers, the compressive strength of the dust aggregates should be used, for which Blum & Schräpler (2004), Güttler et al. (2009) and Seizinger et al. (2012) found values of ~10³ Pa for ϕ = 0.2 and > 10⁴ Pa for ϕ > 0.3 (see Eq. (31) in Seizinger et al. 2012), which increases υ_stick by two orders of magnitude. We thus also considered values of υ_stick 10 to 100 times larger than that given by Eq. (14) – see below.

Two regimes can be distinguished: elastic collision and plastic collision. The separation between the two regimes is given by the yield velocity υ_yield. Thornton & Ning (1998) expressed υ_yield as a function of the contact radius a_yield when yield occurs or of the limiting contact pressure p_yield. However, as these material properties are difficult to obtain, and are in particular unknown for water ice, they suggested to turn to values from impact experiments. In theirs, Weidling et al. (2012) found a value for the yield velocity of SiO₂ aggregates about 5 times larger than their sticking velocity. According to Musiolik et al. (2016), and Yasui et al. (2017), water ice has more favourable sticking properties, and we adopted υ_yield = 10υ_stick. More recently, Kimura et al. (2020) demonstrated that silicates are more prone to sticking together with a higher surface energy. We thus kept a factor of ten between the two velocities for silicates as well and took

$${\upsilon }_{\text{yield }}=10{\upsilon }_{\text{stick }}=84.6{\left(\frac{{\gamma }_{\text{s}}^5{s}^4}{m^3{ℰ}^2}\right)}^{1/6}.$$(15)

Consequently, the coefficient of restitution e can now be defined as follows:

$$e\equiv \frac{{\upsilon }_{\text{rel},\text{f}}}{{\upsilon }_{\text{rel}}},$$(16)

where υ_rel,f is the relative velocity after bounce. After the collision, a portion of the kinetic energy is used for plastic deformation given by (1 − e²)E_kin, and the remaining part serves to separate the grains with the remaining kinetic energy e² E_kin. Depending on the relative velocity, the coefficient of restitution can then be expressed (Thornton & Ning 1998; Garcia 2018) as follows:

$${\upsilon }_{\text{rel }}\le {\upsilon }_{\text{stick }}e=0$$(17)

$${\upsilon }_{\text{stick }}<{\upsilon }_{\text{rel }}\le {\upsilon }_{\text{yield }}e=1$$(18)

$$\begin{array}{c}{\upsilon }_{\text{rel }}\ge {\upsilon }_{\text{yield }}e=[1.2\sqrt{3}\left(1-\frac{1}{6}{\left(\frac{{\upsilon }_{\text{yield }}}{{\upsilon }_{\text{rel }}}\right)}^2\right)\\ {\times {\left(1+2\sqrt{1.2{\left(\frac{{\upsilon }_{\text{rel }}}{{\upsilon }_{\text{yield }}}\right)}^2-0.2}\right)}^{-1/2}-{\left(\frac{{\upsilon }_{\text{stick }}}{{\upsilon }_{\text{rel }}}\right)}^2]}^{1/2}.\end{array}$$(19)

When υ_rel ≫ υ_yield, the coefficient of restitution tends towards zero, as shown in Fig. 2. This indicates that almost all the energy is used in plastic deformation, and the grain is strongly compacted. Beyond a velocity threshold υ_end, the resulting velocity from the remaining kinetic energy e² E_kin becomes so low that after a bounce, the grain is too slow to collide with another grain and can no longer stick. To be consistent with experiments from Shimaki & Arakawa (2012b), we set υ_end = 0.03 υ_rel. Solving Eq. (19) for e = 0.03 with υ_yield = 10υ_stick allowed us to express υ_end as a function of υstick so that

$${\upsilon }_{\text{end }}\approx 2\times {10}^7{v}_{\text{stick }}\approx 2\times {10}^8{\left(\frac{{\gamma }_{\text{s}}^5{s}^4}{m^3{ℰ}^2}\right)}^{1/6}\text{. }$$(20)

Typical values of υ_end are on the order of 10–100 km s⁻¹, much larger than fragmentation thresholds.

Fig. 2 Coefficient of restitution e when υrel ≥ υyield, given by Eq. (19).

Fig. 3 Model comparison for compression pressure. Tatsuuma et al. (2023) connect the regime with small ϕ given by Kataoka et al. (2013) to the large ϕ with the model of Güttler et al. (2010) and Seizinger et al. (2012).

Sticking probability and filling factor

Empirically, Weidling et al. (2012) determined the probability that a grain sticks, or not. As previously mentioned, when υ_stick < υ_rel < υ_end, an aggregate can bounce instead of stick. Naturally, the closer υ_rel is to υ_stick, the higher the chances of sticking, and conversely for υ_end. Thus, Garcia (2018), inspired by Weidling et al. (2012), defined the probability 𝒫 of sticking for a grain based on its velocity according to

$$=\{\begin{array}{c}1,{\upsilon }_{\text{rel }}<{\upsilon }_{\text{stick }},\\ \frac{\log \left(v_{\text{rel }}\right)-\log \left({\upsilon }_{\text{end }}\right)}{\log \left(v_{\text{stick }}\right)-\log \left({\upsilon }_{\text{end }}\right)},{\upsilon }_{\text{stick }}\le {\upsilon }_{\text{rel }}<{\upsilon }_{\text{end }},\\ 0.{\upsilon }_{\text{rel }}\ge {\upsilon }_{\text{end }},\end{array}$$(21)

Depending on the sticking probability, the value of the growth rate will be different. When the grain does not stick, there is no contribution to the growth rate. Thus, it can be expressed as

$${\left(\frac{\text{d}m}{\text{d}t}\right)}_{\text{bounce }}={\left(\frac{\text{d}m}{\text{d}t}\right)}_{\text{grow}}.$$(22)

However, when considering bouncing, a problem arises. It is impossible to write the filling factor of a grain in terms of its mass when the deformation is plastic. Garcia (2018) provided a formula for the final filling factor during bounce based on the initial filling factor and an arbitrary time step t

$${\varphi }_{\text{bounce }}\approx {\varphi }_{\text{i}}{\left(1-\frac{\text{Δ}V}{V_{\text{i}}}\right)}^{-n}\text{, }$$(23)

where n = t/τ_coll and ∆V is the volume change during compression due to plastic deformation. Shimaki & Arakawa (2012a) provided an expression to compute ∆V, which is proportional to the kinetic energy,

$$\text{Δ}V=\frac{\left(1-e^2\right)E_{\text{kin}}}{2{\text{Υd}}_{}},$$(24)

where Υ_d = Υ₀ϕ^𝒩 is the dynamic compression resistance, and Υ₀ is the reference dynamic compression resistance. Mellor (1974) gave Υ₀ = 9.8 MPa and 𝒩 = 4 for water ice. However, Υ, which depends on porosity, is rarely studied because it is difficult to determine, especially for highly porous materials. Therefore, this formula cannot be used directly as Υ remains unknown for most materials that one might want to simulate, such as silicates. However, for numerous bounces, the process can be approximated as static compression. Kataoka et al. (2013) provided a formula for static compression pressure for highly porous grains only, as given in Eq. (A.9), and Güttler et al. (2009, 2010) and Seizinger et al. (2012) provided an expression for values of ϕ > 0.2. Very recently, Tatsuuma et al. (2023) have developed a formula for compression pressure that depends only on E_roll and ϕ, which is valid for both small and large ϕ, shown in Fig. 3,

$$P_{\text{bounce }}=\left(\frac{E_{\text{roll }}}{a_0^3}\right){\left(\frac{1}{\varphi }-\frac{1}{0.74}\right)}^{-3}\text{. }$$(25)

The value ϕ = 0.74 corresponds to the maximum possible filling factor, which is that of either a cubic or a hexagonal close packed arrangement of equal spheres.

In order to illustrate what the porosity evolution of a grain experiencing only growth and bouncing would be, we consider for the time being all values of the filling factor before restricting it again to values above ϕ = 0.3. Figure 4 allows for the comparison of the bounce model from Garcia (2018) with the equation involving Υ, which is valid only for water ice (or a mixture of ice and silicates), and our model where P_bounce replaces Υ in Eq. (24). The difference between both models is observed when the grain is highly compacted, with a filling factor ϕ > 0.5. With the model from Garcia (2018), the grains continue to be efficiently compacted until they reach the maximum possible filling factor. In our model, as grains become highly compacted, it becomes increasingly difficult to compress them. This is due to the fact that as ϕ approaches 1, the value of the compression pressure changes only slightly. The dashed and dot-dashed curves show that considering values of υ_stick 10 or 100 times larger only delays the onset of bouncing but has very little influence on the final value of the filling factor. The light tan polygons trace the parameter space explored in the experiments from Weidling et al. (2012) and Kothe et al. (2013), as estimated from the latter’s Fig. 7.

We now relax our working assumption and return to the situation where bouncing is only effective for ϕ ≥ 0.3, which excludes the grey-shaded area in Fig. 4. Grains experiencing only growth quickly dive below the threshold and never bounce again. A mechanism able to compact grains, such as the one presented in Sect. 2.5, is thus needed for bouncing to reappear. The sizes at which such compacted grains will cross above the threshold depend mainly on the compacting mechanism. Their subsequent evolution depends very little on the value of υstick (see Sect. 2.5).

In conclusion, to compute the filling factor due to bouncing, we used for values larger than 0.3 only

$${\varphi }_{\text{coll \& bounce }}=\{\begin{array}{c}{\varphi }_{\text{coll }},{\upsilon }_{\text{rel }}<{\upsilon }_{\text{yield }},\\ {\varphi }_{\text{coll}}𝒫+(1-𝒫){\varphi }_{\text{bounce }},{\upsilon }_{\text{yield }}\le {\upsilon }_{\text{rel }}<{\upsilon }_{\text{end }},\\ {\varphi }_{\text{bounce }},{\upsilon }_{\text{rel }}\ge {\upsilon }_{\text{end }},\end{array}$$(26)

where ϕ_coll is the filling factor resulting from collisions in one of the regimes described in Appendix A.2, and ϕ_bounce is computed using Eq. (23). The final filling factor is the largest value between ϕ_gas, ϕ_grav, and ϕ_{coll & bounce}. Our model is an improvement over the model from Garcia (2018). It now takes into account the saturation effect of compression for filling factors close to one and can be used with different materials as long as the surface energy, Young’s modulus, and monomer size are known.

Fig. 4 Growth and bouncing of a water ice grain with a0 = 0.2 µm in our standard disc model (see Sect. 3.3) for two different distances. Our model (solid lines) uses Eq. (25), while the model from Garcia (2018, dotted lines) uses the Υd given by Shimaki & Arakawa (2012a). The dashed and dot-dashed lines show the impact of multiplying Vstick from Eq. (14) by 10 or 100. The light tan polygons trace the parameter space probed in experiments by Weidling et al. (2012) and Kothe et al. (2013). The grey shaded area marks the range of filling factors for which bouncing is not effective (see text).

2.4 Fragmentation

The other natural process that appears in the life of grains is fragmentation. When the relative velocity between grains exceeds the fragmentation threshold υ_frag, then the grain fragments instead of growing (Tanaka et al. 1996). The kinetic energy upon impact is such that the grain’s structure cannot absorb it, breaking the bonds between the monomers constituting the aggregate.

A model developed by Kobayashi & Tanaka (2010) and used by Vericel et al. (2021) allows for a gradual fragmentation depending on the value of υ_rel with respect to υ_frag

$${\left(\frac{\text{d}m}{\text{d}t}\right)}_{\text{frag }}=-4\pi \frac{{\upsilon }_{\text{rel}}^3}{{\upsilon }_{\text{rel }}^2+{\upsilon }_{\text{frag }}^2}{\rho }_{\text{d}}{s}^2.$$(27)

In this model, when the relative velocity is close to the threshold, the mass loss is less significant. A fragmenting grain loses half of its mass after a collision time τ_coll (υ_rel = υ_frag) or more (υ_rel > υ_frag). In the case where υ_rel ≫ υ_frag, we recover the scenario described by Gonzalez et al. (2015), where the entire grain fragments after a collision time, independently of υ_rel.

In reality, both experimental studies (Blum & Wurm 2000; Shimaki & Arakawa 2012a; Weidling et al. 2012) and numerical experiments (Dominik & Tielens 1997; Okuzumi et al. 2009, 2012; Kataoka et al. 2013; Krijt et al. 2015; Planes et al. 2021) showed more complex behaviours that are challenging to model. Fragmentation involves a cascade of sizes and depends on parameters such as composition, porosity, shape, impact parameter, mass ratio, and impact velocity. The presented model follows the one introduced by Kobayashi & Tanaka (2010) using a fragmentation definition that relies on several approximations to capture the essential physics according to the relative collision velocity. To determine υ_frag, an estimate of the energy required to break the bond between two monomers E_break is given by Dominik & Tielens (1997),

$$E_{\text{break }}=C_{\text{break }}{F}_{\text{c}}{\delta }_{\text{c}}\text{, }$$(28)

where C_break is a constant that they take equal to 1.8. The term F_c is the critical pulling force between two monomers,

$$F_{\text{c}}=3\pi {\gamma }_{\text{s}}{a}_0.$$(29)

We can then derive the expression for E_break:

$$E_{\text{break }}\simeq C_{\text{break }}48{\left(\frac{{\gamma }_{\text{s}}^5{a}_0^4}{{ℰ}^2}\right)}^{\frac{1}{3}}.$$(30)

To find the total energy, we just need to know the number of monomers and the number of bonds between monomers,

$$E_{\text{frag }}=\kappa N_{\text{tot }}{E}_{\text{break }},$$(31)

where N_tot is the total number of monomers, and κ is a numerical factor representing the number of bonds and depending on the species and porosity. For erosion, κ = 1–2, and for fragmentation, κ = 3–10 (Blum & Wurm 2000; Wada et al. 2007). We can then define υ_frag in terms of E_frag as

$${\upsilon }_{\text{frag }}=2\sqrt{\frac{E_{\text{frag }}}{m}}=2\sqrt{\frac{48\kappa N_{\text{tot }}{C}_{\text{break }}}{m}}{\left(\frac{{\gamma }_{\text{s}}^5{a}_0^4}{{ℰ}^2}\right)}^{\frac{1}{6}}\text{. }$$(32)

These equations depend on γ_s and a₀. An error in estimating these parameters can thus lead to changes in the values of E_frag and υ_frag.

Research on the value of υ_frag according to different types of materials is an active topic, and the community does not yet agree on the values. Different studies have yielded a range of values for different materials like silicates and water, and there is a complex interplay of factors. For example, Blum & Wurm (2008) and Güttler et al. (2010) deduced υ_{frag, Si} ~ 1 m s⁻¹ for silicates. Wada et al. (2009, 2013) found a similar value for silicates, around υ_{frag, Si} ~ 5 m s⁻¹, as well as that for water, ${\upsilon }_{{\text{frag, H}}_2\text{O}}\sim 60-70$ m s⁻¹. However, the method used by Yamamoto et al. (2014) yielded a different value for water, ${\upsilon }_{{\text{frag, H}}_2\text{O}}=56\text{\hspace{0.17em}\hspace{0.17em}m\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}$, based on the link between surface energy (and hence fragmentation threshold) and material melting temperature. The values for silicates differ widely: Yamamoto et al. (2014) found γ_s = 0.3 J m⁻², while Kimura et al. (2020) found γ_s = 0.15 J m⁻². Silicates are stronger than previously thought due to experimental issues with the original estimations of γ_s. Silicates easily absorb water from the air, which reduces the surface energy of the original aggregate, thereby reducing the adhesion ability of the silicate grains. In the case of water ice, Shimaki & Arakawa (2012a) measured ${\upsilon }_{{\text{frag, H}}_2\text{O}}~15\text{\hspace{0.17em}\hspace{0.17em}m\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}$, while collision simulations gave higher values (Wada et al. 2009). The same value was computed by Gonzalez et al. (2015) using the experimental measurements of the energy required to fragment a unit mass of grains (~55 J kg⁻¹). For the purpose of comparison with certain previous studies (Garcia & Gonzalez 2020; Vericel et al. 2021), we adopted the value of ${\upsilon }_{{\text{frag, H}}_2\text{O}}\approx 15\text{\hspace{0.17em}\hspace{0.17em}m\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}$ for water.

However, pure water ice grains are unlikely. Most often, grains are composed of solid materials like silicates, covered by a layer of another more volatile element such as water or organic matter. The fragmentation velocity is primarily determined by the material forming the surface layer, and not the internal material. Therefore, we chose different thresholds for the silicates, as there are still uncertainties about their values, to mimic the effect of a layer of volatiles: υ_frag = 5, 10, or 20 m s⁻¹ for more fragile materials like water or CO, and υ_frag = 40 ms⁻¹ for modelling pure silicate grains or those surrounded by organic materials.

Contrary to Shimaki & Arakawa (2012a), who demonstrated that a mixture of silicates and water is more fragile than pure ice, we assumed the opposite here, based on measurements conducted in recent years (Kimura et al. 2020; San Sebastián et al. 2020). Pure silicate grains are now considered more resistant than water ice grains and we suppose a combination of the two materials is more robust than water ice but less so than silicates. However, this remains to be confirmed experimentally. Finally, one should note that we assumed surface energies and fragmentation thresholds to be independent of monomer size for the values we have chosen.

Fig. 5 Filling factor ϕf /ϕi as a function of υrel. Left: experimental values for different initial filling factors ϕi. Those provided by Ringl et al. (2012) are in orange, and those from Gunkelmann et al. (2016) are in blue, red, and green. Right: comparison between our model given by Eq. (35) and the fits of the curves given by Ringl et al. (2012) and Gunkelmann et al. (2016) in our fragmentation module.

2.5 Compaction during fragmentation

The second improvement we made to the porosity model is to consider grain compaction during fragmentation. The main motivation behind this enhancement is to attempt to explain the high filling factor values for intermediate-sized grains (10 µm-cm) found on comets (Güttler et al. 2019) and inferred from observations through polarimetric measurements (Kataoka et al. 2015, 2016, 2019; Tazaki et al. 2019; Tazaki & Dominik 2022). Since neither growth, static gas compaction, nor bouncing can account for such high factors (ϕ > 0.1) as grains grow by increasing their porosity, we turned to fragmentation to obtain more compact grains. Sirono (2004) showed that the filling factor remains constant during fragmentation. For simplicity, this assumption, which considers that the impact energy is used to break bonds between monomers rather than for internal grain restructuring, was adopted by Garcia (2018).

However, Ringl et al. (2012) and Gunkelmann et al. (2016) found that the filling factor after fragmentation is approximately 1.5 to 2 times larger than the initial filling factor. The remaining energy for a grain to undergo compaction can be computed by subtracting the energy E_frag used to fragment the grain from the kinetic energy E_kin released upon impact. We can rewrite E_frag as the energy E_break required to break bonds between monomers, multiplied by the number of ejected monomers:

$$E_{\text{comp }}=E_{\text{kin }}-E_{\text{frag }}=E_{\text{kin }}-\frac{\left(2m_{\text{i}}-m_{\text{f}}\right)}{m_0}\kappa E_{\text{break }}\text{. }$$(33)

Here, (2m_i – m_f)/m₀ corresponds to the number of monomers that were ejected during fragmentation. In the case of fragmentation, surface monomers are the first to be ejected since they have fewer bonds (Ringl et al. 2012; Gunkelmann et al. 2016). We therefore considered κ = 3 defined by Eq. (31), which has a lower value for fragmentation compared to those proposed by Blum & Wurm (2000) and Wada et al. (2007). In this sense, we modeled the fact that ejected grains are closer to the surface, and the remaining energy is used for rearranging the internal monomers. A value of κ = 10 means more energy is used to break bonds, allowing aggregates to remain at larger sizes before starting to compact at a constant fragmentation threshold. We can compute, similarly to bouncing (Eq. (24)), the volume change after fragmentation ∆V, which is given by

$$\text{Δ}V=\frac{E_{\text{kin }}-\left(2m_{\text{i}}-m_{\text{f}}\right)\kappa E_{\text{break }}{m}_0^{-1}}{2\left(\frac{E_{\text{roll }}}{a_0^3}\right){\left(\frac{1}{\varphi }-\frac{1}{0.74}\right)}^{-3}}.$$(34)

Furthermore, we fitted the data obtained by Ringl et al. (2012) and Gunkelmann et al. (2016), see the left panel of Fig. 5, which shows the ratio ϕ_f /ϕ_i of the final filling factor to the initial one as a function of the relative velocity υ_rel for three values of ϕ_i. It should be noted that Gunkelmann et al. (2016) defined υ_frag as the point when a grain starts to lose a monomer, which is equal to 0.17 m s⁻¹, and differs from our definition as the point where the grain loses half of its mass. We deduced that the ratio ϕ_f/ϕ_i varies as exp (1 − (υ_rel /υ_frag)²). In the right panel of Fig. 5, we present a comparison between our model and the fits to data from Ringl et al. (2012) and Gunkelmann et al. (2016) for two distances to the star, 10 and 100 au. Our model shows good agreement with the fits. Divergence between the model and fits occurs when υ_rel/υ_frag exceeds 1.15 to 1.2. However, grains never reach these values as they fragment before reaching such velocities. One should note the axes are different. On the left panel, the initial filling factors are constant for all values of υ_rel/υ_frag. On the contrary, on the right panel, the fit extracted from Ringl et al. (2012) and Gunkelmann et al. (2016) has been implemented in our growth and fragmentation model and compared to our compaction model. In this case, the value of ϕ₀ varies after each collision. The filling factor after fragmentation is computed in the same way as bouncing, using Eq. (23).

$${\varphi }_{\text{frag–comp }}\approx {\varphi }_{\text{i}}{\left(1-\exp \left(1-{\left(\frac{{\upsilon }_{\text{rel }}}{{\upsilon }_{\text{frag }}}\right)}^2\right)\frac{\text{Δ}V}{V_{\text{i}}}\right)}^{-n},$$(35)

where n represents the number of collisions during a time step. The effect of compaction can be seen in Fig. 6, produced with PAMDEAS with growth, bouncing, fragmentation and radial drift (see Sect. 3.1). In the left panel, where compaction is not taken into account, grains reach a growth-fragmentation equilibrium and stay in the same range of size and filling factor. ϕ remains below 0.3, bouncing thus never occurs. In the right panel, with compaction, the filling factor of fragmenting grains increases while their size decreases. Because ϕ reaches values larger than 0.3, bouncing then appears and must be considered. It prevents grains from growing significantly when compacted. The inset shows that using a value of υ_stick 100 times larger than that given by Eq. (14) – see Sect. 2.3 – has a very limited impact and only slightly lowers the final filling factor.

Fig. 6 Filling factor ϕ as a function of size s for different distances to the star r for silicate grains composed of 0.2 µm monomers in our standard disc model with growth and fragmentation. Left: without compaction. Right: with compaction. The inset is a zoom of the light grey box showing only tracks for R0 = 10 and 100 au. The dash-dotted lines show the impact of multiplying υstick from Eq. (14) by 100.

2.6 Rotational disruption

Rotational disruption was identified as a possible growth barrier in protoplanetary discs by Tatsuuma & Kataoka (2021): highly porous grains can be disrupted by the gas-flow torque when the tensile stress from the centrifugal force exceeds their tensile strength. Michoulier & Gonzalez (2022a) studied the effects of rotational disruption in 1D simulations, in this work we implemented their equations in PHANTOM.

2.7 Summary

The algorithm to compute the final filling factor can be summarised as follows:

Compute ϕcoll, ϕgas, ϕgrav
Compute ϕmin = max(ϕcoll, ϕgas, ϕgrav)
if υrel < υfrag then
    Compute ϕgrow
    if (Grains can bounce) then
        Compute ϕcoll & bounce
        Compute ϕf = max(ϕcoll & bounce, ϕmin)
    else
        Compute ϕf = max(ϕgrow, ϕmin)
  end if
else
  if (Fragmentation with compaction) then
        Compute ϕfrag–comp
        Compute ϕf = max(ϕfrag–comp, ϕmin)
  else
  Compute ϕf = max(ϕi, ϕmin)
  end if
end if
ϕfinal = min(ϕf, 0.74)

3 Numerical simulations

3.1 The 1D code PAMDEAS

PAMDEAS (Porous Aggregate Model and Dust Evolution in protoplAnetary discS) is a one-dimensional code, presented in Michoulier & Gonzalez (2022a), to study the evolution ofporous grains within protoplanetary discs considering different physical processes, such as growth, fragmentation, radial drift or rotational disruption. We added new physics such as aeolian erosion (Rozner et al. 2020; Grishin et al. 2020; Michoulier et al. 2024), bouncing and compaction during fragmentation. This code allowed us to follow the evolution of a given number of particles from a set of initial conditions in a static, vertically isothermal, non self-gravitating gas disc. The gas surface density and temperature profiles are given by power laws of indices p and q, respectively,

$${\text{Σg}}_{}(r)={\text{Σ}}_{\text{g},0}{\left(r/R_0\right)}^{-p}$$(36)

and

$$T_{\text{g}}(r)=T_{\text{g},0}{\left(r/R_0\right)}^{-q},$$(37)

where R₀ is a reference radius. The disc extends from R_in to R_out and has a total mass M_disc. As the gas disc structure is held fixed, the dust-to-gas ratio ε is kept constant and the feedback of dust on gas is neglected. Each grain evolves separately, starting from an initial size equal to the monomer size, and the code tracks its full evolution. Despite its limitations, simulations with PAMDEAS are useful to understand the different stages experienced by a single grain.

Fig. 7 Radial (left panel) and vertical (right panel) profiles of size s obtained with PAMDEAS using our standard disc model after 100 yr. Grains grow from the monomer size at different initial positions. Colour indicates time t.

3.2 The 3D smoothed particle hydrodynamics code PHANTOM

PHANTOM (Price et al. 2018) is a 3D smoothed particle hydrodynamics (SPH; Lucy 1977; Gingold & Monaghan 1977) code for hydrodynamics and magnetohydrodynamics, designed to be efficient. It is public and widely used. In this paper, we introduce a new module in which we implemented the algorithms to treat dust porosity with all the physics presented in Sect. 2.22.6, both in the ‘dust-as-mixture’ (also known as ‘one-fluid’) and ‘dust-as-particles’ (or ‘two-fluid’) formalisms. This development is tightly coupled to the dust growth module implemented and described in detail by Vericel et al. (2021). PHANTOM, which takes into account collective effects and dust vertical settling, and simulates the coupled evolution of gas and dust due to aerodynamic drag (including the back-reaction of dust on gas), produces more realistic dust distributions that can be compared to observations of discs.

We first setup a gas disc with an exponentially decreasing surface density profile:

$$\text{Σ}(r)={\text{Σ}}_0\left(1-\sqrt{\frac{R_{\text{in }}}{r}}\right){\left(\frac{r}{R_{\text{c}}}\right)}^{-p}\exp \left[-{\left(r/R_{\text{c}}\right)}^{2-p}\right]\text{, }$$(38)

where R_c is the cut-off radius. With the exponential tapering, the outer density profile is smooth, rather than being sharply truncated at R_out, allowing better control over how relaxation behaves in the outer region. Furthermore, using a pseudo-relaxed disc as the initial state makes the computation faster as particles do not have to significantly move radially. Simulating the innermost region with sufficient resolution is computationally expensive because it leads to a large dynamic time ratio between the inner and outer edges. As a result, we set the accretion radius, which represents the radius within which a particle is considered accreted, to R_in. We also remove any particle moving further than 1000 au away to prevent completely isolated particles. The gas is locally isothermal, with a temperature profile set by Eq. (37), and non self-gravitating.

Similarly to Price & Laibe (2015), and in contrast to Vericel et al. (2021), we first performed a simulation solely with gas to avoid any spurious behaviour of the dust grains during gas relaxation. We increased the number of particles and the disc mass by about 20% with respect to the target values to compensate for the mass loss due to particle accretion onto the star during relaxation. When the density profile of the outer gas region stabilizes, typically after 7 to 8 orbits at the outer edge, the disc is fully relaxed.

We then added the dust phase with the same spatial distribution as the relaxed gas, with an initially uniform dust-to-gas ratio ε₀. Unlike in PAMDEAS, we did not take the monomer size as the initial grain size. Indeed, if the simulation was initialised with a single user-defined value of the initial size for all grains, the filling factor computed from s with Eq. (C.1) and the mass computed from Eq. (8) would result in less porous, more massive grains in the upper disc layers. This is unphysical. Furthermore grains in the inner region are expected to have grown during the early stages of star formation, as shown in Bate (2022) and Lebreuilly et al. (2023). We thus used an initial state where each SPH particle has a size depending on its location (r, ɀ), following fits to a PAMDEAS simulation of the evolution of grains from monomers for 100 yr. Figure 7 shows the resulting radial and vertical size profiles. The left panel shows that the midplane radial size distribution can be fitted by a power law

$$s(r)=s_0{\left(r/R_0\right)}^{-2}\text{, }$$(39)

and in the right panel, the vertical dependence of the final size is ∝ exp (−ɀ²/H²), where H is the disc scale height, is obtained, here at 100 au. After the initialisation, the evolution variable is the grain mass, and the mass growth rate and filling factor are computed using the equations in Sect. 2.

In this paper, the gas and dust disc is evolved as a single set of SPH particles, using the dust-as-mixture algorithms of Price & Laibe (2015) and Ballabio et al. (2018), based on Laibe & Price (2014).

3.3 Setup

We chose to use a disc model that represents an ‘average disc’ (Williams & Best 2014). The mass of the star was fixed at M_star = 1 M_⊙, and the mass of the disc was M_disc = 0.01 M_⊙. The temperature was set by the choice of the aspect ratio (H/R)₀ = 0.0895 at R₀ = 100 au, with q = 0.5.

For PAMDEAS, we used R_in = 1 au, R_out = 300 au and p = 1. All grains evolved from an initial size equal to the monomer size. We ran the simulations up to t = 1 Myr.

For PHANTOM, we took R_in = 10 au, R_out = 400 au and p = 0.75. This value of p, combined with the exponential tapering, led to a profile whose slope, after relaxation, is similar to that of PAMDEAS between 10 and 300 au. We set the number of particles to 1.2 million, the turbulent viscosity parameter (Shakura & Sunyaev 1973) to α = 5 × 10⁻³ (by setting α_AV = 0.1658) and the initial dust-to-gas ratio ε₀ to a typical value of 1%. We evolved the simulations for 300 000 yr.

For both codes, the size of monomers was set to a₀ = 0.2 µm, in agreement with recent observations. Indeed, Tazaki & Dominik (2022) and Tazaki et al. (2023) have shown that to accurately reproduce the polarization degree with respect to the scattering angle and the SED, the small grains present on the surface of the IM Lupi disc are most likely fractal aggregates with 0.2 µm monomers. Verrios et al. (2022) reached the same conclusion with respect to the dust settling in IM Lupi, requiring porosities of ϕ ≲ 0.1 to match the observations. We chose to use two different species commonly found in discs: water ice and silicates.

For silicate grains, we chose an intrinsic density of ρ_s = 2 700 kg m⁻³ and a surface energy γ_s = 0.2 J m⁻², in agreement with Yamamoto et al. (2014), who estimate γ_s = 0.3 J m⁻². This value is of the same order of magnitude as γ_s = 0.15 J m⁻², found by Kimura et al. (2015, 2020), through experimental measurements with sicastar^® aggregates (micromod Partikeltechnologie GmbH). The Young’s modulus was ɛ = 72 GPa (Yamamoto et al. 2014), which implies, assuming that the critical separation between two monomers δ_c before their separation is of the same order of magnitude as ξ_crit (Chokshi et al. 1993), a value of ξ_crit ≈ 6 A.

For pure water ice, the intrinsic density of monomers was ρ_s = 1 000 kgm⁻³, with a surface energy of γ_s = 0.1 J m⁻². The critical separation is still debated, as are many other properties of dust in protoplanetary discs. We chose a Young’s modulus of ɛ = 9.4 GPa, following Yamamoto et al. (2014), which results in a value of ξ_crit ~ 10 A, close to ξ_crit ~ 8 A used by Wada et al. (2011) and Tatsuuma & Kataoka (2021).

The different values of the fragmentation thresholds for both species are given in Sect. 2.4. In our simulations, silicates serve as the reference material, and υ_frag = 20 m s⁻¹ is the reference threshold. We have conducted numerous simulations with different parameters. For easier reading, each simulation is assigned a name, as listed in Table 1.

4 Results

4.1 Effect of porosity on grain growth

In this section, we compare the evolution of compact (i.e. ϕ = 1) and porous dust, including only growth and fragmentation for now.

4.1.1 1D simulations with PAMDEAS

First, we will examine the effect of porosity using the 1D code PAMDEAS. Here the disc is static, and only the evolution of grains is taken into account. Figure 8 makes it easy to observe the influence of porosity between the GF-Si-comp-Vf20 simulations on the left and GF-Si-a02-Vf20 on the right. The upper left (right) panel shows the size of compact (porous) silicate grains as a function of the distance to the star r. The colour indicates the coupling with the gas through the Stokes number. Each grain grows from the monomer size without drifting because it is coupled to the gas (vertical lines). The grains start drifting inwards when St ~ 5 × 10⁻³, and they continue to grow until they reach the fragmentation threshold (horizontal plateau), maintaining an equilibrium between growth and fragmentation. While the Stokes numbers (computed as per Eq. (5) in Vericel et al. 2021) are similar in both simulations, with St ~ 10⁻¹ when the grains reach the fragmentation threshold, the sizes are drastically different. Compact grains reach about 100 µm, whereas porous grains reach about 3 cm. These porous grains have a filling factor ϕ ~ 2 × 10⁻³, making them ~50 000 times more massive than compact grains¹. Indeed, porous grains are capable of growing to larger sizes due to their larger cross-sectional area (Garcia & Gonzalez 2020). Moreover, porous grains remain well coupled to the gas even at sizes of 100 µm.

The bottom panels of Fig. 8 display the size evolution. Porous and compact grains reach their maximum size in similar timescales. A porous grain initially at 15 au grows to 3 cm in 1000 yr, while a compact grain starting at the same location reaches 100 µm. Since the fragmentation threshold Stokes numbers are similar (St ~ 0.1), the aerodynamic evolution is the same, and hence the grains drift at the same speed. While drifting, the grains reach the inner regions where relative velocities are higher. This explains why the size plateaus have decreased to 8 mm for porous grains and 70 µm for compact grains at a distance of r = 10 au, even though the grains were able to grow to larger sizes at a larger distance.

A similar behaviour can be seen for the comparison between the GF-H2O-comp-Vfl5 and GF-H2O-aO2-Vfl5 simulations on Fig. 9. Porous water ice grains grow to a size of 20 cm, while compact grains only grow to 200 µm. These sizes are larger than for silicate grains because of a lower intrinsic density, which increases the growth-fragmentation equilibrium size. At the fragmentation threshold, porous water ice grains are 5 × 10⁵ times more massive than compact grains.

In all cases, compaction by gas or self-gravity is never reached, as fragmentation maintains the grain sizes below a few metres, a range where gas compaction appears. Porous grains can therefore grow to much larger sizes (from centimetres to decimetres) and masses compared to their compact counterparts (a few hundred micrometres), demonstrating that porosity must be considered in dust evolution.

Fig. 8 Comparison between PAMDEAS simulations GF-Si-comp-Vf20 (left) and GF-Si-aO2-Vf20 (right). The top panels show the size s as a function of distance r, with colour indicating the Stokes number St. The bottom panels show the size s as a function of time t, with colour indicating the Stokes number St as well. The leftmost curve corresponds to the grain initially at 15 au, and the rightmost curve to 300 au.

Fig. 9 Same as the top panels of Fig. 8 for PAMDEAS simulations GF-H2O-comp-Vf15 (left) and GF-H2O-a02-Vf15 (right).

4.1.2 3D global simulations with PHANTOM

We now explore the effects of porosity with PHANTOM. Figure 10 compares the GF-Si-comp-Vf20 simulation (left) to GF-Si-a02-Vf20 (right). The top panels display the dust-to-gas ratio in the meridional plane (r, ɀ) at different simulation times. Grey areas represent the initial dust-to-gas ratio, red regions are dust-enriched, and blue regions are dust-depleted. It can be noticed that in the early times, porous grains settle faster in the midplane compared to compact grains, as they grow rapidly in the inner regions. However, in the outer regions, porous grains, due to their low density, remain more coupled to the gas. They drift and settle more slowly, resulting in lower dust concentration and a thicker dust disc. This effect is evident at times t = 50 or 100 kyr. At the end of the simulation, the compact dust disc is less radially extended (≈ 100 au) and has a higher dust concentration, around ε ≈ 3–4, with a thickness of ~10 au. In contrast, the porous dust disc extends out to 150 au with a dust-to-gas ratio of ε ≈ 0.8–1. Porous grains thus allow the dust in the outer regions to be retained for a longer time while achieving large dust concentrations. Similarly to simulations with growth and fragmentation in Garcia (2018), the disc of porous grains after 300 kyr is thinner than that of compact grains.

The bottom panels of Fig. 10 display the radial size distribution of dust grains. Colour represents the Stokes number St, where red to green grains are strongly coupled to the gas, blue grains are marginally coupled, and purple grains are decoupled. To understand what happens in dust-concentrated regions, particles with ε < 5 × 10⁻³ have been deliberately excluded, representing regions highly depleted in dust where little activity occurs. Porous grains can grow to more substantial sizes, on the order of millimetres or centimetres, with filling factors ϕ ~ 5 × 10⁻³−10⁻²(see Fig. 14, left panel on second row), across a large portion of the disc, up to 200 au. In contrast, compact grains struggle to reach millimetre sizes and the largest ones remain in the 100–500 µm range. This is inconsistent with observations that report the presence of millimetre-sized grains in the midplane. Hence, porosity helps grain growth to more significant sizes and masses, approximately 5 to 5000 times larger, with dust concentrations and disc thicknesses similar to compact grains, and sizes compatible with observations. However, the resulting filling factors (Fig. 14) are too small by a factor of ~ 10–20 compared to observations (see Sect. 5). In both cases, grains do not cross St = 1 and do not decouple from the gas.

Despite special attention to avoid numerical artefacts, some remain, inherent to the evolution model itself. For porous grains at time t = 100 kyr, a region close to 500 au where grains have grown slightly more than those closer in can be observed. This can be explained rather simply. Initially, in the farthest regions (>300 au), there are only monomers. Monomers at 600 au have larger St values than those at 500 au. They drift and settle faster than grains closer to the star can grow, resulting in a more dust-enriched zone, also visible in the top right panel. All grains drift, settle, and grow, and this artefact related to the model disappears without consequence.

The same behaviour is seen for water ice grains, whose evolution is shown in Fig. 11. Compared to silicate grains, the top panels demonstrate that water ice dust discs are thicker (~20 au) and extend radially further (> 400 au). This is due to the lower density of water ice grains, resulting in smaller Stokes numbers for a given position and size. Water ice grains thus settle and drift more slowly than silicate grains. However, dust-to-gas ratios of around ε ≈ 2–3 for compact grains and ε ≈ 0.7–1 for porous grains are still reached. Compact water ice grains struggle to reach millimetre sizes in the midplane (lower panels), although their growth is faster than that of compact silicate grains due to their lower density and larger cross-sectional area for a given mass. They remain relatively small, in the hundreds of microns range. In contrast, porous water ice grains can easily attain centimetre scale sizes. Millimetre- to centimetre-sized grains are found out to 400–500 au. The resulting filling factors (see also Fig. 14, bottom left), around 10⁻³ − 5 × 10⁻³, are also too small compared to observations by a factor 20 to 100. Grains are still maintained at St values around 0.1; none manage to decouple in the inner regions. This is not the case for some grains that reach St = 1 in the outer disc region, where a lower gas density and a higher growth rate enable these St values to be reached. Thus, porosity allows grains to grow more easily, reaching millimetre sizes, and forming dust discs as thin as those composed of compact grains.

Fig. 10 Comparison between PHANTOM simulations GF-Si-comp-Vf20 (left) and GF-Si-a02-Vf20 (right). The top panels show the dust-to-gas ratio in colour in the (r, ɀ) plane. The bottom panels show the radial grain size distribution, with colour representing the Stokes number St.

4.2 Effect of bouncing and compaction during fragmentation

In this section, we investigate the effects of compaction during fragmentation and bouncing at different fragmentation thresholds. Disc observations reveal that larger dust particles in the midplane are not as porous as the values obtained in the absence of compaction (Sect. 4.1). Kataoka et al. (2016), Guidi et al. (2022), and Zhang et al. (2023) find that grains have filling factors on the order of 10% with a product sϕ larger than 100 µm, meaning grains with s > 1 mm.

Fig. 11 Same as Fig. 10 for PHANTOM simulations GF-H2O-comp-Vf15 (left) and GF-H2O-a02-Vf15 (right) using water ice instead of silicate grains.

4.2.1 1D simulations with PAMDEAS

Figure 12 compares grain sizes as a function of distance from the star between simulations without bouncing and compaction during fragmentation (GF) on the left, and those with (GBFc) on the right. The top three rows correspond to simulations with silicates for three fragmentation thresholds: υ_frag,Si = 10, 20, and 40 m s⁻¹, from top to bottom, while the bottom row shows simulations with water ice for υ_frag,H2O = 15 m s⁻¹. For all GBFc simulations, the grains grow to the same sizes and filling factors as the GF simulations. However, when the fragmentation threshold is reached, the grains in the left panels remain in an equilibrium between growth and fragmentation because, during fragmentation, the grains retain their filling factor. In the right panels, aggregates are compacted during fragmentation, the filling factor increases, and sizes become smaller until a new equilibrium is reached. The value of ϕ for compacted grains is not the maximum possible value ϕ_max = 0.74. The maximum value that is reached is instead around 0.5–0.6 as in equilibrium, grain growth tends to decrease ϕ, while compaction by fragmentation or bouncing tends to increase ϕ. According to Fig. 6 and Eq. (A.9), the closer ϕ gets to the maximum value, the harder it is to compact the aggregate. The slow growth of grains thus compensates for compaction towards the maximum value.

For simulations with υ_frag,Si = 10 m s⁻¹, the grains reach sizes of a few millimetres, with a maximum around 3 mm between 20 and 50 au. When compacted, the aggregate sizes then drop to about 40–80 µm at maximum compaction before being accreted by the star. The filling factor reaches a minimum of 10⁻² just before the fragmentation threshold, before being compacted to ϕ = 0.5. The aggregates are compacted very efficiently, as a grain initially growing at 300 au ends up fully compacted during its drift beyond 100 au. Lastly, one can note that the compacted grains drift as fast as those that do not undergo compaction, as indicated by the black points and line. This is due to the fact that the product sϕ during compaction remains roughly the same as when there is no compaction. As St ∝ sϕ, a similar Stokes number results in a similar drift.

With υ_frag,Si = 20 m s⁻¹, grains in the right panel reach sizes of around 3 cm before being compacted. As the fragmentation threshold is higher, the grains have larger Stokes numbers and therefore drift more during compaction. The maximum compaction occurs between 10 and 20 au, resulting in sizes ranging from 100 to 300 µm and filling factors between 0.3 and 0.5.

For υ_frag,Si = 40 m s⁻¹, grains are able to grow almost without encountering the fragmentation barrier. In both cases, the grains start to fragment when reaching sizes between 20 and 50 cm, with filling factors smaller than 3 × 10⁻³. With such a high fragmentation threshold, grains undergo mostly pure growth in the majority of the disc and compaction does not operate. Moreover, due to this free growth, grains in the inner regions are accreted more rapidly since the St values they reach are larger, around 0.5–1 between 10 and 20 au, and over 0.1 in the rest of the disc.

Finally, in the case of water ice (bottom row of Fig. 12), the aggregates reach decimetre sizes in both cases. When compacted, their size is reduced by a factor of 500. The size reached at maximum compaction ranges from 300 to 500 µm, with a filling factor ϕ ~ 0.4. Thus, a behaviour similar to simulations with silicates and υ_frag,Si = 20 m s¹־ is observed, but with larger variations in size and filling factor.

Figure 13 shows the evolution of size and filling factor, with colour representing the distance r. All grains start to grow in the hit & stick regime from the size of a monomer a₀ = 0.2 µm and ϕ = 1. They then transition to a different growth regime (ϕ_{Ep, St<1}) resulting in a change in slope around the micrometre size, depending on the species (see Fig. 6).

In left panels for simulations without compaction, once the fragmentation threshold is reached, the grains remain confined in the ‘tip’ at the bottom right. On the other hand, in the right panels, grain compaction is observed, leading to a decrease in size while the filling factor increases towards unity. With υ_frag,Si = 10 m s¹־, grains are compacted up to ϕ = 0.5 and a minimum size of 40 µm, even at larger distances r. Bouncing doesn’t occur in this case, as fragmentation is highly effective. For υ_frag,Si = 20 m s⁻¹, grains fragment at larger sizes and smaller ϕ. However, once the fragmentation threshold is reached, they are compacted up to ϕ = 0.6 and sizes of 100 µm. Bouncing begins to appear but has no influence on grain evolution and is not visible in the figure. For υ_{frag, Si} = 40 m s¹־, aggregates fragment after growing to sizes of around a decimetre and ϕ ~ 2 × 10⁻³ They then fragment and are compacted to ϕ = 0.7 and sizes of 700 µm.

Finally, in the case of water ice, the same behaviour is observed. Aggregates are compacted until they reach ϕ = 0.4 and s = 300 µm. As mentioned in Sect. 2.3, bouncing is accompanied by growth (but hardly visible in Fig. 13). Similarly to Fig. 6, the characteristic plateau is seen, where grains move from left to right and bottom to top, as they are compacted by bouncing during growth before continuing to drift and be accreted.

Thus, compaction during fragmentation has a significant impact on grain evolution. Taking porosity into account allows grains to overcome the bouncing barrier, which occurs at grain sizes of a few to several tens of microns. Bouncing doesn’t interfere with aggregate growth before reaching the fragmentation threshold, as ϕ < 0.3. When ϕ becomes larger than 0.3 again, grains grow and don’t bounce but fragment. It is necessary for the grains to be compacted and at sufficiently small sizes for bouncing to play a role again before fragmentation.

Fig. 12 Comparison between PAMDEAS simulations with growth and fragmentation (GF-*-aO2-*, left) and simulations with growth, bouncing, and fragmentation with compaction (GBFc-*-aO2-*, right), up to t = 1 Myr. Grain sizes are given as a function of the distance from the star, with colour indicating the filling factor. Simulation parameters are indicated in each row. Black dots and dashed lines mark the time t = 50 kyr, providing a ‘snapshot’ of the size distribution at that moment.

Fig. 13 Comparison between PAMDEAS simulations with growth and fragmentation (GF-*-a02-*, left) and simulations with growth, bouncing, and fragmentation with compaction (GBFc-*-a02-*, right) up to t = 1 Myr. The filling factor ϕ is plotted against grain size, with colour indicating the distance r. Simulation parameters are indicated in each row.

4.2.2 3D global simulations with PHANTOM

Figure 14 compares PHANTOM simulations without (top) and with (bottom) compaction during fragmentation and bouncing. We show the radial grain size distribution, colour-coded with the filling factor, for simulations with silicates and υ_frag,Si = 10, 20 and 40 m s¹־, and with water ice, from top to bottom in each panel (in the same order as in Fig. 12). In all cases, the maximum grain sizes obtained with or without compaction are similar, unlike in the 1D simulations, despite efficient compaction. It should be noted that 1D and 3D simulations cannot be directly compared as the former follow the evolution of single grains in the disc midplane, while the latter deal with a population of grains at various altitudes, whose size and porosity evolution depends on varying disc conditions as they settle and drift. The effect of compaction (bottom panel) can best be seen at t = 100 kyr interior to 200 au where the largest grains are still relatively porous, with ϕ ranging from a few 10⁻² to a few 10⁻³ depending on the simulation, and mostly at large r, while the smallest grains have been compacted to ϕ ~ 0.4 or larger. In simulations without compaction (top panel), the larger grains are also more porous than smaller ones, as expected from their evolution during growth (Sect. 2), but the porosity range is much smaller. At t = 300 Myr, the most compacted grains, in the midplane, are substantially smaller than in simulations without compaction.

Simulations with silicate grains and υ_frag,Si = 10 and 20 m s⁻¹ (first two rows in both panels) are similar: without or with compaction, grains reach sizes of a few millimetres at t = 100 kyr for υ_frag = 10 m s⁻¹, while they grow faster and larger, up to a few centimetres, for υ_{frag, Si} = 20 m s⁻¹. However, after 300 kyr, compacted grains are at most 600–700 µm in size for υ_frag = 10 m s⁻¹ and a few millimetres for υ_{frag, Si} = 20 m s⁻¹. The former case is difficult to reconcile with observations of grains of mm size or larger in protoplanetary discs. When υ_{frag, Si} = 40 m s⁻¹ (third rows), the threshold is high enough for grains to remain mostly in the pure growth regime, similarly to what was seen with PAMDEAS (Sect. 4.2.1), with almost no fragmentation or compaction, except in the very inner disc. Here, grain sizes reach several tens of metres and their Stokes number exceeds unity. This breaks the terminal velocity approximation used in the dust-as-mixture formalism (Youdin & Goodman 2005; Laibe & Price 2014) and those simulations should not be considered valid. Finally, for simulations GF-H2O-a02-Vf15 and GBFc-H2O-a02-Vf15 (bottom rows), the same pattern as for simulations GF-Si-a02-Vf20 and GBFc-Si-a02-Vf20 is observed, with grains reaching slightly larger sizes. The joint evolution of grain size and filling factor is described in Appendix D.

The thickness of the of the dust discs can be examined in Fig. 15, showing the size in the (r, ɀ) plane. Like in other figures, only particles with ε ≥ 5 × 10⁻³ are shown to eliminate dust-depleted regions. The light grey background indicates the gas disc’s thickness. An interesting result can be noted: whether the grains are compacted or not, the thickness of the dust discs is very similar between 10 and 200 au, where dust is most abundant, and is only slightly larger for compacted dust. In all cases, compared to the compact grain discs of simulations GF-Si-comp-Vf20 (Fig. 10) and GF-H2O-comp-Vf15 (Fig. 11), the porous grain discs are just as thin in the inner regions, or even thinner in the early stages for r < 200–300 au. However, the porous grain discs are thicker farther from the star and have larger sizes than their compact grains counterparts.

Finally, the maximum size of compact grains is larger in PHANTOM simulations compared with PAMDEAS by a factor of ~5 (comparing the left panels of Fig. 8 to the left panels of Fig. 10) because dust settling in 3D increases dust density in the midplane, which helps grains to grow to larger sizes. On the other hand, the sizes, filling factors, and St values of porous grains are similar with both codes. When grains are highly porous, their growth is fast and the increase in growth rate due to settling provides little assistance – the limiting factor is the fragmentation threshold. Porous grains that undergo compaction are in a intermediate situation.

Fig. 14 Comparison of the radial grain size distribution, with colour representing the filling factor, between PHANTOM simulations with growth and fragmentation (GF-*-a02-*, top panel) and with growth, bouncing, and fragmentation with compaction (GBFc-*-a02-*, bottom panel). Simulation parameters are indicated in each row. We note that the colorbar range is different from that in Fig. 12.

Fig. 15 Comparison between PHANTOM simulations with growth and fragmentation (GF-*-a02-*, top panel) and simulations with growth, bouncing, and fragmentation with compaction (GBFc-*-a02-*, bottom panel). The grain size s is shown in the meridian plane (r, ɀ), with the light grey background indicating the gas disc thickness. Simulation parameters are indicated in each row.

4.3 Effects of a snow line

Here, we examine the influence of the CO ice line on the GBFc-Si-a02-Vf20 simulation. The sublimation temperature for CO is 20 K, which corresponds to a position of the snow line at r ~ 100 au in our disc model. We set the fragmentation threshold for the outer region, where grains are assumed to be surrounded by CO ice, to υ_frag,ext = 5 m s⁻¹. In the inner region, instead of taking a threshold of 15 m s⁻¹, which would correspond to grains made entirely of water ice, we set υ_frag,int = 20 m s⁻¹ to model silicate grains assumed to be surrounded by a layer of water ice. The outer layer only affects adhesive properties, as discussed in Vericel & Gonzalez (2020). The result from this simulation, GBFc-Si-a02-Vf205, are shown in Fig. 16.

The top panel shows the grain size s in the (r, ɀ) plane, with a cut-off at ε ≥ 5 × 10⁻³, and the light grey gas disc in the background, while the central panel displays the radial grain size distribution coloured by the Stokes number. Both panels show two distinct regions. Interior to 100 au, grains have settled into a very thin disc and formed large aggregates with sizes of several millimetres up to a centimetre. They have St ~ 5 × 10⁻² when they are compacted, mainly between 50 and 100 au, down to a few millimetres. Beyond 100 au, grains have not been able to grow to large sizes due to the lower fragmentation threshold exterior to the snow line, where the balance between growth and fragmentation keeps them at smaller sizes ranging from ~100 µm when they are still porous, down to a few tens of µm once compacted. Their smaller St slows down their settling and drift. The dust disc has thus settled much less compared to simulation GBFc-Si-a02-Vf20 shown in Fig. 15, with a thicker and more extended outer disc.

Finally, the bottom panel shows the filling factor versus the grain size s, with colour indicating the distance r. With the snow line, the interpretation becomes more complex, as grains can fragment and compact in two different regions, leading to compaction even for small sizes in the outer regions. Grains between 100 and 150 au are compacted, and once past the snow line, they have the opportunity to grow again before reaching the inner fragmentation threshold of 20 m s⁻¹. Thus, we observe not one but two vertical columns in the plot due to compaction during fragmentation.

Snow lines in this configuration are an effective way to form relatively compact dust grains with ϕ ~ 0.1 far from the star, as long as the fragmentation threshold is not too high; otherwise, the grains would drift inward faster and would not be compacted efficiently.

Fig. 16 PHANTOM GBFcS-Si-a02-Vf205 simulation. Top: aggregate sizes in the same plane (r, ɀ). Middle: radial grain size distribution, with colour representing the Stokes number St. Bottom: filling factor ϕ plotted against grain size s, with colour indicating the distance r.

4.4 Effects of rotational disruption

We present in Appendix E the results of 3D simulations of rotational disruption, complementing the 1D PAMDEAS simulations of Michoulier & Gonzalez (2022a) and removing some of their limitations. They show that rotational disruption has a negligible effect on the dust evolution in discs.

5 Discussion

We have shown that simulations with porous grains allow for the growth of grains to larger sizes and masses regardless of the species or fragmentation threshold. Hence, as indicated by Garcia & Gonzalez (2020), porosity enables large dust grain formation. Our 1D simulations using PAMDEAS and our 3D simulations with PHANTOM both result in centimetre-sized grains for porous grains when including growth and fragmentation only. With compact or compacted grains due to fragmentation, PAMDEAS tends to underestimate the maximum grain size by a factor of a few. Nevertheless, PAMDEAS efficiently models the dust evolution, with or without considering porosity evolution.

We found that simulations involving compaction contradict the earlier simulations without compaction, which indicated that larger grains were more porous and located in the midplane (Garcia 2018; Garcia & Gonzalez 2020). With compaction, we find that large aggregates in the midplane are compact. In the upper layers, grains are nearly as large but more porous and consequently less massive, undergoing settling. Finally, small grains are either monomers or fractal aggregates composed of a few dozen monomers. This suggests that the size and porosity distribution of dust mainly depends on grain altitude. Our simulations with porous and compacting grains provide better explanations for observations that report both non-porous millimetre-sized grains in the midplane and more porous millimetre-sized grains above and below.

Filling factors of dust grains in discs are generally between 0.1 and 1 according to observations (Guidi et al. 2022; Zhang et al. 2023; Ohashi et al. 2023). Figure 17 summarizes the best models able to match observations of HL Tau (Zhang et al. 2023) and DG Tau (Ohashi et al. 2023). In both studies, the authors assumed either compact (ϕ = 1) of porous (ϕ = 0.1 for Zhang et al. 2023 and 0.2 for Ohashi et al. 2023) grains and fit for the grain sizes for which SED and polarization data are best reproduced. In the inner disc (interior to 40 au), compact grains of a few hundred µm are found to match both discs. However, for porous grains, sizes of 1 mm or larger are preferred for HL Tau, while values from 100 µm closer in up to 1 mm further out work best for DG Tau. For HL Tau, Zhang et al. (2023) found that compact grains of any size cannot reproduce the data between 20 and 60 au. At larger distances, the best fit is found for similar sizes to the inner disc for HL Tau, but for compact grains of 1 cm or porous grains larger than 3 cm for DG Tau. However, the authors of both studies cautioned that the low signal-to-noise ratio in the outer regions makes it difficult to distinguish between different grain populations.

Our simulations provide an explanation for these grain distributions. Thanks to porosity, all grains are capable of rapid growth to reach millimetre sizes that are compacted in the inner regions with consistent filling factors (ϕ ~ 0.1–1). In the outer regions, smaller grains (10–100 µm) are prevalent, with varying porosities (ϕ ~ 0.01–0.5), which is consistent with the observations of HL Tau. The large grains (>3 cm) seen in the DG Tau disc, however, are challenging to reproduce in our simulations. Grains with these sizes are present, but they are more porous with ϕ ~ 10⁻³–10⁻².

Table 2 presents a summary of the classifications of cometary dust and its morphology from the Rosetta and Stardust missions (Güttler et al. 2019). These observations found that the most common grains are represented by the porous particle group, with ϕ in the [0.05,0.9] range, similar to the compacted grains in our simulations, although their size can be much larger. Fractal aggregates are rarer, also in agreement with our simulations. The solid particle group is similar to our monomers. Compaction of grains during their evolution is thus necessary to explain both disc and comet observations.

Some limitations of our work are inherent to the codes themselves (see Sect. 2), in particular for PAMDEAS with its one-dimensional nature assuming a static gas disc, limitations that are overcome with PHANTOM. However, the growth and fragmentation model remains the mono-disperse model (Stepinski & Valageas 1997) in both cases. In particular, it cannot account for collisions between particules of very different sizes or masses. High-mass-ratio collisions have been found to produce stronger compression (Tanaka et al. 2023). Taking them into account would help aggregates to reach their compacted state sooner.

The mono-disperse model nonetheless results in global size distributions that are similar to those obtained with polydisperse models (Gonzalez et al. 2017; Vorobyov et al. 2018). The reader is referred to Laibe et al. (2008); Vericel & Gonzalez (2020); Vericel et al. (2021); Michoulier et al. (2024) for additional discussions of this formalism.

The 3D simulations use the dust-as-mixture formalism in the terminal velocity approximation, which is no longer valid when dust decouples from the gas. Moreover, the dust-as-mixture (one fluid) simulations do not depict the formation of self-induced dust traps, whereas previous dust-as-particles simulations showed this mechanism (Michoulier & Gonzalez 2022b). Future dust-as-particles simulations, albeit significantly more time-consuming, should verify the presence of these traps, as well as demonstrate whether dust can effectively decouple.

The fragmentation model is not entirely consistent. We set the fragmentation threshold to a predefined, fixed value. In practice, the fragmentation threshold should depend on the size and porosity of grains. Garcia (2018) tested a model with a variable threshold with limited success. However, there is no existing model or study that establishes this relationship. A consistent model is needed to better understand the fragmentation of aggregates and the its influence on the evolution of dust grains.

In a similar vein, the difficulty to obtain material properties for different compounds due to the complexity of experimental studies and the large parameter space they would need to cover makes assumptions unavoidable in order to develop a general model. Some of them may be lifted in the future as more data becomes available, for example measurements of the Young’s modulus over a wide range of porosities for several species of astrophysical interest.

Fig. 17 Grains sizes compatible with all the data obtained from SED and polarization observations of HL Tau (Zhang et al. 2023) and DG Tau (Ohashi et al. 2023) for different values of the assumed filling factor.

6 Conclusions

A key missing piece in our understanding of the puzzle of planet formation is how sub-micron-sized dust grains coagulate into planetesimals while avoiding the fragmentation and radial drift barriers. These obstacles must be overcome for dust to survive long enough to form kilometre-sized rocky bodies. While solutions have been proposed to explain planetesimal formation, mechanisms such as the streaming instability (Youdin & Goodman 2005) involve specific conditions, such as the presence of large dust particles settled in the midplane. However, there is no clear answer as to how to achieve these conditions. Growth up to a few microns is relatively straightforward, but reaching millimetre-sized or larger grains, as observed, is difficult because of radial drift.

To address these issues, we explored the evolution of porous grains, including various mechanisms that destroy them. We developed a new physical model and code module to account for growth and porosity evolution, as well as bouncing and fragmentation with compaction, within the mono-disperse approximation. We investigated the influence of porosity on dust evolution using the 1D code PAMDEAS (Michoulier & Gonzalez 2022a) and the 3D SPH code PHANTOM (Price et al. 2018). We performed simulations varying the porosity, compaction, the constituent species of the grains, and the fragmentation threshold. Our conclusions are as follows:

1. We confirm that porosity enables the formation of millimetre- to centimetre-sized grains (Garcia 2018). The main limiting factor is the fragmentation threshold.

2. With growth and fragmentation alone, we found large values for the filling factor (ϕ ~ 10⁻⁴–10⁻³) that are incompatible with observations, which show compact grains in the inner regions.

3. When including compaction during fragmentation and bouncing, we formed grains of several hundred micrometres to a few millimetres, but this time with filling factors compatible with observations (ϕ ∊ [0.1; 1]).

4. The 3D simulations with high turbulent viscosity showed signs of rotational disruption in the intermediate layers of the disc when aggregates are settling in the outer regions, but we found the overall impact of rotational disruption on dust evolution to be small.

Figure 18 summarises the key mechanisms involved in the evolution of grains. Overall, we found that considering porosity, compaction effects, and snow lines can reconcile observations with theoretical models of planet formation. Future work may involve developing a coherent model for the fragmentation threshold, further investigating the impact of snow lines, and using simulations to generate synthetic images to refine our understanding of observations.

Fig. 18 Diagram summarising the different physical processes related to dust seen and studied in this paper.

Acknowledgements

We thank the anonymous referee for their report, which allowed us to improve the clarity of this paper. The authors acknowledge funding from ANR (Agence Nationale de la Recherche) of France under contract number ANR-16-CE31-0013 (Planet-Forming-Disks) and thank the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) for its financial support within the Plan France 2030 of the French government operated by the ANR. This research was partially supported by the Programme National de Physique Stellaire and the Programme National de Planétologie of CNRS (Centre National de la Recherche Scientifique)/INSU (Institut National des Sciences de l’Univers), France. We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements No 210021 and No 823823 (DUSTBUSTERS). DJP is grateful for Australian Research Council Discovery Project funding via DP220103767 and DP240103290. We thank the Australian-French Association for Research and Innovation (AFRAN) for financial support. Figures were made with the Python library matplotlib (Hunter 2007). The data necessary to recreate the Phantom simulations and to reproduce the corresponding figures are available at https://doi.org/10.5281/zenodo.12729815.

Appendix A Equations for porosity evolution during growth

A.1 Hit & stick regime

In the hit & stick regime, the mass of the grain doubles during a collision, and a volume of void V_void is formed. According to Okuzumi et al. (2009), this captured volume V_void corresponds to 99% of the initial volume V_i, thus adding to the initial volumes of the two grains before the collision. The filling factor ϕ ∝m/V is therefore multiplied by a factor of 2/2.99. For an arbitrary number of collisions n, there is a geometric relationship between the initial filling factor ϕ₀ and the final ϕ_f, given by ϕ_f = (2/2.99)ⁿϕ₀ (Garcia & Gonzalez 2020). Moreover, the mass of a grain after the same number of collisions has grown by m_f = 2ⁿm₀, m₀ being the monomer mass. By isolating n and combining the two expressions, we obtained the filling factor in the hit & stick regime (assuming E_kin ≪ E_roll, Garcia & Gonzalez 2020):

$${\varphi }_{\text{h\&s}}={\left(\frac{m}{m_0}\right)}^{\psi }\text{, }$$(A.1)

with ψ = ln (2/2.99) / ln 2 ≈ −0.58, where for the monomer ϕ₀ = 1.

A.2 Collisional compression regime

As both the grain mass and the relative velocity υ_rel increase, the kinetic energy increases. Since υ_rel depends on the Stokes number (Sect. 2.1), which in turn depends on the drag regime, there is an expression for each case (Garcia & Gonzalez 2020). Suyama et al. (2012) provide the volume of a grain after collision Vf in terms of the volume of a grain before collision V_i. This relation can be transformed to incorporate filling factors using the fact that ϕ_f = m_f/V_fρ_s = 2m_i/V_fρ_s (Garcia & Gonzalez 2020):

$${\varphi }_{\text{f}}=\frac{2m_{\text{i}}}{{\rho }_{\text{s}}}{\left(\frac{{(3/5)}^5\left(E_{\text{Kin }}-3b{E}_{\text{roll }}\right)}{N_{\text{tot }}^{5}b{E}_{\text{roll }}{V}_0^{10/3}}+\frac{V_{\text{i}}^{-10/3}}{2^4}\right)}^{3/10}.$$(A.2)

Here, N_tot is the total number of monomers, which is 2m_i/m₀, and V₀ is the volume of a monomer. This model is recursive, meaning that the evolution of porosity is tracked after each collision. However, in the case of global simulations, the integration time step rarely aligns with the time interval between collisions τ_coll. Therefore, Garcia & Gonzalez (2020) transformed this discrete model into a continuous one to depend solely on the grain mass and quantities related to r, enabling its implementation in codes for global simulations. To achieve this, they assumed the approximation E_kin » E_roll.

In the Epstein and Stokes regimes for St < 1, one obtains

$${\varphi }_{\text{Ep},\text{St}<1}={\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)}^{3/8}{\left(\frac{m}{m_0}\right)}^{-1/8},$$(A.3)

$${\varphi }_{\text{St},\text{St}<1}={\left(\frac{A{c}_{\text{g}}{a}_0}{9v_{\text{mol}}\left(2^{1/5}-1\right)}\right)}^{1/3}.$$(A.4)

$$v_{\text{mol}}=\frac{5\sqrt{\pi }}{64}\frac{m_{\text{mol}}{c}_{\text{g}}}{{\rho }_{\text{g}}{\sigma }_{\text{mol}}}$$(A.5)

is the gas molecular kinematic viscosity, where m_mol and σ_mol are the mass and molecular cross section of the H₂ molecule, and A is a factor used to simplify the expressions and is given by

$$A=\frac{243\sqrt{2}\pi }{15625}\frac{\text{Ro\hspace{0.17em}\hspace{0.17em}}\alpha a_0^4{\rho }_{\text{s}}^2{c}_{\text{g}}{\text{ΩK}}_{}}{\rho b{E}_{\text{roll }}},$$(A.6)

where ρ = ρ_g + ρ_d is the total density It is observed that in the Stokes regime with St < 1, the filling factor depends only on the local disc conditions. For regimes with St > 1, collisional compression is less efficient due to the decrease in υ_rel. By neglecting E_kin, it can be derived from Okuzumi et al. (2012) and Garcia & Gonzalez (2020) that ${\varphi }_{\text{f}}\propto m_{\text{f}}^{-1/5}$ for St > 1. Therefore, the filling factors for St > 1 can be expressed as

$${\varphi }_{\text{Ep},\text{St}>1}={\varphi }_{\text{Ep},\text{St}<1}\left(M_4\right){\left(\frac{m}{M_4}\right)}^{-1/5},$$(A.7)

$${\varphi }_{\text{St},\text{St}>1}={\varphi }_{\text{St},\text{St}<1}\left(M_5\right){\left(\frac{m}{M_5}\right)}^{-1/5},$$(A.8)

where M₄ and M₅ are the transition masses between the Epstein and Stokes regimes for St < 1 and St> 1. The transition masses are given in Appendix B. The equations giving ϕ in the various regimes can also be rewritten as a function of the size s instead of the mass m, they are listed in Appendix C. However, the natural variable to describe dust evolution is m.

The filling factor resulting from the hit & stick and collisional compression regimes is noted ϕ_coll.

A.3 Static compression regime

In addition to undergoing compression due to collisions, aggregates can also be statically compacted either by the surrounding gas or by their own gravity. Kataoka et al. (2013) and Garcia & Gonzalez (2020) provided the relationship between static compression pressure and the filling factor:

$$\varphi ={\left(\frac{a_0^3{P}_{\text{comp }}}{E_{\text{roll }}}\right)}^{1/3},$$(A.9)

which is only valid for small filling factors ϕ < 0.1. In the case of compression by gas, the pressure is given by

$$P_{\text{comp }}=\frac{F_{\text{D}}}{\pi s^2}=\frac{m\text{Δ}\upsilon {\text{ΩK}}_{}}{\pi s^2\text{St}},$$(A.10)

which allows one to obtain the expression for the filling factor ϕ_gas (Garcia & Gonzalez 2020):

$${\varphi }_{\text{gas }}={\left(\frac{m_0{a}_0}{\pi E_{\text{roll }}}\frac{\text{Δ}\upsilon {\text{ΩK}}_{}}{\text{St}}\right)}^{3/7}{\left(\frac{m}{m_0}\right)}^{1/7}.$$(A.11)

For self-gravity, the exerted pressure is given by

$$P_{\text{comp }}=\frac{\text{G}{m}^2}{\pi s^4},$$(A.12)

and the expression for the filling factor ϕ_grav is (Garcia & Gonzalez 2020)

$${\varphi }_{\text{grav }}={\left(\frac{\text{G}{m}_0^2}{\pi a_0{E}_{\text{roll }}}\right)}^{3/5}{\left(\frac{m}{m_0}\right)}^{2/5}\text{, }$$(A.13)

where G is the universal gravitational constant. To obtain the final filling factor of a grain with mass m at a distance r from the star, we take the value of ϕ given by the phenomenon that compacts the grain the most, which is the largest value among static compression by gas, static compression by gravity, and compression by collision: ϕ_min = max(ϕ_coll, ϕ_gas, ϕ_grav). All these regimes are presented in Fig. 1 in the (s,ϕ) plane for different distances to the star.

Fig. B.1 Transition masses for our standard disc model and for silicate grains.

A.4 Growth from above the compression limit

The minimum filling factor must be ϕ_min. However, there are cases where the evolution of a grain (e.g. via bouncing) can produce a filling factor larger than ϕ_min, placing it in the space above the curves showed in Fig. 1. To take into account the increase of porosity of grains during growth above the limit, we can compute the filling factor as a function of the initial filling factor ϕ_i, and the final and initial masses m_f and m_i in the case where υ_rel < υ_frag·

$${\varphi }_{\text{grow }}={\varphi }_{\text{i}}{\left(\frac{m_{\text{f}}}{m_{\text{i}}}\right)}^{\psi }\text{\hspace{1em} if }{E}_{\text{kin }}<3b{E}_{\text{roll }},$$(A.14)

$${\varphi }_{\text{grow }}={\varphi }_{\text{i}}{\left(\frac{m_{\text{f}}}{m_{\text{i}}}\right)}^{-1/5}\text{\hspace{1em} otherwise.}$$(A.15)

The resulting filling factor is then the maximum of ϕ_grow and ϕ_min.

Appendix B Transition masses

The transition masses are the masses at which the behaviour of a grain changes. The different expressions of ϕ found for each regime in Sect. 2.2 are independent, the transition masses allowed us to relate them. For the transition masses M₁ and M₂ that connect the hit & stick regime to the Epstein and Stokes regimes with St < 1, we find

$$\frac{M_1}{m_0}={\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)}^{\frac{3/8}{\psi +1/8}},$$(B.1)

$$\frac{M_2}{m_0}={\left(\frac{A{c}_{\text{g}}{a}_0}{9v_{\text{mol}}\left(2^{1/5}-1\right)}\right)}^{1/3\psi }.$$(B.2)

The transition mass M₃ corresponds to the transition between the Epstein and Stokes regimes

$$\frac{M_3}{m_0}={\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)}^3{\left(\frac{A{c}_{\text{g}}{a}_0}{9v_{\text{mol}}\left(2^{1/5}-1\right)}\right)}^{-8/3},$$(B.3)

$$={\left(\frac{M_1}{m_0}\right)}^{8\psi +1}{\left(\frac{M_2}{m_0}\right)}^{-8\psi }$$(B.4)

Finally, the expressions for M₄ and M₅ are obtained by equating the expression of ϕ for the Epstein or Stokes regimes with St < 1 to the expression of ϕ obtained by setting St = 1

$$\frac{M_4}{m_0}={\left(\frac{\rho c_{\text{g}}}{{\rho }_{\text{s}}{a}_0{\text{ΩK}}_{}}\right)}^4{\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)}^{-1},$$(B.5)

$$\frac{M_5}{m_0}={\left(\frac{9v_{\text{mol}}\rho }{2{\rho }_{\text{s}}{a}_0^2{\Omega }_{\text{K}}}\right)}^{3/2}{\left(\frac{A{c}_{\text{g}}{a}_0}{9v_{\text{mol}}\left(2^{1/5}-1\right)}\right)}^{-1/6}.$$(B.6)

Here, we have not considered the other Stokes regimes that are reached for even larger sizes than the linear regime as static grain compression is always dominant compared to these regimes The various transition masses as well as the different expansion and compression regimes are shown in Fig B.1, inspired by Fig A1 of Garcia & Gonzalez (2020) We note that the transitionmass M₂ is always largerthanmass M₁ The Stokes regimes are reached in the inner regions of a disc, at a distance less than 10 au

Aappendix C: Filling factor as a function of size

The equations giving the filling factor ϕ as a function of mass m in Appendix A can be rewritten as a function of the size s as

$${\varphi }_{\text{h\&s}}={\left(\frac{s}{a_0}\right)}^{\frac{3\psi }{1-\psi }},$$(C.1)

$${\varphi }_{\text{Ep}-\text{St}<1}={\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)}^{1/3}{\left(\frac{s}{a_0}\right)}^{-1/3},$$(C.2)

$${\varphi }_{\text{St}-\text{St}<1}={\left(\frac{A{c}_{\text{g}}{a}_0}{9v\left(2^{1/5}-1\right)}\right)}^{1/3},$$(C.3)

$${\varphi }_{\text{Ep}-\text{St}>1}={\left(\frac{\rho c_{\text{g}}}{{\rho }_{\text{S}}{\Omega }_{\text{K}}{a}_0}\frac{A}{2^{3/40}-1}\right)}^{1/4}{\left(\frac{s}{a_0}\right)}^{-1/2}\text{,}$$(C.4)

$${\varphi }_{\text{St}-\text{St}>1}={\left(\frac{\rho c_{\text{g}}{a}_0}{2{\rho }_{\text{s}}{\Omega }_{\text{K}}{a}_0^2}\frac{A}{2^{1/5}-1}\right)}^{1/4}{\left(\frac{s}{a_0}\right)}^{-1/2}.$$(C.5)

Appendix D Co-evolution of grain size and filling factor in PHANTOM simulations with and without compaction

Figure D.1 shows the joint evolution of grain size and filling factor in PHANTOM simulations without (top) and with (bottom) compaction during fragmentation and bouncing, with the same order of simulations as in Fig 14 Without compaction (top panel), the trajectories for all simulations are very similar to the 1D case (Fig. 12), with a spread corresponding to the disc’s vertical extent. When compaction is considered (bottom panels), the compaction trajectory is more vertical compared to that obtained with PAMDEAS, and the bouncing plateau at ϕ ~ 0.4, primarily populated with grains very close to the disc inner edge, is more clearly visible. The overall evolution appears more complex and the s − ϕ plane is more densely populated. A large spread of filling factors can be reached for a given grain size and at a given distance, grains of various s and ϕ can be found, indicating the coexistence of porous and compacted dust.

Fig. D.1 Same as Fig. 14 but showing the filling factor ϕ plotted against grain size, with colour indicating the distance r.

Simulations with silicate grains and υ_{frag, Si} = 10 and 20 m s⁻¹ (first two rows) are again quite similar, with grains reaching larger sizes and porosities before being fragmented and compacted for the latter value. The main difference is that at t = 300 kyr almost all grains are compacted or undergoing compaction for υ_{frag, Si} = 20 m s⁻¹, with very few left on the growth branch, while for 10 m s⁻¹ a large fraction of grains are still growing, mainly in the outer disc. Despite the fact that simulations GF-Si-a02-Vf40 and GBFc-Si-a02-Vf40 (third rows) are not valid, it is still interesting to note that, due to larger sizes being reached, the gas and self-gravity compaction regimes are visible for solids larger than a few m (see Fig. 1). Again, the water ice simulations (bottom rows) are similar to that with silicate grains and υ_{frag, Si} = 20 m s⁻¹ , with larger s and smaller ϕ, except for a fraction of still growing grains at large r.

Appendix E Influence of rotational disruption

Michoulier & Gonzalez (2022a) studied the effect of rotational disruption on the evolution of porous grains using 1D simulations with PAMDEAS. 3D disc simulations taking into account the evolution of gas and dust are necessary to understand in more detail when and where rotational disruption influences dust evolution. To that effect, we performed simulations modelling this phenomenon for silicates with υ_frag,Si = 10 and 20 m s⁻¹, and for water ice with υ_frag,Si = 15 m s⁻¹, listed in the last three rows of Table 1. The grain properties were somewhat different compared to Michoulier & Gonzalez (2022a): the monomer size was set to 0.2 instead of 0.1 µm, and γ_s = 0.2 J m⁻² for silicates. In fact, with these new parameters, the tensile strength S_max should be larger by a factor of two.

Fig. E.1 Comparison between PHANTOM simulations with growth, bouncing, fragmentation with compaction, and without (GBFc-*-a02-*, top panel) and with (GBFcD-*-a02-*, bottom panel) rotational disruption showing the radial grain distribution, with colour representing the dust filling factor. Simulation parameters are indicated in each row.

Simulations with PAMDEAS showed signs of rotational disruption only for very small values of the viscosity parameter α < 5 × 10⁻⁴, lower than that used with PHANTOM. The conditions necessary for its appearance in the midplane (which PAMDEAS simulates) are therefore challenging to achieve.

Figure E.1 compares 3D simulations including growth, bouncing and fragmentation with compaction, without (top) and with (bottom) rotational disruption. It shows the radial grain distribution, with colour representing the dust filling factor ϕ. For simulations GBFc-Si-a02-Vf10 and GBFcD-Si-a02-Vf10 (top rows in both panels), very few differences are observed, which is expected since the fragmentation threshold is low enough for grains to fragment before they can be rotation-ally disrupted (which only happens at higher relative velocities and therefore larger sizes). For simulations GBFc-Si-a02-Vf20 and GBFcD-Si-a02-Vf20 (middle rows), rotational disruption is able to destroy grains during their fall towards the midplane, which slows down their settling and limits dust enrichment in the midplane. Since the dust density is slightly lower, grains do not grow as much as in the case without disruption. This is mostly visible for grains outside of 50 au at t = 100 kyr. After 300 kyr, the differences have been largely erased by the subsequent dust evolution. Finally, for simulations GBFc-H2O-a02-Vf15 and GBFcD-H2O-a02-Vf15 (bottom rows), disruption has a minimal effect on the size and porosity of dust grains, only visible outside of 150 au. Indeed, the growth of water ice grains is very rapid and largely counteracts rotational disruption. It is at t = 300 kyr that the effect becomes more visible since, unlike the simulation without disruption, many small grains are observed between 200 and 500 au.

In general, disruption has little effect on the size and porosity of grains and their evolution in the disc, with fragmentation and compaction largely dominating grain growth and destruction.

References

All Tables

All Figures

	Fig. 1 Filling factor ϕ as a function of size s for different distances from the star R0 for silicate grains composed of 0.2 µm monomers in the case of pure growth in our standard disc model (see Sect. 3.3).
In the text

	Fig. 2 Coefficient of restitution e when υrel ≥ υyield, given by Eq. (19).
In the text

	Fig. 3 Model comparison for compression pressure. Tatsuuma et al. (2023) connect the regime with small ϕ given by Kataoka et al. (2013) to the large ϕ with the model of Güttler et al. (2010) and Seizinger et al. (2012).
In the text

Fig. 4 Growth and bouncing of a water ice grain with a0 = 0.2 µm in our standard disc model (see Sect. 3.3) for two different distances. Our model (solid lines) uses Eq. (25), while the model from Garcia (2018, dotted lines) uses the Υd given by Shimaki & Arakawa (2012a). The dashed and dot-dashed lines show the impact of multiplying Vstick from Eq. (14) by 10 or 100. The light tan polygons trace the parameter space probed in experiments by Weidling et al. (2012) and Kothe et al. (2013). The grey shaded area marks the range of filling factors for which bouncing is not effective (see text).

In the text

Fig. 5 Filling factor ϕf /ϕi as a function of υrel. Left: experimental values for different initial filling factors ϕi. Those provided by Ringl et al. (2012) are in orange, and those from Gunkelmann et al. (2016) are in blue, red, and green. Right: comparison between our model given by Eq. (35) and the fits of the curves given by Ringl et al. (2012) and Gunkelmann et al. (2016) in our fragmentation module.

In the text

Fig. 6 Filling factor ϕ as a function of size s for different distances to the star r for silicate grains composed of 0.2 µm monomers in our standard disc model with growth and fragmentation. Left: without compaction. Right: with compaction. The inset is a zoom of the light grey box showing only tracks for R0 = 10 and 100 au. The dash-dotted lines show the impact of multiplying υstick from Eq. (14) by 100.

In the text

	Fig. 7 Radial (left panel) and vertical (right panel) profiles of size s obtained with PAMDEAS using our standard disc model after 100 yr. Grains grow from the monomer size at different initial positions. Colour indicates time t.
In the text

Fig. 8 Comparison between PAMDEAS simulations GF-Si-comp-Vf20 (left) and GF-Si-aO2-Vf20 (right). The top panels show the size s as a function of distance r, with colour indicating the Stokes number St. The bottom panels show the size s as a function of time t, with colour indicating the Stokes number St as well. The leftmost curve corresponds to the grain initially at 15 au, and the rightmost curve to 300 au.

In the text

	Fig. 9 Same as the top panels of Fig. 8 for PAMDEAS simulations GF-H2O-comp-Vf15 (left) and GF-H2O-a02-Vf15 (right).
In the text

	Fig. 10 Comparison between PHANTOM simulations GF-Si-comp-Vf20 (left) and GF-Si-a02-Vf20 (right). The top panels show the dust-to-gas ratio in colour in the (r, ɀ) plane. The bottom panels show the radial grain size distribution, with colour representing the Stokes number St.
In the text

	Fig. 11 Same as Fig. 10 for PHANTOM simulations GF-H2O-comp-Vf15 (left) and GF-H2O-a02-Vf15 (right) using water ice instead of silicate grains.
In the text

Fig. 12 Comparison between PAMDEAS simulations with growth and fragmentation (GF-*-aO2-*, left) and simulations with growth, bouncing, and fragmentation with compaction (GBFc-*-aO2-*, right), up to t = 1 Myr. Grain sizes are given as a function of the distance from the star, with colour indicating the filling factor. Simulation parameters are indicated in each row. Black dots and dashed lines mark the time t = 50 kyr, providing a ‘snapshot’ of the size distribution at that moment.

In the text

	Fig. 13 Comparison between PAMDEAS simulations with growth and fragmentation (GF-*-a02-*, left) and simulations with growth, bouncing, and fragmentation with compaction (GBFc-*-a02-*, right) up to t = 1 Myr. The filling factor ϕ is plotted against grain size, with colour indicating the distance r. Simulation parameters are indicated in each row.
In the text

	Fig. 14 Comparison of the radial grain size distribution, with colour representing the filling factor, between PHANTOM simulations with growth and fragmentation (GF-*-a02-*, top panel) and with growth, bouncing, and fragmentation with compaction (GBFc-*-a02-*, bottom panel). Simulation parameters are indicated in each row. We note that the colorbar range is different from that in Fig. 12.
In the text

	Fig. 15 Comparison between PHANTOM simulations with growth and fragmentation (GF-*-a02-*, top panel) and simulations with growth, bouncing, and fragmentation with compaction (GBFc-*-a02-*, bottom panel). The grain size s is shown in the meridian plane (r, ɀ), with the light grey background indicating the gas disc thickness. Simulation parameters are indicated in each row.
In the text

	Fig. 16 PHANTOM GBFcS-Si-a02-Vf205 simulation. Top: aggregate sizes in the same plane (r, ɀ). Middle: radial grain size distribution, with colour representing the Stokes number St. Bottom: filling factor ϕ plotted against grain size s, with colour indicating the distance r.
In the text

	Fig. 17 Grains sizes compatible with all the data obtained from SED and polarization observations of HL Tau (Zhang et al. 2023) and DG Tau (Ohashi et al. 2023) for different values of the assumed filling factor.
In the text

	Fig. 18 Diagram summarising the different physical processes related to dust seen and studied in this paper.
In the text

	Fig. B.1 Transition masses for our standard disc model and for silicate grains.
In the text

	Fig. D.1 Same as Fig. 14 but showing the filling factor ϕ plotted against grain size, with colour indicating the distance r.
In the text

	Fig. E.1 Comparison between PHANTOM simulations with growth, bouncing, fragmentation with compaction, and without (GBFc-*-a02-*, top panel) and with (GBFcD-*-a02-*, bottom panel) rotational disruption showing the radial grain distribution, with colour representing the dust filling factor. Simulation parameters are indicated in each row.
In the text