# Grain growth for astrophysics with Discontinuous Galerkin schemes

Maxime Lombart1, Guillaume Laibe1,2
1Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis,-Laval, France.
2Institut Universitaire de France
maxime.lombart<EMAIL_ADDRESS>
(Accepted 2020 November 23. Received 2020 November 19; in original form 2020
October 20.)

###### Abstract

Depending on their sizes, dust grains store more or less charges, catalyse
more or less chemical reactions, intercept more or less photons and stick more
or less efficiently to form embryos of planets. Hence the need for an accurate
treatment of dust coagulation and fragmentation in numerical modelling.
However, existing algorithms for solving the coagulation equation are over-
diffusive in the conditions of 3D simulations. We address this challenge by
developing a high-order solver based on the Discontinuous Galerkin method.
This algorithm conserves mass to machine precision and allows to compute
accurately the growth of dust grains over several orders of magnitude in size
with a very limited number of dust bins.

###### keywords:

methods: numerical — (ISM:) dust, extinction — protoplanetary discs

††pagerange: Grain growth for astrophysics with Discontinuous Galerkin
schemes–Data availability††pubyear: 2020

##  1 Introduction

Solid particles pervade the interstellar medium at all scales. Although they
represent a small amount of its total mass, they deeply influence its
evolution by setting the local chemical, thermal and charge balances. Dust
plays also a key role in the formation of planets, since solid bodies grow
over thirty orders of magnitude in mass to form cores of planets. Spatially
resolved observations of young stellar objects strongly suggest that at least
some planets have to form in less that one million of years (e.g. ALMA
Partnership et al. 2015; Avenhaus et al. 2018; Pinte et al. 2020). Key is to
understand how dust growth can be so efficient. However, planet formation is
an out-of-equilibrium non-linear multi-scales and multi-physics process. For
example, dust grains differentiate from the gas as they settle vertically and
drift radially in the disc (i.e Testi et al. 2014 and references therein).
This creates instabilities which concentrate the solids even more, affecting
the collisional rate of the grains, and thus, their growth or fragmentation.
Since dust dynamics strongly depends on the grain size, growth operates a
strong feed-back on the spatial distribution of the particles.

Figure 1: An illustration of the growth over-diffusion problem: numerical
schemes of order 0 over-estimate the formation of large grains at low
resolution. The plot has been realised with the scheme presented in Kovetz &
Olund (1969) for the case of a constant kernel $K=1$ with $N=15$
logarithmically-spaced dust bins.

Hence, 3D dust/gas simulations that include growth and fragmentation are
compulsory to understand dust evolution during the early stages of planet
formation (e.g. Safronov 1972; Hayashi & Nakagawa 1975; Weidenschilling 1980;
Ohtsuki et al. 1990; Wetherill 1990; Tanaka et al. 1996; Dominik et al. 2007;
Ormel et al. 2007; Birnstiel et al. 2010). The simplest way to formalise the
evolution of a local mass distribution of dust grains is by the mean of the
deterministic mean-field Smoluchowski equation, which assumes binary
collisions (Smoluchowski, 1916). This equation does not have generic analytic
solutions. Integrated non-linearities challenge numerical solvers to obtain
accurate solutions (see Fig. 1). As such, this equation has been thoroughly
studied since a century (e.g. Müller 1928; Schumann 1940; Chandrasekhar 1943;
Melzak 1953; McLeod 1962a; Golovin 1963; Berry 1967; Scott 1968; Trubnikov
1971; Hidy & Brock 1972; Drake 1972; Gillespie 1975b; Silk & White 1978; Silk
& Takahashi 1979; Gelbard et al. 1980; Aldous 1999; Friedlander et al. 2000;
Ramkrishna 2000; Filbet & Laurencot 2004; Jacobson 2005; Pruppacher & Klett
2010), and applied extensively to several fields such aerosols science,
chemistry, meteorology, biology and astrophysics.

It has been shown that classical solvers require a sufficient resolution in
mass to avoid artificial formation of aggregates of large masses (Soong, 1974;
Berry & Reinhardt, 1974; Trautmann & Wanner, 1999; Khain & Pinsky, 2018). This
artificial diffusion may become particularly important when the mass interval
considered is large (Fig. 1). Typically, for planet formation, a few hundreds
of mass bins are required to compute dust growth from interstellar sizes to
pebbles. Usually, this fact is of no importance given current computational
capacities. However, 3D hydrodynamical simulations can hardly handle more than
(a few) ten(s) of mass bins in practice. Compromises have therefore been
performed either by simplifying their growth or their dynamics. However, 1-2 D
hydrodynamical codes integrating the Smoluchowski equation (e.g. Birnstiel et
al. 2010) provide different results compared to 3D hydrodynamical codes with
monodisperse growth models (e.g. Gonzalez et al. 2017), showing the necessity
of a comprehensive approach. This implies to develop a solver which solves
accurately the Smoluchowski equation with a limited number of bins, tractable
by 3D hydrodynamical codes.

Reaching high accuracy with a low number of bins while conserving mass of a
finite interval of mass is a characteristic property of finite volume high-
order solvers, which stem therefore as a natural way to address the growth
over-diffusion problem. In this study, we present a high-order solver for the
Smoluchowski equation based on the Discontinous Galerkin method, following the
pioneering work of Liu et al. (2019). Important properties of the Smoluchowski
equation discussed in the astrophysical context are presented in Sect. 2. The
novel Discontinous Galerkin numerical scheme is presented in Sect. 3. The
performances of the solver regarding the over-diffusion problem are studied in
Sect. 4. Applicability of the algorithm to young stellar objects or in other
astrophysical contexts are discussed in Sect. 5.

##  2 Smoluchowski equation

###  2.1 Short summary

The Smoluchowski equation describes mass conservation for a distribution
aggregates where mass transfers are allowed. This equation exists under a
discret form (monomers forming polymers) or a continuous limit form when mass
quantization becomes negligible (Müller, 1928). The Smoluchowski equation is a
non-linear integro-differential hyperbolic equation that depend on a collision
function called the growth kernel (or kernel) which quantifies the collision
rate between two grains. Explicit solutions exist only for the so-called
constant (Smoluchowski, 1916; Schumann, 1940; Scott, 1968), additive (Golovin,
1963; Scott, 1968) and multiplicative kernels (McLeod, 1962a; Scott, 1968),
implying numerical resolution for physical problems. Among the known
solutions, self-similar solutions are particularly important since they
provide asymptotic behaviour of the mass distribution at large times
(Schumann, 1940; Friedlander & Wang, 1966; Wang, 1966; Menon & Pego, 2004;
Niethammer et al., 2016; Laurençot, 2018). A generic feature of these
solutions is the exponentially fast decay of the solution at large masses.
Gelation, i.e. formation of aggregates of infinite mass form in a finite time
for kernels sustaining explosive growth (Leyvraz & Tschudi, 1981). In
astrophysics, collisions occurs essentially through ballistic impacts
modulated by focusing due to long-range interactions (Safronov, 1972;
Dullemond & Dominik, 2005). Kernels are non-explosive and mass remains
rigorously conserved during the grow process.

###  2.2 Conservative form

Mass conservation for a distribution of growing grains has been originally
formalised by Smoluchowski (1916). Growth is modelled via binary collisions
between spheres having known mean probabilities. The by-products of collisions
are called aggregates or polymers. In Smoluchowski (1916), aggregates are
assumed to also have spherical shapes. Spatial correlations are neglected. The
smallest colliding elements are referred as monomers. For physical systems
involving aggregates made of large numbers of monomers, it is often convenient
to assume continuous mass distributions. The population density of grains
within an elementary mass range $\mathrm{d}m$ is characterised by its number
density $n\\!\left(m\right)$. The continuous Smoluchowski equation is given by

$\displaystyle\frac{\partial n\left(m,t\right)}{\partial t}=$
$\displaystyle\frac{1}{2}\int\limits_{0}^{m}\\!K\\!\left(m-m^{\prime},m^{\prime}\right)n\\!\left(m-m^{\prime},t\right)n\\!\left(m^{\prime},t\right)\mathrm{d}m^{\prime}$
(1)
$\displaystyle-n\\!\left(m,t\right)\int\limits_{0}^{\infty}\\!K\\!\left(m,m^{\prime}\right)n\\!\left(m^{\prime},t\right)\mathrm{d}m^{\prime},$

where $t$ denotes time and $m$ and $m^{\prime}$ the masses of two colliding
polymers. The averaged probabilities of collision are encoded inside the
coagulation kernel $K\left(m,m^{\prime}\right)$, which is a symmetric function
of $m$ and $m^{\prime}$ for binary collisions (see Sect. 2.3). Fig. 2 shows
the physical meaning of the non-linear integro-differential equation Eq. 1.
The number of grains encompassed within a given interval of masses varies
since i) binary collisions of aggregates of appropriate masses can increase
this population (first term of the right-hand side of Eq. 1), but ii) those
grains may themselves collide with other grains to form larger aggregates
(second term of the right-hand side of Eq. 1). This equation can be put under
a convenient dimensionless form by introducing (Scott, 1968; Drake, 1972)

