# Compaction during fragmentation and bouncing produces realistic dust grain porosities in protoplanetary discs

Stéphane Michoulier, Jean-François Gonzalez, Daniel J. Price

► To cite this version:

Stéphane Michoulier, Jean-François Gonzalez, Daniel J. Price. Compaction during fragmentation and bouncing produces realistic dust grain porosities in protoplanetary discs. Astronomy and Astrophysics - A&A, In press, 10.1051/0004-6361/202449719. [hal-04626481](https://hal.science/hal-04626481)

HAL Id: hal-04626481

https://hal.science/hal-04626481

Submitted on 26 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Compaction during fragmentation and bouncing
produces realistic dust grain porosities in protoplanetary
discs

Stéphane Michoulier, Jean-François Gonzalez, Daniel J. Price

► To cite this version:

Stéphane Michoulier, Jean-François Gonzalez, Daniel J. Price. Compaction during fragmentation and bouncing produces realistic dust grain porosities in protoplanetary discs. Astronomy & Astrophysics - A&A, 2024, 688, pp.A31. 10.1051/0004-6361/202449719. hal-04626481v2

HAL Id: hal-04626481

https://hal.science/hal-04626481v2

Submitted on 30 Jul 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Compaction during fragmentation and bouncing produces realistic dust grain porosities in protoplanetary discs

Stéphane Michoulier<fnref n="1" />, Jean-François Gonzalez<fnref n="1" />, and Daniel J. Price<fnref n="2" />

Received 26 February 2024; accepted 21 June 2024

## ABSTRACT

*Context.* In protoplanetary discs, micron-sized dust grows to form millimetre- to centimetre-sized pebbles but encounters several barriers during its evolution. Collisional fragmentation and radial drift impede further dust growth to planetesimal size. Fluffy grains have been hypothesised to solve these problems. While porosity leads to faster grain growth, the implied porosity values obtained from previous simulations were larger than suggested by observations.

*Aims.* In this paper, we study the influence of porosity on dust evolution, taking into account growth, bouncing, fragmentation, compaction, rotational disruption, and snow lines, in order to understand their impact on dust evolution.

*Methods.* We developed a module for porosity evolution for the 3D smoothed particle hydrodynamics code PHANTOM that accounts for dust growth and fragmentation. This mono-disperse model is integrated into both a 1D code and the 3D code to capture the overall evolution of dust and gas.

*Results.* We show that porosity helps dust growth and leads to the formation of larger solids than when considering compact grains, as predicted by previous work. Our simulations taking into account compaction during fragmentation show that large millimetre grains are still formed but are ten to 100 times more compact. Thus, millimetre sizes with typical filling factors of ~ 0.1 match the values measured on comets or via polarimetric observations of protoplanetary discs.

**Key words.** methods: numerical — planets and satellites: formation — protoplanetary discs

## 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 de-

pend 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); Blum & Wurm (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

<footnote n="1">¹ Universite Claude Bernard Lyon 1, CRAL UMR5574, ENS de Lyon, CNRS, Villeurbanne, F-69622, France
e-mail: jean-francois.gonzalez@ens-lyon.fr</footnote>

<footnote n="2">² School of Physics and Astronomy, Monash University, Vic. 3800, Australia</footnote>
# Compaction during fragmentation and bouncing produces realistic dust grain porosities in protoplanetary discs

Stéphane Michoulier<fnref n="1" />, Jean-François Gonzalez<fnref n="1" />, and Daniel J. Price<fnref n="2" />

Received 24 February 2024 / Accepted 21 June 2024

## ABSTRACT

*Context.* In protoplanetary discs, micron-sized dust grows to form millimetre- to centimetre-sized pebbles but encounters several barriers during its evolution. Collisional fragmentation and radial drift impede further dust growth to planetesimal size. Fluffy grains have been hypothesised to solve these problems. While porosity leads to faster grain growth, the implied porosity values obtained from previous simulations were larger than suggested by observations.

*Aims.* In this paper, we study the influence of porosity on dust evolution, taking into account growth, bouncing, fragmentation, compaction, rotational disruption, and snow lines, in order to understand their impact on dust evolution.

*Methods.* We developed a module for porosity evolution for the 3D smoothed particle hydrodynamics code PHANTOM that accounts for dust growth and fragmentation. This mono-disperse model is integrated into both a 1D code and the 3D code to capture the overall evolution of dust and gas.

*Results.* We show that porosity helps dust growth and leads to the formation of larger solids than when considering compact grains, as predicted by previous work. Our simulations taking into account compaction during fragmentation show that large millimetre grains are still formed but are ten to 100 times more compact. Thus, millimetre sizes with typical filling factors of ~0.1 match the values measured on comets or via polarimetric observations of protoplanetary discs.

**Key words.** methods: numerical – planets and satellites: formation – protoplanetary disks

## 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

<footnote n="1">¹ Universite Claude Bernard Lyon 1, CRAL UMR5574, ENS de Lyon, CNRS, Villeurbanne 69622, France
e-mail: jean-francois.gonzalez@ens-lyon.fr</footnote>

<footnote n="2">² School of Physics and Astronomy, Monash University, Vic. 3800, Australia</footnote>
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; Suyama et al. 2012; Okuzumi et al. 2012; Kataoka et al. 2013a; 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; Okuzumi et al. 2009; Suyama et al. 2012; Okuzumi et al. 2012; Kataoka et al. 2013a) 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; Suyama et al. 2008; Wada et al. 2009; Suyama et al. 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 $\tau_{coll}$

$$ \left(\frac{dm}{dt}\right)_{\text{grow}} \approx \frac{m}{\tau_{\text{coll}}}. \qquad (1) $$

The characteristic time can be expressed in terms of the dust number density $n_d$, collision cross-sectional area ($\sigma_d = 4\pi s^2$), and relative velocity $v_{rel}$ during collision as $\tau_{coll} = (n_d\sigma v_{rel})^{-1}$. The mass change rate can then be computed as

$$ \left(\frac{dm}{dt}\right)_{\text{grow}} = 4\pi\rho_d s^2 v_{\text{rel}}, \qquad (2) $$

where $\rho_d$ is the dust volume density, while the relative velocity $v_{rel}$ can be expressed as

$$ v_{\text{rel}} = \sqrt{2v_t \frac{\sqrt{Sc} - 1}{Sc}}. \qquad (3) $$

The turbulent velocity $v_t$ is given by

$$ v_t = \sqrt{2^{1/2}R_o \alpha c_g}, \qquad (4) $$

where $R_o$ is the Rossby number, which is a constant equal to 3.0 (Stepinski & Valageas 1997), $\alpha$ 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

$$ Sc = (1 + St) \sqrt{1 + \frac{\Delta v^2}{v_t^2}}, \qquad (5) $$

where $\Delta v$ 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

$$ St = t_s \Omega_K, \qquad (6) $$

where $t_s$ is the stopping time, that is, the time needed for a grain to reach the gas velocity, and $\Omega_K$ the Keplerian frequency. More details about the implementation and the growth model are provided in Vericel et al. (2021).
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; Schäfer 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 $\tau_{coll}$

$$ \left(\frac{dm}{dt}\right)_{\text{grow}} \approx \frac{m}{\tau_{\text{coll}}}. \qquad (1) $$

The characteristic time can be expressed in terms of the dust number density $n_d$, collision cross-sectional area ($\sigma_d = 4\pi s^2$), and relative velocity $v_{rel}$ during collision as $\tau_{coll} = (n_d\sigma v_{rel})^{-1}$. The mass change rate can then be computed as

$$ \left(\frac{dm}{dt}\right)_{\text{grow}} = 4\pi\rho_d s^2 v_{\text{rel}}, \qquad (2) $$

where $\rho_d$ is the dust volume density, while the relative velocity $v_{rel}$ can be expressed as

$$ v_{\text{rel}} = \sqrt{2v_t \frac{\sqrt{Sc} - 1}{Sc}}. \qquad (3) $$

The turbulent velocity $v_t$ is given by

$$ v_t = \sqrt{2^{1/2}R_o \alpha c_g}, \qquad (4) $$

where $R_o$ is the Rossby number, which is a constant equal to 3.0 (Stepinski & Valageas 1997), $\alpha$ 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

$$ Sc = (1 + St) \sqrt{1 + \frac{\Delta v^2}{v_t^2}}, \qquad (5) $$

where $\Delta v$ 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

$$ St = t_s \Omega_K, \qquad (6) $$

