11institutetext: Universite Claude Bernard Lyon 1, CRAL UMR5574, ENS de Lyon,
CNRS, Villeurbanne, F-69622, France
11email<EMAIL_ADDRESS>22institutetext: School of Physics
and Astronomy, Monash University, Vic. 3800, Australia

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

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

(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 develop a module for porosity evolution for the 3D Smoothed
Particle Hydrodynamics (SPH) 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 10 to 100 times more
compact. Thus, mm sizes with typical filling factors of $\sim 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 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 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) show 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 the rotation disruption (Tatsuuma &
Kataoka, 2021), where the gas flow torques the grains, causing them to rotate.
When the centrifugal force exceeds the maximum tensile strength of the grains,
they shatter into fragments. Tatsuuma & Kataoka (2021) and Michoulier &
Gonzalez (2022a) studied this process for porous aggregates in protoplanetary
discs.

All of these barriers theoretically hinder dust growth. In order to understand
how the barriers are bypassed in nature to form planetesimals from micron-
sized grains, various solutions have been proposed. Some of them involve
creating dust traps with local pressure maxima, such as vortices (Barge &
Sommeria, 1995; Meheut et al., 2012), the 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; Drażkowska et al., 2014; Draż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, like 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
concentrates dust in filamentary structures, which can ultimately collapse
gravitationally.

One solution that has gained recent 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; Draż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, contributing to planetesimal formation
through coagulation (Garcia, 2018). Laboratory experiments studying growth of
highly porous aggregates are restricted to low velocity (Blum & Schräpler,
2004), 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).

Güttler et al. (2019), by looking at dust in comets, 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 show that porosity must be considered to explain
spectral energy distributions (SEDs, Kataoka et al., 2016; Zhang et al.,
2023), polarisation (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 smaller grains on disc surfaces are identified as fractal or porous
aggregates (Tazaki et al., 2023).

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
limitations of our models, in Sect. 5 and 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) assume 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 adopt 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_{\mathrm{coll}}$

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

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

$\left(\frac{\mathrm{d}m}{\mathrm{d}t}\right)_{\mathrm{grow}}=4\pi\rho_{\mathrm{d}}s^{2}v_{\mathrm{rel}},$
(2)

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

$v_{\mathrm{rel}}=\sqrt{2}v_{\mathrm{t}}\frac{\sqrt{\mathrm{Sc}-1}}{\mathrm{Sc}}.$
(3)

The turbulent velocity $v_{\mathrm{t}}$ is given by

$v_{\mathrm{t}}=\sqrt{2^{1/2}\mathrm{Ro}\,\alpha}\,c_{\mathrm{g}},$ (4)

$\mathrm{Ro}$ is the Rossby number, which is a constant equal to 3 (Stepinski
& Valageas, 1997), $\alpha$ the turbulent viscosity parameter (Shakura &
Sunyaev, 1973) and $c_{\mathrm{g}}$ the sound speed. The Schmidt number
$\mathrm{Sc}$ of a dust grain measures the grain’s coupling to the vortex

$\mathrm{Sc}=\left(1+\mathrm{St}\right)\sqrt{1+\frac{\Delta
v^{2}}{v_{\mathrm{t}}^{2}}},$ (5)

where $\Delta v$ is the difference in velocity between grains due to friction
with the gas. St measures the coupling between gas and dust, and is expressed
as

$\mathrm{St}=t_{\mathrm{s}}\Omega_{\mathrm{K}}.$ (6)