$\left\\{\begin{aligned} &x\equiv m/m_{0},\,y\equiv
m^{\prime}/m_{0},\,\mathcal{K}(x,y)=K(m,m^{\prime})/K_{0},\\\
&\tau=(K_{0}N_{0})t,\,f(x,\tau)=m_{0}\,n(m,t)/N_{0}.\end{aligned}\right.$ (2)

$N_{0}$ is the initial total number density of particles, $m_{0}$ is the
initial mean mass of the particles and $K_{0}$ is a normalising constant with
dimensions $[\mathrm{length}]^{3}/\mathrm{time}$. We adopt the variables $x$
and $\tau$ for sake of clarity and homogeneity with the existing literature
(e.g. Friedlander et al. 2000, and references therein; Jacobson 2005, and
references therein). $x$ denotes therefore masses. Eq. 1 transforms into

$\displaystyle\frac{\partial f(x,\tau)}{\partial\tau}=$
$\displaystyle\frac{1}{2}\int\limits_{0}^{x}\\!\mathcal{K}(y,x-y)f(y,\tau)f(x-y,\tau)\mathrm{d}y$
(3)
$\displaystyle-f(x,\tau)\int\limits_{0}^{\infty}\\!\mathcal{K}(y,x)f(y,\tau)\mathrm{d}y.$

Eq. 3 is physically ill-posed, since the probability to form aggregates of
mass larger than the initial mass of the system may be non-zero. Recently,
Tanaka et al. (1996) have shown that Eq. 3 can be equivalently written under
the conservative form

$\left\\{\begin{aligned} &\frac{\partial
g\left(x,\tau\right)}{\partial\tau}+\frac{\partial
F_{\mathrm{coag}}\left[g\right]\left(x,\tau\right)}{\partial x}=0\\\
&F_{\mathrm{coag}}\left[g\right]\left(x,\tau\right)=\int\limits_{0}^{x}\\!\\!\int\limits_{x-u}^{\infty}\mathcal{K}\left(u,v\right)g\left(u,\tau\right)\frac{g\left(v,\tau\right)}{v}\mathrm{d}u\mathrm{d}v,\end{aligned}\right.$
(4)

where $g\left(x,\tau\right)\equiv xf\left(x,\tau\right)$ is the mass density
of polymers per unit mass, and
$F_{\mathrm{coag}}\left[g\right]\left(x,\tau\right)$ is the flux of mass
density across the mass $x$ triggered by coagulation (Filbet & Laurencot,
2004). Under this conservative form, the infinite upper bound of the second
integral in $F_{\mathrm{coag}}$ can simply be replaced by $x_{\rm max}-u$.
This prevents the formation of aggregates of masses larger than $x_{\rm max}$
by settling the passing-through mass flux to be rigorously zero.

###  2.3 Kernels

Physically, the coagulation kernel is defined according to

$K\\!\left(m,m^{\prime}\right)\equiv\beta\\!\left(m,m^{\prime},\Delta
v\right)\Delta
v\\!\left(m,m^{\prime}\right)\sigma\\!\left(m,m^{\prime}\right),$ (5)

where $\Delta v$ is the mean relative velocity between two aggregates of
masses $m$ and $m^{\prime}$, $\sigma$ is the mean effective cross section of
collision and $\beta$ denotes the mean sticking probability of the grains. The
coagulation kernel encodes the microphysics of collisions inside $\beta$,
$\sigma$ and $\Delta v$, those parameters depending a priori on the sizes of
the colliding grains, or the kinetic and thermodynamical parameters of an
eventual surrounding flow. A kernel of particular importance for physical
problems is the Ballistic kernel (Table 1). In this case, $\sigma$ corresponds
simply to the geometric cross-section of the grains (focusing effects due to
electrostatic or gravitational forces being neglected), and $\beta$ and
$\Delta v$ are treated as constants (which may be a relevant approximation at
least over moderate ranges of masses). Coagulation kernel can also be seen as
mathematical objects useful to study the properties of the Smoluchowski
equation under various conditions or to derive explicit analytic solutions.
The expression of the four kernels discussed in this work is given in Table 1.

Figure 2: Illustration of the Smoluchowski equation Eq. 1. Polymers of mass
$m_{i}$ are represented in orange. The green and blue polymers have masses
lower than $m_{i}$. Creation (resp. growth) of polymers of mass $m_{i}$
increases (resp. decreases) its number density.

###  2.4 Analytic solutions

Explicit analytic solutions exist in the case of simple kernels and specific
initial conditions. We review these solutions hereafter since they will be
used in Sect. 4 to benchmark the numerical algorithms.

####  2.4.1 Constant kernel

For the constant kernel $\mathcal{K}\\!\left(x,y\right)=1$ and the initial
condition $f\left(x,0\right)=\exp\left(-x\right)$, the solution of Eq. 3 is
(Müller, 1928; Schumann, 1940; Melzak, 1957; Rajagopal, 1959; Scott, 1968;
Silk & Takahashi, 1979)

$\left\\{\begin{aligned}
&f_{1}(\tau)\equiv\frac{4}{(2+\tau)^{2}},\,f_{2}(\tau)\equiv\frac{\tau}{2+\tau},\\\
&f(x,\tau)=f_{1}(\tau)\exp\left(-\left\\{1-f_{2}(\tau)\right\\}x\right).\\\
\end{aligned}\right.$ (6)

Physically, a constant kernel $\mathcal{K}=1$ implies that the frequency of
collisions between two aggregates is independent of their size.

Kernel | $\mathcal{K}(x,y)$
---|---
Size-independent | $1$
Sum | $x+y$
Product | $xy$
Ballistic | $\pi\left(x^{1/3}+y^{1/3}\right)^{2}\Delta v$
Table 1: Functional form of the different coagulation kernels $\mathcal{K}$
considered in this study.

####  2.4.2 Additive kernel

The solution for the additive kernel $\mathcal{K}(x,y)=x+y$ with the initial
condition $f_{0}(x,0)=\exp(-x)$ has been derived by Golovin (1963). Scott
(1968) extended the derivation for a general initial condition. For an initial
condition under the form $f(x,0)=\exp(-x)$, the solution of Eq. 3 is

$\left\\{\begin{aligned} &T\equiv 1-\exp(-\tau),\\\
&f(x,\tau)=\frac{\left(1-T\right)\exp\left(-x\left\\{1+T\right\\}\right)}{xT^{1/2}}I_{1}\\!\left(2xT^{1/2}\right),\\\
\end{aligned}\right.$ (7)

where $I_{1}$ is the modified Bessel function of first kind. Physically, the
additive kernel implies that the frequency of collisions increases according
to the size of the grains. Large aggregates form faster compared to case of a
constant kernel, leading to broader dust distributions at large masses. The
asymptotic tail presents therefore a smoother decay compared to the case
$\mathcal{K}=1$.

####  2.4.3 Multiplicative kernel

Originally, McLeod (1962b) derived a solution for the multiplicative kernel
$\mathcal{K}(x,y)=xy$ with the initial condition $f_{0}(x,0)=x^{-1}\exp(-x)$
only for a small interval of time. The general solution for this problem was
later found by Ernst et al. (1984)

$\left\\{\begin{aligned} &T\equiv\left\\{\begin{aligned}
&1+\tau\quad\mathrm{if}\,\,\,\,\tau\leq 1\\\
&2\tau^{1/2}\quad\mathrm{otherwise}\end{aligned},\right.\\\
&f(x,\tau)=\frac{\exp\left(-Tx\right)I_{1}\\!\left(2x\tau^{1/2}\right)}{x^{2}\tau^{1/2}}.\\\
\end{aligned}\right.$ (8)

The multiplicative kernel is a typical kernel to study the occurrence of
gelation, since at $\tau=1$, aggregates with infinite masses form and mass
conservation is mathematically no longer satisfied. Physically, the
multiplicative kernel means an explosive increase of the collisional
frequencies with respect to grain sizes. Massive grains form faster compared
to the case of the additive kernel. In the same time, the mass density of
small grains decreases quickly.

###  2.5 Numerical methods

No known analytic solutions exist for the Smoluchowski coagulation equation
with physical kernels, implying numerical resolution. Various numerical
schemes have been developed for this purpose. Two classes of algorithms have
been developed. A first class of solvers consists of Monte-Carlo simulations
(e.g. Gillespie 1975a; Liffman 1992; Smith & Matsoukas 1998; Lee & Matsoukas
2000; Debry et al. 2003; Sheng & Shen 2006; Ormel et al. 2007; Zsom &
Dullemond 2008). Although convenient, these methods have two principal
drawbacks. Firstly, a large number of particles is required to ensure
appropriate accuracy of the number density distribution $f$. Secondly, the
scheme is not deterministic and simulations can be reproduced only in a
statistical sense, which is not satisfying when interfacing with
hydrodynamics. A second class of solvers consist of deterministic algorithms.
These methods have been summarised in Kostoglou & Karabelas (1994); Kumar &
Ramkrishna (1996); Ramkrishna (2000); Pruppacher & Klett (2010); Khain &
Pinsky (2018). A short but comprehensive summary is given hereafter.

####  2.5.1 Method of moments

The method of moments seems to be the first numerical method proposed to solve
the Smoluchowski equation (Hulburt & Katz, 1964). A system of ordinary
differential equations is written over the $k$th moments
$M_{k}\equiv\int_{0}^{\infty}x^{k}f(x,\tau)\mathrm{d}x$ of the number density
function. Approximations either for the reconstruction of $f$ (Hulburt & Katz,
1964) or for the derivation of fractional moments (Estrada & Cuzzi, 2008) are
then required to close this system of ordinary differential equations. The
Standard Moment Method (SMM) requires an analytical integration of the kernel.
To avoid this difficulty, Quadrature Moment Methods (QMM), where integrals are
approximated by Gaussian quadrature methods, have been developed. Solutions of
moments can be used directly to derive the total number of particles $M_{0}$,
the total mass $M_{1}$ or other physical quantity such as dust opacities
(Marchisio et al., 2003; Estrada & Cuzzi, 2008). Number densities $f$ are
reconstructed using polynomials (Pruppacher & Klett, 1980; Piskunov & Golubev,
2002).

####  2.5.2 Point-based methods

The number density function $f$ is sampled over a mass grid. The main
difficulty lies in representing the continuous distribution $f$ as accurately
as possible using the values of $f$ at the sampling points. Different
algorithms have been developed using this approach:

##### Interpolation method

This method was developed by Berry (1967); Berry & Reinhardt (1974). The
continuous Smoluchowski equation is written in terms of $g(x,\tau)\equiv
xf(x,\tau)$, the mass density function. The mass interval is discretised using
a logarithmic grid. A system of ordinary differential equations is derived
with respect to the variable $g$ evaluated on the grid points. Gain and loss
terms are evaluated separately, and integrals are calculated by using high-
order Lagrangian interpolations. Middleton & Brock (1976); Suck & Brock (1979)
improved this method by using Simpson’s rules for the integrals and cubic
splines interpolations.

##### Method of orthogonal collocation

The method of weighted residuals (Finlayson, 1972) is a general method for
obtaining numerical solutions to differential equations. The unknown solution
is tested over a set of weight functions and is adapted to give the best
approximated solution to the differential equation. The Smoluchowski equation
is multiplied by the weight function $\phi$ and integrated over all the mass
domain to form the residual

$\displaystyle R\equiv\int_{0}^{\infty}$ $\displaystyle\left(\frac{\partial
f(x,\tau)}{\partial\tau}\right.-\int_{0}^{x}\mathcal{K}(x-y,y)f(x-y,\tau)f(x,\tau)\mathrm{d}y$
(9)
$\displaystyle\left.+\int_{0}^{\infty}\mathcal{K}(x,y)f(x,\tau)f(y,\tau)\mathrm{d}y\right)\phi(x)\mathrm{d}x=0.$

The number density $f$ is approximated by polynomials. The collocation method
corresponds to the case where $\phi(x)=\delta(x-x_{0})$. The coagulation
equation is evaluated at the collocation points $x_{0}$. This gives a set of
ordinary differential equations equal to the degree of freedom of the
polynomials used. Integrals are usually performed using Gaussian quadrature
rules (Eyre et al., 1988).

Figure 3: Illustration of the pair interaction methods. A particle of mass
$x_{n+l}=x_{n}+x_{l}$ forms from collision between particles of masses $x_{l}$
and $x_{n}$. The resulting mass $x_{n+l}$ is distributed onto adjacent bins,
generating numerical over-diffusion towards large masses.

##### Pair interaction methods

Numerical integration of the Smoluchowski equation consists of summing
contributions of pairwise collisions between all grid points of different
masses. For non-regular mass samplings, aggregates do usually not have masses
corresponding to an existing grid point. To ensure mass conservation, the mass
of the aggregate is distributed over the two relevant adjacent grid points
(Fig. 3). The first pair-interaction solver has been developed by Kovetz &
Olund (1969). In this algorithm, a system of ordinary differential equations
is obtained over the quantities
$N(x_{i})=\int_{a_{i}}^{b_{i}}\\!f(x)\mathrm{d}x$ where $x_{i}$ denotes the
mass of individual particles of the $i$-th point, and
$a_{i}\equiv(x_{i+1}-x_{i})/2$ and $b_{i}\equiv(x_{i}-x_{i-1})/2$. In
practice, logarithmic grids are used to cover wide ranges of masses. In the
context of planet formation, widely used solvers follow this approach (e.g.
Brauer et al. 2008; Birnstiel et al. 2010). The principal drawback of this
method is that redistribution of mass towards large grains tend to over-
predict the number of large aggregates, triggering artificial formation of
large bodies (Fig. 1). A large number of grid points is therefore required to
avoid an artificial broadening of number density of particles $f$ (Berry &
Reinhardt, 1974; Soong, 1974; Khain & Pinsky, 2018). Moreover, a sufficient
number of grid points is also needed to avoid difficulties related to
collisions that form aggregates of masses larger than the largest mass point.
Jacobson (2005) extended the Kovetz & Olund (1969) algorithm by distributing
the mass between grid points and writing the scheme in a semi-implicit form.
This solver ensures mass conservation to machine precision. Bott (1998);
Simmel et al. (2002); Wang et al. (2007) developed also binary-pairs
interaction methods. Mass is advected towards adjacent grid points by a mass
flux expressed with a high-order scheme. These methods do not introduce a
significant numerical broadening. Other methods have been developed by
Hounslow et al. (1988); Lister et al. (1995) where four binary interaction
mechanisms of gain and loss of particles are considered to deal correctly the
rate of change of particle and mass.

####  2.5.3 Finite element methods

In these methods, the continuous mass distribution is discretised over a
finite number of mass elements (intervals, cells, bins).

##### Moments with finite elements

The first finite element scheme for coagulation was developed by Bleck (1970)
by discretising mass distributions over logarithmic bins. $f$ is approximated
by its moment of order zero over each bin to obtain a system of ordinary
differential equations. Over-diffusion for large grains is observed with this
piecewise constant approximation. A change of variable $x\rightarrow x^{-3}$
is operated to reduce diffusivity at large masses. The method of Soong (1974)
follows Bleck (1970). The Smoluchowski equation is written in terms of mass
density distributions $g$ and approximated by piecewise exponential functions.
This allows to reduce drastically the diffusive effect at large masses.
Gelbard et al. (1980); Landgrebe & Pratsinis (1990) proposed a similar method,
where the Smoluchowski equation is decomposed over bins of indices $j$ in
terms of $Q_{j}=\int_{I_{j}}xf(x,\tau)\mathrm{d}x$. A precise account of gain
and loss of particles in terms of fluxes of $Q$ is performed. Trautmann &
Wanner (1999) extends the work of Gelbard et al. (1980), also finding
numerical diffusion when using piecewise constant approximation, and
addressing it by using piecewise exponential approximations. Another moment
method that involves polynomial approximations for the first two moments
$M_{0}$ and $M_{1}$ of $f$ has been proposed by (Enukashvily, 1980; Kumar &
Ramkrishna, 1996; Tzivion et al., 1999).

##### Discontinuous Galerkin method

The discontinuous Galerkin method is a weighted residual method where the
weight $\phi(x)$ consists of orthogonal polynomials (Lagrange polynomials,
Legendre polynomials, cubic splines). The numerical solution of the
Smoluchowski equation is decomposed on each bin over this basis and a system
of ordinary differential equations is obtained for the coefficients (e.g.
Pilinis 1990; Erasmus et al. 1994; Mahoney & Ramkrishna 2002, see Sect. 3).
Generally, the integrals are performed by Gaussian quadrature rules (Gelbard &
Seinfeld, 1978; Rigopoulos & Jones, 2003; Sandu, 2006).

####  2.5.4 Finite element schemes in the conservative form

The conservative form Eq. 4 has been exploited for numerical simulations only
lately. Filbet & Laurencot (2004) derived a finite volume scheme of order zero
where volume integrals over flux divergences are replaced by flux terms at the
interfaces by the mean of the divergence theorem. This scheme conserves mass
exactly and has been further extended by (Filbet, 2008; Bourgade & Filbet,
2008; Forestier-Coste & Mancini, 2012). The mass interval can be sampled
uniformly or non-uniformly. Finite volume schemes of higher orders solving for
the conservative form have been investigated recently (Gabriel & Tine, 2010;
Liu et al., 2019). Gabriel & Tine (2010) used WENO reconstruction (Jiang &
Peng, 2000) to approximate the coagulation flux at interfaces. Liu et al.
(2019) developed a numerical scheme based on the discontinuous Galerkin
method. This method provides the further advantage to choose the order of the
scheme in a flexible manner. Integrals are calculated using Gaussian
Quadrature rules, which implies sub-sampling of the mass intervals.

Figure 4: Sketch of the discontinuous Galerkin method. In each cell, the
solution is approximated by high-order polynomials $k$ to increase accuracy.

###  2.6 Requirements from hydrodynamical simulations

Densities must remain strictly positive and total mass conserved rigorously to
ensure the stability of hydrodynamical simulations. These two properties are
genuinely ensured by finite volume methods based on the conservative form Eq.
4. The double-integral formulation allows to simply quench the formation of
aggregates with unphysical masses, by setting for the integral bound the
maximum mass allowed. These constrains may not always be satisfied with simple
integral formulations.

On the other hand, observational constrains on young stellar objects are
essentially provided by high-contrast spectro-polarimetry at infrared
wavelengths (SPHERE/VLT, GPI, Subaru/HiCIAO) and millimetre interferometry
(ALMA). These observations probe (sub)micron-to-millimetre-in-size dust
distributions in discs, which corresponds to 4 orders of magnitude in size,
i.e. 12 orders of magnitude in mass for compact grains. With current
computational capacities, 3D dust/gas simulations of dusty discs can handle
$\sim$10-20 dust species simultaneously (e.g. PHANTOM, Price et al. 2018 or
RAMSES, Lebreuilly et al. 2020). The global accuracy of second-order
hydrodynamical solvers is of order $\sim 10^{-3}$. We aim therefore to design
a versatile algorithm for coagulation of accuracy $\sim 10^{-3}$ with $\sim
15$ dust bins distributed over 12 orders of magnitude in mass that allows
tractable simulations. We therefore face the issue of over-diffusion
associated to piecewise constant reconstructions with few mass bins, and high-
order schemes appear as a natural way to overcome this difficulty. It is much
preferable for hydrodynamics to handle a fix grid of sizes, to avoid
interpolations when updating forces. We seek therefore for a growth algorithm
that works efficiently with a fixed grid.

Additionally, we seek for an algorithm which allows for convergence analysis
in 3D hydrodynamical simulations. As explained above, multiplying the number
of dust bins provides prohibitive computational costs. Instead, the order of
the scheme may be varied, should it be parametrised in a flexible manner. This
requirement tends to favour Discontinuous Galerkin schemes with respect to
WENO schemes, although they provide in theory equivalent accuracies. Compared
to regular Galerkin schemes, discontinuous Galerkin solvers decompose the
solution over several mass bins. This helps to better capture the exponential
decay of the solution at large masses and avoid over-diffusion biases. For
these reasons, we have chosen to focus on the Discontinuous Galerkin method to
solve for the Smoluchowski equation in astrophysical contexts, an approach
recently pioneered by Liu et al. (2019).

Monofluid dust/gas hydrodynamical solvers provide a natural architecture to
include a coagulation equation. Indeed, relative drifts between grains of
different sizes are genuinely computed, eventually in the terminal velocity
approximation (e.g. Laibe & Price 2014; Hutchison et al. 2018; Lebreuilly et
al. 2019). Monofluid formalism also ensures exact conservation of momentum,
i.e. no thrust due to mass transfers propel the mixture. Sub-grid fluctuations
should be prescribed by an accurate model that describes local turbulence or
Brownian motion.

##  3 Discontinuous Galerkin algorithm

###  3.1 Discontinuous Galerkin method

The discontinuous Galerkin method is presented for the general scalar
hyperbolic conservative equation

$\left\\{\begin{aligned} &\frac{\partial
g(x,\tau)}{\partial\tau}+\frac{\partial F[g](x,\tau)}{\partial x}=0,\\\
&(x,\tau)\in\mathbb{R}_{+},\end{aligned}\right.$ (10)

where $g$ is a density of a conservative quantity and $F\\!\left[g\right]$ the
associated flux.

Let partition the domain of interest
$[x_{\mathrm{min}},x_{\mathrm{max}}]\in\mathbb{R}$ in $N$ subintervals
(alternatively, cells or bins), not necessarily of equal sizes. Each cell is
defined by $I_{j}=(x_{j-1/2},x_{j+1/2}],j\in[\\![1,N]\\!]$. The size of the
$j$-th cell is defined as $h_{j}=x_{j+1/2}-x_{j-1/2}$. The cell is centred
around the position $x_{j}=\left(x_{j+1/2}+x_{j-1/2}\right)/2$. We define
$\mathcal{V}^{k}$ the space of polynomials of degree $k$ in each cell $I_{j}$

$\mathcal{V}^{k}=\left\\{v:v|_{I_{j}}\in
P^{k}\left(I_{j}\right),j\in[\\![1,N]\\!]\right\\}.$ (11)

We denote $g_{j}\in\mathcal{V}^{k}$ the approximate solution of $g$ in the bin
$I_{j}$. The terminology discontinuous Galerkin (DG) comes from the fact that
in $\mathcal{V}^{k}$, the functions are allowed to have jumps at the
interfaces $x_{j+1/2}$. One obtains a weak formulation of Eq. 10 by
multiplying by a test function $\phi\in\mathcal{V}^{k}$, integrating over
$I_{j}$ and finally integrating by parts (Cockburn & Shu, 1989)

$\displaystyle\int_{I_{j}}\frac{\partial g_{j}}{\partial t}\phi\mathrm{d}x$
$\displaystyle-\int_{I_{j}}F[g]\left(x,t\right)\frac{\partial\phi}{\partial
x}\mathrm{d}x$ (12)
$\displaystyle+F[g]\left(x_{j+1/2},t\right)\phi(x_{j+1/2})$
$\displaystyle-F[g]\left(x_{j-1/2},t\right)\phi(x_{j-1/2})=0.$

Eq. 12 allows to fix unequivocally the degrees of freedom of the function
$g_{j}$. The residual of Eq. 10 on bin $I_{j}$ is defined as

$\displaystyle R_{j}\equiv\int_{I_{j}}\frac{\partial g_{j}}{\partial
t}\phi\mathrm{d}x$
$\displaystyle-\int_{I_{j}}F[g]\left(x,t\right)\frac{\partial\phi}{\partial
x}\mathrm{d}x$ (13)
$\displaystyle+F[g]\left(x_{j+1/2},t\right)\phi(x_{j+1/2})$
$\displaystyle-F[g]\left(x_{j-1/2},t\right)\phi(x_{j-1/2}).$

DG schemes consist of choosing a local orthogonal polynomials basis on $I_{j}$
to replace the test function and to approximate the solution. Residuals
$R_{j}$ are therefore null in the sense of orthogonalisation on the basis. In
practice, Legendre polynomials are used (Cockburn & Shu, 1989). We denote
hereafter the $i$-th Legendre polynomial by $\phi_{i}\left(\xi\right)$, where
$\xi\in[-1,1]$. Polynomial functions $\phi_{i}\left(\xi\right)$ are orthogonal
in $L^{2}\left([-1,1]\right)$ with respect to the inner product with weight
unity. Fig. 4 shows a sketch of the DG method. In each cell, the function $g$
is approximated by Legendre polynomials. The accuracy of the approximation
increases with respect to the order of the polynomials. The approximation of
$g$ in cell $I_{j}$ writes

$\displaystyle\forall x\in I_{j},\;g(x)\approx
g_{j}\left(x,t\right)=\sum_{i=0}^{k}g_{j}^{i}\left(t\right)\phi_{i}(\xi_{j}\left(x\right)),$
(14) $\displaystyle
g_{j}\left(x,t\right)=\bm{g}^{T}_{j}(t)\cdot\bm{\phi}(\xi_{j}(x)),\;\mathrm{with}\;\bm{g}_{j}=\begin{bmatrix}g_{j}^{0}\\\
\vdots\\\
g_{j}^{k}\end{bmatrix}\mathrm{and}\;\bm{\phi}=\begin{bmatrix}\phi_{0}\\\
\vdots\\\ \phi_{k}\end{bmatrix},$

where $g_{j}^{i}$ is the component of $g_{j}$ on the Legendre polynomials
basis. The function
$\xi_{j}\left(x\right)\equiv\frac{2}{h_{j}}\left(x-x_{j}\right)$ is used to
map the interval $I_{j}$ onto the interval $[-1,1]$. Normalising the Legendre
basis gives

$\int\limits_{-1}^{1}\bm{\phi}(\xi)\bm{\phi}^{T}(\xi)\mathrm{d}\xi=d_{i}\delta_{ik}\;\mathrm{with}\;d_{i}\equiv\frac{2}{2i+1},$
(15)

where $d_{i}$ is the coefficient of normalisation. By combining Eqs. 12, 14
and 15 one obtains

$\displaystyle\frac{\mathrm{d}\bm{g}_{j}\left(t\right)}{\mathrm{d}t}=\bm{L}[g]\;\mathrm{with}$
(16) $\displaystyle\begin{aligned}
\bm{L}[g]\equiv&\frac{2}{h_{j}}\begin{bmatrix}1/d_{0}&&\\\ &\ddots&\\\
&&1/d_{i}\end{bmatrix}\\\ &\begin{aligned}
\Bigg{(}\Bigg{.}\int_{I_{j}}&F\left[g\right]\left(x,t\right)\partial_{x}\bm{\phi}\left(\xi_{j}\left(x\right)\right)\mathrm{d}x\\\
&\begin{aligned}
-\bigg{[}\bigg{.}&F\left[g\right]\left(x_{j+1/2},t\right)\bm{\phi}\left(\xi_{j}\left(x_{j+1/2}\right)\right)\\\
&\Bigg{.}\bigg{.}-F\left[g\right]\left(x_{j-1/2},t\right)\bm{\phi}\left(\xi_{j}\left(x_{j-1/2}\right)\right)\bigg{]}\Bigg{)},\end{aligned}\end{aligned}\end{aligned}$

where $\bm{L}$ is the operator that results from applying the DG procedure to
Eq. 10 with a Legendre polynomials basis. With the procedure described above,
the original system of partial differential equations (PDE) Eq. 12 is
transformed into a system of ordinary differential equations (ODE) Eq. 16 onto
the coefficients $g_{j}^{i}(t)$. The initial condition $g_{j}\left(x,0\right)$
is generated by the piecewise $L^{2}$ projection of an initial mass density
distribution $g_{0}(x)$ on each bin, i.e.

$\displaystyle\forall j\in[\\![1,N]\\!],$ (17)
$\displaystyle\int_{I_{j}}\left(g_{j}\left(x,0\right)-g_{0}\left(x\right)\right)\bm{\phi}^{T}(\xi_{j}(x))\mathrm{d}x=\bm{0}.$

Orthogonality of Legendre polynomials ensures

$\displaystyle\int_{I_{j}}g_{j}\bm{\phi}^{T}\mathrm{d}x$
$\displaystyle=\frac{h_{j}}{2}\int_{-1}^{1}\bm{\phi}(\xi)\bm{\phi}^{T}(\xi)\mathrm{d}\xi\bm{g}_{j}(t)$
(18)
$\displaystyle=\frac{h_{j}}{2}\mathrm{diag}[d_{0},...,d_{k}]\bm{g}_{j}(t).$

Then, the components of $\bm{g}_{j}$ are given by

$\displaystyle\forall j\in[\\![1,N]\\!],\forall i\in[\\![0,k]\\!],$ (19)
$\displaystyle
g_{j}^{i}(0)=\frac{2}{h_{j}d_{i}}\int\limits_{-1}^{1}g_{0}\left(\frac{h_{j}}{2}\xi_{j}+x_{j}\right)\phi_{i}(\xi_{j})\mathrm{d}\xi_{j}.$

Hence, the DG method consists in solving the following Cauchy problem

$\left\\{\begin{aligned} &\forall j\in[\\![1,N]\\!],\forall
i\in[\\![0,k]\\!],\\\
&\frac{\mathrm{d}\bm{g}_{j}\left(t\right)}{\mathrm{d}t}=\bm{L}[g],\\\
&g_{j}^{i}(0)=\frac{2}{h_{j}d_{i}}\int\limits_{-1}^{1}g_{0}\left(\frac{h_{j}}{2}\xi_{j}+x_{j}\right)\phi_{i}(\xi_{j})\mathrm{d}\xi_{j},\end{aligned}\right.$
(20)

where $\bm{L}$ is detailed in Eq. 16.

###  3.2 Evaluation of the flux

####  3.2.1 Regularised flux

The continuous Smoluchowski coagulation Eq. 3 is defined over an unbounded
interval of masses $x\in\mathbb{R}_{+}$. Before applying the DG procedure, Eq.
3 is restrained to a physical mass interval. Moreover, growth from a gaseous
reservoir is excluded, meaning that $x>0$. The mass interval is therefore
reduced to the interval $[x_{\mathrm{min}}>0,x_{\mathrm{max}}<+\infty]$
(Filbet & Laurencot, 2004; Liu et al., 2019). The coagulation flux can be
truncated according to two procedures (Filbet & Laurencot, 2004). On the one
hand

$\displaystyle
F_{\mathrm{coag}}^{\mathrm{c}}\left[g\right]\left(x,\tau\right)=$ (21)
$\displaystyle\qquad\int\limits_{x_{\mathrm{min}}}^{x}\\!\int\limits_{x-u+x_{\mathrm{min}}}^{x_{\mathrm{max}}-u+x_{\mathrm{min}}}\mathcal{K}\left(u,v\right)g\left(u,\tau\right)\frac{g\left(v,\tau\right)}{v}\mathrm{d}v\mathrm{d}u,$

where $F_{\mathrm{coag}}^{\mathrm{c}}$ is the conservative flux, meaning that
no particle of mass larger than $x_{\mathrm{max}}$ is allowed to form. On the
other hand

$F_{\mathrm{coag}}^{\mathrm{nc}}\left[g\right]\left(x,\tau\right)=\int\limits_{x_{\mathrm{min}}}^{x}\\!\int\limits_{x-u+x_{\mathrm{min}}}^{x_{\mathrm{max}}}\mathcal{K}\left(u,v\right)g\left(u,\tau\right)\frac{g\left(v,\tau\right)}{v}\mathrm{d}v\mathrm{d}u,$
(22)

where $F_{\mathrm{coag}}^{\mathrm{nc}}$ is the non-conservative flux which
allows formation of particles of mass $x>x_{\mathrm{max}}$.
$F_{\mathrm{coag}}^{\mathrm{c}}$ is useful in realistic simulations of growth,
whereas $F_{\mathrm{coag}}^{\mathrm{nc}}$ should be used to compare numerical
solution to analytic solutions of Eq. 1.

####  3.2.2 Method for evaluating the flux

A crucial difference between this scheme and usual DG solvers is that the
coagulation flux $F_{\mathrm{coag}}^{\mathrm{nc}}$ is non local. The
evaluation of the numerical flux $F_{\mathrm{coag}}^{\mathrm{nc}}[g]$ at the
interface $x_{j+1/2}$ depends on the evaluation of $g_{j}$ in all cells.
Mathematically, $F_{\mathrm{coag}}^{\mathrm{nc}}$ is a double integral of a
product of polynomials. Then the flux is a continuous function of mass $x$. At
the interface $x_{j+1/2}$, the numerical flux reduces to
$F_{\mathrm{coag}}^{\mathrm{nc}}\left[g\right]=F_{\mathrm{coag}}^{\mathrm{nc}}\left[g\right]\left(x_{j+1/2},t\right)$.
In usual DG solvers, the numerical flux is a discontinuous function and must
be reconstructed at the interfaces (e.g Cockburn & Shu 1989; Zhang & Shu
2010).

The principal difficulty lies in carefully evaluating the flux at interfaces.
This relies on handling the numerical integration of the polynomials $g_{j}$
in every relevant cell. Liu et al. (2019) uses a Gaussian quadrature method
with a Legendre polynomials basis to approximate the flux. The lower bound of
the inner integral $x-u$ does usually not correspond to a grid point. To
accurately perform the Gauss quadrature, some grid elements must be sub-
divided, increasing drastically the cost of the numerical procedure,
especially for high-order polynomials. To avoid prohibitive computational
costs due to cell oversampling, we take advantage of the polynomial
approximation by calculating integrals analytically. This requires integrable
kernels, which is the case for the four kernels presented in this study. This
approach maintains a reasonable computational cost by not multiplying the
number of sampling points. This also avoid to add errors due to the numerical
integration and to approximate kernels by piecewise constant functions.

####  3.2.3 Mathematical procedure

To integrate analytically the numerical flux, let define the function
$\tilde{g}$ that approximates the function $g$ over the entire mass interval

$\displaystyle\forall x\in[x_{\mathrm{min}},x_{\mathrm{max}}],$ (23)
$\displaystyle\tilde{g}\left(x,\tau\right)\equiv$
$\displaystyle\sum_{l=1}^{N}\sum_{i=0}^{k}g_{l}^{i}\left(\tau\right)\phi_{i}(\xi_{l}(x))[\theta(x-x_{l-1/2})-\theta(x-x_{l+1/2})].$

We assume that the kernel function is explicitly integrable and can be written
as $\mathcal{K}(u,v)=\mathcal{K}_{1}(u)\mathcal{K}_{2}(v)$, which is
effectively the case for the three simple kernels and the ballistic kernel
(see Sect. 2.3). For instance, the additive kernel writes
$\mathcal{K}_{\mathrm{kadd}}(u,v)=u+v=\mathcal{K}_{1}^{1}(u)\mathcal{K}_{2}^{1}(v)+\mathcal{K}_{1}^{2}(u)\mathcal{K}_{2}^{2}(v)$,
where $\mathcal{K}_{1}^{1}(u)=u$, $\mathcal{K}_{2}^{1}(v)=1$,
$\mathcal{K}_{1}^{2}(u)=1$ and $\mathcal{K}_{2}^{2}(v)=v$. The numerical flux
is split in two terms. The numerical flux writes

$\displaystyle
F_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g}](x,t)=\sum_{l^{\prime}=1}^{N}\sum_{i^{\prime}=0}^{k}\sum_{l=1}^{N}\sum_{i=0}^{k}g_{l^{\prime}}^{i^{\prime}}(t)g_{l}^{i}(t)$
(24) $\displaystyle\begin{aligned}
\int\limits_{x_{\mathrm{min}}}^{x}\int\limits_{x-u+x_{\mathrm{min}}}^{x_{\mathrm{max}}}&\frac{\mathcal{K}(u,v)}{v}\\\
&\phi_{i^{\prime}}(\xi_{l^{\prime}}(u))[\theta(u-x_{l^{\prime}-1/2})-\theta(u-x_{l^{\prime}+1/2})]\\\
&\phi_{i}(\xi_{l}(v))[\theta(v-x_{l-1/2})-\theta(v-x_{l+1/2})]\mathrm{d}v\mathrm{d}u.\end{aligned}$

In the DG Eq. 12, the numerical flux is evaluated on grid points $x_{j+1/2}$
and $x_{j-1/2}$ with $j\in[\\![1,N]\\!]$. $k$ is the order of the Legendre
polynomials to approximate the solution. Therefore,
$F_{\mathrm{coag}}^{\mathrm{nc}}$ depends on $j$ and $k$. The flux is sampled
over a $2$D array $(N,k+1)$ in order to use vectorial operations to reduce the
computational time. The numerical flux is

$\left\\{\begin{aligned} &F_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g}](x,t)=\\\
&\sum_{l^{\prime}=1}^{N}\sum_{i^{\prime}=0}^{k}\sum_{l=1}^{N}\sum_{i=0}^{k}g_{l^{\prime}}^{i^{\prime}}(t)g_{l}^{i}(t)T(x,x_{\mathrm{min}},x_{\mathrm{max}},i^{\prime},i,l^{\prime},l),\\\
&T(x,x_{\mathrm{min}},x_{\mathrm{max}},i^{\prime},i,l^{\prime},l)=\\\
&\begin{aligned}
&\int\limits_{x_{\mathrm{min}}}^{x}\mathcal{K}_{1}(u)\phi_{i^{\prime}}(\xi_{l^{\prime}}(u))[\theta(u-x_{l^{\prime}-1/2})-\theta(u-x_{l^{\prime}+1/2})]\\\
&\int\limits_{x-u+x_{\mathrm{min}}}^{x_{\mathrm{max}}}\frac{\mathcal{K}_{2}(v)}{v}\phi_{i}(\xi_{l}(v))\\\
&\qquad\qquad[\theta(v-x_{l-1/2})-\theta(v-x_{l+1/2})]\mathrm{d}v\mathrm{d}u.\end{aligned}\end{aligned}\right.$
(25)

A priori, the boundaries for the intervals of integration can be arbitrarily
large. We therefore rescale these intervals to avoid any numerical issues
related to large numbers when calculating the terms $T$ in the variables
$\xi_{l^{\prime}}$ and $\xi_{l}$. To avoid critical typos, the term $T$ is
derived with Mathematica by starting with the inner integral on $\xi_{l}$ and
then the integral on $\xi_{l^{\prime}}$. Further details about the derivation
of the algorithm are given in supplementary material on GitHub (see Data
Availability Sect. Data availability) for reproducibility. These integrals do
not commute. The high-order solver is written in Fortran. Reducing the number
of integrals is key to avoid numerical issues with differences of large
numbers. For this purpose, the expression of $T$ is split in several terms
provided on GitHub (see Data Availability Sect. Data availability). For
robustness, all these integrals are calculated with Mathematica. The
Mathematica function FortranForm is used to translate integral expressions to
Fortran. For large expressions, it is necessary to split them with the
function MonomialList. The scheme to evaluate
$T(x,x_{\mathrm{min}},x_{\mathrm{max}},i^{\prime},i,l^{\prime},l)$ in Fortran
is given on GitHub (see Data Availability Sect. Data availability).

A $4\mathrm{D}$ array with element
$T(x,x_{\mathrm{min}},x_{\mathrm{max}},i^{\prime},i,l^{\prime},l)$ and a
$4\mathrm{D}$ array with element $g_{l^{\prime}}^{i^{\prime}}(t)g_{l}^{i}(t)$
are computed. The element $(j,k)$ of the $2$D array corresponding to the flux
is obtained by multiplying these two $4\mathrm{D}$ arrays and summing over of
all elements. $F_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g}]$ is then evaluated
in $x_{j-1/2}$ and $x_{j+1/2}$ for all $j$.

###  3.3 Evaluation of the integral of the flux

Let denote $\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}$ the term of Eq. 12
corresponding to the integral of the numerical flux.
$\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}$ writes

$\left\\{\begin{aligned}
&\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g},j,k](t)=\\\
&\sum_{l^{\prime}=1}^{N}\sum_{i^{\prime}=0}^{k}\sum_{l=1}^{N}\sum_{i=0}^{k}\,g_{l^{\prime}}^{i^{\prime}}(t)\,g_{l}^{i}(t)\mathcal{T}\left(x_{\mathrm{min}},x_{\mathrm{max}},j,k,i^{\prime},i,l^{\prime},l\right)\\\
&\begin{aligned}
&\mathcal{T}\left(x_{\mathrm{min}},x_{\mathrm{max}},j,k,i^{\prime},i,l^{\prime},l\right)\equiv\\\
&\int\limits_{I_{j}}\int\limits_{x_{\mathrm{min}}}^{x}\int\limits_{x-u+x_{\mathrm{min}}}^{x_{\mathrm{max}}}\frac{\mathcal{K}(u,v)}{v}\,\partial_{x}\phi_{k}(\xi_{j}(x))\\\
&\qquad\phi_{i^{\prime}}(\xi_{l^{\prime}}(u))[\theta(u-x_{l^{\prime}-1/2})-\theta(u-x_{l^{\prime}+1/2})]\\\
&\qquad\phi_{i}(\xi_{l}(v))[\theta(v-x_{l-1/2})-\theta(v-x_{l+1/2})]\mathrm{d}v\,\mathrm{d}u\,\mathrm{d}x.\end{aligned}\end{aligned}\right.$
(26)