where $t_s$ is the stopping time, that is, the time needed for a grain to reach the gas velocity, and $\Omega_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), Suyama et al. (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_0$, size $a_0$, and intrinsic density $\rho_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 $\phi$ related to the porosity $p$ as the ratio between the volume occupied by the monomers $V_{\text{mat}}$ and the volume of the aggregate $V$:

$$ \phi = \frac{V_{\text{mat}}}{V} \quad \text{and} \quad p + \phi = 1. \qquad (7) $$

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

$$ m = \rho_s \phi \frac{4\pi}{3} s^3. \qquad (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 $v_{\text{rel}}$. In the referential of the centre of mass, the kinetic energy is expressed as

$$ E_{\text{kin}} = \frac{m}{4} v_{\text{rel}}^2, \qquad (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_{\text{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_s a_0 \xi_{\text{crit}}, \qquad (10) $$

where $\gamma_s$ is the surface energy (energy per unit area) of a monomer, $\xi_{\text{crit}}$ is the critical rolling displacement, and $a_0$ is the monomer size. The value of $\xi_{\text{crit}}$ is still poorly constrained. According to Dominik & Tielens (1997), the critical separation $\delta_c$ between two monomers before they separate,

$$ \delta_c = \left( \frac{27\pi^2 \gamma_s^2 a_0}{2 \varepsilon^2} \right)^{\frac{1}{3}}, \qquad (11) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Filling factor $\phi$ as a function of size *s* for different distances from the star $R_0$ for silicate grains composed of 0.2 µm monomers in the case of pure growth in our standard disc model (see Sect. 3.3).</figcaption></figure>

is of the same order of magnitude as $\xi_{\text{crit}}$ (Chokshi et al. 1993). Thus, we can rewrite Eq. (10) as

$$ E_{\text{roll}} = \left( 2916\pi^8 \frac{\gamma_s^5 a_0^4}{\varepsilon^2} \right)^{\frac{1}{3}} \approx 302 \left( \frac{\gamma_s^5 a_0^4}{\varepsilon^2} \right)^{\frac{1}{3}}, \qquad (12) $$

with $\varepsilon$ 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_{\text{kin}} < 3b E_{\text{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 $\phi$ of a grain, starting as a compact ($\phi = 1$) monomer ($s = a_0 = 0.2$ µm) as its size *s* increases, for different distances to the star.

## 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 $\phi = 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.
## 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_0$, size $a_0$, and intrinsic density $\rho_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 $\phi$ related to the porosity $p$ as the ratio between the volume occupied by the monomers $V_{\text{mat}}$ and the volume of the aggregate $V$:

$$ \phi = \frac{V_{\text{mat}}}{V} \quad \text{and} \quad p + \phi = 1. \qquad (7) $$

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

$$ m = \rho_s \phi \frac{4\pi}{3} s^3. \qquad (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 $v_{\text{rel}}$. In the referential of the centre of mass, the kinetic energy is expressed as

$$ E_{\text{kin}} = \frac{m}{4} v_{\text{rel}}^2, \qquad (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_{\text{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_s a_0 \xi_{\text{crit}}, \qquad (10) $$

where $\gamma_s$ is the surface energy (energy per unit area) of a monomer, $\xi_{\text{crit}}$ is the critical rolling displacement, and $a_0$ is the monomer size. The value of $\xi_{\text{crit}}$ is still poorly constrained. According to Dominik & Tielens (1997), the critical separation $\delta_c$ between two monomers before they separate,

$$ \delta_c = \left( \frac{27\pi^2}{2} \frac{\gamma_s a_0}{\mathcal{E}^2} \right)^{\frac{1}{3}}, \qquad (11) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Filling factor $\phi$ as a function of size *s* for different distances from the star $R_0$ for silicate grains composed of 0.2 µm monomers in the case of pure growth in our standard disc model (see Sect. 3.3).</figcaption></figure>

is of the same order of magnitude as $\xi_{\text{crit}}$ (Chokshi et al. 1993). Thus, we can rewrite Eq. (10) as

$$ E_{\text{roll}} = \left( 2916\pi^8 \frac{\gamma_s^5 a_0^4}{\mathcal{E}^2} \right)^{\frac{1}{3}} \approx 302 \left( \frac{\gamma_s^5 a_0^4}{\mathcal{E}^2} \right)^{\frac{1}{3}}, \qquad (12) $$

with $\mathcal{E}$ 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_{\text{kin}} < 3b E_{\text{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 $\phi$ of a grain, starting as a compact ($\phi = 1$) monomer ($s = a_0 = 0.2$ µm) as its size *s* increases, for different distances to the star.

## 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 $\phi = 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
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 $v_{\text{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 $v_{\text{stick}}$ for the collision of two identical grains

$$v_{\text{stick}} = 4.23 \left( \frac{\gamma_s^5 s^4}{m^3 \varepsilon^2} \right)^{1/6} . \quad (13)$$

Dominik & Tielens (1997) provided another formula for $v_{\text{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 $\varepsilon$ of water ice varies with porosity by about an order of magnitude. Its value as a function of $\phi$ is not known for silicates. Since $v_{\text{stick}} \propto \varepsilon^{-1/3}$, by using the value of $\varepsilon$ for monomers and retaining the factor of two, a reasonable compromise is achieved. Therefore, we kept

$$\begin{aligned} v_{\text{stick}} = & 8.46 \left( \frac{\gamma_s^5 s^4}{m^3 \varepsilon^2} \right)^{1/6} \\ & 2.82 \, \text{mm s}^{-1} \left( \frac{\rho_s}{1000 \, \text{kg m}^{-3}} \right)^{-1/2} \left( \frac{\gamma_s}{0.1 \, \text{J m}^{-2}} \right)^{5/6} \\ & \times \left( \frac{\varepsilon}{10 \, \text{GPa}} \right)^{-1/3} \left( \frac{s}{1 \, \text{mm}} \right)^{-5/6} \phi^{-1/2}, \end{aligned} \quad (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 $v_{\text{stick}} = 2.89 \, \text{mm s}^{-1}$ and $5.26 \, \text{mm s}^{-1}$, respectively, for millimeter-sized grains with $\phi = 0.3$. Weidling et al. (2012) argued that instead of $\gamma_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 $v_{\text{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 $v_{\text{stick}}$ of several to tens of cm s$^{-1}$. 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^3$ Pa for $\phi = 0.2$ and $\ge 10^4$ Pa for $\phi > 0.3$ (see Eq. 31 in Seizinger et al. 2012), which increases $v_{\text{stick}}$ by two orders of magnitude. We thus also considered values of $v_{\text{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 $v_{\text{yield}}$. Thornton & Ning (1998) expressed $v_{\text{yield}}$ as a function of the contact radius $a_{\text{yield}}$ when yield occurs or of the limiting contact pressure $p_{\text{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$_2$ 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 $v_{\text{yield}} = 10v_{\text{stick}}$. More recently, Kimura et al.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Coefficient of restitution $e$ when $v_{\text{rel}} \ge v_{\text{yield}}$, given by Eq. (19).</figcaption></figure>

(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

$$v_{\text{yield}} = 10v_{\text{stick}} = 84.6 \left( \frac{\gamma_s^5 s^4}{m^3 \varepsilon^2} \right)^{1/6} . \quad (15)$$

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

$$e \equiv \frac{v_{\text{rel,f}}}{v_{\text{rel}}}, \quad (16)$$

where $v_{\text{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^2)E_{\text{kin}}$, and the remaining part serves to separate the grains with the remaining kinetic energy $e^2E_{\text{kin}}$. Depending on the relative velocity, the coefficient of restitution can then be expressed (Thornton & Ning 1998; Garcia 2018) as follows:

$$v_{\text{rel}} \le v_{\text{stick}} \qquad e = 0 \quad (17)$$

$$v_{\text{stick}} < v_{\text{rel}} \le v_{\text{yield}} \qquad e = 1 \quad (18)$$

$$v_{\text{rel}} \ge v_{\text{yield}} \qquad e = \left[ 1.2 \sqrt{3} \left( 1 - \frac{1}{6} \left( \frac{v_{\text{yield}}}{v_{\text{rel}}} \right)^2 \right) \times \left( 1 + 2 \sqrt{1.2 \left( \frac{v_{\text{rel}}}{v_{\text{yield}}} \right)^2 - 0.2} \right)^{-1/2} - \left( \frac{v_{\text{stick}}}{v_{\text{rel}}} \right)^2 \right]^{1/2} . \quad (19)$$

When $v_{\text{rel}} \gg v_{\text{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 $v_{\text{end}}$, the resulting velocity from the remaining kinetic energy $e^2E_{\text{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 $v_{\text{end}} = 0.03v_{\text{rel}}$. Solving Eq. (19) for $e = 0.03$ with $v_{\text{yield}} = 10v_{\text{stick}}$ allowed us to express $v_{\text{end}}$ as a function of $v_{\text{stick}}$ so that

$$v_{\text{end}} \approx 2 \times 10^7 v_{\text{stick}} \approx 2 \times 10^8 \left( \frac{\gamma_s^5 s^4}{m^3 \varepsilon^2} \right)^{1/6} . \quad (20)$$

Typical values of $v_{\text{end}}$ are on the order of 10 – 100 km s$^{-1}$, much larger than fragmentation thresholds.
wave (Thornton & Ning 1998). Several velocities can be distinguished. The first is the sticking velocity $v_{\text{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 $v_{\text{stick}}$ for the collision of two identical grains

$$v_{\text{stick}} = 4.23 \left( \frac{\gamma_s^5 s^4}{m^3 \mathcal{E}^2} \right)^{1/6} . \qquad (13)$$

Dominik & Tielens (1997) provided another formula for $v_{\text{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 $\mathcal{E}$ of water ice varies with porosity by about an order of magnitude. Its value as a function of $\phi$ is not known for silicates. Since $v_{\text{stick}} \propto \mathcal{E}^{-1/3}$, by using the value of $\mathcal{E}$ for monomers and retaining the factor of two, a reasonable compromise is achieved. Therefore, we kept

$$\begin{aligned} v_{\text{stick}} &= 8.46 \left( \frac{\gamma_s^5 s^4}{m^3 \mathcal{E}^2} \right)^{1/6} \\ &= 2.82 \, \text{mm s}^{-1} \left( \frac{\rho_s}{1000 \, \text{kg m}^{-3}} \right)^{-1/2} \left( \frac{\gamma_s}{0.1 \, \text{J m}^{-2}} \right)^{5/6} \\ &\quad \times \left( \frac{\mathcal{E}}{10 \, \text{GPa}} \right)^{-1/3} \left( \frac{s}{1 \, \text{mm}} \right)^{-5/6} \phi^{-1/2}, \end{aligned} \qquad (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 $v_{\text{stick}} = 2.89 \, \text{mm s}^{-1}$ and $5.26 \, \text{mm s}^{-1}$, respectively, for millimeter-sized grains with $\phi = 0.3$. Weidling et al. (2012) argued that instead of $\gamma_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 $v_{\text{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 $v_{\text{stick}}$ of several tens of cm s$^{-1}$. 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 $\sim 10^3$ Pa for $\phi = 0.2$ and $\ge 10^4$ Pa for $\phi > 0.3$ (see Eq. (31) in Seizinger et al. 2012), which increases $v_{\text{stick}}$ by two orders of magnitude. We thus also considered values of $v_{\text{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 $v_{\text{yield}}$. Thornton & Ning (1998) expressed $v_{\text{yield}}$ as a function of the contact radius $a_{\text{yield}}$ when yield occurs or of the limiting contact pressure $p_{\text{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$_2$ 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 $v_{\text{yield}} = 10v_{\text{stick}}$. More recently, Kimura et al. (2020) demonstrated that silicates are more prone to sticking

<figure><img src="image_3.png" /><figcaption>Fig. 2. Coefficient of restitution $e$ when $v_{\text{rel}} \ge v_{\text{yield}}$, given by Eq. (19).</figcaption></figure>

together with a higher surface energy. We thus kept a factor of ten between the two velocities for silicates as well and took

$$v_{\text{yield}} = 10v_{\text{stick}} = 84.6 \left( \frac{\gamma_s^5 s^4}{m^3 \mathcal{E}^2} \right)^{1/6} . \qquad (15)$$

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

$$e \equiv \frac{v_{\text{rel,f}}}{v_{\text{rel}}}, \qquad (16)$$

where $v_{\text{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^2)E_{\text{kin}}$, and the remaining part serves to separate the grains with the remaining kinetic energy $e^2E_{\text{kin}}$. Depending on the relative velocity, the coefficient of restitution can then be expressed (Thornton & Ning 1998; Garcia 2018) as follows:

$$\begin{aligned} v_{\text{rel}} \le v_{\text{stick}} &\qquad e = 0 \\ v_{\text{stick}} < v_{\text{rel}} \le v_{\text{yield}} &\qquad e = 1 \end{aligned} \qquad (17) \qquad (18)$$

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

When $v_{\text{rel}} \gg v_{\text{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 $v_{\text{end}}$, the resulting velocity from the remaining kinetic energy $e^2E_{\text{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 $v_{\text{end}} = 0.03v_{\text{rel}}$. Solving Eq. (19) for $e = 0.03$ with $v_{\text{yield}} = 10v_{\text{stick}}$ allowed us to express $v_{\text{end}}$ as a function of $v_{\text{stick}}$ so that

$$v_{\text{end}} \approx 2 \times 10^7 v_{\text{stick}} \approx 2 \times 10^8 \left( \frac{\gamma_s^5 s^4}{m^3 \mathcal{E}^2} \right)^{1/6} . \qquad (20)$$

Typical values of $v_{\text{end}}$ are on the order of 10–100 km s$^{-1}$, much larger than fragmentation thresholds.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Model comparison for compression pressure. Tatsuuma et al. (2023) connect the regime with small $\phi$ given by Kataoka et al. (2013b) to the large $\phi$ with the model of Güttler et al. (2010) and Seizinger et al. (2012).</figcaption></figure>

## Sticking probability and filling factor

Empirically, Weidling et al. (2012) determined the probability that a grain sticks, or not. As previously mentioned, when $v_{\text{stick}} < v_{\text{rel}} < v_{\text{end}}$, an aggregate can bounce instead of stick. Naturally, the closer $v_{\text{rel}}$ is to $v_{\text{stick}}$, the higher the chances of sticking, and conversely for $v_{\text{end}}$. Thus, Garcia (2018), inspired by Weidling et al. (2012), defined the probability $P$ of sticking for a grain based on its velocity according to

$$ P = \begin{cases} 1, & v_{\text{rel}} < v_{\text{stick}}, \\ \frac{\log(v_{\text{rel}}) - \log(v_{\text{end}})}{\log(v_{\text{stick}}) - \log(v_{\text{end}})}, & v_{\text{stick}} \le v_{\text{rel}} < v_{\text{end}}, \\ 0, & v_{\text{rel}} \ge v_{\text{end}}, \end{cases} \quad (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{dm}{dt}\right)_{\text{bounce}} = P\left(\frac{dm}{dt}\right)_{\text{grow}}. \quad (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$

$$ \phi_{\text{bounce}} \approx \phi_i \left(1 - \frac{\Delta V}{V_i}\right)^{-n}, \quad (23) $$

where $n = t/\tau_{\text{coll}}$ and $\Delta V$ is the volume change during compression due to plastic deformation. Shimaki & Arakawa (2012a) provided an expression to compute $\Delta V$, which is proportional to the kinetic energy,

$$ \Delta V = \frac{(1 - e^2)E_{\text{kin}}}{2\Upsilon_d}, \quad (24) $$

where $\Upsilon_d = \Upsilon_0\phi^n$ is the dynamic compression resistance, and $\Upsilon_0$ is the reference dynamic compression resistance. Mellor

<figure><img src="image_5.png" /><figcaption>Fig. 4. Growth and bouncing of a water ice grain with $a_0 = 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 $\Upsilon_d$ given by Shimaki & Arakawa (2012a). The dashed and dot-dashed lines show the impact of multiplying $v_{\text{stick}}$ 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).</figcaption></figure>

(1974) gave $\Upsilon_0 = 9.8$ MPa and $N = 4$ for water ice. However, $\Upsilon$, 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 $\Upsilon$ 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. (2013b) 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 $\phi > 0.2$. Very recently, Tatsuuma et al. (2023) have developed a formula for compression pressure that depends only on $E_{\text{roll}}$ and $\phi$, which is valid for both small and large $\phi$, shown in Fig. 3,

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

The value $\phi = 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 $\phi = 0.3$. Figure 4 allows for the comparison of the bounce model from Garcia (2018) with the equation involving $\Upsilon$, which is valid only for water ice (or a mixture of ice and silicates), and our model where $P_{\text{bounce}}$ replaces $\Upsilon$ in Eq. (24). The difference between both models is observed when the grain is highly compacted, with a filling factor $\phi > 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 $\phi$ approaches 1, the value of the compression pressure changes only slightly. The dashed and dot-dashed curves
<figure><img src="image_4.png" /><figcaption>Fig. 3. Model comparison for compression pressure. Tatsuuma et al. (2023) connect the regime with small $\phi$ given by Kataoka et al. (2013) to the large $\phi$ with the model of Güttler et al. (2010) and Seizinger et al. (2012).</figcaption></figure>

### Sticking probability and filling factor

Empirically, Weidling et al. (2012) determined the probability that a grain sticks, or not. As previously mentioned, when $v_{\text{stick}} < v_{\text{rel}} < v_{\text{end}}$, an aggregate can bounce instead of stick. Naturally, the closer $v_{\text{rel}}$ is to $v_{\text{stick}}$, the higher the chances of sticking, and conversely for $v_{\text{end}}$. Thus, Garcia (2018), inspired by Weidling et al. (2012), defined the probability $\mathcal{P}$ of sticking for a grain based on its velocity according to

$$ \mathcal{P} = \begin{cases} 1, & v_{\text{rel}} < v_{\text{stick}}, \\ \frac{\log(v_{\text{rel}}) - \log(v_{\text{end}})}{\log(v_{\text{stick}}) - \log(v_{\text{end}})}, & v_{\text{stick}} \le v_{\text{rel}} < v_{\text{end}}, \\ 0, & v_{\text{rel}} \ge v_{\text{end}}, \end{cases} \quad (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{dm}{dt}\right)_{\text{bounce}} = \mathcal{P}\left(\frac{dm}{dt}\right)_{\text{grow}}. \quad (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$

$$ \phi_{\text{bounce}} \approx \phi_i \left(1 - \frac{\Delta V}{V_i}\right)^{-n}, \quad (23) $$

where $n = t/\tau_{\text{coll}}$ and $\Delta V$ is the volume change during compression due to plastic deformation. Shimaki & Arakawa (2012a) provided an expression to compute $\Delta V$, which is proportional to the kinetic energy,

$$ \Delta V = \frac{(1 - e^2)E_{\text{kin}}}{2\Upsilon_d}, \quad (24) $$

where $\Upsilon_d = \Upsilon_0\phi^N$ is the dynamic compression resistance, and $\Upsilon_0$ is the reference dynamic compression resistance.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Growth and bouncing of a water ice grain with $a_0 = 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 $\Upsilon_d$ given by Shimaki & Arakawa (2012a). The dashed and dot-dashed lines show the impact of multiplying $v_{\text{stick}}$ 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).</figcaption></figure>

Mellor (1974) gave $\Upsilon_0 = 9.8$ MPa and $N = 4$ for water ice. However, $\Upsilon$, 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 $\Upsilon$ 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 $\phi > 0.2$. Very recently, Tatsuuma et al. (2023) have developed a formula for compression pressure that depends only on $E_{\text{roll}}$ and $\phi$, which is valid for both small and large $\phi$, shown in Fig. 3,

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

The value $\phi = 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 $\phi = 0.3$. Figure 4 allows for the comparison of the bounce model from Garcia (2018) with the equation involving $\Upsilon$, which is valid only for water ice (or a mixture of ice and silicates), and our model where $P_{\text{bounce}}$ replaces $\Upsilon$ in Eq. (24). The difference between both models is observed when the grain is highly compacted, with a filling factor $\phi > 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 $\phi$ approaches 1, the value of the compression pressure changes only slightly. The dashed and dot-dashed curves
show that considering values of $v_{\text{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 figure 7.

We now relax our working assumption and return to the situation where bouncing is only effective for $\phi \ge 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 $v_{\text{stick}}$ (see Sect. 2.5).

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

$$ \phi_{\text{coll & bounce}} = \begin{cases} \phi_{\text{coll}}, & v_{\text{rel}} < v_{\text{yield}}, \\ \phi_{\text{coll}}\mathcal{P} + (1-\mathcal{P})\phi_{\text{bounce}}, & v_{\text{yield}} \le v_{\text{rel}} < v_{\text{end}}, \\ \phi_{\text{bounce}}, & v_{\text{rel}} \ge v_{\text{end}}, \end{cases} \quad (26) $$

where $\phi_{\text{coll}}$ is the filling factor resulting from collisions in one of the regimes described in Appendix A.2, and $\phi_{\text{bounce}}$ is computed using Eq. (23). The final filling factor is the largest value between $\phi_{\text{gas}}$, $\phi_{\text{grav}}$, and $\phi_{\text{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.

## 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 $v_{\text{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 $v_{\text{rel}}$ with respect to $v_{\text{frag}}$

$$ \left(\frac{dm}{dt}\right)_{\text{frag}} = -4\pi \frac{v_{\text{rel}}^3}{v_{\text{rel}}^2 + v_{\text{frag}}^2} \rho_d s^2. \quad (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 $\tau_{\text{coll}}$ ($v_{\text{rel}} = v_{\text{frag}}$) or more ($v_{\text{rel}} > v_{\text{frag}}$). In the case where $v_{\text{rel}} \gg v_{\text{frag}}$, we recover the scenario described by Gonzalez et al. (2015), where the entire grain fragments after a collision time, independently of $v_{\text{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. 2013b; 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 $v_{\text{frag}}$, an estimate of the energy required to break the bond between two monomers $E_{\text{break}}$ is given by Dominik & Tielens (1997),

$$ E_{\text{break}} = C_{\text{break}} F_c \delta_c, \quad (28) $$

where $C_{\text{break}}$ is a constant that they take equal to 1.8. The term $F_c$ is the critical pulling force between two monomers,

$$ F_c = 3\pi\gamma_s a_0. \quad (29) $$

We can then derive the expression for $E_{\text{break}}$:

$$ E_{\text{break}} \approx C_{\text{break}} 48 \left( \frac{\gamma_s^5 a_0^4}{\varepsilon^2} \right)^{\frac{1}{3}}. \quad (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}}, \quad (31) $$

where $N_{\text{tot}}$ is the total number of monomers, and $\kappa$ is a numerical factor representing the number of bonds and depending on the species and porosity. For erosion, $\kappa = 1 - 2$, and for fragmentation, $\kappa = 3 - 10$ (Blum & Wurm 2000; Wada et al. 2007). We can then define $v_{\text{frag}}$ in terms of $E_{\text{frag}}$ as

$$ v_{\text{frag}} = 2 \sqrt{\frac{E_{\text{frag}}}{m}} = 2 \sqrt{\frac{48 \kappa N_{\text{tot}} C_{\text{break}}}{m} \left( \frac{\gamma_s^5 a_0^4}{\varepsilon^2} \right)^{\frac{1}{6}}}. \quad (32) $$

These equations depend on $\gamma_s$ and $a_0$. An error in estimating these parameters can thus lead to changes in the values of $E_{\text{frag}}$ and $v_{\text{frag}}$.

Research on the value of $v_{\text{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 $v_{\text{frag, Si}} \sim 1 \text{ m s}^{-1}$ for silicates. Wada et al. (2009, 2013) found a similar value for silicates, around $v_{\text{frag, Si}} \sim 5 \text{ m s}^{-1}$, as well as that for water, $v_{\text{frag, H}_2\text{O}} \sim 60 - 70 \text{ m s}^{-1}$. However, the method used by Yamamoto et al. (2014) yielded a different value for water, $v_{\text{frag, H}_2\text{O}} = 56 \text{ m 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 $\gamma_s = 0.3 \text{ J/m}^2$, while Kimura et al. (2020) found $\gamma_s = 0.15 \text{ J/m}^2$. Silicates are stronger than previously thought due to experimental issues with the original estimations of $\gamma_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 $v_{\text{frag, H}_2\text{O}} \sim 15 \text{ m 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 $v_{\text{frag, H}_2\text{O}} \approx 15 \text{ m s}^{-1}$ for water.

However, pure water ice grains are unlikely. Most often, grains are composed of solid materials like silicates, covered by
show that considering values of $v_{\text{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 $\phi \ge 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 $v_{\text{stick}}$ (see Sect. 2.5).

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

$$ \phi_{\text{coll & bounce}} = \begin{cases} \phi_{\text{coll}}, & v_{\text{rel}} < v_{\text{yield}}, \\ \phi_{\text{coll}}\mathcal{P} + (1 - \mathcal{P})\phi_{\text{bounce}}, & v_{\text{yield}} \le v_{\text{rel}} < v_{\text{end}}, \\ \phi_{\text{bounce}}, & v_{\text{rel}} \ge v_{\text{end}}, \end{cases} \qquad (26) $$

where $\phi_{\text{coll}}$ is the filling factor resulting from collisions in one of the regimes described in Appendix A.2, and $\phi_{\text{bounce}}$ is computed using Eq. (23). The final filling factor is the largest value between $\phi_{\text{gas}}$, $\phi_{\text{grav}}$, and $\phi_{\text{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.

## 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 $v_{\text{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 $v_{\text{rel}}$ with respect to $v_{\text{frag}}$

$$ \left(\frac{dm}{dt}\right)_{\text{frag}} = -4\pi \frac{v_{\text{rel}}^3}{v_{\text{rel}}^2 + v_{\text{frag}}^2} \rho_{\text{d}} s^2. \qquad (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 $\tau_{\text{coll}}$ ($v_{\text{rel}} = v_{\text{frag}}$) or more ($v_{\text{rel}} > v_{\text{frag}}$). In the case where $v_{\text{rel}} \gg v_{\text{frag}}$, we recover the scenario described by Gonzalez et al. (2015), where the entire grain fragments after a collision time, independently of $v_{\text{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 $v_{\text{frag}}$, an estimate of the energy required to break the bond between two monomers $E_{\text{break}}$ is given by Dominik & Tielens (1997),

$$ E_{\text{break}} = C_{\text{break}} F_{c} \delta_{c}, \qquad (28) $$

where $C_{\text{break}}$ is a constant that they take equal to 1.8. The term $F_{c}$ is the critical pulling force between two monomers,

$$ F_{c} = 3\pi\gamma_{s}a_{0}. \qquad (29) $$

We can then derive the expression for $E_{\text{break}}$:

$$ E_{\text{break}} \approx C_{\text{break}} 48 \left( \frac{\gamma_{s}^{5} a_{0}^{4}}{\varepsilon^{2}} \right)^{\frac{1}{3}}. \qquad (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}}, \qquad (31) $$

where $N_{\text{tot}}$ is the total number of monomers, and $\kappa$ is a numerical factor representing the number of bonds and depending on the species and porosity. For erosion, $\kappa = 1-2$, and for fragmentation, $\kappa = 3-10$ (Blum & Wurm 2000; Wada et al. 2007). We can then define $v_{\text{frag}}$ in terms of $E_{\text{frag}}$ as

$$ v_{\text{frag}} = 2 \sqrt{\frac{E_{\text{frag}}}{m}} = 2 \sqrt{\frac{48 \kappa N_{\text{tot}} C_{\text{break}}}{m}} \left( \frac{\gamma_{s}^{5} a_{0}^{4}}{\varepsilon^{2}} \right)^{\frac{1}{6}}. \qquad (32) $$

These equations depend on $\gamma_s$ and $a_0$. An error in estimating these parameters can thus lead to changes in the values of $E_{\text{frag}}$ and $v_{\text{frag}}$.

Research on the value of $v_{\text{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 $v_{\text{frag, Si}} \sim 1 \text{ m s}^{-1}$ for silicates. Wada et al. (2009, 2013) found a similar value for silicates, around $v_{\text{frag, Si}} \sim 5 \text{ m s}^{-1}$, as well as that for water, $v_{\text{frag, H}_2\text{O}} \sim 60-70 \text{ m s}^{-1}$. However, the method used by Yamamoto et al. (2014) yielded a different value for water, $v_{\text{frag, H}_2\text{O}} = 56 \text{ m 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 $\gamma_s = 0.3 \text{ J m}^{-2}$, while Kimura et al. (2020) found $\gamma_s = 0.15 \text{ J m}^{-2}$. Silicates are stronger than previously thought due to experimental issues with the original estimations of $\gamma_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 $v_{\text{frag, H}_2\text{O}} \sim 15 \text{ m 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$^{-1}$). For the purpose of comparison with certain previous studies (Garcia & Gonzalez 2020; Vericel et al. 2021), we adopted the value of $v_{\text{frag, H}_2\text{O}} \approx 15 \text{ m s}^{-1}$ for water.

However, pure water ice grains are unlikely. Most often, grains are composed of solid materials like silicates, covered by
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Filling factor $\phi_f/\phi_i$ as a function of $v_{rel}$. **Left:** Experimental values for different initial filling factors $\phi_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.</figcaption></figure>

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: $v_{frag} = 5, 10, \text{ or } 20 \text{ m s}^{-1}$ for more fragile materials like water or CO, and $v_{frag} = 40 \text{ m s}^{-1}$ 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.

## 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; Tazaki et al. 2019; Kataoka et al. 2019; Tazaki & Dominik 2022). Since neither growth, static gas compaction, nor bouncing can account for such high factors ($\phi > 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_{comp} = E_{kin} - E_{frag} = E_{kin} - \frac{(2m_i - m_f)\kappa E_{break}}{m_0} \quad (33)$$

Here, $(2m_i - m_f)/m_0$ 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 $\kappa = 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 $\kappa = 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 $\Delta V$, which is given by

$$\Delta V = \frac{E_{kin} - (2m_i - m_f)\kappa E_{break} m_0^{-1}}{2 \left(\frac{E_{roll}}{a_0^3}\right) \left(\frac{1}{\phi} - \frac{1}{0.74}\right)^{-3}} \quad (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 $\phi_f/\phi_i$ of the final filling factor to the initial one as a function of the relative velocity $v_{rel}$ for three values of $\phi_i$. It should be noted that Gunkelmann et al. (2016) defined $v_{frag}$ as the point when a grain starts to lose a monomer, which is equal to $0.17 \text{ m s}^{-1}$, and differs from our definition as the point where the grain loses half of its mass. We deduced that the ratio $\phi_f/\phi_i$ varies as $\exp(1 - (v_{rel}/v_{frag})^2)$. 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 $v_{rel}/v_{frag}$ exceeds 1.15 to 1.2. However, grains never
<figure><img src="image_6.png" /><figcaption>Fig. 5. Filling factor $\phi_f/\phi_i$ as a function of $v_{rel}$. Left: experimental values for different initial filling factors $\phi_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.</figcaption></figure>

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: $v_{frag} = 5, 10, \text{ or } 20 \text{ m s}^{-1}$ for more fragile materials like water or CO, and $v_{frag} = 40 \text{ m s}^{-1}$ 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.

## 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 ($\phi > 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_{comp} = E_{kin} - E_{frag} = E_{kin} - \frac{(2m_i - m_f)}{m_0} \kappa E_{break}. \quad (33)$$

Here, $(2m_i - m_f)/m_0$ 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 $\kappa = 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 $\kappa = 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 $\Delta V$, which is given by

$$\Delta V = \frac{E_{kin} - (2m_i - m_f)\kappa E_{break}m_0^{-1}}{2\left(\frac{E_{roll}}{a_0^3}\right)\left(\frac{1}{\phi} - \frac{1}{0.74}\right)^{-3}}. \quad (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 $\phi_f/\phi_i$ of the final filling factor to the initial one as a function of the relative velocity $v_{rel}$ for three values of $\phi_i$. It should be noted that Gunkelmann et al. (2016) defined $v_{frag}$ as the point when a grain starts to lose a monomer, which is equal to 0.17 m s$^{-1}$, and differs from our definition as the point where the grain loses half of its mass. We deduced that the ratio $\phi_f/\phi_i$ varies as $\exp(1 - (v_{rel}/v_{frag})^2)$. 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
<figure><img src="image_7.png" /><figcaption>Fig. 6. Filling factor $\phi$ 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 $R_0 = 10$ and 100 au. The dash-dotted lines show the impact of multiplying $v_{\text{stick}}$ from Eq. (14) by 100.</figcaption></figure>

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 $v_{\text{rel}}/v_{\text{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 $\phi_0$ varies after each collision. The filling factor after fragmentation is computed in the same way as bouncing, using Eq. (23).

$$ \phi_{\text{frag-comp}} \approx \phi_i \left( 1 - \exp \left( 1 - \left( \frac{v_{\text{rel}}}{v_{\text{frag}}} \right)^2 \frac{\Delta V}{V_i} \right) \right)^{-n}, \quad (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. $\phi$ 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 $\phi$ 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 $v_{\text{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.

## 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 $\phi_{\text{coll}}$, $\phi_{\text{gas}}$, $\phi_{\text{grav}}$

Compute $\phi_{\text{min}} = \max(\phi_{\text{coll}}, \phi_{\text{gas}}, \phi_{\text{grav}})$

if $v_{\text{rel}} < v_{\text{frag}}$ then

Compute $\phi_{\text{grow}}$

if (Grains can bounce) then

Compute $\phi_{\text{coll & bounce}}$

Compute $\phi_f = \max(\phi_{\text{coll & bounce}}, \phi_{\text{min}})$

else

Compute $\phi_f = \max(\phi_{\text{grow}}, \phi_{\text{min}})$

end if

else

if (Fragmentation with compaction) then

Compute $\phi_{\text{frag-comp}}$

Compute $\phi_f = \max(\phi_{\text{frag-comp}}, \phi_{\text{min}})$

else

Compute $\phi_f = \max(\phi_i, \phi_{\text{min}})$

end if

end if

$\phi_{\text{final}} = \min(\phi_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 of porous 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,

$$ \Sigma_g(r) = \Sigma_g, \quad 0(r/R_0)^{-p} \quad (36) $$

and

$$ T_g(r) = T_g, \quad 0(r/R_0)^{-q}, \quad (37) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Filling factor $\phi$ 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 $R_0 = 10$ and 100 au. The dash-dotted lines show the impact of multiplying $v_{\text{stick}}$ from Eq. (14) by 100.</figcaption></figure>

occurs when $v_{\text{rel}}/v_{\text{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 $v_{\text{rel}}/v_{\text{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 $\phi_0$ varies after each collision. The filling factor after fragmentation is computed in the same way as bouncing, using Eq. (23).

$$ \phi_{\text{frag-comp}} \approx \phi_i \left(1 - \exp\left(1 - \left(\frac{v_{\text{rel}}}{v_{\text{frag}}}\right)^2\right) \frac{\Delta V}{V_i}\right)^{-n}, \quad (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. $\phi$ 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 $\phi$ 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 $v_{\text{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.

## 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 $\phi_{\text{coll}}$, $\phi_{\text{gas}}$, $\phi_{\text{grav}}$
Compute $\phi_{\text{min}} = \max(\phi_{\text{coll}}, \phi_{\text{gas}}, \phi_{\text{grav}})$
**if** $v_{\text{rel}} < v_{\text{frag}}$ **then**
  Compute $\phi_{\text{grow}}$
  **if** (Grains can bounce) **then**
    Compute $\phi_{\text{coll & bounce}}$
    Compute $\phi_{\text{f}} = \max(\phi_{\text{coll & bounce}}, \phi_{\text{min}})$
  **else**
    Compute $\phi_{\text{f}} = \max(\phi_{\text{grow}}, \phi_{\text{min}})$
  **end if**
**else**
  **if** (Fragmentation with compaction) **then**
    Compute $\phi_{\text{frag-comp}}$
    Compute $\phi_{\text{f}} = \max(\phi_{\text{frag-comp}}, \phi_{\text{min}})$
  **else**
    Compute $\phi_{\text{f}} = \max(\phi_{\text{i}}, \phi_{\text{min}})$
  **end if**
**end if**
$\phi_{\text{final}} = \min(\phi_{\text{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 of porous 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
<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

where $R_0$ 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 $\epsilon$ 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.

## 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.2–2.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:

$$ \Sigma(r) = \Sigma_0 \left(1 - \sqrt{\frac{R_{in}}{r}}\right) \left(\frac{r}{R_c}\right)^{-p} \exp \left[-(r/R_c)^{2-p}\right], \quad (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 $\epsilon_0$. 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,z)$, 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(r/R_0)^{-2}, \quad (39) $$

and in the right panel, the vertical dependence of the final size is $\propto \exp(-z^2/H^2)$, 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).
<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

p and q, respectively,

$$ \Sigma_g(r) = \Sigma_g, \quad 0(r/R_0)^{-p} \qquad (36) $$

and

$$ T_g(r) = T_g, \quad 0(r/R_0)^{-q}, \qquad (37) $$

where $R_0$ 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 $\epsilon$ 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.

## 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.2–2.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:

$$ \Sigma(r) = \Sigma_0 \left(1 - \sqrt{\frac{R_{in}}{r}}\right) \left(\frac{r}{R_c}\right)^{-p} \exp \left[-(r/R_c)^{2-p}\right], \qquad (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 $\epsilon_0$. 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, z), 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(r/R_0)^{-2}, \qquad (39) $$

and in the right panel, the vertical dependence of the final size is $\propto \exp(-z^2/H^2)$, 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.
<table><caption>**Table 1.** Simulations carried out with both PAMDEAS and PHANTOM.</caption><thead><tr><th>Name</th><th>species</th><th>bouncing & compaction</th><th>disruption</th><th>a<sub>0</sub> (μm)</th><th>v<sub>frag</sub> (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>GF-Si-comp-Vf10</td><td>Si</td><td>no</td><td>no</td><td>comp</td><td>10</td></tr><tr><td>GF-Si-comp-Vf20</td><td>Si</td><td>no</td><td>no</td><td>comp</td><td>20</td></tr><tr><td>GF-Si-comp-Vf40</td><td>Si</td><td>no</td><td>no</td><td>comp</td><td>40</td></tr><tr><td>GF-Si-a02-Vf10</td><td>Si</td><td>no</td><td>no</td><td>0.2</td><td>10</td></tr><tr><td>GF-Si-a02-Vf20</td><td>Si</td><td>no</td><td>no</td><td>0.2</td><td>20</td></tr><tr><td>GF-Si-a02-Vf40</td><td>Si</td><td>no</td><td>no</td><td>0.2</td><td>40</td></tr><tr><td>GBFc-Si-a02-Vf10</td><td>Si</td><td>yes</td><td>no</td><td>0.2</td><td>10</td></tr><tr><td>GBFc-Si-a02-Vf20</td><td>Si</td><td>yes</td><td>no</td><td>0.2</td><td>20</td></tr><tr><td>GBFc-Si-a02-Vf40</td><td>Si</td><td>yes</td><td>no</td><td>0.2</td><td>40</td></tr><tr><td>GBFcS-Si-a02-Vf205*</td><td>Si</td><td>yes</td><td>no</td><td>0.2</td><td>20</td></tr><tr><td>GF-H2O-comp-Vf15</td><td>H2O</td><td>no</td><td>no</td><td>comp</td><td>15</td></tr><tr><td>GF-H2O-a02-Vf15</td><td>H2O</td><td>no</td><td>no</td><td>0.2</td><td>15</td></tr><tr><td>GBFc-H2O-a02-Vf15</td><td>H2O</td><td>yes</td><td>no</td><td>0.2</td><td>15</td></tr><tr><td>GBFcD-Si-a02-Vf10</td><td>Si</td><td>yes</td><td>yes</td><td>0.2</td><td>10</td></tr><tr><td>GBFcD-Si-a02-Vf20</td><td>Si</td><td>yes</td><td>yes</td><td>0.2</td><td>20</td></tr><tr><td>GBFcD-H20-a02-Vf15</td><td>H2O</td><td>yes</td><td>yes</td><td>0.2</td><td>15</td></tr></tbody></table>

**Notes.** The ‘G’ stands for growth, ‘F’ for fragmentation, ‘B’ for bounce, ‘c’ for compaction, and ‘D’ for rotational disruption. The term ‘comp’ in the column for the monomer size $a_0$ denotes compact grains, with porosity evolution not taken into account. ‘D’ simulations were run only with PHANTOM. (*) This simulation includes a snow line, indicated by ‘S’, with an inner threshold of 20 m s$^{-1}$ and an outer threshold of 5 m s$^{-1}$ to model the CO snow line ($T = 20$ K, approximately ~100 au).

### 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 = 0.0895$ at $R_0 = 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 $\alpha = 5 \times 10^{-3}$ (by setting $\alpha_{AV} = 0.1658$) and the initial dust-to-gas ratio $\epsilon_0$ 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 = 0.2$ µm, in agreement with recent observations. Indeed, Tazaki & Domnik (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 $\phi \le 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 $\rho_s = 2700$ kg m$^{-3}$ and a surface energy $\gamma_s = 0.2$ J m$^{-2}$, in agreement with Yamamoto et al. (2014), who estimate $\gamma_s = 0.3$ J m$^{-2}$. This value is of the same order of magnitude as $\gamma_s = 0.15$ J m$^{-2}$, found by Kimura et al. (2015) and Kimura et al. (2020), through experimental measurements with sicastar® aggregates (micromod Partikeltechnologie GmbH). The Young’s modulus was $\mathcal{E} = 72$ GPa (Yamamoto et al. 2014), which implies, assuming that the critical separation between two monomers $\delta_c$ before their separation is of the same order of magnitude as $\xi_{crit}$ (Chokshi et al. 1993), a value of $\xi_{crit} \approx 6$ Å.

For pure water ice, the intrinsic density of monomers was $\rho_s = 1000$ kg m$^{-3}$, with a surface energy of $\gamma_s = 0.1$ J m$^{-2}$. The critical separation is still debated, as are many other properties of dust in protoplanetary discs. We chose a Young’s modulus of $\mathcal{E} = 9.4$ GPa, following Yamamoto et al. (2014), which results in a value of $\xi_{crit} \sim 10$ Å, close to $\xi_{crit} \sim 8$ Å 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 $v_{frag} = 20$ m s$^{-1}$ 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. $\phi = 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 \sim 5 \times 10^{-3}$, and they continue to grow until they reach the fragmentation threshold (horizontal plateau), maintaining an equilibrium between growth and fragmentation. While
<table><caption>**Table 1.** Simulations carried out with both PAMDEAS and PHANTOM.</caption><thead><tr><th>Name</th><th>Species</th><th>Bouncing & compaction</th><th>Disruption</th><th>a<sub>0</sub> (μm)</th><th>v<sub>frag</sub> (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>GF-Si-comp-Vf10</td><td>Si</td><td>No</td><td>No</td><td>Comp</td><td>10</td></tr><tr><td>GF-Si-comp-Vf20</td><td>Si</td><td>No</td><td>No</td><td>Comp</td><td>20</td></tr><tr><td>GF-Si-comp-Vf40</td><td>Si</td><td>No</td><td>No</td><td>Comp</td><td>40</td></tr><tr><td>GF-Si-a02-Vf10</td><td>Si</td><td>No</td><td>No</td><td>0.2</td><td>10</td></tr><tr><td>GF-Si-a02-Vf20</td><td>Si</td><td>No</td><td>No</td><td>0.2</td><td>20</td></tr><tr><td>GF-Si-a02-Vf40</td><td>Si</td><td>No</td><td>No</td><td>0.2</td><td>40</td></tr><tr><td>GBFc-Si-a02-Vf10</td><td>Si</td><td>Yes</td><td>No</td><td>0.2</td><td>10</td></tr><tr><td>GBFc-Si-a02-Vf20</td><td>Si</td><td>Yes</td><td>No</td><td>0.2</td><td>20</td></tr><tr><td>GBFc-Si-a02-Vf40</td><td>Si</td><td>Yes</td><td>No</td><td>0.2</td><td>40</td></tr><tr><td>GBFcS-Si-a02-Vf205<sup>(*)</sup></td><td>Si</td><td>Yes</td><td>No</td><td>0.2</td><td>20</td></tr><tr><td>GF-H2O-comp-Vf15</td><td>H<sub>2</sub>O</td><td>No</td><td>No</td><td>Comp</td><td>15</td></tr><tr><td>GF-H2O-a02-Vf15</td><td>H<sub>2</sub>O</td><td>No</td><td>No</td><td>0.2</td><td>15</td></tr><tr><td>GBFc-H2O-a02-Vf15</td><td>H<sub>2</sub>O</td><td>Yes</td><td>No</td><td>0.2</td><td>15</td></tr><tr><td>GBFcD-Si-a02-Vf10</td><td>Si</td><td>Yes</td><td>Yes</td><td>0.2</td><td>10</td></tr><tr><td>GBFcD-Si-a02-Vf20</td><td>Si</td><td>Yes</td><td>Yes</td><td>0.2</td><td>20</td></tr><tr><td>GBFcD-H20-a02-Vf15</td><td>H<sub>2</sub>O</td><td>Yes</td><td>Yes</td><td>0.2</td><td>15</td></tr></tbody></table>

**Notes.** The ‘G’ stands for growth, ‘F’ for fragmentation, ‘B’ for bounce, ‘c’ for compaction, and ‘D’ for rotational disruption. The term ‘comp’ in the column for the monomer size $a_0$ denotes compact grains, with porosity evolution not taken into account. ‘D’ simulations were run only with PHANTOM. (*) This simulation includes a snow line, indicated by ‘S’, with an inner threshold of 20 m s$^{-1}$ and an outer threshold of 5 m s$^{-1}$ to model the CO snow line ($T = 20$ K, approximately ~100 au).

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_{\text{star}} = 1 M_{\odot}$, and the mass of the disc was $M_{\text{disc}} = 0.01 M_{\odot}$. The temperature was set by the choice of the aspect ratio $(H/R)_0 = 0.0895$ at $R_0 = 100$ au, with $q = 0.5$.

For PAMDEAS, we used $R_{\text{in}} = 1$ au, $R_{\text{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_{\text{in}} = 10$ au, $R_{\text{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 $\alpha = 5 \times 10^{-3}$ (by setting $\alpha_{\text{AV}} = 0.1658$) and the initial dust-to-gas ratio $\epsilon_0$ 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 = 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 $\phi \le 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 $\rho_s = 2700$ kg m$^{-3}$ and a surface energy $\gamma_s = 0.2$ J m$^{-2}$, in agreement with Yamamoto et al. (2014), who estimate $\gamma_s = 0.3$ J m$^{-2}$. This value is of the same order of magnitude as $\gamma_s = 0.15$ J m$^{-2}$, found

by Kimura et al. (2015, 2020), through experimental measurements with sicastar® aggregates (micromod Partikeltechnologie GmbH). The Young’s modulus was $\epsilon = 72$ GPa (Yamamoto et al. 2014), which implies, assuming that the critical separation between two monomers $\delta_c$ before their separation is of the same order of magnitude as $\xi_{\text{crit}}$ (Chokshi et al. 1993), a value of $\xi_{\text{crit}} \approx 6$ Å.

For pure water ice, the intrinsic density of monomers was $\rho_s = 1000$ kg m$^{-3}$, with a surface energy of $\gamma_s = 0.1$ J m$^{-2}$. The critical separation is still debated, as are many other properties of dust in protoplanetary discs. We chose a Young’s modulus of $\epsilon = 9.4$ GPa, following Yamamoto et al. (2014), which results in a value of $\xi_{\text{crit}} \sim 10$ Å, close to $\xi_{\text{crit}} \sim 8$ Å 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 $v_{\text{frag}} = 20$ m s$^{-1}$ 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. $\phi = 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
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison between PAMDEAS simulations GF-Si-comp-Vf20 (left) and GF-Si-a02-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.</figcaption></figure>

the Stokes numbers (computed as per Eq. 5 in Vericel et al. 2021) are similar in both simulations, with St ~ 10⁻<fnref n="1" /> 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<fnref n="1" />. 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 years, 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-Vf15 and GF-H2O-a02-Vf15 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

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Same as the top panels of Fig. 8 for PAMDEAS simulations GF-H2O-comp-Vf15 (left) and GF-H2O-a02-Vf15 (right).</figcaption></figure>

decimetres) and masses compared to their compact counterparts (a few hundred micrometres), demonstrating that porosity must be considered in dust evolution.

### 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, z) 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.

<footnote n="1">¹ To estimate the mass ratio between porous and compact grains, we can write the mass ratio as $m_{\text{porous}}/m_{\text{compact}} = \phi(s_{\text{porous}}/s_{\text{compact}})^3$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison between PAMDEAS simulations GF-Si-comp-Vf20 (**left**) and GF-Si-a02-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.</figcaption></figure>

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 \sim 5 \times 10^{-3}$, 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 \sim 10^{-1}$ 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 $\phi \sim 2 \times 10^{-3}$, making them ~50 000 times more massive than compact grains<fnref n="1" />. 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 \sim 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-Vf15 and GF-H2O-a02-Vf15 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

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Same as the top panels of Fig. 8 for PAMDEAS simulations GF-H₂O-comp-Vf15 (**left**) and GF-H₂O-a02-Vf15 (**right**).</figcaption></figure>

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 \times 10^5$ 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.

### 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, z$) 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 ($\approx 100$ au) and has a higher dust concentration, around $\epsilon \approx 3-4$, with a thickness of $\sim 10$ au. In contrast, the porous dust disc extends out to 150 au with a dust-to-gas ratio of $\epsilon \approx 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 $\epsilon < 5 \times 10^{-3}$ 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 $\phi \sim 5 \times 10^{-3}-10^{-2}$ (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

<footnote n="1">¹ To estimate the mass ratio between porous and compact grains, we can write the mass ratio as $m_{porous}/m_{compact} = \phi(s_{porous}/s_{compact})^3$.</footnote>
<figure><img src="image_11.png" /><figcaption>**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,z) plane. The bottom panels show the radial grain size distribution, with colour representing the Stokes number St.</figcaption></figure>

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 $\epsilon \approx 2-3$ for compact grains and $\epsilon \approx 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 millions 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^{-3} - 5 \times 10^{-3}$, 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.

## 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\phi$ larger than 100 µm, meaning grains with $s > 1$ mm.

### 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 sili-
<figure><img src="image_11.png" /><figcaption>**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, z) plane. The bottom panels show the radial grain size distribution, with colour representing the Stokes number St.</figcaption></figure>

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 $\epsilon \approx 2–3$ for compact grains and $\epsilon \approx 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^{-3} – 5 \times 10^{-3}$, 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.

## 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\phi$ larger than 100 µm, meaning grains with $s > 1$ mm.
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

cates for three fragmentation thresholds: $v_{frag, Si} = 10, 20, \text{ and } 40 \text{ m s}^{-1}$, from top to bottom, while the bottom row shows simulations with water ice for $v_{frag, H2O} = 15 \text{ m s}^{-1}$. 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 $\phi$ for compacted grains is not the maximum possible value $\phi_{max} = 0.74$. The maximum value that is reached is instead around 0.5 – 0.6 as in equilibrium, grain growth tends to decrease $\phi$, while compaction by fragmentation or bouncing tends to increase $\phi$. According to Fig. 6 and Eq. (A.9), the closer $\phi$ 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 $v_{frag, Si} = 10 \text{ m s}^{-1}$, 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 $\mu$m at maximum compaction before being accreted by the star. The filling factor reaches a minimum of $10^{-2}$ just before the fragmentation threshold, before being compacted to $\phi = 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\phi$ during compaction remains roughly the same as when there is no compaction. As $\text{St} \propto s\phi$, a similar Stokes number results in a similar drift.

With $v_{frag, Si} = 20 \text{ m s}^{-1}$, 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 $\mu$m and filling factors between 0.3 and 0.5.

For $v_{frag, Si} = 40 \text{ m s}^{-1}$, 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 \times 10^{-3}$. 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 $\mu$m, with a filling factor $\phi \sim 0.4$. Thus, a behaviour similar to simulations with silicates and $v_{frag, Si} = 20 \text{ m s}^{-1}$ is observed, but with larger variations in size and filling factor.
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

### 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: $v_{frag, Si} = 10, 20,$ and $40 \text{ m s}^{-1}$, from top to bottom, while the bottom row shows simulations with water ice for $v_{frag, H2O} = 15 \text{ m s}^{-1}$. 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 $\phi$ for compacted grains is not the maximum possible value $\phi_{max} = 0.74$. The maximum value that is reached is instead around 0.5–0.6 as in equilibrium, grain growth tends to decrease $\phi$, while compaction by fragmentation or bouncing tends to increase $\phi$. According to Fig. 6 and Eq. (A.9), the closer $\phi$ 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 $v_{frag, Si} = 10 \text{ m s}^{-1}$, 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^{-2}$

just before the fragmentation threshold, before being compacted to $\phi = 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\phi$ during compaction remains roughly the same as when there is no compaction. As $\text{St} \propto s\phi$, a similar Stokes number results in a similar drift.

With $v_{frag, Si} = 20 \text{ m s}^{-1}$, 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 $v_{frag, Si} = 40 \text{ m s}^{-1}$, 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 \times 10^{-3}$. 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 $\phi \sim 0.4$. Thus, a behaviour similar to simulations
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** 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. 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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**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 $\phi$ is plotted against grain size, with colour indicating the distance $r$. Simulation parameters are indicated in each row.</figcaption></figure>

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 = 0.2~\mu$m and $\phi = 1$. They then transition to a different growth regime ($\phi_{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 $v_{frag, Si} = 10~m~s^{-1}$, grains are compacted up to $\phi = 0.5$ and a minimum size of 40 $\mu$m, even at larger distances $r$. Bouncing doesn’t occur in this case, as fragmentation is highly effective. For $v_{frag, Si} = 20~m~s^{-1}$, grains fragment at larger sizes and smaller $\phi$. However, once the fragmentation threshold is reached, they are compacted up to $\phi = 0.6$ and sizes of 100 $\mu$m. Bouncing begins to appear but has no influence on grain evolution and is not visible in the figure. For $v_{frag, Si} = 40~m~s^{-1}$, aggregates fragment after growing to sizes of around a decimetre and $\phi \sim 2 \times 10^{-3}$. They then fragment and are compacted to $\phi = 0.7$ and sizes of 700 $\mu$m.

Finally, in the case of water ice, the same behaviour is observed. Aggregates are compacted until they reach $\phi = 0.4$ and $s = 300~\mu$m. As mentioned in Sect. 2.3, bouncing is accompa-

nied 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 $\phi < 0.3$. When $\phi$ 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.

### 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 $v_{frag, Si} = 10, 20$ and $40~m~s^{-1}$, 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
<figure><img src="image_13.png" /></figure>

**Fig. 12.** 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. 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.

with silicates and $v_{frag, Si} = 20 \text{ m s}^{-1}$ 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 = 0.2 \ \mu\text{m}$ and $\phi = 1$. They then transition to a different growth regime ($\phi_{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 $v_{frag, Si} = 10 \text{ m s}^{-1}$, grains are compacted up to $\phi = 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 $v_{frag, Si} = 20 \text{ m s}^{-1}$, grains fragment at larger sizes and smaller $\phi$. However, once the fragmentation threshold is reached, they are compacted up to $\phi = 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 $v_{frag, Si} = 40 \text{ m s}^{-1}$, aggregates fragment after growing to sizes of around a decimetre and $\phi \sim 2 \times 10^{-3}$. They then fragment and are compacted to $\phi = 0.7$ and sizes of 700 µm.

<figure><img src="image_14.png" /></figure>

**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 $\phi$ is plotted against grain size, with colour indicating the distance *r*. Simulation parameters are indicated in each row.

Finally, in the case of water ice, the same behaviour is observed. Aggregates are compacted until they reach $\phi = 0.4$ and $s = 300 \ \mu\text{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 $\phi < 0.3$. When $\phi$ 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.

### 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 $v_{frag, Si} = 10, 20$ and $40 \text{ m s}^{-1}$, and with water ice, from top to bottom in each panel (in the same order as in Fig. 12). In all cases, the
<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

the former follow the evolution of single grains in the disc mid-plane, 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 $\phi$ ranging from a few $10^{-2}$ to a few $10^{-3}$ depending on the simulation, and mostly at large $r$, while the smallest grains have been compacted to $\phi \sim 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 $v_{frag, Si} = 10$ and $20 \text{ m s}^{-1}$ (first two rows in both panels) are similar: without or with compaction, grains reach sizes of a few millimetres at $t = 100$ kyr for $v_{frag} = 10 \text{ m s}^{-1}$, while they grow faster and larger, up to a few centimetres, for $v_{frag, Si} = 20 \text{ m s}^{-1}$. However, after 300 kyr, compacted grains are at most 600-700 $\mu$m in size for $v_{frag} = 10 \text{ m s}^{-1}$ and a few millimetres for $v_{frag, Si} = 20 \text{ m s}^{-1}$. The former case is difficult to reconcile with observations of grains of mm size or larger in protoplanetary discs. When $v_{frag, Si} = 40 \text{ m s}^{-1}$ (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
<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

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 *v*<sub>frag, Si</sub> = 10 and 20 m s<sup>-1</sup> (first two rows in both panels) are similar: without or with compaction, grains reach sizes of a few millimetres at *t* = 100 kyr for *v*<sub>frag</sub> = 10 m s<sup>-1</sup>, while they grow faster and larger, up to a few centimetres, for *v*<sub>frag, Si</sub> = 20 m s<sup>-1</sup>. However, after 300 kyr, compacted grains are at most 600–700 µm in size for *v*<sub>frag</sub> = 10 m s<sup>-1</sup> and a few millimetres for *v*<sub>frag, Si</sub> = 20 m s<sup>-1</sup>. The former case is difficult to reconcile with observations of grains of mm size or larger in protoplanetary discs. When *v*<sub>frag, Si</sub> = 40 m s<sup>-1</sup> (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
<figure><img src="image_16.png" /><figcaption>**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*, *z*), with the light grey background indicating the gas disc thickness. Simulation parameters are indicated in each row.</figcaption></figure>

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 dust discs can be examined in Fig. 15, showing the size in the (*r*, *z*) plane. Like in other figures, only particles with $\epsilon \ge 5 \times 10^{-3}$ 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 mid-plane, 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 an intermediate situation.

### 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 $v_{frag, ext} = 5 \text{ m s}^{-1}$. 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 $v_{frag, int} = 20 \text{ m s}^{-1}$ 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*, *z*) plane, with a cut-off at $\epsilon \ge 5 \times 10^{-3}$, and the light grey gas disc in the back-
<figure><img src="image_16.png" /><figcaption>**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*, *z*), with the light grey background indicating the gas disc thickness. Simulation parameters are indicated in each row.</figcaption></figure>

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 dust discs can be examined in Fig. 15, showing the size in the (*r*, *z*) plane. Like in other figures, only particles with $\epsilon \ge 5 \times 10^{-3}$ 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.

### 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 $v_{frag, ext} = 5 \text{ m s}^{-1}$. 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 $v_{frag, int} = 20 \text{ m s}^{-1}$ to model silicate grains assumed to be surrounded by a layer of water ice. The outer layer only affects adhesive properties, as discussed
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** PHANTOM GBFcS-Si-a02-Vf205 simulation. **Top:** Aggregate sizes in the same plane (r, z). **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.</figcaption></figure>

ground, while the central panel displays the radial grain size dis-
tribution 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 mil-
limetres 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 fragmen-
tation 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 com-
paction 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.

4.4. Effects of rotational disruption

We present in Appendix E the results of 3D simulations of ro-
tational 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 Gar-
cia & Gonzalez (2020), porosity enables large dust grain forma-
tion. 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 evolu-
tion, 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 (Garc-
cia 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 conse-
quently 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 dis-
tribution of dust mainly depends on grain altitude. Our simula-
tions with porous and compacting grains provide better expla-
nations 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 as-
sumed 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 repro-
duced. 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
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** PHANTOM GBFcS-Si-a02-Vf205 simulation. *Top*: aggregate sizes in the same plane (*r*, *z*). *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*.</figcaption></figure>

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*, *z*) plane, with a cut-off at *ε* ≥ 5 × 10<sup>-3</sup>, 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<sup>-2</sup> 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<sup>-1</sup>. 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.

### 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
<figure><img src="image_18.png" /><figcaption>**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.</figcaption></figure>

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 ($\phi \sim 0.1 - 1$). In the outer regions, smaller grains (10–100 µm) are prevalent, with varying porosities ($\phi \sim 0.01 - 0.5$), which is consistent with the observations of HL Tau. The large grains ($\ge 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 $\phi \sim 10^{-3} - 10^{-2}$.

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 $\phi$ 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 particles 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

<table><caption>**Table 2.** Summary of the classifications of cometary grains from the Rosetta and Stardust missions, adapted from Güttler et al. (2019).</caption><thead><tr><th>Group</th><th>s</th><th>φ</th><th>Contribution</th></tr></thead><tbody><tr><td>Porous</td><td>1 µm – 1 m</td><td>0.05 to 0.9</td><td>Dominates the size distribution</td></tr><tr><td>Fractal</td><td>1 µm – 10 mm</td><td>&lt; 0.05</td><td>Low fraction</td></tr><tr><td>Solid</td><td>0.1 µm – 0.5 mm</td><td>&gt; 0.9</td><td>Very common</td></tr></tbody></table>

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.

## 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:
<figure><img src="image_18.png" /><figcaption>**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.</figcaption></figure>

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 ($\phi \sim 0.1-1$). In the outer regions, smaller grains (10–100 µm) are prevalent, with varying porosities ($\phi \sim 0.01-0.5$), which is consistent with the observations of HL Tau. The large grains ($\ge 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 $\phi \sim 10^{-3}-10^{-2}$.

Table 2 presents a summary of the classifications of cometary dust and its morphology from the Rosetta and Star-dust missions (Güttler et al. 2019). These observations found that the most common grains are represented by the porous particle group, with $\phi$ 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 particles 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.

<table><caption>**Table 2.** Summary of the classifications of cometary grains from the Rosetta and Star-dust missions, adapted from Güttler et al. (2019).</caption><thead><tr><th>Group</th><th>s</th><th>φ</th><th>Contribution</th></tr></thead><tbody><tr><td>Porous</td><td>1 µm–1 m</td><td>0.05 to 0.9</td><td>Dominates the size distribution</td></tr><tr><td>Fractal</td><td>1 µm–10 mm</td><td>&lt;0.05</td><td>Low fraction</td></tr><tr><td>Solid</td><td>0.1 µm–0.5 mm</td><td>&gt;0.9</td><td>Very common</td></tr></tbody></table>

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 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.

## 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
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Diagram summarising the different physical processes related to dust seen and studied in this paper.</figcaption></figure>

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 ($\phi \sim 10^{-4}-10^{-3}$) 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 ($\phi \in [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.

*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 Stelлаire 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).

## References

Arakawa, S., Okuzumi, S., Tatsuuma, M., et al. 2023, ApJ, 951, L16

Auffinger, J. & Laibe, G. 2018, MNRAS, 473, 796

Ballabio, G., Dipierro, G., Veronesi, B., et al. 2018, MNRAS, 477, 2766

Barge, P. & Sommeria, J. 1995, A&A, 295, L1

Barrière-Fouchet, L., Gonzalez, J.-F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Bate, M. R. 2022, MNRAS, 514, 2145

Blum, J. & Schräpler, R. 2004, Phys. Rev. Lett., 93, 115503

Blum, J. & Wurm, G. 2000, Icarus, 143, 138

Blum, J. & Wurm, G. 2008, ARA&A, 46, 21

Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1

Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43

Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806

Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102

Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647

Drążkowska, J. & Alibert, Y. 2017, A&A, 608, A92

Drążkowska, J., Windmark, F., & Dullemond, C. P. 2014, A&A, 567, A38

Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971

Fouchet, L., Gonzalez, J. F., & Maddison, S. T. 2010, A&A, 518, A16

Fouchet, L., Maddison, S. T., Gonzalez, J. F., & Murray, J. R. 2007, A&A, 474, 1037

Garcia, A. 2018, PhD thesis, Université de Lyon, available at https://theses.hal.science/tel-01977317/document, in French

Garcia, A. J. L. & Gonzalez, J.-F. 2020, MNRAS, 493, 1788

Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 375

Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Gonzalez, J.-F., Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015, Planet. Space Sci., 116, 48

Gonzalez, J. F., Pinte, C., Maddison, S. T., Ménard, F., & Fouchet, L. 2012, A&A, 547, A58
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Diagram summarising the different physical processes related to dust seen and studied in this paper.</figcaption></figure>

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 ($\phi \sim 10^{-4}-10^{-3}$) 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 micrometers to a few millimetres, but this time with filling factors compatible with observations ($\phi \in [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.

**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 Stelлаire 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.12729015.

## References

Arakawa, S., Okuzumi, S., Tatsuuma, M., et al. 2023, ApJ, 951, L16

Auffinger, J., & Laibe, G. 2018, MNRAS, 473, 796

Ballabio, G., Dipierro, G., Veronesi, B., et al. 2018, MNRAS, 477, 2766

Barge, P., & Sommeria, J. 1995, A&A, 295, L1

Barrière-Fouchet, L., Gonzalez, J.-F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Bate, M. R. 2022, MNRAS, 514, 2145

Blum, J., & Schräpler, R. 2004, Phys. Rev. Lett., 93, 115503

Blum, J., & Wurm, G. 2000, Icarus, 143, 138

Blum, J., & Wurm, G. 2008, ARA&A, 46, 21

Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1

Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43

Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806

Cuzzi, J. N., Dobrovołskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102

Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647

Drążkowska, J., & Alibert, Y. 2017, A&A, 608, A92

Drążkowska, J., Windmark, F., & Dullemond, C. P. 2014, A&A, 567, A38

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

Fouchet, L., Maddison, S. T., Gonzalez, J. F., & Murray, J. R. 2007, A&A, 474, 1037

Fouchet, L., Gonzalez, J. F., & Maddison, S. T. 2010, A&A, 518, A16

Garcia, A. 2018, PhD thesis, Université de Lyon, France

Garcia, A. J. L., & Gonzalez, J.-F. 2020, MNRAS, 493, 1788

Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375

Gonzalez, J. F., Pinte, C., Maddison, S. T., Ménard, F., & Fouchet, L. 2012, A&A, 547, A58

Gonzalez, J.-F., Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015, Planet. Space Sci., 116, 48

Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Grishin, E., Rozner, M., & Perets, H. B. 2020, ApJ, 898, L13

Guidi, G., Isella, A., Testi, L., et al. 2022, A&A, 664, A137

Gunkelmann, N., Ringl, C., & Urbassek, H. M. 2016, A&A, 589, A30

Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130