where $t_{\mathrm{s}}$ is the stopping time, that is the time needed for a
grain to reach the gas velocity, and $\Omega_{\mathrm{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) use 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 consider as compact
spheres with mass $m_{0}$, size $a_{0}$, and intrinsic density
$\rho_{\mathrm{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 consider 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 define the filling factor $\phi$, related to the porosity $\mathfrak{p}$,
as the ratio between the volume occupied by the monomers $V_{\mathrm{mat}}$
and the volume of the aggregate $V$

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

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

$m=\rho_{\mathrm{s}}\phi\frac{4\pi}{3}s^{3}.$ (8)

Two energies can be associated with an aggregate. The first is the kinetic
energy upon impact when two identical grains with mass $m$ collide with a
relative velocity $v_{\mathrm{rel}}$. In the referential of the centre of
mass, the kinetic energy is expressed as

$E_{\mathrm{kin}}=\frac{m}{4}\,v_{\mathrm{rel}}^{2},$ (9)

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

The second energy is the rolling energy $E_{\mathrm{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_{\mathrm{roll}}=6\pi^{2}\gamma_{\mathrm{s}}a_{0}\xi_{\mathrm{crit}},$ (10)

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

$\delta_{\mathrm{c}}=\left(\frac{27\pi^{2}}{2}\frac{\gamma_{\mathrm{s}}^{2}a_{0}}{\mathcal{E}^{2}}\right)^{\frac{1}{3}},$
(11)

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

$E_{\mathrm{roll}}=\left(2916\pi^{8}\frac{\gamma_{\mathrm{s}}^{5}a_{0}^{4}}{\mathcal{E}^{2}}\right)^{\frac{1}{3}}\approx
302\left(\frac{\gamma_{\mathrm{s}}^{5}a_{0}^{4}}{\mathcal{E}^{2}}\right)^{\frac{1}{3}},$
(12)

with $\mathcal{E}$ the 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_{\mathrm{kin}}<3b\,E_{\mathrm{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\ \mu$m) as its size $s$ increases, for different distances to
the star.

Figure 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\ \mu$m
monomers in the case of pure growth in our standard disc model (see Sect.
3.3).

###  2.3 Bouncing

In addition to growth, grains can also bounce. We have improved upon the model
initially developed by Garcia (2018). Wada et al. (2011), Shimaki & Arakawa
(2012a), and Arakawa et al. (2023) have numerically and experimentally shown
that grains are no longer capable of bouncing when the filling factor drops
below $\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 wave (Thornton & Ning, 1998).
Several velocities can be distinguished. The first is the sticking velocity
$v_{\mathrm{stick}}$, corresponding to the velocity at which grains no longer
stick together systematically. Johnson et al. (1971) and Thornton & Ning
(1998) provide the sticking velocity $v_{\mathrm{stick}}$ for the collision of
two identical grains

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

Dominik & Tielens (1997) provide another formula for $v_{\mathrm{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); Shimaki & Arakawa
(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_{\mathrm{stick}}\propto\mathcal{E}^{-1/3}$, by using the value of
$\mathcal{E}$ for monomers and retaining the factor of 2, a reasonable
compromise is achieved. Therefore, we keep

$\displaystyle v_{\mathrm{stick}}$
$\displaystyle=8.46\left(\frac{\gamma_{\mathrm{s}}^{5}s^{4}}{m^{3}\mathcal{E}^{2}}\right)^{1/6}$
(14) $\displaystyle=2.82\
\mathrm{mm\,s}^{-1}\left(\frac{\rho_{\mathrm{s}}}{1000\
\mathrm{kg\,m}^{-3}}\right)^{-1/2}\left(\frac{\gamma_{\mathrm{s}}}{0.1\
\mathrm{J\,m}^{-2}}\right)^{5/6}$
$\displaystyle\quad\times\left(\frac{\mathcal{E}}{10\
\mathrm{GPa}}\right)^{-1/3}\left(\frac{s}{1\
\mathrm{mm}}\right)^{-5/6}\phi^{-1/2},$

where we have used Eq. (8) to obtain a dependence on size and filling factor.
With the material properties we adopt for silicate and water ice grains (see
Sect. 3.3), we obtained $v_{\mathrm{stick}}=2.89\leavevmode\nobreak\ $mm s-1
and 5.26 mm s-1, respectively, for mm-sized grains with $\phi=0.3$. Weidling
et al. (2012) argued that instead of $\gamma_{\mathrm{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_{\mathrm{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_{\mathrm{stick}}$ of several to tens of
cm s-1. Additionally, other studies suggest 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
$\gtrsim 10^{4}$ Pa for $\phi>0.3$ (see equation 31 in Seizinger et al.,
2012), which increases $v_{\mathrm{stick}}$ by two orders of magnitude. We
thus also consider values of $v_{\mathrm{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_{\mathrm{yield}}$. Thornton & Ning (1998) express $v_{\mathrm{yield}}$ as a
function of the contact radius $a_{\mathrm{yield}}$ when yield occurs or of
the limiting contact pressure $p_{\mathrm{yield}}$. However, as these material
properties are difficult to obtain, and are in particular unknown for water
ice, they suggest to turn to values from impact experiments. In theirs,
Weidling et al. (2012) find a value for the yield velocity of SiO2 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_{\mathrm{yield}}=10v_{\mathrm{stick}}$. More
recently, Kimura et al. (2020) demonstrated that silicates are more prone to
sticking together with a higher surface energy. We thus kept a factor of 10
between the two velocities for silicates as well, and take

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

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

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

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

$\displaystyle v_{\mathrm{rel}}$ $\displaystyle\leq
v_{\mathrm{stick}}\quad\quad e=0$ (17) $\displaystyle
v_{\mathrm{stick}}<v_{\mathrm{rel}}$ $\displaystyle\leq
v_{\mathrm{yield}}\quad\quad e=1$ (18) $\displaystyle v_{\mathrm{rel}}$
$\displaystyle\geq v_{\mathrm{yield}}\quad\quad
e=\left[1.2\sqrt{3}\left(1-\frac{1}{6}\left(\frac{v_{\mathrm{yield}}}{v_{\mathrm{rel}}}\right)^{2}\right)\right.\times$
$\displaystyle\left.\left(1+2\sqrt{1.2\left(\frac{v_{\mathrm{rel}}}{v_{\mathrm{yield}}}\right)^{2}-0.2}\right)^{-1/2}-\left(\frac{v_{\mathrm{stick}}}{v_{\mathrm{rel}}}\right)^{2}\right]^{1/2}$
(19)

Figure 2: Coefficient of restitution $e$ when $v_{\mathrm{rel}}\geq
v_{\mathrm{yield}}$, given by Eq. (2.3).

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

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

Typical values of $v_{\mathrm{end}}$ are of the order of $10-100$ km s-1, much
larger than fragmentation thresholds.

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

#### Sticking probability and filling factor

Empirically, Weidling et al. (2012) determined the probability that a grain
sticks, or not. As previously mentioned, when
$v_{\mathrm{stick}}<v_{\mathrm{rel}}<v_{\mathrm{end}}$, an aggregate can
bounce instead of stick. Naturally, the closer $v_{\mathrm{rel}}$ is to
$v_{\mathrm{stick}}$, the higher the chances of sticking, and conversely for
$v_{\mathrm{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}=\left\\{\begin{array}[]{l l}\displaystyle
1,&\;v_{\mathrm{rel}}<v_{\mathrm{stick}},\\\ \\\
\displaystyle\frac{\log{\left(v_{\mathrm{rel}}\right)}-\log{\left(v_{\mathrm{end}}\right)}}{\log{\left(v_{\mathrm{stick}}\right)}-\log{\left(v_{\mathrm{end}}\right)}},&\;v_{\mathrm{stick}}\leq
v_{\mathrm{rel}}<v_{\mathrm{end}},\\\ \\\ \displaystyle
0.&\;v_{\mathrm{rel}}\geq v_{\mathrm{end}},\end{array}\right.$ (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{\mathrm{d}m}{\mathrm{d}t}\right)_{\mathrm{bounce}}=\mathcal{P}\left(\frac{\mathrm{d}m}{\mathrm{d}t}\right)_{\mathrm{grow}}.$
(22)

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

$\phi_{\mathrm{bounce}}\approx\phi_{\mathrm{i}}\left(1-\frac{\Delta
V}{V_{\mathrm{i}}}\right)^{-n},$ (23)

where $n=t/\tau_{\mathrm{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_{\mathrm{kin}}}{2\varUpsilon_{\mathrm{d}}},$ (24)

where $\varUpsilon_{\mathrm{d}}=\varUpsilon_{0}\phi^{\mathcal{N}}$ is the
dynamic compression resistance, and $\varUpsilon_{0}$ is the reference dynamic
compression resistance. Mellor (1974) give $\varUpsilon_{0}=9.8$ MPa and
$\mathcal{N}=4$ for water ice. However, $\varUpsilon$, 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 $\varUpsilon$ 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) provides a formula
for static compression pressure for highly porous grains only, as given in Eq.
(48), and Güttler et al. (2009, 2010) and Seizinger et al. (2012) provide 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_{\mathrm{roll}}$ and $\phi$ which is valid for both small and large $\phi$,
shown on Fig. 3,

$P_{\mathrm{bounce}}=\left(\frac{E_{\mathrm{roll}}}{a_{0}^{3}}\right)\left(\frac{1}{\phi}-\frac{1}{0.74}\right)^{-3}.$
(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.

Figure 4: Growth and bouncing of a water ice grain with $a_{0}=0.2\ \mu$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 $\varUpsilon_{d}$ given by Shimaki & Arakawa (2012a). The dashed and
dot-dashed lines show the impact of multiplying $v_{\mathrm{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).

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 $\varUpsilon$, which is valid only for water ice (or a
mixture of ice and silicates), and our model where $P_{\mathrm{bounce}}$
replaces $\varUpsilon$ 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_{\mathrm{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\geq 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 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_{\mathrm{stick}}$ (see Sect. 2.5).

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

$\phi_{\mathrm{coll\,\&\,bounce}}=\left\\{\begin{array}[]{l
l}\displaystyle\phi_{\mathrm{coll}},&\;v_{\mathrm{rel}}<v_{\mathrm{yield}},\\\
\\\
\displaystyle\phi_{\mathrm{coll}}\mathcal{P}+(1-\mathcal{P})\phi_{\mathrm{bounce}},&\;v_{\mathrm{yield}}\leq
v_{\mathrm{rel}}<v_{\mathrm{end}},\\\ \\\
\displaystyle\phi_{\mathrm{bounce}},&\;v_{\mathrm{rel}}\geq
v_{\mathrm{end}},\end{array}\right.$ (26)

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

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

In this model, when the relative velocity is close to the threshold, the mass
loss is less significant. A fragmenting grain loses half of its mass after a
collision time $\tau_{\mathrm{coll}}$ ($v_{\mathrm{rel}}=v_{\mathrm{frag}}$)
or more ($v_{\mathrm{rel}}>v_{\mathrm{frag}}$). In the case where
$v_{\mathrm{rel}}\gg v_{\mathrm{frag}}$, we recover the scenario described by
Gonzalez et al. (2015), where the entire grain fragments after a collision
time, independently of $v_{\mathrm{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) show more complex behaviours that are challenging to
model. Fragmentation involves a cascade of sizes and depends on arameters such
as composition, porosity, shape, impact parameter, mass ratio, impact
velocity, etc. 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_{\mathrm{frag}}$, an estimate of the
energy required to break the bond between two monomers $E_{\mathrm{break}}$ is
given by Dominik & Tielens (1997),

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

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

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

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

$E_{\mathrm{break}}\simeq
C_{\mathrm{break}}48\left(\frac{\gamma_{\mathrm{s}}^{5}a_{0}^{4}}{\mathcal{E}^{2}}\right)^{\frac{1}{3}}.$
(30)

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

$E_{\mathrm{frag}}=\kappa N_{\mathrm{tot}}E_{\mathrm{break}},$ (31)

where $N_{\mathrm{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_{\mathrm{frag}}$
in terms of $E_{\mathrm{frag}}$ as

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

These equations depend on $\gamma_{\mathrm{s}}$ and $a_{0}$. An error in
estimating these parameters can thus lead to changes in the values of
$E_{\mathrm{frag}}$ and $v_{\mathrm{frag}}$.

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

Figure 5: Filling factor $\phi_{\mathrm{f}}/\phi_{\mathrm{i}}$ as a function
of $v_{\mathrm{rel}}$. Left: Experimental values for different initial filling
factors $\phi_{\mathrm{i}}$ provided by Ringl et al. (2012) in orange and
Gunkelmann et al. (2016) in blue, red and green. Right: Comparison between our
model given by Eq. (35), and fits of the curves given by Ringl et al. (2012)
and Gunkelmann et al. (2016) in our fragmentation module.

However, pure water ice grains are unlikely. Most often, grains are composed
of solid materials like silicates, covered by a layer of another more volatile
element such as water or organic matter. The fragmentation velocity is
primarily determined by the material forming the surface layer, and not the
internal material. Therefore, we choose different thresholds for the
silicates, as there are still uncertainties about their values, to mimic the
effect of a layer of volatiles: $v_{\mathrm{frag}}=5$, 10, or $20\
\mathrm{m\,s^{-1}}$ for more fragile materials like water or CO, and
$v_{\mathrm{frag}}=40\ \mathrm{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 assume 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 $\mu$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_{\mathrm{frag}}$ used
to fragment the grain from the kinetic energy $E_{\mathrm{kin}}$ released upon
impact. We can rewrite $E_{\mathrm{frag}}$ as the energy $E_{\mathrm{break}}$
required to break bonds between monomers, multiplied by the number of ejected
monomers

$E_{\mathrm{comp}}=E_{\mathrm{kin}}-E_{\mathrm{frag}}=E_{\mathrm{kin}}-\frac{(2m_{\mathrm{i}}-m_{\mathrm{f}})}{m_{0}}\kappa
E_{\mathrm{break}}.$ (33)

Here, $(2m_{\mathrm{i}}-m_{\mathrm{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 will therefore consider
$\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 model 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_{\mathrm{kin}}-(2m_{\mathrm{i}}-m_{\mathrm{f}})\kappa
E_{\mathrm{break}}m_{0}^{-1}}{2\displaystyle\left(\frac{E_{\mathrm{roll}}}{a_{0}^{3}}\right)\left(\frac{1}{\phi}-\frac{1}{0.74}\right)^{-3}}.$
(34)

Figure 6: Filling factor $\phi$ as a function of size $s$ for different
distances to the star $r$ for silicate grains composed of $0.2\ \mu$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_{\mathrm{stick}}$ from Eq. (14) by 100.

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_{\mathrm{f}}/\phi_{\mathrm{i}}$ of the final filling factor to the
initial one as a function of the relative velocity $v_{\mathrm{rel}}$ for
three values of $\phi_{\mathrm{i}}$. It should be noted that Gunkelmann et al.
(2016) define $v_{\mathrm{frag}}$ as the point when a grain starts to lose a
monomer, which is equal to 0.17 $\mathrm{m\,s^{-1}}$, and differs from our
definition as the point where the grain loses half of its mass. We deduce that
the ratio $\phi_{\mathrm{f}}/\phi_{\mathrm{i}}$ varies as
$\exp{\left(1-\left(v_{\mathrm{rel}}/v_{\mathrm{frag}}\right)^{2}\right)}$. 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_{\mathrm{rel}}/v_{\mathrm{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_{\mathrm{rel}}/v_{\mathrm{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_{\mathrm{frag-
comp}}\approx\phi_{\mathrm{i}}\left(1-\exp\left(1-\left(\frac{v_{\mathrm{rel}}}{v_{\mathrm{frag}}}\right)^{2}\right)\frac{\Delta
V}{V_{\mathrm{i}}}\right)^{-n},$ (35)

where $n$ represents the number of collisions during a time step. The effect
of compaction can be seen in Fig. 6, produced with Pamdeas with growth,
bouncing, fragmentation and radial drift (see Sect. 3.1). In the left panel,
where compaction is not taken into account, grains reach a growth-
fragmentation equilibrium and stay in the same range of size and filling
factor. $\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_{\mathrm{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_{\mathrm{coll}}$, $\phi_{\mathrm{gas}}$, $\phi_{\mathrm{grav}}$

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

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

Compute $\phi_{\mathrm{grow}}$

if (Grains can bounce) then

Compute $\phi_{\mathrm{coll\,\&\,bounce}}$

Compute $\phi_{\mathrm{f}}=\max(\phi_{\mathrm{coll\,\&\,bounce}}$,
$\phi_{\mathrm{min}})$

else

Compute $\phi_{\mathrm{f}}=\max(\phi_{\mathrm{grow}}$, $\phi_{\mathrm{min}})$

end if

else

if (Fragmentation with compaction) then

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

Compute $\phi_{\mathrm{f}}=\max(\phi_{\mathrm{frag-comp}}$,
$\phi_{\mathrm{min}})$

else

Compute $\phi_{\mathrm{f}}=\max(\phi_{\mathrm{i}},\phi_{\mathrm{min}})$

end if

end if

$\phi_{\mathrm{final}}=\min(\phi_{\mathrm{f}},0.74)$

##  3 Numerical simulations

###  3.1 The 1D code Pamdeas

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

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 allows 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_{\mathrm{g}}(r)=\Sigma_{\mathrm{g,\ 0}}(r/R_{0})^{-p}$ (36)

and

$T_{\mathrm{g}}(r)=T_{\mathrm{g,\ 0}}(r/R_{0})^{-q},$ (37)

where $R_{0}$ is a reference radius. The disc extends from $R_{\mathrm{in}}$
to $R_{\mathrm{out}}$ and has a total mass $M_{\mathrm{disc}}$. As the gas
disc structure is held fixed, the dust-to-gas ratio $\varepsilon$ 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 SPH 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_{\mathrm{in}}}{r}}\right)\left(\frac{r}{R_{\mathrm{c}}}\right)^{-p}\exp{\left[-\left(r/R_{\mathrm{c}}\right)^{2-p}\right]},$
(38)

where $R_{\mathrm{c}}$ is the cut-off radius. With the exponential tapering,
the outer density profile is smooth, rather than being sharply truncated at
$R_{\mathrm{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_{\mathrm{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 perform a simulation solely with gas to avoid any spurious behaviour
of the dust grains during gas relaxation. We increase 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 add the dust phase with the same spatial distribution as the relaxed
gas, with an initially uniform dust-to-gas ratio $\varepsilon_{0}$. Unlike in
Pamdeas, we do not take the monomer size as the initial grain size. Indeed, if
the simulation was initialized with a single, user-defined, value of the
initial size for all grains, the filling factor computed from $s$ with Eq.
(61) 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 use 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},$ (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 initialization, the evolution variable is
the grain mass, and the mass growth rate and filling factor are computed using
the equations in Sect. 2.

Table 1: List of simulations carried out with both Pamdeas and Phantom. ‘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. *This simulation includes a snow line, indicated by ‘S’, with an inner threshold of $20\ \mathrm{\mathrm{m\,s^{-1}}}$ and an outer threshold of $5\ \mathrm{\mathrm{m\,s^{-1}}}$ to model the CO snow line ($T=20$ K, approximately $\sim$ 100 au). ‘D’ simulations were run only with Phantom. Name | species | bouncing $\&$ compaction | disruption |  $a_{0}$ ($\mu$m) |  $v_{\mathrm{frag}}$ ($\mathrm{m\,s^{-1}}$)
---|---|---|---|---|---
GF-Si-comp-Vf10 | Si | no | no | comp | 10
GF-Si-comp-Vf20 | Si | no | no | comp | 20
GF-Si-comp-Vf40 | Si | no | no | comp | 40
GF-Si-a02-Vf10 | Si | no | no | 0.2 | 10
GF-Si-a02-Vf20 | Si | no | no | 0.2 | 20
GF-Si-a02-Vf40 | Si | no | no | 0.2 | 40
GBFc-Si-a02-Vf10 | Si | yes | no | 0.2 | 10
GBFc-Si-a02-Vf20 | Si | yes | no | 0.2 | 20
GBFc-Si-a02-Vf40 | Si | yes | no | 0.2 | 40
GBFcS-Si-a02-Vf205* | Si | yes | no | 0.2 | 20
GF-H2O-comp-Vf15 | H2O | no | no | comp | 15
GF-H2O-a02-Vf15 | H2O | no | no | 0.2 | 15
GBFc-H2O-a02-Vf15 | H2O | yes | no | 0.2 | 15
GBFcD-Si-a02-Vf10 | Si | yes | yes | 0.2 | 10
GBFcD-Si-a02-Vf20 | Si | yes | yes | 0.2 | 20
GBFcD-H20-a02-Vf15 | H2O | yes | yes | 0.2 | 15

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 have chosen to use a disc model which represents an ‘average disc’
(Williams & Best, 2014). The mass of the star is fixed at
$M_{\mathrm{star}}=1\ \mathrm{M_{\odot}}$, and the mass of the disc is
$M_{\mathrm{disc}}=0.01\ \mathrm{M_{\odot}}$. The temperature is 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 use $R_{\mathrm{in}}=1$ au, $R_{\mathrm{out}}=300$ au and
$p=1$. All grains evolve from an initial size equal to the monomer size. We
run the simulations up to $t=1$ Myr.

For Phantom, we take $R_{\mathrm{in}}=10$ au, $R_{\mathrm{out}}=400$ au and
$p=0.75$. This value of $p$, combined with the exponential tapering, leads 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_{\rm AV}=0.1658$) and the initial dust-to-gas ratio
$\varepsilon_{0}$ to a typical value of 1%. We evolve the simulations for
300,000 yr.

For both codes, the size of monomers is set to $a_{0}=0.2\ \mu$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 $\mu$m
monomers. Verrios et al. (2022) reached the same conclusion with respect to
the dust settling in IM Lupi, requiring porosities of $\phi\lesssim 0.1$ to
match the observations. We have chosen to use two different species commonly
found in discs: water ice and silicates.

For silicate grains, we have chosen an intrinsic density of
$\rho_{\mathrm{s}}=2\,700\ \mathrm{kg\,m^{-3}}$ and a surface energy
$\gamma_{\mathrm{s}}=0.2\ \mathrm{J\,m^{-2}}$, in agreement with Yamamoto et
al. (2014), who estimate $\gamma_{\mathrm{s}}=0.3\ \mathrm{J\,m^{-2}}$. This
value is of the same order of magnitude as $\gamma_{\mathrm{s}}=0.15\
\mathrm{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 is $\mathcal{E}=72$ GPa
(Yamamoto et al., 2014), which implies, assuming that the critical separation
between two monomers $\delta_{\mathrm{c}}$ before their separation is of the
same order of magnitude as $\xi_{\mathrm{crit}}$ (Chokshi et al., 1993), a
value of $\xi_{\mathrm{crit}}\approx 6$ Å.

For pure water ice, the intrinsic density of monomers is
$\rho_{\mathrm{s}}=1\,000\ \mathrm{kg\,m^{-3}}$, with a surface energy of
$\gamma_{\mathrm{s}}=0.1\ \mathrm{J\,m^{-2}}$. The critical separation is
still debated, as are many other properties of dust in protoplanetary discs.
We have chosen a Young’s modulus of $\mathcal{E}=9.4$ GPa, following Yamamoto
et al. (2014), which results in a value of $\xi_{\mathrm{crit}}\sim 10$ Å,
close to $\xi_{\mathrm{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_{\mathrm{frag}}=20\ \mathrm{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

Figure 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 to 300 au.

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 $\mathrm{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 $\mathrm{St}\sim 10^{-1}$ when the grains reach the
fragmentation threshold, the sizes are drastically different. Compact grains
reach about 100 $\mu$m, whereas porous grains reach about 3 cm. These porous
grains have a filling factor $\phi\sim 2\times 10^{-3}$, making them
$\sim$50,000 times more massive than compact grains111To estimate the mass
ratio between porous and compact grains, we can write the mass ratio as
$m_{\mathrm{porous}}/m_{\mathrm{compact}}=\phi(s_{\mathrm{porous}}/s_{\mathrm{compact}})^{3}$..
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 $\mu$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 $\mu$m. Since the fragmentation threshold
Stokes numbers are similar ($\mathrm{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
$\mu$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 $\mu$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\times 10^{5}$ times
more massive than compact grains.

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

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 hundreds of micrometres), demonstrating that porosity must
be considered in dust evolution.

####  4.1.2 3D global Simulations with Phantom

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

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. Gray 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\ \mathrm{kyr}$. At the end of the simulation, the compact dust
disc is less radially extended ($\approx 100\ \mathrm{au}$) and has a higher
dust concentration, around $\varepsilon\approx 3-4$, with a thickness of $\sim
10\ \mathrm{au}$. In contrast, the porous dust disc extends out to 150
$\mathrm{au}$ with a dust-to-gas ratio of $\varepsilon\approx 0.8-1$. Porous
grains thus allow retaining dust in the outer regions 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 $\varepsilon<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 $\mathrm{au}$. In contrast, compact grains struggle to reach millimetre
sizes and the largest ones remain in the 100-500 $\mathrm{\mu 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 $\sim 10-20$ compared to observations (see Sect.
5). In both cases, grains do not cross $\mathrm{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\
\mathrm{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 $\mathrm{au}$), there are only
monomers. Monomers at 600 $\mathrm{au}$ have larger St values than those at
500 $\mathrm{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 ($\sim 20\ \mathrm{au}$) and extend radially
further (¿ 400 $\mathrm{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 $\varepsilon\approx 2-3$ for compact
grains and $\varepsilon\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 $\mathrm{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 $\mathrm{St}=1$ in the outer disc region, where lower gas density
and a higher growth rate enable reaching these St values. 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 $\mu$m,
meaning grains with $s>1$ mm.

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

####  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_{\mathrm{frag,\>Si}}=10$, 20, and 40 $\mathrm{m\,s^{-1}}$, from top to
bottom, while the bottom row shows simulations with water ice for
$v_{\mathrm{frag,\>H2O}}=15\ \mathrm{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_{\mathrm{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. (48), 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_{\mathrm{frag,\>Si}}=10\ \mathrm{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
$\mathrm{St}\propto s\phi$, a similar Stokes number results in a similar
drift.

With $v_{\mathrm{frag,\>Si}}=20\ \mathrm{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 (St) 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_{\mathrm{frag,\>Si}}=40\ \mathrm{\mathrm{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 $\mathrm{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_{\mathrm{frag,\>Si}}=20\
\mathrm{\mathrm{m\,s^{-1}}}$ is observed, but with larger variations in size
and filling factor.

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

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_{\mathrm{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_{\mathrm{frag,\>Si}}=10\ \mathrm{\mathrm{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_{\mathrm{frag,\>Si}}=20\ \mathrm{\mathrm{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_{\mathrm{frag,\>Si}}=40\
\mathrm{\mathrm{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 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.

Figure 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. Note that the colorbar range is
different from that in Fig. 12. Figure 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 gray background indicating the gas disc
thickness. Simulation parameters are indicated in each row.

####  4.2.2 3D global Simulations with Phantom

Figure 14 compares Phantom simulations without (top) and with (bottom)
compaction during fragmentation and bouncing. We show the radial grain size
distribution, colour-coded with the filling factor, for simulations with
silicates and $v_{\mathrm{frag,\>Si}}=10$, 20 and $40\ \mathrm{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 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 $\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_{\mathrm{frag,\>Si}}=10$ and 20
$\mathrm{m\,s^{-1}}$ (first two rows in both panels) are similar: without or
with compaction, grains reach sizes of a few mm at $t=100$ kyr for
$v_{\mathrm{frag}}=10$ $\mathrm{m\,s^{-1}}$, while they grow faster and
larger, up to a few cm, for $v_{\mathrm{frag,\>Si}}=20$ $\mathrm{m\,s^{-1}}$.
However, after 300 kyr, compacted grains are at most 600-700 $\mu$m in size
for $v_{\mathrm{frag}}=10$ $\mathrm{m\,s^{-1}}$ and a few mm for
$v_{\mathrm{frag,\>Si}}=20$ $\mathrm{m\,s^{-1}}$. The former case is difficult
to reconcile with observations of grains of mm size or larger in
protoplanetary discs. When $v_{\mathrm{frag,\>Si}}=40$ $\mathrm{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 approximation used in the dust-as-
mixture formalism (Youdin & Goodman, 2005; Laibe & Price, 2014) and those
simulations should not be considered valid. Finally, for simulations
GF-H2O-a02-Vf15 and GBFc-H2O-a02-Vf15 (bottom rows), the same pattern as for
simulations GF-Si-a02-Vf20 and GBFc-Si-a02-Vf20 is observed, with grains
reaching slightly larger sizes. The joint evolution of grain size and filling
factor is described in Appendix D.

The thickness of the of the dust discs can be examined in Fig. 15, showing the
size in the ($r$, $z$) plane. Like in other figures, only particles with
$\varepsilon\geq 5\times 10^{-3}$ are shown to eliminate dust-depleted
regions. The light gray 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 $\sim 5$ (comparing the left panels of
Figure 8 to the left panels of 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\sim 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_{\mathrm{frag,\>ext}}=5\ \mathrm{m\,s^{-1}}$. In
the inner region, instead of taking a threshold of 15 $\mathrm{m\,s^{-1}}$,
which would correspond to grains made entirely of water ice, we set
$v_{\mathrm{frag,\>int}}=20\ \mathrm{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.

Figure 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 $\phi$ plotted
against grain size $s$, with colour indicating the distance $r$.

The top panel shows the grain size $s$ in the ($r$, $z$) plane, with a cut-off
at $\varepsilon\geq 5\times 10^{-3}$, and the light gray 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 $\mathrm{St}\sim 5\times 10^{-2}$ when they are compacted, mainly
between 50 and 100 au, down to a few mm. 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 $\sim 100$ $\mu$m when they are still
porous, down to a few tens of $\mu$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 $\mathrm{m\,s^{-1}}$. 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 $\phi\sim 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 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 cm-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.

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

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 assume either compact ($\phi=1$) of porous ($\phi=0.1$ for Zhang
et al. 2023 and 0.2 for Ohashi et al. 2023) grains and fit for the grain sizes
for which SED and polarization data are best reproduced. In the inner disc
(interior to 40 au), compact grains of a few hundred $\mu$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 $\mu$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 caution 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\ \mu$m) are
prevalent, with varying porosities ($\phi\sim 0.01-0.5$), which is consistent
with the observations of HL Tau. The large grains ($\gtrsim 3$ cm) seen in the
DG Tau disc are however challenging to reproduce in our simulations. However,
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 particules of very different sizes or masses. High-mass-ratio
collisions have been found to produce stronger compression (Tanaka et al.,
2023). Taking them into account would help aggregates to reach their compacted
state sooner. The mono-disperse model nonetheless results in global size
distributions that are similar to those obtained with polydisperse models
(Gonzalez et al., 2017; Vorobyov et al., 2018). The reader is referred to
Laibe et al. (2008); Vericel & Gonzalez (2020); Vericel et al. (2021);
Michoulier et al. (2024) for additional discussions of this formalism.

Table 2: Summary of the classifications of cometary grains from the Rosetta and Stardust missions, adapted from Güttler et al. (2019). Group | $s$ | $\phi$ | Contribution
---|---|---|---
Porous |  1 $\mu$m – 1 m |  $0.05$ to $0.9$ |  Dominates the size distribution
Fractal | 1 $\mu$m – 10 mm | $<0.05$ | Low fraction
Solid | 0.1 $\mu$m – 0.5 mm | $>0.9$ | Very common
Figure 18: Diagram summarizing the different physical processes related to
dust seen and studied in this paper.

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 on 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 porosity, compaction, the
constituent species of the grains, and the fragmentation threshold. Our
conclusions are as follows:

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

With growth and fragmentation alone, we find 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. 3.

When including compaction during fragmentation and bouncing, we formed grains
of several hundred $\mu$m to a few mm, but this time with filling factors
compatible with observations ($\phi\in[0.1;1]$).

  4. 4.

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 Stellaire and the Programme
National de Planétologie of CNRS (Centre National de la Recherche
Scientifique)/INSU (Institut National des Sciences de l’Univers), France. We
gratefully acknowledge support from the PSMN (Pôle Scientifique de
Modélisation Numérique) of the ENS de Lyon for the computing resources. This
project has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreements No
210021 and No 823823 (DUSTBUSTERS). DJP is grateful for Australian Research
Council Discovery Project funding via DP220103767 and DP240103290. We thank
the Australian-French Association for Research and Innovation (AFRAN) for
financial support. Figures were made with the Python library matplotlib
(Hunter, 2007).

## References

  * Arakawa et al. (2023) Arakawa, S., Okuzumi, S., Tatsuuma, M., et al. 2023, ApJ, 951, L16
  * Auffinger & Laibe (2018) Auffinger, J. & Laibe, G. 2018, MNRAS, 473, 796
  * Ballabio et al. (2018) Ballabio, G., Dipierro, G., Veronesi, B., et al. 2018, MNRAS, 477, 2766
  * Barge & Sommeria (1995) Barge, P. & Sommeria, J. 1995, A&A, 295, L1
  * Barrière-Fouchet et al. (2005) Barrière-Fouchet, L., Gonzalez, J.-F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185
  * Bate (2022) Bate, M. R. 2022, MNRAS, 514, 2145
  * Blum & Schräpler (2004) Blum, J. & Schräpler, R. 2004, Phys. Rev. Lett., 93, 115503
  * Blum & Wurm (2000) Blum, J. & Wurm, G. 2000, Icarus, 143, 138
  * Blum & Wurm (2008) Blum, J. & Wurm, G. 2008, ARA&A, 46, 21
  * Brauer et al. (2008) Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Chokshi et al. (1993) Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806
  * Cuzzi et al. (1993) Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
  * Dominik & Tielens (1997) Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647
  * Drażkowska & Alibert (2017) Drażkowska, J. & Alibert, Y. 2017, A&A, 608, A92
  * Drażkowska et al. (2014) Drażkowska, J., Windmark, F., & Dullemond, C. P. 2014, A&A, 567, A38
  * Dullemond & Dominik (2005) Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971
  * Fouchet et al. (2010) Fouchet, L., Gonzalez, J. F., & Maddison, S. T. 2010, A&A, 518, A16
  * Fouchet et al. (2007) Fouchet, L., Maddison, S. T., Gonzalez, J. F., & Murray, J. R. 2007, A&A, 474, 1037
  * Garcia (2018) Garcia, A. 2018, Phd thesis, Université de Lyon, available at https://theses.hal.science/tel-01977317/document, in French
  * Garcia & Gonzalez (2020) Garcia, A. J. L. & Gonzalez, J.-F. 2020, MNRAS, 493, 1788
  * Gingold & Monaghan (1977) Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 375
  * Gonzalez et al. (2017) Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984
  * Gonzalez et al. (2015) Gonzalez, J.-F., Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015, Planet. Space Sci., 116, 48
  * Gonzalez et al. (2012) Gonzalez, J. F., Pinte, C., Maddison, S. T., Ménard, F., & Fouchet, L. 2012, A&A, 547, A58
  * Grishin et al. (2020) Grishin, E., Rozner, M., & Perets, H. B. 2020, ApJ, 898, L13
  * Guidi et al. (2022) Guidi, G., Isella, A., Testi, L., et al. 2022, A&A, 664, A137
  * Gunkelmann et al. (2016) Gunkelmann, N., Ringl, C., & Urbassek, H. M. 2016, A&A, 589, A30
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
  * Güttler et al. (2009) Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130
  * Güttler et al. (2019) Güttler, C., Mannel, T., Rotundi, A., et al. 2019, A&A, 630, A24
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Hyodo et al. (2019) Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90
  * Jacquet et al. (2011) Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591
  * Johansen & Youdin (2007) Johansen, A. & Youdin, A. 2007, ApJ, 662, 627
  * Johnson et al. (1971) Johnson, K. L., Kendall, K., & Roberts, A. D. 1971, Proceedings of the Royal Society of London Series A, 324, 301
  * Kataoka et al. (2016) Kataoka, A., Muto, T., Momose, M., Tsukagoshi, T., & Dullemond, C. P. 2016, ApJ, 820, 54
  * Kataoka et al. (2015) Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78
  * Kataoka et al. (2019) Kataoka, A., Okuzumi, S., & Tazaki, R. 2019, ApJ, 874, L6
  * Kataoka et al. (2013a) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013a, A&A, 557, L4
  * Kataoka et al. (2013b) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013b, A&A, 554, A4
  * Kimura et al. (2015) Kimura, H., Wada, K., Senshu, H., & Kobayashi, H. 2015, ApJ, 812, 67
  * Kimura et al. (2020) Kimura, H., Wada, K., Yoshida, F., et al. 2020, MNRAS, 496, 1667
  * Klahr & Bodenheimer (2003) Klahr, H. H. & Bodenheimer, P. 2003, ApJ, 582, 869
  * Kobayashi & Tanaka (2010) Kobayashi, H. & Tanaka, H. 2010, Icarus, 206, 735
  * Kothe et al. (2013) Kothe, S., Blum, J., Weidling, R., & Güttler, C. 2013, Icarus, 225, 75
  * Kretke & Lin (2007) Kretke, K. A. & Lin, D. N. C. 2007, ApJ, 664, L55
  * Krijt et al. (2015) Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2015, A&A, 574, A83
  * Laibe et al. (2008) Laibe, G., Gonzalez, J. F., Fouchet, L., & Maddison, S. T. 2008, A&A, 487, 265
  * Laibe & Price (2014) Laibe, G. & Price, D. J. 2014, MNRAS, 440, 2136
  * Lebreuilly et al. (2023) Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326
  * Li et al. (2019) Li, R., Youdin, A. N., & Simon, J. B. 2019, ApJ, 885, 69
  * Loren-Aguilar & Bate (2015) Loren-Aguilar, P. & Bate, M. R. 2015, MNRAS, 453, L78
  * Lucy (1977) Lucy, L. B. 1977, AJ, 82, 1013
  * Meheut et al. (2012) Meheut, H., Keppens, R., Casse, F., & Benz, W. 2012, A&A, 542, A9
  * Mellor (1974) Mellor, M. 1974, A Review of Basic Snow Mechanics (U.S. Army Cold Regions Research and Engineering Laboratory)
  * Michoulier & Gonzalez (2022a) Michoulier, S. & Gonzalez, J.-F. 2022a, MNRAS, 517, 3064
  * Michoulier & Gonzalez (2022b) Michoulier, S. & Gonzalez, J. F. 2022b, in SF2A-2022: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. Eds.: J. Richard, ed. J. Richard, A. Siebert, E. Lagadec, N. Lagarde, O. Venot, J. Malzac, J. B. Marquette, M. N’Diaye, & B. Briot, 323–324
  * Michoulier et al. (2024) Michoulier, S., Gonzalez, J.-F., Grishin, E., & Petetin, C. 2024, A&A, 686, A32
  * Musiolik et al. (2016) Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016, ApJ, 827, 63
  * Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
  * Ohashi et al. (2023) Ohashi, S., Momose, M., Kataoka, A., et al. 2023, ApJ, 954, 110
  * Okuzumi et al. (2012) Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
  * Okuzumi et al. (2009) Okuzumi, S., Tanaka, H., & Sakagami, M.-a. 2009, ApJ, 707, 1247
  * Ormel et al. (2007) Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215
  * Paardekooper & Mellema (2004) Paardekooper, S. J. & Mellema, G. 2004, A&A, 425, L9
  * Paraskov et al. (2006) Paraskov, G. B., Wurm, G., & Krauss, O. 2006, ApJ, 648, 1219
  * Pinte et al. (2019) Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nature Astronomy, 3, 1109
  * Planes et al. (2021) Planes, M. B., Millán, E. N., Urbassek, H. M., & Bringa, E. M. 2021, MNRAS, 503, 1717
  * Price & Laibe (2015) Price, D. J. & Laibe, G. 2015, MNRAS, 451, 813
  * Price et al. (2018) Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, Publications of the Astronomical Society of Australia, 35, e031
  * Ringl et al. (2012) Ringl, C., Bringa, E. M., Bertoldi, D. S., & Urbassek, H. M. 2012, ApJ, 752, 151
  * Ringl et al. (2012) Ringl, C., Bringa, E. M., Bertoldi, D. S., & Urbassek, H. M. 2012, ApJ, 752, 151
  * Rozner et al. (2020) Rozner, M., Grishin, E., & Perets, H. B. 2020, MNRAS, 496, 4827
  * San Sebastián et al. (2020) San Sebastián, I. L., Dolff, A., Blum, J., Parisi, M. G., & Kothe, S. 2020, MNRAS, 497, 2418
  * Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69
  * Schaffer et al. (2021) Schaffer, N., Johansen, A., & Lambrechts, M. 2021, A&A, 653, A14
  * Schräpler & Blum (2011) Schräpler, R. & Blum, J. 2011, ApJ, 734, 108
  * Seizinger et al. (2012) Seizinger, A., Speith, R., & Kley, W. 2012, A&A, 541, A59
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shimaki & Arakawa (2012a) Shimaki, Y. & Arakawa, M. 2012a, Icarus, 218, 737
  * Shimaki & Arakawa (2012b) Shimaki, Y. & Arakawa, M. 2012b, Icarus, 221, 310
  * Shimaki & Arakawa (2021) Shimaki, Y. & Arakawa, M. 2021, Icarus, 369, 114646
  * Sirono (2004) Sirono, S.-i. 2004, Icarus, 167, 431
  * Stepinski & Valageas (1997) Stepinski, T. F. & Valageas, P. 1997, A&A, 319, 1007
  * Suyama et al. (2008) Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310
  * Suyama et al. (2012) Suyama, T., Wada, K., Tanaka, H., & Okuzumi, S. 2012, ApJ, 753, 115
  * Tanaka et al. (2023) Tanaka, H., Anayama, R., & Tazaki, R. 2023, ApJ, 945, 68
  * Tanaka et al. (1996) Tanaka, H., Inaba, S., & Nakazawa, K. 1996, Icarus, 123, 450
  * Tatsuuma & Kataoka (2021) Tatsuuma, M. & Kataoka, A. 2021, ApJ, 913, 132
  * Tatsuuma et al. (2023) Tatsuuma, M., Kataoka, A., Okuzumi, S., & Tanaka, H. 2023, ApJ, 953, 6
  * Tazaki & Dominik (2022) Tazaki, R. & Dominik, C. 2022, A&A, 663, A57
  * Tazaki et al. (2023) Tazaki, R., Ginski, C., & Dominik, C. 2023, ApJ, 944, L43
  * Tazaki et al. (2019) Tazaki, R., Tanaka, H., Kataoka, A., Okuzumi, S., & Muto, T. 2019, ApJ, 885, 52
  * Thornton & Ning (1998) Thornton, C. & Ning, Z. 1998, Powder Technology, 99, 154
  * Vericel & Gonzalez (2020) Vericel, A. & Gonzalez, J.-F. 2020, MNRAS, 492, 210
  * Vericel et al. (2021) Vericel, A., Gonzalez, J.-F., Price, D. J., Laibe, G., & Pinte, C. 2021, MNRAS, 507, 2318
  * Verrios et al. (2022) Verrios, H. J., Price, D. J., Pinte, C., Hilder, T., & Calcino, J. 2022, ApJ, 934, L11
  * Vorobyov et al. (2018) Vorobyov, E. I., Akimkin, V., Stoyanovskaya, O., Pavlyuchenkov, Y., & Liu, H. B. 2018, A&A, 614, A98
  * Wada et al. (2013) Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62
  * Wada et al. (2007) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320
  * Wada et al. (2009) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490
  * Wada et al. (2011) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Weidenschilling & Cuzzi (1993) Weidenschilling, S. J. & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine, 1031
  * Weidling et al. (2012) Weidling, R., Güttler, C., & Blum, J. 2012, Icarus, 218, 688
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Williams & Best (2014) Williams, J. P. & Best, W. M. J. 2014, ApJ, 788, 59
  * Windmark et al. (2012) Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73
  * Yamamoto et al. (2014) Yamamoto, T., Kadono, T., & Wada, K. 2014, ApJ, 783, L36
  * Yang et al. (2017) Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
  * Yasui et al. (2017) Yasui, M., Schulson, E. M., & Renshaw, C. E. 2017, Journal of Geophysical Research (Solid Earth), 122, 6014
  * Youdin & Johansen (2007) Youdin, A. & Johansen, A. 2007, ApJ, 662, 613
  * Youdin & Goodman (2005) Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459
  * Zhang et al. (2023) Zhang, S., Zhu, Z., Ueda, T., et al. 2023, ApJ, 953, 96
  * Zhu et al. (2014) Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122
  * Zhu & Yang (2021) Zhu, Z. & Yang, C.-C. 2021, MNRAS, 501, 467
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57

##  Appendix A Equations for porosity evolution during growth

###  A.1 Hit & Stick Regime

In the hit & stick regime, the mass of the grain doubles during a collision,
and a volume of void $V_{\mathrm{void}}$ is formed. According to Okuzumi et
al. (2009), this captured volume $V_{\mathrm{void}}$ corresponds to $99\%$ of
the initial volume $V_{\mathrm{i}}$, thus adding to the initial volumes of the
two grains before the collision. The filling factor $\phi\propto m/V$ is
therefore multiplied by a factor of 2/2.99. For an arbitrary number of
collisions $n$, there is a geometric relationship between the initial filling
factor $\phi_{0}$ and the final $\phi_{\mathrm{f}}$, given by
$\phi_{\mathrm{f}}=(2/2.99)^{n}\phi_{0}$ (Garcia & Gonzalez 2020). Moreover,
the mass of a grain after the same number of collisions has grown by
$m_{\mathrm{f}}=2^{n}m_{0}$, $m_{0}$ being the monomer mass. By isolating $n$
and combining the two expressions, we obtain the filling factor in the hit &
stick regime (assuming $E_{\mathrm{kin}}\ll E_{\mathrm{roll}}$, Garcia &
Gonzalez 2020)

$\phi_{\mathrm{h\&s}}=\left(\frac{m}{m_{0}}\right)^{\psi},$ (40)

with $\psi=\ln\left({2/2.99}\right)/\ln{2}\simeq-0.58$, where for the monomer
$\phi_{0}=1$.

###  A.2 Collisional compression regime

As both the grain mass and the relative velocity $v_{\mathrm{rel}}$ increase,
the kinetic energy increases. Since $v_{\mathrm{rel}}$ depends on the Stokes
number (Sect. 2.1), which in turn depends on the drag regime, there is an
expression for each case (Garcia & Gonzalez 2020). Suyama et al. (2012)
provide the volume of a grain after collision $V_{\mathrm{f}}$ in terms of the
volume of a grain before collision $V_{\mathrm{i}}$. This relation can be
transformed to incorporate filling factors using the fact that
$\phi_{\mathrm{f}}=m_{\mathrm{f}}/V_{\mathrm{f}}\rho_{\mathrm{s}}=2m_{\mathrm{i}}/V_{\mathrm{f}}\rho_{\mathrm{s}}$
(Garcia & Gonzalez 2020)

$\phi_{\mathrm{f}}=\frac{2m_{\mathrm{i}}}{\rho_{\mathrm{s}}}\left(\frac{(3/5)^{5}\left(E_{\mathrm{Kin}}-3b\,E_{\mathrm{roll}}\right)}{N_{\mathrm{tot}}^{5}\,b\,E_{\mathrm{roll}}V_{0}^{10/3}}+\frac{V_{\mathrm{i}}^{-10/3}}{2^{4}}\right)^{3/10}.$
(41)

Here, $N_{\mathrm{tot}}$ is the total number of monomers, which is
$2m_{\mathrm{i}}/m_{0}$, and $V_{0}$ is the volume of a monomer. This model is
recursive, meaning that the evolution of porosity is tracked after each
collision. However, in the case of global simulations, the integration time
step rarely aligns with the time interval between collisions
$\tau_{\mathrm{coll}}$. Therefore, Garcia & Gonzalez (2020) transformed this
discrete model into a continuous one to depend solely on the grain mass and
quantities related to $r$, enabling its implementation in codes for global
simulations. To achieve this, they assumed the approximation
$E_{\mathrm{kin}}\gg E_{\mathrm{roll}}$.

In the Epstein and Stokes regimes for $\mathrm{St}<1$, one obtains

$\displaystyle\phi_{\mathrm{Ep,St<1}}$
$\displaystyle=\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)^{3/8}\left(\frac{m}{m_{0}}\right)^{-1/8},$
(42) $\displaystyle\phi_{\mathrm{St,St<1}}$
$\displaystyle=\left(\frac{Ac_{\mathrm{g}}a_{0}}{9\nu_{\mathrm{mol}}\left(2^{1/5}-1\right)}\right)^{1/3}.$
(43)
$\nu_{\mathrm{mol}}=\frac{5\sqrt{\pi}}{64}\frac{m_{\mathrm{mol}}c_{\mathrm{g}}}{\rho_{\mathrm{g}}\sigma_{\mathrm{mol}}},$
(44)

is the gas molecular kinematic viscosity, where $m_{\mathrm{mol}}$ and
$\sigma_{\mathrm{mol}}$ are the mass and molecular cross section of the H2
molecule, and $A$ is a factor used to simplify the expressions and is given by

$A=\frac{243\sqrt{2}\pi}{15625}\frac{\mathrm{Ro}\,\alpha\,a_{0}^{4}\,\rho_{\mathrm{s}}^{2}c_{\mathrm{g}}\Omega_{\mathrm{K}}}{\rho\,b\,E_{\mathrm{roll}}},$
(45)

where $\rho=\rho_{\mathrm{g}}+\rho_{\mathrm{d}}$ is the total density. It is
observed that in the Stokes regime with $\mathrm{St}<1$, the filling factor
depends only on the local disc conditions. For regimes with $\mathrm{St}>1$,
collisional compression is less efficient due to the decrease in
$v_{\mathrm{rel}}$. By neglecting $E_{\mathrm{kin}}$, it can be derived from
Okuzumi et al. (2012) and Garcia & Gonzalez (2020) that
$\phi_{\mathrm{f}}\propto m_{\mathrm{f}}^{-1/5}$ for $\mathrm{St}>1$.
Therefore, the filling factors for $\mathrm{St}>1$ can be expressed as

$\displaystyle\phi_{\mathrm{Ep,St>1}}$
$\displaystyle=\phi_{\mathrm{Ep,St<1}}(M_{4})\left(\frac{m}{M_{4}}\right)^{-1/5},$
(46) $\displaystyle\phi_{\mathrm{St,St>1}}$
$\displaystyle=\phi_{\mathrm{St,St<1}}(M_{5})\left(\frac{m}{M_{5}}\right)^{-1/5},$
(47)

where $M_{4}$ and $M_{5}$ are the transition masses between the Epstein and
Stokes regimes for $\mathrm{St}<1$ and $\mathrm{St}>1$. The transition masses
are given in Appendix B. The equations giving $\phi$ in the various regimes
can also be rewritten as a function of the size $s$ instead of the mass $m$,
they are listed in Appendix C. However, the natural variable to describe dust
evolution is $m$.

The filling factor resulting from the hit & stick and collisional compression
regimes is noted $\phi_{\mathrm{coll}}$.

###  A.3 Static compression regime

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

$\phi=\left(\frac{a_{0}^{3}P_{\mathrm{comp}}}{E_{\mathrm{roll}}}\right)^{1/3},$
(48)

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

$P_{\mathrm{comp}}=\frac{F_{\mathrm{D}}}{\pi s^{2}}=\frac{m\Delta
v\Omega_{\mathrm{K}}}{\pi s^{2}\mathrm{St}},$ (49)

which allows to obtain the expression for the filling factor
$\phi_{\mathrm{gas}}$ (Garcia & Gonzalez 2020):

$\phi_{\mathrm{gas}}=\left(\frac{m_{0}\,a_{0}}{\pi
E_{\mathrm{roll}}}\frac{\Delta
v\,\Omega_{\mathrm{K}}}{\mathrm{St}}\right)^{3/7}\left(\frac{m}{m_{0}}\right)^{1/7}.$
(50)

For self-gravity, the exerted pressure is given by

$P_{\mathrm{comp}}=\frac{\mathrm{G}m^{2}}{\pi s^{4}},$ (51)

and the expression for the filling factor $\phi_{\mathrm{grav}}$ is (Garcia &
Gonzalez 2020):

$\phi_{\mathrm{grav}}=\left(\frac{\mathrm{G}m_{0}^{2}}{\pi
a_{0}E_{\mathrm{roll}}}\right)^{3/5}\left(\frac{m}{m_{0}}\right)^{2/5},$ (52)

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

###  A.4 Growth from above the compression limit

The minimum filling factor must be $\phi_{\mathrm{min}}$. However, there are
cases where the evolution of a grain (e.g. via bouncing) can produce a filling
factor larger than $\phi_{\mathrm{min}}$, placing it in the space above the
curves showed in Fig. 1. To take into account the increase of porosity of
grains during growth above the limit, we can compute the filling factor as a
function of the initial filling factor $\phi_{\mathrm{i}}$, and the final and
initial masses $m_{\mathrm{f}}$ and $m_{\mathrm{i}}$ in the case where
$v_{\mathrm{rel}}<v_{\mathrm{frag}}$.

$\displaystyle\phi_{\mathrm{grow}}$
$\displaystyle=\phi_{\mathrm{i}}\left(\frac{m_{\mathrm{f}}}{m_{\mathrm{i}}}\right)^{\psi}\quad\text{if
}E_{\mathrm{kin}}<3bE_{\mathrm{roll}},$ (53)
$\displaystyle\phi_{\mathrm{grow}}$
$\displaystyle=\phi_{\mathrm{i}}\left(\frac{m_{\mathrm{f}}}{m_{\mathrm{i}}}\right)^{-1/5}\>\text{otherwise}.$
(54)

The resulting filling factor is then the maximum of $\phi_{\mathrm{grow}}$ and
$\phi_{\mathrm{min}}$.

##  Appendix B Transition masses

Figure 19: Transition masses for our standard disc model and for silicate
grains.

The transition masses are the masses at which the behavior of a grain changes.
The different expressions of $\phi$ found for each regime in Sect. 2.2 are
independent, the transition masses allow us to relate them. For the transition
masses $M_{1}$ and $M_{2}$ that connect the hit & stick regime to the Epstein
and Stokes regimes with $\mathrm{St}<1$, we find

$\displaystyle\frac{M_{1}}{m_{0}}$
$\displaystyle=\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)^{\frac{3/8}{\psi+1/8}},$
(55) $\displaystyle\frac{M_{2}}{m_{0}}$
$\displaystyle=\left(\frac{Ac_{\mathrm{g}}a_{0}}{9\nu_{\mathrm{mol}}\left(2^{1/5}-1\right)}\right)^{1/3\psi}.$
(56)

The transition mass $M_{3}$ corresponds to the transition between the Epstein
and Stokes regimes

$\displaystyle\frac{M_{3}}{m_{0}}$
$\displaystyle=\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)^{3}\left(\frac{Ac_{\mathrm{g}}a_{0}}{9\nu_{\mathrm{mol}}\left(2^{1/5}-1\right)}\right)^{-8/3},$
(57)
$\displaystyle=\left(\frac{M_{1}}{m_{0}}\right)^{8\psi+1}\left(\frac{M_{2}}{m_{0}}\right)^{-8\psi}$
(58)

Finally, the expressions for $M_{4}$ and $M_{5}$ are obtained by equating the
expression of $\phi$ for the Epstein or Stokes regimes with $\mathrm{St}<1$ to
the expression of $\phi$ obtained by setting $\mathrm{St}=1$.

$\displaystyle\frac{M_{4}}{m_{0}}$ $\displaystyle=\left(\frac{\rho
c_{\mathrm{g}}}{\rho_{\mathrm{s}}a_{0}\Omega_{\mathrm{K}}}\right)^{4}\left(\frac{A}{2\left(2^{3/40}-1\right)}\right)^{-1},$
(59) $\displaystyle\frac{M_{5}}{m_{0}}$
$\displaystyle=\left(\frac{9\nu_{\mathrm{mol}}\rho}{2\rho_{\mathrm{s}}a_{0}^{2}\Omega_{\mathrm{K}}}\right)^{3/2}\left(\frac{Ac_{\mathrm{g}}a_{0}}{9\nu_{\mathrm{mol}}\left(2^{1/5}-1\right)}\right)^{-1/6}.$
(60)

Here, we have not considered the other Stokes regimes that are reached for
even larger sizes than the linear regime as static grain compression is always
dominant compared to these regimes. The various transition masses as well as
the different expansion and compression regimes are shown in Fig. 19, inspired
by Fig. A1 of Garcia & Gonzalez (2020). Note that the transition mass $M_{2}$
is always larger than mass $M_{1}$. The Stokes regimes are reached in the
inner regions of a disc, at a distance less than 10 au.

##  Appendix C Filling factor as a function of size

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

$\phi_{\mathrm{h\&s}}=\left(\frac{s}{a_{0}}\right)^{\frac{3\,\psi}{1-\psi}},$
(61)

$\displaystyle\phi_{\mathrm{Ep-St}<1}$
$\displaystyle=\left(\frac{A}{2(2^{3/40}-1)}\right)^{1/3}\left(\frac{s}{a_{0}}\right)^{-1/3},$
(62) $\displaystyle\phi_{\mathrm{St-St}<1}$
$\displaystyle=\left(\frac{A\,c_{\mathrm{g}}\,a_{0}}{9\nu(2^{1/5}-1)}\right)^{1/3},$
(63)

$\displaystyle\phi_{\mathrm{Ep-St}>1}$
$\displaystyle=\left(\frac{\rho\,c_{\mathrm{g}}}{\rho_{\mathrm{s}}\,\Omega_{\mathrm{K}}\,a_{0}}\frac{A}{2^{3/40}-1}\right)^{1/4}\left(\frac{s}{a_{0}}\right)^{-1/2},$
(64) $\displaystyle\phi_{\mathrm{St-St}>1}$
$\displaystyle=\left(\frac{\rho\,c_{\mathrm{g}}\,a_{0}}{2\rho_{\mathrm{s}}\,\Omega_{\mathrm{K}}\,a_{0}^{2}}\,\,\frac{A\,\,}{2^{1/5}-1}\right)^{1/4}\left(\frac{s}{a_{0}}\right)^{-1/2}.$
(65)

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

Figure 20: Same as Fig. 14 but showing the filling factor $\phi$ plotted
against grain size, with colour indicating the distance $r$.

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

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

##  Appendix E Influence of rotational disruption

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

Michoulier & Gonzalez (2022a) studied the effect of rotational disruption on
the evolution of porous grains using 1D simulations with Pamdeas. 3D disc
simulations taking into account the evolution of gas and dust are necessary to
understand in more detail when and where rotational disruption influences dust
evolution. To that effect, we performed simulations modelling this phenomenon
for silicates with $v_{\mathrm{frag,\>Si}}=10$ and 20 $\mathrm{m\,s^{-1}}$,
and for water ice with $v_{\mathrm{frag,\>Si}}=15$ $\mathrm{m\,s^{-1}}$,
listed in the last three rows of Table 1. The grain properties are somewhat
different compared to Michoulier & Gonzalez (2022a): the monomer size is set
to 0.2 instead of 0.1 $\mu$m, and $\gamma_{\mathrm{s}}=0.2\
\mathrm{J\,m^{-2}}$ for silicates. In fact, with these new parameters, the
tensile strength $S_{\mathrm{max}}$ should be larger by a factor of two.

Simulations with Pamdeas showed signs of rotational disruption only for very
small values of the viscosity parameter $\alpha<5\times 10^{-4}$, lower than
that used with Phantom. The conditions necessary for its appearance in the
midplane (which Pamdeas simulates) are therefore challenging to achieve.

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

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