$\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g}]$ is evaluated similarly
to the flux. A triple integral is derived with Mathematica with the changes of
variables

$\xi_{l}\equiv\frac{2}{h_{l}}\left(v-x_{l}\right),\,\xi_{l^{\prime}}\equiv\frac{2}{h_{l^{\prime}}}\left(u-x_{l^{\prime}}\right),\,\xi_{j}\equiv\frac{2}{h_{j}}\left(x-x_{j}\right).$
(27)

To derive tractable equations for the integrals involving Heaviside
distributions, we start to compute integrals over the variable $\xi_{l}$, then
calculating the integral over $\xi_{l^{\prime}}$ and finally, over $x$. The
details of the calculations and the scheme in Fortran to evaluate
$\mathcal{T}\left(x_{\mathrm{min}},x_{\mathrm{max}},j,k,i^{\prime},i,l^{\prime},l\right)$
are given in supplementary material on GitHub (see Data Availability Sect.
Data availability) for completeness.
$\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}$ is computed as a product of 4D
arrays similarly to $F_{\mathrm{coag}}^{\mathrm{nc}}$. Accuracy on $T$ and
$\mathcal{T}$ depends only the quality of the polynomial approximation of $g$
by $\tilde{g}$, since the integrals corresponding to
$F_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g}]$ and
$\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}[\tilde{g}]$ are calculated
analytically.

###  3.4 Slope limiter

For most of astrophysical kernels, the solution of the Smoluchowski
coagulation equation has been mathematically shown to decays with an
exponential tail in at large masses (Schumann, 1940; Menon & Pego, 2004). This
part is challenging to approximate with polynomials, and numerical estimates
$g_{j}$ of $g$ in the bin $I_{j}$ can lead to negative values, which is not
acceptable physically.

To preserve the positivity of solution, the requirement $g_{j}(x,t)\geq 0$ for
$x\in I_{j}$ needs to be enforced. The idea is to use a scaling limiter which
controls the maximum/minimum of the reconstructed polynomials (Liu & Osher,
1996; Zhang & Shu, 2010; Liu et al., 2019). This is achieved by a
reconstruction step based on cell averaging. Let us consider the polynomials
$g_{j}(x)$ of order $k$ that approximates $g(x)$ on $I_{j}$. Let denote $m$
and $M$ two positive reals $M_{j}\equiv\underset{x\in
I_{j}}{\mathrm{max}}\,g_{j}(x)$, $m_{j}\equiv\underset{x\in
I_{j}}{\mathrm{min}}\,g_{j}(x)$ and define the scaled polynomials

$\displaystyle
p_{j}\left(x\right)\equiv\gamma_{j}\left(g_{j}(x)-\overline{g}_{j}\right)+\overline{g}_{j},$
(28)
$\displaystyle\gamma_{j}=\mathrm{min}\left\\{\left|\frac{M-\overline{g}_{j}}{M_{j}-\overline{g}_{j}}\right|,\left|\frac{m-\overline{g}_{j}}{m_{j}-\overline{g}_{j}}\right|,1\right\\}.$

where $\overline{g}_{j}$ refers to the cell average of $g$ in $I_{j}$

$\overline{g}_{j}\equiv\frac{1}{h_{j}}\int_{I_{j}}g_{j}(x,t)\mathrm{d}x.$ (29)

For all $j$, we assume $\overline{g}_{j}\in[m,M]$. $p_{j}(x)$ is a polynomial
of order $k$ such as $\overline{p}_{j}=\overline{g}_{j}$. Liu & Osher (1996)
proved that $\forall x\in I_{j},\;p_{j}(x)\in[m,M]$. This scaling limiter
allows to build a maximum-principle-satisfying DG scheme, in the sense that
the numerical solution never goes out of the range $[m,M]$. The main
difficulty is to ensure the property $\overline{g}_{j}\in[m,M]$ during the
evolution without loosing high accuracy.

In the DG scheme given by Eq. 16, polynomials $g_{j}(x)$ are replaced by the
scaled polynomials $p_{j}(x)$ such as

$\displaystyle p_{j}\left(x\right)$
$\displaystyle=\gamma_{j}\left(g_{j}(x)-\overline{g}_{j}\right)+\overline{g}_{j}$
(30)
$\displaystyle=\sum_{i=0}^{k}\gamma_{j}g_{j}^{i}(t)\phi_{1,i}(\xi_{j}(x))+\sum_{i=0}^{k}g_{j}^{i}(t)\phi_{2,i}(\xi_{j}(x))$

with

$\displaystyle\phi_{1,i}(\xi_{j}(x))\equiv\left(\phi_{i}(\xi_{j}(x))-\frac{1}{2}\int_{I_{j}}\phi_{i}(\xi_{j}(x))\mathrm{d}x\right),$
(31)
$\displaystyle\phi_{2,i}(\xi_{j}(x))\equiv\frac{1}{2}\int_{I_{j}}\phi_{i}(\xi_{j}(x))\mathrm{d}x.$

Replacing $g_{j}$ by $p_{j}$ in Eq. 25 gives four terms for the function $T$:
$T_{11}[\phi_{1,i^{\prime}}\phi_{1,i}]$,
$T_{12}[\phi_{1,i^{\prime}}\phi_{2,i}]$,
$T_{21}[\phi_{2,i^{\prime}}\phi_{1,i}]$ and
$T_{22}[\phi_{2,i^{\prime}}\phi_{2,i}]$. For each term, a corresponding
coefficient $g_{l^{\prime},i^{\prime}}(t)g_{l,i}(t)$ is associated, namely
$\gamma_{l^{\prime}}g_{l^{\prime},i^{\prime}}(t)g_{l,i}(t)$,
$\gamma_{l^{\prime}}g_{l^{\prime},i^{\prime}}(t)g_{l,i}(t)$,
$\gamma_{l}g_{l^{\prime},i^{\prime}}(t)g_{l,i}(t)$ and
$g_{l^{\prime},i^{\prime}}(t)g_{l,i}(t)$ (no $\gamma$ in the last term).
$F_{\mathrm{coag}}^{\mathrm{nc}}$ is evaluated by summing over those four
terms. The same procedure is applied for
$\mathcal{F}_{\mathrm{coag}}^{\mathrm{nc}}$. Therefore, the positivity of
$\tilde{g}$ is ensured in each cell.

###  3.5 High-order time stepping

####  3.5.1 CFL condition

Forward Euler discretisation of Eq. 12 gives

$\displaystyle\overline{g}_{j}^{n+1}=$ (32)
$\displaystyle\overline{g}_{j}^{n}-\frac{\Delta t}{\Delta
x_{j}}\left[F_{\mathrm{coag}}^{\mathrm{nc}}\left[g_{j}\right]\left(x_{j+1/2},t\right)-F_{\mathrm{coag}}^{\mathrm{nc}}\left[g_{j}\right]\left(x_{j-1/2},t\right)\right],$

for the $n$-th time step. The Courant-Friedrichs-Lewy condition (CFL) of the
scheme is chosen to guarantee the positivity of the cell average
$\overline{g}_{j}^{n+1}>0$ (Filbet & Laurencot, 2004), i.e.

$\Delta t<\frac{\Delta
x_{j}\overline{g}_{j}^{n}}{|F_{\mathrm{coag}}^{\mathrm{nc}}\left[g_{j}\right]\left(x_{j+1/2},t\right)-F_{\mathrm{coag}}^{\mathrm{nc}}\left[g_{j}\right]\left(x_{j-1/2},t\right)|}.$
(33)

This CFL condition associated with the slope limiter (see Sect. 3.4) ensures
the positivity of the global scheme. The CFL condition is initially dominated
by small grains and softens as grains grow.

####  3.5.2 Strong Stability Preserving Runge-Kutta method

In Eq. 4, the spatial derivative $\partial_{x}F_{\mathrm{coag}}[g]$ is
approximated by the nonlinearly stable operator $-L[g]$ given in Eq. 16. For
hyperbolic conservation laws, nonlinear stability is characterised by the
total variation diminishing (TVD) semi-norm

$TV\left(g\right)\equiv\sum_{j}|\overline{g}_{j+1}-\overline{g}_{j}|.$ (34)

The spatial discretisation $-L[g]$ has the property that the total variation
of the numerical solution does not increase for a forward Euler integration

$g^{n+1}=g^{n}+\Delta tL[g],\;\Delta t\leq\Delta t_{\mathrm{FE}},$ (35)

when $\Delta t_{\mathrm{FE}}$ the CFL condition determined in Eq. 33, i.e.
$TV\left(g^{n+1}\right)\leq TV\left(g^{n}\right)$. TVD property can be
generalised to high-order time discretisation with a Strong Stability
Preserving (SSP) scheme (Shu & Osher, 1988; Gottlieb et al., 2001; Zhang &
Shu, 2010; Liu et al., 2019). The method is SSP if the following condition
holds

$TV\left(g^{n+1}\right)\leq TV\left(g^{n}\right),$ (36)

and the timestep satisfies

$\Delta t_{\mathrm{SSP}}\leq c\Delta t_{\mathrm{FE}},$ (37)

where $c$ is a positive coefficient. Stability arguments are based on convex
decomposition of high-order methods in term of the first-order Euler elements.
This ensures that SSP preserves high-order accuracy in time for any convex
functional (e.g. $TV$). In practice, errors are dominated by mass
discretisation. We use a SSP Runge-Kutta (SSPRK) third-order method (Gottlieb
et al., 2009; Zhang & Shu, 2010; Liu et al., 2019) which writes, with $c=1$,

$\displaystyle\bm{g}_{j}^{(1)}$ $\displaystyle=\bm{g}_{j}^{n}+\Delta
t_{\mathrm{SSP}}\bm{L}[g_{j}^{n}],$ (38) $\displaystyle\bm{g}_{j}^{(2)}$
$\displaystyle=\frac{3}{4}\bm{g}_{j}^{n}+\frac{1}{4}\left(\bm{g}_{j}^{(1)}+\Delta
t_{\mathrm{SSP}}\bm{L}[g_{j}^{(1)}]\right),$ $\displaystyle\bm{g}_{j}^{n+1}$
$\displaystyle=\frac{1}{3}\bm{g}_{j}^{n}+\frac{2}{3}\left(\bm{g}_{j}^{(2)}+\Delta
t_{\mathrm{SSP}}\bm{L}[g_{j}^{(2)}]\right).$

This SSPRK third-order method ensures that $\overline{g}_{j}\in[m,M]$ for
$(m,M)\in\mathbb{R}_{+}$ at all times. Hence, under a suitable CFL condition,
SSP high-order time discretisation preserves the property
$\overline{g}_{j}\in[m,M]$ of the DG scheme and the linear scaling presented
in Sect. 3.4 satisfies a maximum principle.

###  3.6 Algorithm flowchart

Associating SSPRK with a DG scheme provides overall an high-order scheme that
maintains overall a uniform high-order accuracy of the solution (Zhang & Shu,
2010; Liu et al., 2019). We use the SSPRK of third order given by Eq. 38.
Splitting the algorithm into the following steps ensures positivity:

  1. 1.

Initialisation: From the initial data $g_{0}(x)$,

    1. (a)

generate $\forall j\in[\\![1,N]\\!],\;g_{j}(x,0)\in\mathcal{V}^{k}$ by
piecewise $L^{2}$ projection and get the components on Legendre basis Eq. 19,

    2. (b)

define $[m,M]$ for which $\overline{g}_{j}(x,0)\in[m,M]$,

    3. (c)

replace $g_{j}$ by $p_{j}$,

  2. 2.

Evolution: Use the scheme Eq. 38 to compute $\forall j\in[\\![1,N]\\!],\forall
i\in[\\![1,k]\\!],\;(g_{j}^{i})^{n+1}$,

  3. 3.

Reconstruction: Use Eq. 30 to reconstruct $p_{j}(x,t)$.

##  4 Numerical tests

The high-order solver presented in Sect. 3 is benchmarked against the
analytical solutions presented in Sect. 2.4, similarly to Liu et al. (2019).
Accuracy tests are performed with a small number of bins, consistently with
hydrodynamical requirements.

###  4.1 Error measurements

Numerical simulations are carried out to i) investigate the experimental order
of convergence (EOC, Kumar et al. 2014; Liu et al. 2019 ) , and ii) determine
the efficiency of the algorithm. Relative errors are measured using a
continuous norm and a discrete norm. The $L^{1}$ norm is a natural choice for
equations of conservation. The continuous $L^{1}$ norm can be approximated by
using a high order Gaussian quadrature rule

$\displaystyle\left\|f\right\|_{1}$
$\displaystyle\equiv\int_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}|f(x)|\mathrm{d}x$
(39)
$\displaystyle=\sum_{j=1}^{N}\int_{I_{j}}|f(x)|\mathrm{d}x\approx\sum_{j=1}^{N}\frac{h_{j}}{2}\sum_{\alpha=1}^{R}\omega_{\alpha}|f(x_{j}^{\alpha})|,$

where $N$ is the number of bins, $h_{j}$ is the size of bin $I_{j}$,
$\omega_{\alpha}$ are the weights and $x_{j}^{\alpha}$ are the corresponding
Gauss points in cell $I_{j}$. We use $R=16$ for sufficient accuracy. The
numerical error $e_{\mathrm{c},N}$ is measured with the continuous $L^{1}$
norm as

$e_{\mathrm{c},N}(\tau)\equiv\sum_{j=1}^{N}\frac{h_{j}}{2}\sum_{\alpha=1}^{R}\omega_{\alpha}|g_{j}(x_{j}^{\alpha},\tau)-g(x_{j}^{\alpha},\tau)|,$
(40)

where $g$ and $g_{j}$ are the analytic and the numerical solutions of the
Smoluchowski equation. Eq. 40 is computed with Mathematica using $16$ digits
for sufficient precision. The discrete $L^{1}$ norm is defined by evaluating
$g_{j}$ and $g$ at the geometric mean
$\hat{x}_{j}\equiv\sqrt{x_{j-1/2}x_{j+1/2}}$ of the bin $I_{j}$. The numerical
error measured with this discrete $L^{1}$ norm is

$e_{\mathrm{d},N}(\tau)\equiv\sum_{j=1}^{N}h_{j}|g_{j}(\hat{x}_{j},\tau)-g(\hat{x}_{j},\tau)|.$
(41)

We follow Liu et al. (2019) to calculate the experimental order of convergence
(EOC)

$\mathrm{EOC}\equiv\frac{\ln\left(\frac{e_{N}(\tau)}{e_{2N}(\tau)}\right)}{\ln(2)},$
(42)

where $e_{N}$ is the error evaluated for $N$ cells and $e_{2N}$ for $2N$
cells. For the calculation of the EOC, the numerical errors are calculated at
time $\tau=0.01$ for the order of convergence of the DG scheme not to be
altered by time stepping errors. The moments of the numerical solution are
defined according to

$\displaystyle M_{p,N}\left(\tau\right)$
$\displaystyle=\int\limits_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}x^{p-1}\tilde{g}(x,\tau)\mathrm{d}x$
(43) $\displaystyle=\sum_{j=1}^{N}\int_{I_{j}}x^{p-1}g_{j}(x,\tau)\mathrm{d}x$
$\displaystyle=\sum_{j=1}^{N}\sum_{i=0}^{k}g_{j}^{i}(\tau)\int_{I_{j}}x^{p-1}\phi_{i}\left(\xi_{j}(x)\right)\mathrm{d}x.$

The total mass of the system writes

$M_{1,N}(\tau)=\sum_{j=1}^{N}\sum_{i=0}^{k}g_{j}^{i}(\tau)\frac{h_{j}}{2}\underbrace{\int\limits_{-1}^{1}\phi_{i}\left(\xi_{j}\right)\mathrm{d}\xi_{j}}_{=\delta_{00}=2}=\sum_{j=1}^{N}h_{j}g_{j}^{0}(\tau).$
(44)

Absolute errors on moments are given by

$e_{M_{p,N}}(\tau)\equiv\frac{|M_{p,N}(\tau)-M_{p}(\tau)|}{M_{p}(\tau)},$ (45)

where $M_{p}(\tau)$ is the moment of order $p$ at time $\tau$ for the exact
solution. In usual convergence tests, errors are normalised with respect to
the number of degrees of freedom of the algorithm. This is not the case here,
since we compare absolute gains for the purpose interfacing it with an
hydrodynamical solver.

###  4.2 Practical implementation of the tests

Numerical tests are performed by comparing numerical solutions the constant,
additive and multiplicative kernels to the solutions given in Eqs. 6, 7 and 8.
Solutions are integrated over the intervals $x\in[10^{-3},10^{6}]$ for the
constant and the additive kernels, and $x\in[10^{-3},10^{3}]$ for the
multiplicative kernel. Tests are performed with Fortran, errors are calculated
with Mathematica at machine precision. Quadruple precision is required for the
additive kernel with $k=2$, and for all kernels with $k=3$. The results are
shown for Legendre polynomials of order $k=0,1,2,3$. Above order $3$,
numerical errors due to arithmetics of large numbers are not negligible
anymore. A safety coefficient of $1/2$ is applied on the CFL condition, i.e.
the coagulation time-step used in practice is
$\mathrm{d}\tau_{\mathrm{coag}}=1/2\,\mathrm{d}\tau_{\mathrm{CFL}}$. Initial
conditions are set to satisfy the analytic solution at initial time $\tau=0$.
The analytical and numerical solutions are compared when particles of large
masses are formed at final times $\tau$ that depend on the kernels.
Simulations are performed by dividing $\tau$ into constant dumps of value
$\mathrm{d}\tau$ (300 for the constant and the additive kernels, 10000 for the
multiplicative kernel). Each dump is subdivided in several coagulation steps
satisfying the CFL condition. The analytical derivation of the coagulation
flux allows the algorithm to be efficient, i.e. to reach desired accuracy with
a low computational time. To quantify efficiency, the computational time is
compared to the one obtained with the scheme of Liu et al. (2019) with a
number of Gauss points $Q=k+1$ on a simulation in double precision with $N=20$
bins, $k=1$ for the additive kernel and $k=2$ for the constant and
multiplicative kernels. The Liu scheme is implemented by following the
description of Liu et al. (2019) step-by-step, without additional
optimisations. Simulations are performed in sequential on an Intel Core i$7$
$2.8\mathrm{GHz}$. We use the gfortran v9.2.0 compiler. Such a comparison is
delicate to perform and interpret, since it is implementation-dependant.
Should the number of Gauss points in the Liu algorithm be increased to better
approximate the integral terms calculated here analytically, this may result
in an increase of computational time by several orders of magnitudes, giving
the false impression that the Liu algorithm is not performant. Hence the
choice $Q=k+1$. Qualitatively, our scheme is more effective by a factor of
several unities for same precision and without requiring sub-binning, except
for the additive kernel for which the Liu scheme exhibits serendipitous super-
convergence (Liu et al., 2019).

###  4.3 Constant kernel

####  4.3.1 Positivity and mass conservation

Fig. 6 shows the numerical solutions obtained for $N=20$ bins, varying the
order of the polynomials $k$. The analytical and numerical solutions are
compared at time $\tau=30000$. As expected, the solution remains positive, as
a result from combining the slope limiter (see Sect. 3.4) and the SSP Runge-
Kutta time stepping (see Sect. 3.5.2). The piecewise linear solution ($k=1$)
appears curved due to the logarithmic scale of the $x$-axis. Fig. 5 shows the
numerical absolute error $e_{M_{1},N}$ on the moment $M_{1,N}$ for $N=20$ bins
from $\tau=0$ to $\tau=30000$. The total mass remains conserved to machine
precision until $\tau=10^{4}$.

Figure 5: Test case, constant kernel: evolution of the numerical absolute
error $e_{M_{1},N}$ on the moment $M_{1,N}$ for $N=20$ bins. The divergence at
long times is explained by accumulation of errors due to numerical diffusion
for even orders $k=0$ and $k=2$. Total mass is conserved at machine precision
until $\tau=10^{4}$.

####  4.3.2 Accuracy of the numerical solution

As expected, the accuracy of the numerical solution improves with the order of
the scheme. Fig. 7 shows the numerical solution obtained at $\tau=30000$ (note
the 16 orders of magnitude in mass on the $y$ axis in log). The major part of
the total mass of the system is located around the maximum of the curve. Fig.
7 shows that around this maximum, schemes of order $k=1,2,3$ provide errors of
order $\sim 0.1-1\%$ when $k=0$ generates errors of order $\sim 30\%$. Fig. 7
also shows that numerical diffusion is drastically reduced in the exponential
tail as the order of the scheme increases, since a gain of a factor $\sim 100$
is obtained with order $3$ compared to order $0$.

Figure 6: Test case, constant kernel: the numerical solution $g_{j}(x,\tau)$
is plotted for $N=20$ bins and $k=0,1,2,3$ from $\tau=0$ to $\tau=30000$, and
compared to the analytic solution $g(x,\tau)$. Vertical grey lines delimit the
bins. The accuracy improves for larger values of $k$. Order $3$ approximates
the bump where the major part of the mass is concentrated with accuracy of
order $\sim 0.1\%$. Figure 7: Test case, constant kernel: numerical solution
$g_{j}(x,\tau$) evaluated with the geometric mean $\hat{x}_{j}$ over each bin
$I_{j}$. At the location of the maximum, orders $k=1,2,3$ achieve an absolute
error of $\sim 0.1-1\%$, to be compared with $30\%$ obtained with $k=0$.
Accuracy in the exponential tail is improved by a factor $100$ with $k=3$
compared to $k=0$.

####  4.3.3 Convergence analysis

Numerical errors introduced in Sect. 4.1 are shown on Fig. 8 at $\tau=0.01$.
$e_{\mathrm{c},N}$ and $e_{\mathrm{d},N}$ are plotted as a functions of the
number of bins per decade $N_{\mathrm{bin}/\mathrm{dec}}$, to infer the EOC
independently from the global mass interval. With the continuous $L^{1}$ norm,
the EOC is of order $k+1$ on a geometric grid, similarly to Liu et al. (2019).
With the discrete $L^{1}$ norm, the EOC is of order $k+2$ for odd polynomials,
and $k+1$ for even polynomials. We recover second order of convergence
(EOC=$2$) for the finite volume scheme with $k=0$ that was predicted by Filbet
& Laurencot (2004). Fig. 8 shows that the expected accuracy of order $\sim
0.1\%$ on $e_{d,N}$ is achieved with more than $10$ bins/decade for orders $0$
and $1$, with $\sim 9$ bins/decade for order $2$ and with $\sim 5$ bins/decade
for order $3$. Accuracy of order $\sim 1\%$ is achieved with $\sim 9$
bins/decade for orders $0$ and $1$, with $\sim 5$ bins/decade for order $2$,
and with $\sim 2$ bins/decade for order $3$.

Figure 8: Test case, constant kernel: the continuous $L^{1}$ error $e_{c,N}$
and the discrete $L^{1}$ error $e_{d,N}$ are plotted as functions of the
number of bins per decade. With $e_{c,N}$, the experimental order of
convergence is EOC = $k+1$. With $e_{d,N}$, EOC = $k+1$ for polynomials of odd
orders and EOC = $k+2$ for polynomials of even orders. The DG scheme achieves
on $e_{d,N}$ an accuracy of $0.1\%$ with more than $10$ bins/decade for
$k=0,1$, with $\sim 9$ bins/decade for $k=2$ and with $\sim 5$ bins/decade for
$k=3$. An accuracy of $1\%$ is achieved with $\sim 9$ bins/decade for $k=0,1$,
with $\sim 5$ bins/decade for $k=2$ and $\sim 2$ bins/decade for $k=3$.

####  4.3.4 Stability in time

Time evolution of the numerical errors $e_{c,N}$ and $e_{d,N}$ are shown in
Fig. 9. The results are shown for $N=20$ bins for $k=0,1,2,3$ at time
$\tau=30000$, when particles of large masses have formed. We verify that
$e_{c,N}$ and $e_{d,N}$ remain bounded.

Figure 9: Test case, constant kernel: numerical errors $e_{c,N}$ with the
$L^{1}$ continuous norm, $e_{d,N}$ with the discrete $L^{1}$ norm. All these
errors are calculated for $N=20$. Errors remain bounded at large times.

####  4.3.5 Computational efficiency

Fig. 10 shows that similar accuracies are obtained with this scheme and the
scheme described in Liu et al. (2019). Computational time is compared on a
simulation with $N=20$ bins, $k=2$ and a final time $\tau=30000$ after $\sim
10^{3}$ timesteps. The computational time for the Liu et al. (2019) scheme is
around 16 seconds (real time). The computational time for this scheme is
around 4 seconds (real time). An improvement of factor 4 is therefore achieved
for the computational time by estimating integrals analytically.

Figure 10: Test case, constant kernel: comparison with the scheme of Liu et
al. (2019). Similar accuracies are reached, but being $\sim 4\times$ more
effective due to numerical integration.

###  4.4 Additive kernel

####  4.4.1 Positivity and mass conservation

Fig. 11 shows numerical solutions obtained for $N=20$ bins and $k=0,1,2,3$ at
time $\tau=3$. The numerical solutions remains positive as grains grow. Fig.
12 shows the evolution of the numerical absolute error $e_{M_{1},N}$ on the
first moment $M_{1,N}$. The total mass remains conserved to machine precision
until $\tau=1$.

Figure 11: Test case, additive kernel: the numerical solution $g_{j}(x,\tau)$
is plotted for $N=20$ bins and $k=0,1,2,3$ from $\tau=0$ to $\tau=3$, and
compared to the analytic solution $g(x,\tau)$. Vertical grey lines delimit the
bins. The accuracy improves for larger values of $k$. Order $3$ approximates
the bump where the major part of the mass is concentrated with accuracy of
order $\sim 0.1\%$.  Figure 12: Test case, additive kernel: evolution of the
numerical absolute error $e_{M_{1},N}$ on the moment $M_{1,N}$ for $N=20$
bins. The divergence at long times is explained by accumulation of errors due
to numerical diffusion for orders $k=0$, $k=2$ and $k=3$. Total mass is
conserved at machine precision until $\tau=1$.

####  4.4.2 Accuracy of the numerical solution

Fig. 13 shows numerical solutions obtained at $\tau=3$ on a logarithmic scale.
Fig. 13 reveals a strong numerical diffusion for order $0$. Numerical errors
are indeed integrated and diffused extremely efficiently towards large masses
by the additive kernel. In this case, the mass density for large-masses
particles is over-estimated by several orders of magnitude. High-order schemes
reduce this numerical diffusion as expected. Fig. 13 shows that around the
maximum, schemes of order $k=1,2,3$ provide errors of order $\sim 0.1-1\%$
when $k=0$ generates errors of order $\sim 10\%$. Numerical diffusion is
reduced in the exponential tail as the order of the scheme increases, up to
reaching a gain of a factor $\sim 10000$ with order $3$ compared to order $0$.

Figure 13:  Test case, additive kernel: At the location of the maximum, orders
$k=1,2,3$ achieve an absolute error of $\sim 0.1-1\%$, to be compared with
$10\%$ obtained with $k=0$. Accuracy in the exponential tail is improved by a
factor $10000$ by $k=3$ compared to $k=0$.

####  4.4.3 Convergence analysis

Numerical errors are shown on Fig. 14 at $\tau=0.01$. Accuracy of order $\sim
0.1\%$ on $e_{d,N}$ errors are achieved with more than $10$ bins/decade for
order $0$ and $1$, with $\sim 9$ bins/decade for orders $2$ and $3$. Accuracy
of order $\sim 1\%$ is achieved with $\sim 9$ bins/decade for orders $0$ and
$1$, with $\sim 5$ bins/decade for order $2$ and $\sim 2$ bins/decade for
order $3$.

Figure 14: Test case, additive kernel: similar to Fig. 8. The DG scheme
achieves on $e_{d,N}$ an accuracy of order $0.1\%$ with more than $10$
bins/decade for $k=0,1$, with $\sim 5$ bins/decade for $k=2,3$. An accuracy of
order $1\%$ is achieved with $\sim 9$ bins/decade for $k=0,1$, with $\sim 5$
bins/decade for $k=2$ and with $\sim 2$ bins/decade for $k=3$.

####  4.4.4 Stability in time

Evolution of the numerical errors $e_{c,N}$ and $e_{d,N}$ are shown in Fig.
15. The results are shown for $N=20$ bins for $k=0,1,2,3$ at $\tau=3$, when
particles with large masses have formed. At order $0$, $e_{c,N}$ (resp.
$e_{d,N}$) increases significantly after $\tau\approx 5\cdot 10^{-1}$ (resp.
$\tau\approx 10^{-1}$). On the contrary, $e_{c,N}$ and $e_{d,N}$ remain
bounded for longer times at orders $1$, $2$ and $3$.

Figure 15: Test case, additive kernel: numerical errors $e_{c,N}$ with the
$L^{1}$ continuous norm, $e_{d,N}$ with the discrete $L^{1}$ norm. All these
errors are calculated for $N=20$. Errors remain bounded at large times for
orders $k=1,2,3$.

####  4.4.5 Computational efficiency

Computational time is compared to Liu et al. (2019) on a simulation with
$N=20$ bins, $k=1$ and a final time $\tau=3$. Fig. 10 shows similar accuracy
for both schemes. The computational time for the Liu et al. (2019) scheme is
around 3 seconds (real time) for a number of Gauss quadrature points $Q=2$.
The computational time for this scheme is 1 second, providing an improvement
by a factor 3. Fig. 16 also shows that for the additive kernel, the Liu scheme
with $Q=2$ is counter-intuitively more accurate than for $Q=16$ and the DG
scheme. This result can be explained by a serendipitous compensation of errors
when approximating the integrals with a Gauss quadrature of low order.

Figure 16: Test case, additive kernel: comparison with the scheme of Liu et
al. (2019). Unexpected accuracy occurs for integral estimates with $Q=2$ Gauss
points due to serendipitous error compensations. Our algorithm is $\sim
3\times$ more effective due to analytical integration compared to the Lui
scheme with $Q=2$.

###  4.5 Multiplicative kernel

####  4.5.1 Positivity and mass conservation

Fig. 17 shows the numerical solutions obtained for $N=20$ bins and $k=0,1,2,3$
after $\tau=100$. The numerical solutions remain positive as grain grow. Fig.
18 shows the evolution of $e_{M_{1},N}$. Total mass remains conserved to
machine precision until $\tau<1$. At $\tau=1$, gelation occurs, particles with
infinite mass are formed (McLeod, 1962b; Ernst et al., 1984; Filbet &
Laurencot, 2004) and total mass is no longer conserved anymore.

Figure 17: Test case, multiplicative kernel: numerical solution
$g_{j}(x,\tau)$ is plotted for $N=20$ bins for $k=0,1,2,3$ from $\tau=0$ to
$\tau=100$ and compared to the analytic solution $g(x,\tau)$. Vertical grey
lines delimit the bins. Accuracy of order $\sim 0.1\%$ is achieved at all
orders.

Figure 18:  Test case, multiplicative kernel: evolution of the numerical
absolute error $e_{M_{1},N}$ on the moment $M_{1,N}$ for $N=20$ bins. Mass is
conserved anymore when gelation occurs at $\tau=1$.

####  4.5.2 Accuracy of the numerical solution

Fig. 19 shows the numerical solution for the multiplicative kernel at
$\tau=100$. Accuracy of order $\sim 0.1\%$ is obtained at all orders, even
$k=0$. Physically, growth is effective enough for advection in the mass space
to be more efficient than numerical diffusion.

Figure 19: Test case, multiplicative kernel: Accuracy of order $\sim 0.1\%$ is
achieved at any order. Growth is so efficient than is overtakes numerical
diffusion.

####  4.5.3 Convergence analysis

Numerical errors are shown on Fig. 20 at $\tau=0.01$. Accuracy of order $\sim
0.1\%$ on $e_{d,N}$ errors are achieved with $\sim 15$ bins/decade for orders
$0$ and $1$, with $\sim 7$ bins/decade for order $2$, and with $\sim 4$
bins/decade for order $3$. Accuracy of order $\sim 1\%$ is achieved with $\sim
7$ bins/decade for orders $0$ and $1$, with $\sim 2$ bins/decade for order
$2$, and with $\sim 1$ bins/decade for order $3$.

Figure 20: Test case, multiplicative kernel: the continuous $L^{1}$ error
$e_{c,N}$ and the discrete $L^{1}$ error $e_{d,N}$ are plotted as functions of
the number of bins per decade. With $e_{c,N}$, the experimental order of
convergence is EOC = $k+1$. With $e_{d,N}$, EOC = $k+1$ for polynomials of odd
orders and EOC = $k+2$ for polynomials of even orders. The DG scheme achieves
on $e_{d,N}$ an accuracy of $0.1\%$ with $\sim 15$ bins/decade for $k=0,1$,
with $\sim 7$ bins/decade for $k=2$ and with $\sim 4$ bins/decade for $k=3$ .
An accuracy of $1\%$ is achieved with $\sim 7$ bins/decade for $k=0,1$, with
$\sim 2$ bins/decade for $k=2$ and with $\sim 1$ bins/decade for $k=3$.

####  4.5.4 Stability in time

The evolution of the numerical errors $e_{c,N}$ and $e_{d,N}$ are shown in
Fig. 21. The results are shown for $N=20$ bins fo $k=0,1,2,3$ at time
$\tau=100$, when particles with large masses have formed. We observe that
$e_{c,N}$ and $e_{d,N}$ remain bounded, even after the occurence of gelation
at $\tau=1$.

Figure 21: Test case, multiplicative kernel: numerical errors $e_{c,N}$ with
the $L^{1}$ continuous norm, $e_{d,N}$ with the discrete $L^{1}$ norm. All
these errors are calculated for $N=20$. Errors remain bounded at large times.

####  4.5.5 Computational efficiency

Fig. 22 shows similar accuracies for the Liu et al. (2019) scheme and our
implementation. With $k=2$, the computational time for the Liu et al. (2019)
scheme is around 8 minutes for a number of Gauss quadrature points $Q=3$. The
computational time is for this scheme 1 minute and 40 seconds, providing an
improvement by a factor 5.

Figure 22: Test case, multiplicative kernel: comparison between the numerical
solutions provided by this scheme and the scheme of Liu et al. (2019). Similar
accuracies are reached, but being $\sim 5\times$ more effective due to
analytical integration.

##  5 Discussion

The Discontinuous Galerkin scheme presented in Sect. 3 involves polynomials of
high-order, implying issues with differences of large real numbers. Order
$k=3$ appears as a maximum limit for the order of the scheme in its current
form in practice. So far, the ratio betwen $10^{6}$ coagulation time-steps and
one hydrodynamical time-step with PHANTOM using $10^{6}$ SPH particles is of
order $\sim 10-100$. We are confident the we can reach a one-to-one ratio by
i) taking advantage of more ingenious time-stepping (e.g. Carrillo & Goudon
2004; Goudon et al. 2013, ii) adopt a more relevant choice for the basis (e.g.
Soong 1974) and iii) use GPU parallelisation, since calls to the coagulation
solver by the hydrodynamical code are independent. These strategies to further
gain accuracy and computational efficiency will be tested in a next future.

The most relevant kernel for astrophysics is the ballistic kernel (Sect. 2.3).
Large-scale values of $\Delta v$ are provided by 2D piecewise constant
functions from hydrodynamic codes. In discs, the $\Delta v$ function
encompasses radial drift, vertical settling and turbulence at large scales.
The ballistic kernel splits in three terms

$\displaystyle\mathcal{K}_{b}(u,v)$
$\displaystyle=\pi(u^{2/3}+2u^{1/3}v^{1/3}+v^{2/3})\Delta v(u,v)$ (46)
$\displaystyle=\mathcal{K}_{b,1}(u,v)+\mathcal{K}_{b,2}(u,v)+\mathcal{K}_{b,3}(u,v),$

$\mathcal{K}_{b,1}(u,v)\equiv\pi u^{2/3}\Delta v(u,v)$,
$\mathcal{K}_{b,2}(u,v)\equiv\pi 2u^{1/3}v^{1/3}\Delta v(u,v)$ and
$\mathcal{K}_{b,3}(u,v)\equiv\pi v^{2/3}\Delta v(u,v)$. The numerical flux is
also split in three terms that are evaluated analytically. Models of
differential velocities are also used to model sub-grid small-scale values of
$\Delta v$ (Brownian motion, dusty turbulence at small scales). Shall these
kernels not be integrable, we will estimate them with an appropriate
interpolation.

Moreover, the algorithm presented above solves for the Smoluchowski equation
with pure growth. Although fragmentation plays a key role in regulating the
number of small grains and preventing the formation of large bodies, it has
not being included in the solver yet. The algorithm presented in Sect. 3 has
been designed to incorporate fragmentation genuinely by adding the extra
fragmentation flux (Paul & Kumar, 2018)

$\displaystyle F_{\mathrm{frag}}\left[g\right]\left(x,\tau\right)\equiv$ (47)
$\displaystyle\int\limits_{0}^{\infty}\int\limits_{x}^{\infty}\int\limits_{0}^{x}\frac{w}{yz}b(w,y,z)\mathcal{K}(y,z)g(y,\tau)g(z,\tau)\mathrm{d}w\mathrm{d}y\mathrm{d}z,$

similarly e.g. to Birnstiel et al. (2010). The kernel $\mathcal{K}$ provides
the fragmentation rate between two particles of masses $x$ and $y$. The
function $b$ is the breakage rate related to the formation of a particle of
mass $x$ from particles of mass $y$ and $w$. Known functional forms of the
fragmentation kernel should authorise direct analytic integrations, similarly
to the derivations performed in Sect. 3.2. For peculiar regimes, fragmentation
kernels can alternatively be interpolated. Astrophysical mass distributions
are expected to be dominated by large grains. Hence, the CFL condition for
fragmentation should be similar to the one for growth (Vericel & Gonzalez,
2020). If so, numerical integration will be performed explicitly. If not,
implicit time-stepping can be implemented in a manageable way since the number
of dust bins has been kept minimal with analytic integrations (i.e. linear
algebra with $\sim 15\times 15$ matrices).

Eq. 1 restrains dust interactions to binary collisions between aggregates of
spherical shapes. Multiple collisions are not expected to play a critical role
in astrophysics, since dust volume densities are extremely low. On the other
hand, dust aggregates are expected to be porous or have fractal structures. In
particular, small bodies that have not been recompacted by collisions are
expected to be fluffy. Eq. 1 also reduces probability distributions of
velocities to their mean values. This approximation may quench grain growth
occurring through rare collisional events, e.g. between large bodies having
low relative velocities (Windmark et al., 2012; Garaud et al., 2013). Finally,
growth is in essence stochastic, but fluctuations of the solution can not be
computed with Eq. 1. This is not critical, those being hardly constrained by
observations. Although the solver presented in Sect. 3 can not be used
directly to treat the additional physical processes described above, the
method could be adapted to do so. Lastly, extending Eq. 1 to multiple
compositions, without or with change of states has been done in other
communities. This comes to the cost of multiplying the number of variables by
the number of materials considered. The algorithm presented in Sect. 3 is a
first step towards reducing the number of dust bins to allow for solving for
multiple compositions in 3D. This would have strong implications for planet
formation, e.g. by handling snow lines consistently and providing constrains
for meteoritic data.

##  6 Conclusion

We have presented an high-order algorithm that solves accurately the
coagulation equation with a limited number of dust bins ($\sim 15$).
Specifically:

  1. 1.

Mass is conserved to machine precision for astrophysical kernels,

  2. 2.

Positivity is guaranteed by combining an appropriate slope-limiter to a Total
Variation Diminishing time-stepping,

  3. 3.

Creating aggregates of masses larger that the mass of the reservoir is
mathematically excluded by a control of the growth flux,

  4. 4.

Errors of order $0.1-1\%$ are achieved by high-order discretisation in time
and space that can be modulated for convergence purpose. They shall not
dominate the error budget over hydrodynamics,

  5. 5.

Combining a low number of bins and analytic integrations allows manageable
costs in memory and time,

  6. 6.

Additional physics should be implementable in a versatile way.

The next step consists of performing 3D hydrodynamical simulations of star and
planet formation with accurate dust growth. The design of the algorithm allows
to implement additional processes such as fragmentation in a genuine way. This
solver encourages the reduction of CO2 emissions related to computational
astrophysics.

## Acknowledgements

GL acknowledges funding from the ERC CoG project PODCAST No 864965. This
project has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
823823. This project was partly supported by the IDEXLyon project (contract
nANR-16-IDEX-0005) under the auspices University of Lyon. We acknowledge
financial support from the national programs (PNP, PNPS, PCMI) of CNRS/INSU,
CEA, and CNES, France. We used Mathematica (Wolfram Research, 14). We thank L.
Tine, E. Déléage, D. Price, D. Mentiplay, T. Guillet, S. Charnoz, R. Teyssier
and the anonymous referee for useful comments and discussions.

## Data availability

The data and supplementary material underlying this article are available in
the repository ”growth” on GitHub at https://github.com/mlombart/growth.git.
Figures can be reproduced following the file README.md. The repository
contains data and Python scripts used to generate figures.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership et al., 2015, The Astrophysical Journal, 808, L3
  * Aldous (1999) Aldous D. J., 1999, Bernoulli, 5, 3
  * Avenhaus et al. (2018) Avenhaus H., et al., 2018, The Astrophysical Journal, 863, 44
  * Berry (1967) Berry E. X., 1967, Journal of Atmospheric Sciences, 24, 688
  * Berry & Reinhardt (1974) Berry E. X., Reinhardt R. L., 1974, Journal of Atmospheric Sciences, 31, 1814
  * Birnstiel et al. (2010) Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79
  * Bleck (1970) Bleck R., 1970, J. Geophys. Res., 75, 5165
  * Bott (1998) Bott A., 1998, Journal of Atmospheric Sciences, 55, 2284
  * Bourgade & Filbet (2008) Bourgade J.-P., Filbet F., 2008, Mathematics of Computation, 77, 851
  * Brauer et al. (2008) Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859
  * Carrillo & Goudon (2004) Carrillo J. A., Goudon T., 2004, J. Sci. Comput., 20, 69–113
  * Chandrasekhar (1943) Chandrasekhar S., 1943, Reviews of Modern Physics, 15, 1
  * Cockburn & Shu (1989) Cockburn B., Shu C.-W., 1989, Mathematics of Computation, 52, 411
  * Debry et al. (2003) Debry E., Sportisse B., Jourdain B., 2003, Journal of Computational Physics, 184, 649
  * Dominik et al. (2007) Dominik C., Blum J., Cuzzi J. N., Wurm G., 2007, Protostars and Planets V, pp 783–800
  * Drake (1972) Drake R. L., 1972, Journal of the Atmospheric Sciences, 29, 537
  * Dullemond & Dominik (2005) Dullemond C. P., Dominik C., 2005, A&A, 434, 971
  * Enukashvily (1980) Enukashvily I. M., 1980, Journal of Atmospheric Sciences, 37, 2521
  * Erasmus et al. (1994) Erasmus L., Eyre D., Everson R., 1994, Computers & Chemical Engineering, 18, 775
  * Ernst et al. (1984) Ernst M. H., Ziff R. M., Hendriks E. M., 1984, Journal of Colloid and Interface Science, 97, 266
  * Estrada & Cuzzi (2008) Estrada P. R., Cuzzi J. N., 2008, ApJ, 682, 515
  * Eyre et al. (1988) Eyre D., Wright C. J., Reutert G., 1988, Journal of Computational Physics, 78, 288
  * Filbet (2008) Filbet F., 2008, AIP Conference Proceedings, 1048, 921
  * Filbet & Laurencot (2004) Filbet F., Laurencot P., 2004, SIAM Journal on Scientific Computing, 25, 2004
  * Finlayson (1972) Finlayson B., 1972, The Method of Weighted Residuals and Variational Principles, Mathematics in Science and Engineering
  * Forestier-Coste & Mancini (2012) Forestier-Coste L., Mancini S., 2012, SIAM Journal on Scientific Computing, 34, B840
  * Friedlander & Wang (1966) Friedlander S. K., Wang C. S., 1966, Journal of Colloid and Interface Science, 22, 126
  * Friedlander et al. (2000) Friedlander S. K., et al., 2000, Smoke, dust, and haze. Vol. 198, Oxford university press New York
  * Gabriel & Tine (2010) Gabriel P., Tine L. M., 2010, ESAIM: Proc., 30, 53
  * Garaud et al. (2013) Garaud P., Meru F., Galvagni M., Olczak C., 2013, ApJ, 764, 146
  * Gelbard & Seinfeld (1978) Gelbard F., Seinfeld J. H., 1978, Journal of Computational Physics, 28, 357
  * Gelbard et al. (1980) Gelbard F., Tambour Y., Seinfeld J. H., 1980, Journal of Colloid and Interface Science, 76, 541
  * Gillespie (1975a) Gillespie D. T., 1975a, Journal of Atmospheric Sciences, 32, 600
  * Gillespie (1975b) Gillespie D. T., 1975b, Journal of Atmospheric Sciences, 32, 1977
  * Golovin (1963) Golovin A., 1963, Izv. Geophys. Ser, 5, 482
  * Gonzalez et al. (2017) Gonzalez J. F., Laibe G., Maddison S. T., 2017, MNRAS, 467, 1984
  * Gottlieb et al. (2001) Gottlieb S., Shu C.-W., Tadmor E., 2001, SIAM Review, 43, 89
  * Gottlieb et al. (2009) Gottlieb S., Ketcheson D. I., Shu C.-W., 2009, J. Sci. Comput., 38, 251–289
  * Goudon et al. (2013) Goudon T., Lagoutière F., Tine L. M., 2013, Mathematical Models and Methods in Applied Sciences, 23, 1177
  * Hayashi & Nakagawa (1975) Hayashi C., Nakagawa Y., 1975, Progress of Theoretical Physics, 54, 93
  * Hidy & Brock (1972) Hidy G., Brock J., 1972, Topics in Current Aerosol Research. Part 2. Pergamon Press
  * Hounslow et al. (1988) Hounslow M. J., Ryall R. L., Marshall V. R., 1988, AIChE Journal, 34, 1821
  * Hulburt & Katz (1964) Hulburt H., Katz S., 1964, Chemical Engineering Science, 19, 555
  * Hutchison et al. (2018) Hutchison M., Price D. J., Laibe G., 2018, MNRAS, 476, 2186
  * Jacobson (2005) Jacobson M. Z., 2005, Fundamentals of Atmospheric Modeling
  * Jiang & Peng (2000) Jiang G.-S., Peng D., 2000, SIAM Journal on Scientific Computing, 21, 2126
  * Khain & Pinsky (2018) Khain A. P., Pinsky M., 2018, Physical processes in clouds and cloud modeling. Cambridge University Press
  * Kostoglou & Karabelas (1994) Kostoglou M., Karabelas A. J., 1994, Journal of Colloid and Interface Science, 163, 420
  * Kovetz & Olund (1969) Kovetz A., Olund B., 1969, Journal of the Atmospheric Sciences, 26, 1060
  * Kumar & Ramkrishna (1996) Kumar S., Ramkrishna D., 1996, Chemical Engineering Science, 51, 1311
  * Kumar et al. (2014) Kumar R., Kumar J., Warnecke G., 2014, Kinetic & Related Models, 7, 713
  * Laibe & Price (2014) Laibe G., Price D. J., 2014, MNRAS, 440, 2136
  * Landgrebe & Pratsinis (1990) Landgrebe J. D., Pratsinis S. E., 1990, Journal of Colloid and Interface Science, 139, 63
  * Laurençot (2018) Laurençot P., 2018, Journal of Statistical Physics,
  * Lebreuilly et al. (2019) Lebreuilly U., Commerçon B., Laibe G., 2019, A&A, 626, A96
  * Lebreuilly et al. (2020) Lebreuilly U., Commerçon B., Laibe G., 2020, A&A, 641, A112
  * Lee & Matsoukas (2000) Lee K., Matsoukas T., 2000, Powder Technology, 110, 82
  * Leyvraz & Tschudi (1981) Leyvraz F., Tschudi H. R., 1981, Journal of Physics A Mathematical General, 14, 3389
  * Liffman (1992) Liffman K., 1992, Journal of Computational Physics, 100, 116
  * Lister et al. (1995) Lister J. D., Smit D. J., Hounslow M. J., 1995, AIChE Journal, 41, 591
  * Liu & Osher (1996) Liu X.-D., Osher S., 1996, SIAM Journal on Numerical Analysis, 33, 760
  * Liu et al. (2019) Liu H., Gröpler R., Warnecke G., 2019, SIAM Journal on Scientific Computing, 41, B448
  * Mahoney & Ramkrishna (2002) Mahoney A. W., Ramkrishna D., 2002, Chemical Engineering Science, 57, 1107
  * Marchisio et al. (2003) Marchisio D. L., Vigil R. D., Fox R. O., 2003, Journal of Colloid and Interface Science, 258, 322
  * McLeod (1962a) McLeod J. B., 1962a, Proc. of the London Mathematical Society, s3-14, 445
  * McLeod (1962b) McLeod J. B., 1962b, Proc. of the London Mathematical Society, s3-14, 445
  * Melzak (1953) Melzak Z. A., 1953, Quarterly of Applied Mathematics, 11, 231
  * Melzak (1957) Melzak Z., 1957, Transactions of the American Mathematical Society, 85, 547
  * Menon & Pego (2004) Menon G., Pego R. L., 2004, Communications on Pure and Applied Mathematics, 57, 1197
  * Middleton & Brock (1976) Middleton P., Brock J., 1976, Journal of Colloid and Interface Science, 54, 249
  * Müller (1928) Müller H., 1928, Fortschrittsberichte über Kolloide und Polymere, 27, 223
  * Niethammer et al. (2016) Niethammer B., Throm S., Velázquez J. J. L., 2016, Journal of Statistical Physics, 164, 399
  * Ohtsuki et al. (1990) Ohtsuki K., Nakagawa Y., Nakazawa K., 1990, Icarus, 83, 205
  * Ormel et al. (2007) Ormel C. W., Spaans M., Tielens A. G. G. M., 2007, A&A, 461, 215
  * Paul & Kumar (2018) Paul J., Kumar J., 2018, Mathematical Methods in the Applied Sciences, 41, 2715
  * Pilinis (1990) Pilinis C., 1990, Atmospheric Environment, 24, 1923
  * Pinte et al. (2020) Pinte C., et al., 2020, ApJ, 890, L9
  * Piskunov & Golubev (2002) Piskunov V. N., Golubev A. I., 2002, Journal of Aerosol Science, 33, 51
  * Price et al. (2018) Price D. J., et al., 2018, PASA, 35, e031
  * Pruppacher & Klett (1980) Pruppacher H. R., Klett J. D., 1980, Nature, 284, 88
  * Pruppacher & Klett (2010) Pruppacher H. R., Klett J. D., 2010, in , Microphysics of Clouds and Precipitation. Springer, pp 10–73
  * Rajagopal (1959) Rajagopal E., 1959, Kolloid-Zeitschrift, 167, 17
  * Ramkrishna (2000) Ramkrishna D., 2000, Population balances: Theory and applications to particulate systems in engineering. Elsevier
  * Rigopoulos & Jones (2003) Rigopoulos S., Jones A. G., 2003, AIChE Journal, 49, 1127
  * Safronov (1972) Safronov V. S., 1972, Evolution of the protoplanetary cloud and formation of the earth and planets.
  * Sandu (2006) Sandu A., 2006, Aerosol Science Technology, 40, 261
  * Schumann (1940) Schumann T. E. W., 1940, Quarterly Journal of the Royal Meteorological Society, 66, 195
  * Scott (1968) Scott W. T., 1968, Journal of Atmospheric Sciences, 25, 54
  * Sheng & Shen (2006) Sheng C., Shen X., 2006, Journal of Aerosol Science, 37, 16
  * Shu & Osher (1988) Shu C.-W., Osher S., 1988, Journal of Computational Physics, 77, 439
  * Silk & Takahashi (1979) Silk J., Takahashi T., 1979, ApJ, 229, 242
  * Silk & White (1978) Silk J., White S. D., 1978, ApJ, 223, L59
  * Simmel et al. (2002) Simmel M., Trautmann T., Tetzlaff G., 2002, Atmospheric Research, 61, 135
  * Smith & Matsoukas (1998) Smith M., Matsoukas T., 1998, Chemical Engineering Science, 53, 1777
  * Smoluchowski (1916) Smoluchowski M. V., 1916, Zeitschrift fur Physik, 17, 557
  * Soong (1974) Soong S.-T., 1974, Journal of Atmospheric Sciences, 31, 1262
  * Suck & Brock (1979) Suck S. H., Brock J. R., 1979, Journal of Aerosol Science, 10, 581
  * Tanaka et al. (1996) Tanaka H., Inaba S., Nakazawa K., 1996, Icarus, 123, 450
  * Testi et al. (2014) Testi L., et al., 2014, Protostars and Planets VI, p. 339
  * Trautmann & Wanner (1999) Trautmann T., Wanner C., 1999, Atmospheric Environment, 33, 1631
  * Trubnikov (1971) Trubnikov B. A., 1971, Soviet Physics Doklady, 16, 124
  * Tzivion et al. (1999) Tzivion S., Reisin T. G., Levin Z., 1999, Journal of Computational Physics, 148, 527
  * Vericel & Gonzalez (2020) Vericel A., Gonzalez J.-F., 2020, MNRAS, 492, 210
  * Wang (1966) Wang C.-S., 1966, PhD thesis, California Institute of Technology.
  * Wang et al. (2007) Wang L.-P., Xue Y., Grabowski W. W., 2007, Journal of Computational Physics, 226, 59
  * Weidenschilling (1980) Weidenschilling S. J., 1980, Icarus, 44, 172
  * Wetherill (1990) Wetherill G. W., 1990, Icarus, 88, 336
  * Windmark et al. (2012) Windmark F., Birnstiel T., Güttler C., Blum J., Dullemond C. P., Henning T., 2012, A&A, 540, A73
  * Wolfram Research (14) Wolfram Research I., 14, Mathematica, Version 14.0, https://www.wolfram.com/mathematica
  * Zhang & Shu (2010) Zhang X., Shu C.-W., 2010, Journal of Computational Physics, 229, 3091
  * Zsom & Dullemond (2008) Zsom A., Dullemond C. P., 2008, A&A, 489, 931