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

HAL
open science

# General non-linear fragmentation with discontinuous Galerkin methods

Maxime Lombart, Charles-Edouard Bréhier, Mark Hutchison, Yueh-Ning Lee

► To cite this version:

Maxime Lombart, Charles-Edouard Bréhier, Mark Hutchison, Yueh-Ning Lee. General non-linear fragmentation with discontinuous Galerkin methods. Monthly Notices of the Royal Astronomical Society, 2024, 533, pp.4410-4434. 10.1093/mnras/stae2039. insu-04843505

HAL Id: insu-04843505

https://insu.hal.science/insu-04843505v1

Submitted on 17 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# General non-linear fragmentation with discontinuous Galerkin methods

Maxime Lombart Ⓟ,<fnref n="1" />,<fnref n="2" />⋆ Charles-Edouard Bréhier Ⓟ,<fnref n="3" /> Mark Hutchison Ⓟ,<fnref n="4" />,<fnref n="5" /> and Yueh-Ning Lee Ⓟ,<fnref n="1" />

Accepted 2024 August 23. Received 2024 August 6; in original form 2024 April 17

## ABSTRACT

Dust grains play a significant role in several astrophysical processes, including gas/dust dynamics, chemical reactions, and radiative transfer. Replenishment of small-grain populations is mainly governed by fragmentation during pair-wise collisions between grains. The wide spectrum of fragmentation outcomes, from complete disruption to erosion and/or mass transfer, can be modelled by the general non-linear fragmentation equation. Efficiently solving this equation is crucial for an accurate treatment of the dust fragmentation in numerical modelling. However, similar to dust coagulation, numerical errors in current fragmentation algorithms employed in astrophysics are dominated by the numerical overdiffusion problem – particularly in three-dimensional hydrodynamic simulations where the discrete resolution of the mass-density distribution tends to be highly limited. With this in mind, we have derived the first conservative form of the general non-linear fragmentation with a mass flux highlighting the mass transfer phenomenon. Then, to address cases of limited mass density resolution, we applied a high-order discontinuous Galerkin scheme to efficiently solve the conservative fragmentation equation with a reduced number of dust bins. An accuracy of 0.1–1 per cent is reached with 20 dust bins spanning a mass range of 9 orders of magnitude.

**Key words:** methods: numerical – (ISM:) dust, extinction.

## 1 INTRODUCTION

Dust grains play a fundamental role at all scales in astrophysics. Dust surfaces act as an efficient catalyst for the formation of H₂ (Cazaux & Tielens 2004), which impacts the star formation rate in galaxies (Yamasawa et al. 2011; Chen et al. 2018). Dust grains absorb, scatter and reemit stellar light, thereby governing the thermal balance between heating and cooling in star forming regions and protoplanetary discs (Mckee & Ostriker 2007; Andrews 2020). Gas and dust dynamics are closely linked through drag forces (Testi et al. 2014; Lesur et al. 2023a). During star formation, large grains decouple dynamically from the gas and concentrate in regions of high gas density (Lebreuilly et al. 2020). In protoplanetary discs, the decoupling of large grains from the gas leads to momentum transfer and radial drift towards pressure maxima (Weidenschilling 1977; Lesur et al. 2023a). When dust becomes sufficiently concentrated in discs, it produces a ‘backreaction’ on the gas which is central for the development of the streaming instability and the eventual formation of planetesimals (Youdin & Goodman 2005; Gonzalez, Laibe & Maddison 2017; Squire & Hopkins 2020; Lesur et al. 2023a, and references therein). Importantly, the efficiency of all these physical processes depends on the grain-size distribution and how it evolves in time. Thus, accurate dust modelling is a much needed feature of modern astrophysical simulations.

Because the grain-size distribution plays such a key role in the above phenomenon, many studies have been devoted to understanding how poly-disperse distributions evolve in time. Evolutionary changes can be dynamical (e.g. advection and diffusion) or collisional (e.g. coagulation and fragmentation) in nature. For purposes of this study, we will focus exclusively on fragmentation. Fragmentation is typically modelled in one of three ways: (i) spontaneous breaking driven by an external force, such as radiative force (Hoang 2019; Hirashita & Hoang 2020), modelled by the *linear fragmentation* equation; (ii) collision between two grains where only one grain fragments, modelled by the *non-linear fragmentation* equation (Kostoglou & Karabelas 2000; Banasiak, Lamb & Laurençot 2019; Lombart, Hutchison & Lee 2022); (iii) collision between two grains where both grains can fragment, modelled by the *general non-linear fragmentation* equation (Safronov 1972; Blum 2006; Hirashita & Yan 2009). Note the third model is a generalization of the second one, which produces a size distribution of solids that can be generically parametrized from experiments. Two main outcomes observed in laboratory experiments are the complete destruction of the two grains or a partial destruction with mass transfer (Güttler et al. 2010; Bukhari Syed et al. 2017; Blum 2018). The phenomenon of mass transfer is important since, during a fragmentation event, a grain can increase in mass – even at high impact velocities. For example, several studies have demonstrated that mass transfer during fragmentation can overcome the so-called bouncing and fragmentation barriers (Windmark et al. 2012; Garaud et al. 2013). Only the general non-linear fragmentation equation can effectively model fragmentation with mass transfer, making it essential for a

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1">¹Department of Earth Sciences, National Taiwan Normal University, Taipei 11677, Taiwan</footnote>

<footnote n="2">²Université Paris Cité, Université Paris-Saclay, CEA, AIM, F-91191, Gif-sur-Yvette, France.</footnote>

<footnote n="3">³Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, 64013 Pau, France</footnote>

<footnote n="4">⁴Universitäts-Sternwarte, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germany</footnote>

<footnote n="5">⁵Hochschule für angewandte Wissenschaften München, Lothstraße 34, D-80335 München, Germany</footnote>

<footnote>* E-mail: maxime.lombart@cea.fr</footnote>
full treatment of the problem. The general non-linear fragmentation equation is formalized within the framework of the Smoluchowski-like equation by a mean-field approach (Safronov 1972; Gillespie & List 1978; Feingold, Tzivion (Tzitzvashvili) & Leviv 1988; Blum 2006; Hirashita & Yan 2009; Jacobson 2011; Banasiak et al. 2019). The model considers spherical grains of the same composition. Because the general non-linear fragmentation equation does not have a generic analytical solution, astrophysical problems inevitably require numerical solutions.

Unfortunately, incorporating accurate numerical solutions to the fragmentation (and/or coagulation), together with tracking the dust-size distribution in three-dimensional (3D) multiphysics hydrodynamics simulations, is currently out of reach. So far only 3D hydrodynamics simulations with a mono-disperse model of dust coagulation and fragmentation have been performed for protoplanetary discs (Vericel et al. 2021). Meanwhile, poly-disperse models of dust coagulation and fragmentation have been confined to one-dimensional (1D) or two-dimensional (2D) hydrodynamic simulations (Suttner & Yorke 2001; Brauer, Dullemond & Henning 2008; Dražkowska et al. 2019; Kobayashi & Tanaka 2021; Lebreuilly et al. 2022; Stammler & Birnstiel 2022; Tu, Li & Lam 2022; Robinson, Booth & Owen 2024) primarily due to the near ubiquitous use of piecewise constant functions to model the dust-size distribution. Piecewise constant functions are known to suffer from a numerical overdiffusion for insufficient mass grid resolution (Grabowski 2022; Birnstiel 2024). This overdiffusion stems from the difficulty to handle the complexity of the integrodifferential and the non-linear properties of the coagulation and the fragmentation equations with only a low-order approximation of the continuous dust-size distribution. Moreover, these algorithms account for each pair-wise collision by redistributing the mass over the mass grid to conserve the total mass. The combinatorics treatment of the collisions requires the need of a large number of mass bins to keep high accuracy. Therefore, these algorithms need a high resolution of the mass grid (more than 100 mass bins) to accurately follow the evolution of the dust-size distribution. However, 3D hydrodynamic codes, such as RAMSES (Teyssier 2002) or PHANTOM (Price et al. 2018), can only handle a few tens of dust mass bins for a multiphysics gas and dust simulation. Further generation of Exascale code will reach further performance, but the energy, the carbon impact and the computational costs of the simulation should remain acceptable. Current algorithms for poly-disperse dust evolution and 3D hydrodynamic codes are simply incompatible. In contrast, numerical schemes using high-order approximations of the continuous dust-size distribution are able to efficiently solve the coagulation and the fragmentation equations with a reduced number of mass bins. This study contains two main parts: (i) the derivation of the flux of fragmentation in mass space, which encompasses all the combinatorics of the pair-wise collisions, in order to obtain the general non-linear fragmentation equation in conservative form and (ii) the application of a high-order solver based on the discontinuous Galerkin (DG) method, inspired by the recent works of Liu, Gröpler & Warnecke (2019) and Lombart et al. (2022). The DG method is shown to efficiently solve the general non-linear fragmentation equation with a reduced number of mass bins while still maintaining high accuracy.

The paper is structured as follows. Properties of the general non-linear fragmentation equation and derivation of its conservative form are presented in Section 2. The DG method applied to the conservative form is presented in Section 3. In Section 4, we analyse the numerical performance of our algorithm on some test cases, including how it copes with the numerical diffusion problem.

Section 5 presents the applicability of the DG method to treat dust fragmentation in astrophysics.

## 2 GENERAL NON-LINEAR FRAGMENTATION

The fragmentation process resulting from the collision of two grains of arbitrary size is described by the general non-linear fragmentation equation, also known as the ‘collision-induced breakup’ or the ‘stochastic breakage’ equation in atmospheric science and mathematics communities (Safronov 1972; Gillespie & List 1978; Feingold et al. 1988; Blum 2006; Giri & Laurençot 2021a). The general non-linear fragmentation model is a natural extension of the non-linear fragmentation model that describes how a small grain fragments after collision with a large grain (Cheng & Redner 1990; Kostoglou & Karabelas 2000; Ernst & Pagonabarraga 2007; Lombart et al. 2022). This fragmentation model is described by a non-linear partial integrodifferential hyperbolic equation that depends on two functions: (i) the fragmentation kernel, and (ii) the distribution function of fragments. This general non-linear fragmentation model was initially formalized in astrophysics by Safronov (1972) and in atmospheric science by List & Gillespie (1976) to study the evolution of the drop size distribution in clouds. Currently in the literature, only one analytical solution exists for a constant fragmentation kernel and a specific form of the distribution of fragments (Feingold et al. 1988). Several recent mathematical works have been dedicated to the study of the general non-linear fragmentation equation, including proving the existence and uniqueness of mass-conserving solutions for a large class of collision kernels and fragment distribution functions (Giri & Laurençot 2021a, b). However, no exact solutions exist for the ballistic collision kernel studied in astrophysics (Safronov 1972; Tanaka, Inaba & Nakazawa 1996; Dullemond & Dominik 2005; Kobayashi & Tanaka 2010; Stammler & Birnstiel 2022) and in atmospheric science (Hu & Srivastava 1995; McFarquhar 2004; Prat & Barros 2007; Jacobson 2011; Khain & Pinsky 2018; Grabowski 2022). Therefore, numerical solutions are required for the ballistic collision kernel. The detailed expression of the ballistic kernel is given in Section 2.3. In this work, numerical results are compared to the exact solution obtained for the constant collision kernel and a distribution of fragments taken from Feingold et al. (1988, Section 4.3) or a multiplicative collision kernel with a power-law distribution of fragments (Section 4.4).

### 2.1 Collision outcomes for fragmentation

In astrophysics, the outcome of two colliding dust grains leading to the formation of fragments is described by one of three scenarios, illustrated in Fig. 1 (Güttler et al. 2010; Windmark et al. 2012; Blum 2018; Birnstiel 2024). First, destructive fragmentation accounts for cases where both grains fragment totally or partially. The second scenario describes mass transfer events where the resulting fragments come from only one body. This scenario has been observed in experiments when a small grain collides with a larger one (Bukhari Syed et al. 2017). During the collision, the smaller grain fragments and leaves a portion of its mass stuck to the larger (intact) grain, causing the mass of the larger grain to increase. The third scenario is the same as the second, but now with cratering of the larger grain. For sufficiently high-velocities, the small grain both transfers mass to and excavates mass from the large grain. Since the non-linear fragmentation equation (Kostoglou & Karabelas 2000; Banasiak et al. 2019; Lombart et al. 2022) cannot account for mass transfer or the fragmentation of both grains, the above scenarios can only be modelled by the general non-linear fragmentation equation. This is of
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustration of the collision outcomes leading to the formation of fragments. Depending on the differential velocity of the two colliding grains, three outcomes are possible: destructive fragmentation, mass transfer, and mass transfer with cratering.</figcaption></figure>

particular interest to the planetary science community because mass transfer may be an avenue to circumvent the fragmentation barrier (Windmark et al. 2012; Garaud et al. 2013) and provide larger seeds for planetesimal formation.

## 2.2 Original rate equation

We begin by formulating a rate equation for fragmentation that describes the gains and losses in number density as a function of grain mass. An alternative form is given later in Section 2.5 that has improved mass conservation properties. Hereafter, we will distinguish between these two formulations using the labels ‘original’ and ‘alternative’, respectively.

We consider a volume composed of a large number of uniformly distributed grains undergoing fragmentation through collisions. The number of grains is considered sufficiently large to be treated in a statistical sense, while the even spatial distribution of the grains permits the exclusion of their motion from consideration. Moreover, colliding grains and fragments are considered to be spherical.

We denote $K(m', m'')$ to be the fragmentation kernel that encodes the collision frequency of grains leading to fragmentation with dimension $[length]^3[time]^{-1}$. The function $K$ is symmetric in variables $m$ and $m'$. Therefore, the mean number of collisions per unit time and unit volume between the grains of mass in the ranges $[m', m' + dm']$ and $[m'', m'' + dm'']$ is given by

$$K(m', m'')n(m', t)n(m'', t)dm'dm'' \text{ (collision rate)}, \quad (1)$$

where $m'$ and $m''$ are the initial masses of the colliding grains and $n(m', t)$ and $n(m'', t)$ are the number densities function per unit mass of grains in mass ranges $[m', m' + dm']$ and $[m'', m'' + dm'']$, respectively. We denote by $\tilde{b}(m; m', m'')$ the distribution of fragments of mass $m$ produced by the collision of grains of mass $m'$ and $m''$. The only physical constraint on the function $\tilde{b}$ is that the mass of a fragment cannot exceed the total mass of the colliding grains; therefore, one has $b(m; m', m'') = 0$ if $m > m' + m''$. In the derivations that follow, this condition is represented by the operator $\mathbb{1}$, defined as

$$\mathbb{1}_{m'+m'' \ge m} = \begin{cases} 1 & \text{if } m' + m'' \ge m, \\ 0 & \text{otherwise.} \end{cases} \quad (2)$$

Therefore, the distribution of fragments writes $\mathbb{1}_{m'+m'' \ge m}\tilde{b}(m; m', m'')$. The term $\tilde{b}(m; m', m'')dm$ is the mean number of fragments of mass whose masses reside in the range $[m, m + dm]$ produced in a collision of one pair of grains of mass $m'$ and $m''$. The dimension of $b$ is $[mass]^{-1}$. The formation rate of grains of mass $[m, m + dm]$ by fragmentation of larger grains of

masses $[m', m' + dm']$ and $[m'', m'' + dm'']$ is equal to the product of the collision rate (equation 1) with the mean number of fragments produced by the collision, $\tilde{b}(m; m', m'')dm$. The formation rate of grains within mass range $[m, m + dm]$ is obtained by considering all collisions, and is expressed as

$$\left[ \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{m'+m'' \ge m} \tilde{b}(m; m', m'') K(m', m'') \times n(m', t)n(m'', t) dm'dm'' \right] dm \text{ (formation rate)}, \quad (3)$$

where the factor $1/2$ prevents grain pairs from being double counted. The loss rate of grains of mass $[m, m + dm]$ by collision with all other grains is given by

$$\left[ \int_{0}^{\infty} K(m, m') n(m, t) n(m', t) dm' \right] dm \text{ (loss rate).} \quad (4)$$

Therefore, the rate equation for the number density of grains of mass $[m, m + dm]$ at time $t$ is obtained by balancing the formation and loss rates

$$\frac{\partial(n(m, t)dm)}{\partial t} = \left[ \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{m'+m'' \ge m} \tilde{b}(m; m', m'') K(m', m'') \times n(m', t)n(m'', t) dm'dm'' \right] - \left[ \int_{0}^{\infty} K(m, m') n(m, t) n(m', t) dm' \right] dm. \quad (5)$$

By dividing by $dm$, we obtain the continuous general non-linear fragmentation equation originally formalized in a mean-field approach by List & Gillespie (1976) and Gillespie & List (1978):

$$\frac{\partial n(m, t)}{\partial t} = \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{m'+m'' \ge m} \tilde{b}(m; m', m'') K(m', m'') \times n(m', t)n(m'', t) dm'dm'' - n(m, t) \int_{0}^{\infty} K(m, m') n(m', t) dm'. \quad (6)$$

To obtain a dimensionless general non-linear fragmentation equation, we first define $n_0(m) = n(m, 0)$ to be the initial number density per mass unit. The total mass density $M$, the total number density of particles $N_0$ and the mean mass of the initial distribution $m_0$ can then be written as

$$M \equiv \int_{0}^{\infty} mn_0(m)dm, \quad N_0 \equiv \int_{0}^{\infty} n_0(m)dm, \quad m_0 \equiv \frac{M}{N_0}. \quad (7)$$

Then, we define the dimensionless variables and functions as follows:

$$\begin{cases} x \equiv \frac{m}{m_0}, \quad y \equiv \frac{m'}{m_0}, \quad z \equiv \frac{m''}{m_0}, \\ K(x, y) = \frac{K(m, m')}{K(m_0, m_0)}, \quad \tau = K(m_0, m_0)N_0t, \\ f(x, \tau) = \frac{m_0}{N_0} n(m, t), \quad b(x; y, z) = m_0 \tilde{b}(m; m', m''), \end{cases} \quad (8)$$
where $K(m_0, m_0)$ is a constant with dimensions of the kernel $[length]^3[time]^{-1}$. To be consistent with the existing literature in mathematics (Banasiak et al. 2019; Giri & Laurençot 2021a), we use the variables $x, \tau$, and $f$ for the dimensionless mass, time, and number density, respectively, to write the general non-linear fragmentation equation in the following dimensionless form:

$$ \frac{\partial f(x, \tau)}{\partial \tau} = \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} \int_{1_{y+z \ge x}} \int_{1_{y+z \le x}} b(x; y, z) K(y, z) f(y, \tau) f(z, \tau) dy dz - f(x, \tau) \int_{0}^{\infty} K(x, y) f(y, \tau) dy. \quad (9) $$

Unless otherwise noted, we will use those dimensionless variables for the remainder of the paper.

## 2.3 Fragmentation kernels

The fragmentation kernel describes the collision rate per unit volume of two grain masses. The expression of the kernel is determined by the mechanism driving the collisions between grains (e.g. turbulence, radial drift, vertical settling). The kernel depends on the mass of the grains, and on the properties of any surrounding gas, such as temperature, pressure, and the characteristics of the flow field. In astrophysics, the collision between grains is modelled by the ballistic kernel,

$$ K(x, y) = E_{\text{coll}}(x, y) P_{\text{frag}}(x, y, \Delta v) \sigma(x, y) \Delta v(x, y), \quad (10) $$

where $\sigma(x, y)$ is the geometric cross-section, $\Delta v$ is the mean relative velocity between two grains of mass $x$ and $y$, while $P_{\text{frag}}$ denotes the probability that the two colliding grains fragment. The collision efficiency, $E_{\text{coll}}$, describes the probability of two grains embedded in a flow field colliding. For instance, a large grain surrounded by gas in the Stokes regime has a small probability of collision with a small grain dragged along the gas stream lines. The limiting case for impact between large and small grains is referred to as the grazing collision trajectory in atmospheric science for droplets (Pruppacher & Klett 2010; Wang 2013; Khain & Pinsky 2018) or grazing impact in astrophysics (Paszun & Dominik 2009; Wada et al. 2009). The standard definition of the collision efficiency is (Paszun & Dominik 2009; Wada et al. 2009; Pruppacher & Klett 2010)

$$ E_{\text{coll}}(s_x, s_y) = \frac{p_{\text{impact}}(s_x, s_y)^2}{(s_x + s_y)^2}, \quad (11) $$

where $s_x$ and $s_y$ are the radius of the grains of masses $x$ and $y$. $p_{\text{impact}}(s_x, s_y)$ is the impact parameter and depends on the size of the two colliding grains. For a collision event, the impact parameter is defined as the projected distance between the centers of mass of the grains in the perpendicular direction to the collision velocity (Wada et al. 2009). By considering spherical grains, we can directly link $E_{\text{coll}}(s_x, s_y)$ and $E_{\text{coll}}(x, y)$, as was done by Pinsky, Khain & Shapiro (2001) by solving the equations of motion of small grains around large grains. In Paszun & Dominik (2009) and Wada et al. (2009), the impact parameter is evaluated as a percentage of the maximum impact parameter $p_{\text{impact,max}} = s_x + s_y$. A value of $E_{\text{coll}}(x, y)$ close to zero means that the small grain follows the gas stream lines. This formalism has been adapted to study the collision between planetesimals and small grains (Guillot, Ida & Ormel 2014; Visser & Ormel 2016). However, for purposes of this study, we will consider only head-on collisions, $E_{\text{coll}}(x, y) = 1$.

The fragmentation probability $P_{\text{frag}}$ is usually defined as

$$ P_{\text{frag}} = \begin{cases} 1 & \text{if } \Delta v > \Delta v_{\text{th}}, \\ 0 & \text{if } \Delta v < \Delta v_{\text{th}}, \end{cases} \quad (12) $$

where $\Delta v$ is the differential velocity between the two colliding grains and $\Delta v_{\text{th}}$ is the threshold differential velocity for which fragmentation occurs. The value of the threshold velocity is determined from experiments (Güttler et al. 2010; Blum 2018) or from theoretical works (Jones, Tielens & Hollenbach 1996; Ormel et al. 2009). In this study, because we only consider fragmentation, we set $P_{\text{frag}} = 1$.

## 2.4 Distribution of fragments

The evolution of the number density $f(x, \tau)$ in equation (6) depends on two physical parameters, the collision kernel $K$ and the distribution function of fragments $b(x; y, z)$ resulting from a collision between two grains of mass $y$ and $z$ (Gillespie & List 1978; Feingold et al. 1988). The function $b$ is symmetric in the mass variables of the two colliding grains

$$ b(x; y, z) = b(x; z, y). \quad (13) $$

Because we do not consider sublimation, the fragmentation process must satisfy the following two mass conservation constraints. First, the mass of a fragment cannot exceed the total mass of the colliding grains:

$$ b(x; y, z) = 0 \text{ if } x > y + z. \quad (14) $$

Secondly, the total mass of the fragments must be equal to the total mass of the colliding grains

$$ \int_{0}^{y+z} x b(x; y, z) dx = y + z. \quad (15) $$

Note that these constraints on the distribution of fragments still allows for the mass transfer phenomenon to occur in equation (6). The number of fragments produced for each collision is defined by

$$ N_{\text{frag}}(y, z) \equiv \int_{0}^{y+z} b(x; y, z) dx. \quad (16) $$

By definition, a fragmentation event produces at least two fragments, therefore $N_{\text{frag}} \ge 2$. The extreme case of $N_{\text{frag}} = 2$ occurs when a small grain breaks in two and one piece is absorbed by the larger grain.

## 2.5 Alternative rate equation

In atmospheric science, the distribution of fragments from droplets collisions is obtained by experiment (Low & List 1982). However, the usual formulation of the distribution of fragment does not strictly satisfy the local mass conservation constraint in equation (15) (Brown 1986; Feingold et al. 1988). To overcome this problem, an alternative rate equation has been formalized by List & Gillespie (1976) and Gillespie & List (1978) where the loss rate of droplets (or, in astrophysics, dust grains) is written with the distribution function of fragments. The loss term in equation (6) is derived by counting all the collisions between droplets of mass $[x, x + dx]$ and any other droplet. In Gillespie & List (1978), the loss rate is defined as the destruction rate of droplets of mass $x$ by fragmentation. We consider the fragmentation of a pair of droplets of masses $(x, y)$. The total mass of fragments created by the broken pair is given by

$$ \int_{0}^{x+y} z b(z; x, y) dz. \quad (17) $$
Therefore, the mass of fragments created by all the broken pairs of droplets of masses (x, y) per unit time is obtained by multiplying equation (17) with the mean number of collisions per unit time and per unit volume between the droplets of mass [x, x + dx] and [y, y + dy] in equation (1), i.e.

$$K(x, y)f(x, \tau)f(y, \tau)dydx \times \int_{0}^{x+y} zb(z; x, y)dz. \quad (18)$$

By the mass conservation, equation (18) equals the mass of all the broken pairs of droplets per unit time. Then, dividing by the mass of one pair x + y, we obtain the number of all broken pairs of masses [x, x + dx] and [y, y + dy] per unit time, i.e.

$$\frac{1}{x+y} \left[ K(x, y)f(x, \tau)f(y, \tau)dydx \times \int_{0}^{x+y} zb(z; x, y)dz \right]. \quad (19)$$

Equivalently, equation (19) gives the number of droplets of mass [x, x + dx] which break per unit time. We obtain the total rate of loss of droplets of mass [x, x + dx] by integrating over y,

$$\left[ \int_{0}^{\infty} K(x, y)f(x, \tau)f(y, \tau) \frac{1}{x+y} \int_{0}^{x+y} zb(z; x, y)dzdy \right] dx, \quad (20)$$

and the alternative rate equation writes

$$\frac{\partial f(x, \tau)}{\partial t} = \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{y+z \ge x} b(x; y, z) K(y, z) f(y, \tau) f(z, \tau) dydz - f(x, \tau) \int_{0}^{\infty} \frac{K(x, y) f(y, \tau)}{x+y} \int_{0}^{x+y} zb(z; x, y) dzdy. \quad (21)$$

The mass conservation with the original equation equation (9) is ensured by the local mass conservation equation (15). However, the alternative rate equation is of particular interest because Feingold et al. (1988) proved that the mass is conserved for any choice of collision kernel $K$ and distribution of fragments $b$, because the function $b$ appears in both term on the right-hand side of equation (21). It means that the total mass is conserved even if locally the mass is not strictly conserved for each breakup event (Hu & Srivastava 1995). Note that equation (9) and equation (21) are two different models to describe the general non-linear fragmentation equation. In the case where the distribution of fragments respects exactly the local mass conservation, we obtain equation (9) from equation (21).

## 2.6 Conservative form

### 2.6.1 Original equation

The fragmentation process is mass conserving, meaning no mass is lost during the process. In order to solve the general non-linear fragmentation equation by using robust numerical schemes, such as finite volume methods or DG methods, which conserve the total mass at machine precision, it is required to derive the conservative form of equation (9), as an hyperbolic conservation law. To our knowledge, there have been no studies that have derived the conservative form of the general non-linear fragmentation equation (equation 9), which we do now. Note that the calculations are formal and the rigorous justification is beyond the scope of this work.

The aim is to find the expression of a mass flux, $F_{frag}$, such that

$$\frac{\partial g(x, \tau)}{\partial \tau} + \frac{\partial F_{frag}[g](x, \tau)}{\partial x} = 0, \quad (22)$$

where $g(x, \tau) \equiv xf(x, \tau)$ is the mass density of grains per unit mass. The term $F_{frag}[g](x, \tau)$ describes a combined mass flux at time $\tau$ crossing mass $x$ resulting from all collisions involving grain pairs with a total mass greater than $x$ and which fragment to produce grains with mass lower than $x$. A relation for the mass flux follows directly from integrating equation (22) with respect to $x$,

$$F_{frag}[g](x, \tau) = - \int_{0}^{x} \frac{\partial g(x', \tau)}{\partial \tau} dx', \quad (23)$$

where the constant $F_{frag}[g](0, \tau)$ is equal to zero since the process of grain nucleation from gas phase to small nano-particles is neglected.

Multiplying equation (9) by $x$ and integrating over $x \in [0, \infty)$, the mass flux is written as

$$F_{frag}[g](x, \tau) = -\frac{1}{2} \int_{0}^{x} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{y+z \ge x'} x' b(x'; y, z) K(y, z) \frac{g(y, \tau)g(z, \tau)}{yz} dydz dx' + \int_{0}^{x} \int_{0}^{\infty} K(y, z) g(z, \tau) \frac{g(y, \tau)}{y} dydz. \quad (24)$$

A non-zero mass flux across mass $x$ only has physical meaning if the mass of the colliding grains is greater than $x$. We can make this explicit by introducing the operator $\mathbb{1}_{y+z \ge x}$ to the right-hand side of equation (24),

$$F_{frag}[g](x, \tau) = -\frac{1}{2} \int_{0}^{x} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{y+z \ge x} x' b(x'; y, z) K(y, z) \frac{g(y, \tau)g(z, \tau)}{yz} dydz dx' + \int_{0}^{x} \int_{0}^{\infty} \mathbb{1}_{y+z \ge x} K(y, z) g(z, \tau) \frac{g(y, \tau)}{y} dydz, \quad (25)$$

noting that the operator $\mathbb{1}_{y+z \ge x}$ imposes $\mathbb{1}_{y+z \ge x'} = 1$, since $x' \in [0, x]$. Thus, the conservative form of equation (9) is equation (22) with the flux given in equation (25). Inserting this expression for the flux into equation (22) and applying the Leibniz integral rule correctly reduces to equation (9), as we would expect. To our knowledge, this is the first time that $F_{frag}[g](x, \tau)$ has been derived.

The first term on the right-hand side of equation (25) describes the flux of mass density across the mass $x$ by fragmentation of grains of masses $y$ and $z$ with $y + z \ge x$ producing fragments of mass lower than $x$, as illustrated in Fig. 2 with the term $F_1(x, \tau)$. The second term describes the flux of mass through $x$ by fragmentation with mass transfer of grains of masses $y \le x$ and $z$ with $y + z \ge x$ producing fragments of mass greater than $x$, as shown in Fig. 2 with the term $F_2(x, \tau)$. The fragmentation flux is a balance between (i) the negative flux of fragments coming from the break-up of larger grains and (ii) the positive flux of smaller grains gaining mass through mass transfer. In the second term, the gain of mass, due to the mass transfer can be interpreted as a form of coagulation. Indeed the flux of mass in the coagulation process (Tanaka et al. 1996; Filbet & Laurençot 2004; Liu et al. 2019; Lombart & Laibe 2021) is mathematically similar to the second term. The conservative form of the general non-linear fragmentation equation, with the flux of mass equation
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Illustration of the fragmentation flux in equation (25). The flux is a balance between two terms. The first term describes the production of particles with mass $x_i < x$ from the collision of two particles with total mass $y+z \ge x$. This first term is called $F_1$. The second term, $F_2$, describes the formation of particles with mass greater than $x$ due to the mass transfer process from the collision of two particles with total mass $y+z \ge x$ but with one particle of mass $y \le x$. The fragmentation flux is a balance between the production of small grains and the production of larger grain due to the mass transfer phenomenon.</figcaption></figure>

(25), directly highlights the mass transfer phenomenon presented in Section 2.1.

### 2.6.2 Alternative equation

The alternative rate equation in equation (21) can also be written in a conservative form with an associated mass flux

$$ 
\left\{
\begin{aligned}
& \frac{\partial g(x, \tau)}{\partial \tau} + \frac{\partial F_{\text{frag,alt}}[g](x, \tau)}{\partial x} = 0, \\
& F_{\text{frag,alt}}[g](x, \tau) = \\
& \quad -\frac{1}{2} \int_0^x \int_0^\infty \int_0^\infty x' b(x'; y, z) \mathbb{1}_{y+z \ge x} \mathcal{K}(y, z) \\
& \quad \times \frac{g(y, \tau)g(z, \tau)}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \\
& \quad + \int_0^x \int_0^\infty \int_0^{z+y} \mathbb{1}_{y+z \ge x} x' b(x'; z, y) \mathcal{K}(z, y) \\
& \quad \times \frac{g(z, \tau)g(y, \tau)}{y(z+y)} \mathrm{d}x' \mathrm{d}y \mathrm{d}z.
\end{aligned}
\right.
\quad (26) $$

The method to derive the flux $F_{\text{frag,alt}}$ is similar to Section 2.6.1. The alternative rate equation in equation (21) is only used for benchmarking the DG method (Section 4) since an analytical solution is derived for the constant collision kernel and a specific distribution of fragments which does not strictly satisfy the local mass conservation in equation (15).

## 2.7 Analytical solution of equation (21)

The only analytical solution to the general non-linear fragmentation equation that we could find in the literature comes from Feingold et al. (1988) in the atmospheric community, who derived an exact solution for equation (9) assuming a constant kernel, $K_0$, and the following fragment distribution function:

$$ \tilde{b}(m; m', m'') = \tilde{\gamma}^2(m' + m'')e^{-\tilde{\gamma}m}, \quad (27) $$

where $K_0$ is a constant with dimension $[\text{length}]^3[\text{time}]^{-1}$ and $\tilde{\gamma}$ is a positive parameter characterizing the physics of the fragmentation process. This distribution of fragments physically represents

the case where the number of fragments depends on the mass of the two colliding droplets $m'$ and $m''$ and the shape of the distribution decreases exponentially with $m$. It means that for one collision a large number of small droplets will be generated. The decaying exponential term in equation (27) allows the distribution of fragment to approximatively respect the physical constraint $\forall m > m' + m''$, $\tilde{b}(m; m'm'') \approx 0$. However, equation (27) does not fully satisfy the local mass conservation condition in equation (15)

$$ 
\begin{aligned}
& \int_0^{m'+m''} m \tilde{b}(m; m', m'') dm \\
& = (m' + m'') \left[ 1 - \underbrace{e^{-\tilde{\gamma}(m'+m'')}(1 + \tilde{\gamma}(m' + m''))}_{\epsilon(\tilde{\gamma}, m', m'')} \right] \\
& < m' + m'',
\end{aligned}
\quad (28) $$

since here $\epsilon(\tilde{\gamma}, m', m'')$ is small but non-zero. Feingold et al. (1988) found that $\epsilon(\tilde{\gamma}, m', m'')$ is a few orders of magnitude less than unity whenever the number of fragments is greater than 10. This can easily be shown for a given mass range and number of fragments by evaluating the integral in equation (16), numerically solving for $\tilde{\gamma}$, and evaluating $\epsilon$. This discussion will be important later for a numerical test that requires a value of $\tilde{\gamma}$ (see Section 4.3).

Let us now show the steps for the derivation of the analytical solution. The dimensionless kernel and distribution of fragments write

$$ \mathcal{K} = 1, \quad b(x; y, z) = \gamma^2(y+z)e^{-\gamma x}, \quad (29) $$

where $\gamma \equiv m_0 \tilde{\gamma}$. As the local mass conservation is not exactly respected, we substitute the distribution of fragments and constant kernel into the alternative rate equation equation (21)

$$ 
\begin{aligned}
\frac{\partial f(x, \tau)}{\partial \tau} = & \frac{\gamma^2}{2} \int_0^\infty \int_0^\infty (y+z) e^{-\gamma x} f(y, \tau) f(z, \tau) \mathrm{d}y \mathrm{d}z \\
& - f(x, \tau) \gamma^2 \int_0^\infty \int_0^\infty f(y, \tau) z e^{-\gamma z} \mathrm{d}z \mathrm{d}y,
\end{aligned}
\quad (30) $$

where the operator $\mathbb{1}_{y+z \ge x}$ is contained in the expression of $b$. Indeed, the function $b$, approximatively respects the physical constraint $\forall x > y + z$, $b(x; y, z) \approx 0$.

By continuing the development, we obtain

$$ \frac{\partial f(x, \tau)}{\partial \tau} = \gamma^2 e^{-\gamma x} \mathcal{M} \mathcal{N}(\tau) - f(x, \tau) \mathcal{N}(\tau), \quad (31) $$

where

$$ \mathcal{M} \equiv \int_0^\infty x f(x, \tau) \mathrm{d}x = 1, \quad \mathcal{N}(\tau) \equiv \int_0^\infty f(x, \tau) \mathrm{d}x, \quad (32) $$

are the dimensionless total mass density and total number of particles. Then, by integrating equation (31) over $x$, we obtain an evolution equation for $\mathcal{N}(\tau)$

$$ \frac{\mathrm{d}\mathcal{N}(\tau)}{\mathrm{d}\tau} = \gamma\mathcal{N}(\tau) - \mathcal{N}(\tau)^2. \quad (33) $$

The solution to equation (33) writes

$$ \mathcal{N}(\tau) = \frac{e^{\gamma\tau}}{1 + \frac{e^{\gamma\tau - 1}}{\gamma}}, \quad (34) $$

where $\mathcal{N}(0) = 1$. We remark that, for a given time, $\mathcal{N}(\tau)$ increases with $\gamma$, a positive integer that characterizes the production of
fragments. Finally, by injecting equation (34) into equation (31) and solving the differential equation, we obtain the exact solution of equation (30):

$$f(x, \tau) = \frac{f(x, 0) + \gamma(e^{\gamma\tau} - 1)e^{-\gamma x}}{1 + \frac{e^{\gamma\tau} - 1}{\gamma}}, \quad (35)$$

where $f(x, 0)$ is the initial number density. Note that the exact solution only retains a memory of the initial condition for a finite time since, as $t \to \infty$, the solution converges to the distribution of fragments

$$f(x, \tau) = \gamma^2 e^{-\gamma x}. \quad (36)$$

These characteristics of the exact solution suggest that the solutions to the general non-linear fragmentation equation depend on the initial condition and the evolution is controlled by the distribution of fragments (Feingold et al. 1988). This evolution is also observed for the multiplicative kernel and a power-law distribution of fragments as shown in Section 4.4.

## 3 DG METHOD

The DG method (Cockburn & Shu 1989; Zhang & Shu 2010; Liu et al. 2019) is an efficient numerical method to solve the non-linear fragmentation equation (Lombart et al. 2022) and the Smoluchowski coagulation equation (Lombart & Laibe 2021) with a reduced number of mass bins. This method is therefore well adapted to treat the coagulation and the fragmentation processes in 3D hydrodynamic simulations. Our objective is to extend the work from Lombart et al. (2022) by applying the DG method to the general non-linear fragmentation equation in equation (9). A complete description of the DG method, in the astrophysical context, is given in Lombart & Laibe (2021). Here, we outline only the principal steps for the general non-linear fragmentation equation.

### 3.1 Summary

The unbounded mass interval in the fragmentation equation is reduced to a physical mass range $[x_{min} > 0, x_{max} < \infty]$ and divided into $N$ bins. Each bin is defined by $I_j = [x_{j-1/2}, x_{j+1/2}]$ for $j \in [1, N]$. The size and the centre position of each bin $j$ are given, respectively, by $h_j = x_{j+1/2} - x_{j-1/2}$ and $x_j = (x_{j+1/2} + x_{j-1/2})/2$. In each bin $j$, the unknown mass density function $g$ is approximated by polynomials of order $k$ defined as a linear combination of Legendre polynomials $\phi_i$ with $i \in [0, k]$,

$$\forall x \in I_j, \quad g(x, \tau) \approx g_j(x, \tau) = \sum_{i=0}^{k} g_i^j(\tau) \phi_i(\xi_j(x)), \quad (37)$$

where the function $\xi$ maps the bin interval $I_j$ into the reference interval $[-1, 1]$ where standard Legendre polynomials are defined. Equation (22) is multiplied by the Legendre polynomials basis function vector $\phi(\xi_j(x)) = [\phi_0(\xi_j(x)), \dots, \phi_k(\xi_j(x))]^T$, then integrated over each bin $I_j$. Therefore, the DG method determines the evolution of the components $g_j^i$ by solving the following equation

$$\int_{I_j} \frac{\partial g_j}{\partial \tau} \phi(\xi_j(x)) dx - \underbrace{\int_{I_j} F_{\text{frag}}[g](x, \tau) \frac{\partial \phi(\xi_j(x))}{\partial x} dx}_{\text{integral of the flux, see Section 3.1.3}} + F_{\text{frag}}[g](x_{j+1/2}, \tau) \phi(\xi_j(x_{j+1/2})) - F_{\text{frag}}[g](x_{j-1/2}, \tau) \phi(\xi_j(x_{j-1/2})) = 0, \quad (38)$$

where the term ‘integral of the flux’ is obtained through integration by parts.

### 3.1.1 Description of the flux

The fragmentation flux, given by equation (25), is truncated into the physically relevant mass range $[x_{min}, x_{max}]$ by replacing 0 and $\infty$ by $x_{min}$ and $x_{max}$, respectively,

$$F_{\text{frag,nc}}[g](x, \tau) = -\frac{1}{2} \int_{x_{\text{min}}}^{x} \int_{x_{\text{min}}}^{x_{\text{max}}} \int_{x_{\text{min}}}^{x_{\text{max}}} x' b(x'; y, z) \mathbb{1}_{y+z \ge x} \\ \times \mathcal{K}(y, z) \frac{g(y, \tau) g(z, \tau)}{yz} dy dz dx' \\ + \int_{x_{\text{min}}}^{x} \int_{x_{\text{min}}}^{x_{\text{max}}} \mathbb{1}_{y+z \ge x} \mathcal{K}(z, y) g(z, \tau) \frac{g(y, \tau)}{y} dy dz, \quad (39)$$

which is a first proposition for the expression of the flux on the physical mass range, but not satisfactory because the flux is not conserved if mass flows out through $x_{max}$ due to mass transfer, i.e. $F_{\text{frag,nc}}[g](x_{max}, \tau) \neq 0$. For that reason, this flux is denoted $F_{\text{frag,nc}}$ for ‘non-conservative’ flux. Note that no mass flows out through $x_{min}$, $F_{\text{frag,c}}[g](x_{min}, \tau) = 0$, meaning that no grains of mass lower than $x_{min}$ are produced. Therefore, it is necessary to modify the flux in order to conserve the mass for which $F_{\text{frag,nc}}[g](x, \tau)|_{x=x_{min}-x_{max}} = 0$. The term $F_{\text{frag,c}}$ stands for the ‘conservative’ flux. To prevent the formation of particles of mass $x \ge x_{max}$ due to the phenomenon of mass transfer, it is sufficient that the total mass of the two colliding particles is lower than $x_{max} + x_{min}$. For instance, if $y + z = x_{max} + x_{min}$, the possible values are $y \le x_{max}$ and $z \le x_{max}$, with the limiting cases being $y = x_{max}$ and $z = x_{min}$, or vice versa. Moreover, the operator $\mathbb{1}_{y+z \ge x}$ has to be changed into $\mathbb{1}_{y+z \ge x+x_{min}}$ since the total mass of the colliding grains is always greater than twice the minimum grain mass. In summary, the condition $x + x_{min} \le y + z \le x_{max} + x_{min}$ must be verified when the range of grain mass is limited to $[x_{min}, x_{max}]$. By applying these modifications to equation (39), we obtain the conservative truncation of the flux

$$F_{\text{frag,c}}[g](x, \tau) = -\frac{1}{2} \int_{x_{\text{min}}}^{x} \int_{x_{\text{min}}}^{x_{\text{max}}} \int_{x_{\text{min}}}^{x_{\text{max}}} x' b(x'; y, z) \mathbb{1}_{y+z \ge x+x_{\text{min}}} \mathbb{1}_{x_{\text{max}}+x_{\text{min}} \ge y+z} \\ \times \mathcal{K}(y, z) \frac{g(y, \tau) g(z, \tau)}{yz} dy dz dx' \\ + \int_{x_{\text{min}}}^{x} \int_{x_{\text{min}}}^{x_{\text{max}}} \mathbb{1}_{y+z \ge x+x_{\text{min}}} \mathbb{1}_{x_{\text{max}}+x_{\text{min}} \ge y+z} \\ \times \mathcal{K}(z, y) g(z, \tau) \frac{g(y, \tau)}{y} dy dz. \quad (40)$$

For numerical purpose, it is necessary to extend the operators $\mathbb{1}$ in the limits of integrals. The first term on the right-hand side of equation (40) splits into two terms by comparing $z$ to $x$. Then, the operators $\mathbb{1}$ are applied on the variable $y$ for all terms to
give

$$
\begin{aligned}
F_{\text{frag,c}}[g](x, \tau) = & -\frac{1}{2} \int_{x_{\min}}^{x} \int_{x_{\min}}^{x} \int_{x-z+x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \mathcal{K}(y, z) \\
& \times \frac{g(y, \tau) g(z, \tau)}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \\
& -\frac{1}{2} \int_{x_{\min}}^{x} \int_{x}^{x_{\max}} \int_{x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \mathcal{K}(y, z) \\
& \times \frac{g(y, \tau) g(z, \tau)}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \\
& + \int_{x_{\min}}^{x} \int_{x-z+x_{\min}}^{x_{\max}-z+x_{\min}} \mathcal{K}(z, y) g(z, \tau) \frac{g(y, \tau)}{y} \mathrm{d}y \mathrm{d}z.
\end{aligned}
\quad (41)
$$

We directly observe that, this new expression of the flux ensures mass conservation of the system evolving in a finite mass range.

### 3.1.2 Evaluation of the flux

The general non-linear fragmentation equation (equation 9) belongs to the family of non-local partial differential equations. The evolution of the number density function $f$ depends on the evaluation of the product of the number density function over all the mass range, similar to the coagulation and non-linear fragmentation equations (Liu et al. 2019; Lombart & Laibe 2021; Lombart et al. 2022). In equation (41), the evaluation of the flux at the interface $x_{j-1/2}$ depends on the evaluation of $g_j$ in all bins, due to the double integral of the mass density function $g$. The approximation of $g$ is a non-continuous function due to the DG method. However, the flux $F_{\text{frag,c}}$, in equation (41), is a continuous function of mass across interfaces. This important characteristic of the flux differs from the usual DG solvers applied on local partial differential equations for which the flux is discontinuous and must be reconstructed at interfaces (Cockburn & Shu 1989; Zhang & Shu 2010; Guillet et al. 2019).

We assume that the distribution of fragments and the collision kernel are integrable, which is true for the analytical solution in Section 2.7. The numerical flux, equation (41), is analytically integrated over the mass variables by approximating $g$ with $g_j$ in bin $I_j$ (equation 37). The numerical flux at interface $x_{j-1/2}$ is

$$
\begin{aligned}
F_{\text{frag,c}}[g](x_{j-1/2}, \tau) = & -\frac{1}{2} \sum_{u=1}^{j-1} \sum_{l'=1}^{j-1} \int_{I_u} \int_{I_{l'}} \int_{x_{j-1/2}-z+x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \mathcal{K}(y, z) \\
& \times \frac{g(y, \tau) g_{l'}(z, \tau)}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \\
& -\frac{1}{2} \sum_{u=1}^{j-1} \sum_{l'=j}^{N} \int_{I_u} \int_{I_{l'}} \int_{x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \mathcal{K}(y, z) \\
& \times \frac{g(y, \tau) g_{l'}(z, \tau)}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \\
& + \sum_{l'=1}^{j-1} \int_{I_{l'}} \int_{x_{j-1/2}-z-z+x_{\min}}^{x_{\max}-z+x_{\min}} \mathcal{K}(z, y) g_{l'}(z, \tau) \frac{g(y, \tau)}{y} \mathrm{d}y \mathrm{d}z.
\end{aligned}
\quad (42)
$$

To be able to analytically calculate the integrals over $y$, it is necessary to approximate the function $g(y, \tau)$ over the entire mass range, with

the following approximation:

$$
\begin{aligned}
& \forall y \in [x_{\min}, x_{\max}], \\
& g(y, \tau) \approx \sum_{l=1}^{N} \sum_{i=0}^{k} g_l^i(\tau) \phi_i(\xi_l(y)) [H(y - x_{l-1/2}) - H(y - x_{l+1/2})],
\end{aligned}
\quad (43)
$$

where $H$ is the Heaviside function. Therefore, equation (42) writes

$$
\begin{aligned}
F_{\text{frag,c}}[g](x_{j-1/2}, \tau) = & -\frac{1}{2} \sum_{l'=1}^{j-1} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times T_{\text{frag},1}(x_{\max}, x_{\min}, j, l', l, i', i) \\
& -\frac{1}{2} \sum_{l'=j}^{N} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times T_{\text{frag},2}(x_{\max}, x_{\min}, j, l', l, i', i) \\
& + \sum_{l'=1}^{j-1} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) T_{\text{coag}}(x_{\max}, x_{\min}, j, l', l, i', i),
\end{aligned}
\quad (44)
$$

where the details of $T_{\text{frag},1}$, $T_{\text{frag},2}$, and $T_{\text{coag}}$ are given in Appendix A. These terms are calculated analytically with MATHEMATICA before being translated into FORTRAN and C++. The algorithm is written in FORTRAN/C++ and tested against the MATHEMATICA version for accuracy. $T_{\text{frag},1}$, $T_{\text{frag},2}$, and $T_{\text{coag}}$ are pre-computed once at the beginning of the algorithm, since they only depend on the chosen mass grid. This significantly improves the performance of the time solver, similar to Lombart & Laibe (2021) and Lombart et al. (2022). In practice, the three terms are stored in arrays with dimensions corresponding to the number of indices they contain. Then, the subarrays for index $j$ are multiplied by $g_{l'}^{i'}(\tau)g_{l}^{i}(\tau)$ and summed over all elements to obtain the three terms in the right-hand side of equation (44). The process is repeated for all $j$ to obtain $F_{\text{frag,c}}[g](x_{j-1/2}, \tau)$ for all $x_{j-1/2}$.

The evaluation of the flux assumes that the collision kernel is a 2D continuous function of mass, which is not always the case for physical problems. For instance, the kernel might depend on time through various physical quantities, such as gas temperature. The implementation of the physical collision kernel in Section 2.3 requires the use of the differential velocities between grains given by the 3D hydrodynamic code. The differential velocity term is a 2D piecewise constant function. One approach to couple the DG scheme with the hydrodynamic solver is to compute the integrals with the continuous cross-section, and then multiply by the 2D array for the differential velocity. This approximation of the physical kernel is given by

$$
\begin{aligned}
K(x, y) &= \pi \sigma(x, y) \Delta v(x, y) \approx K_{\text{approx}}(x, y), \\
K_{\text{approx}}(x, y) &\equiv \\
& \pi \sigma(x, y) \sum_{l'=1}^{N} \sum_{l=1}^{N} \Delta v_{l', l} \mathbf{1}_{x_{l'-1/2} < x < x_{l'+1/2}} \mathbf{1}_{x_{l-1/2} < y < x_{l+1/2}}.
\end{aligned}
\quad (45)
$$

In lieu of working with data obtained from 3D hydrodynamic simulations, we can test the errors from the approximation above by using the well-known differential velocity relation for Brownian motion (Dullemond & Dominik 2005; Brauer et al. 2008), which in
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: Surface plot of the Brownian kernel using the differential velocity relation in equation (46) (blue). The orange patched surface is the approximation of the Brownian kernel using equation (45). The red line indicates the location of the cross-sectional slice found in the right panel. Right: Cross-sectional slice of the Brownian kernel (blue line), with approximations from equation (45) (orange line) and a 2D cubic spline interpolation (green dashed line) overplotted. Vertical grey lines represent the boundaries of the bins.</figcaption></figure>

dimensionless form is written as

$$ \Delta v_{\text{Br}}(x, y) = \sqrt{\frac{1}{x} + \frac{1}{y}}. \qquad (46) $$

The Brownian collision kernel is $K_{\text{Br}}(x, y) = \sigma(x, y)\Delta v_{\text{Br}}(x, y)$. To quantify the errors, we define the continuous 2D absolute L¹ error

$$ e_{2\text{D,c}} = \frac{\sum_{j,i=1}^{N} \int_{I_j} \int_{I_i} |K_{\text{approx}}(x, y) - K_{\text{Br}}(x, y)| \, dx dy}{\sum_{j,i=1}^{N} \int_{I_j} \int_{I_i} K_{\text{Br}}(x, y) \, dx dy}, \qquad (47) $$

where the continuous integrals are evaluated with the Gauss-Legendre quadrature using 16 Gauss points. Fig. 3 shows the performance of the kernel approximation in equation (45). The error $e_{2\text{D,c}}$ generated by this approximation is ~ 3 per cent.

While this approximation provides a simple way of coupling to the differential velocities of hydrodynamic solvers, it is important to remember that, in practice, the differential velocities will be time dependent. Consequently, the physical kernel will also be time dependent and the pre-computed terms in equations (A1)-(A3) have to include the updated array $\Delta v_{i,j}$, e.g. by multiplication with the 2D array $\Delta v_{i,j}$:

$$ T_{\text{frag},1} \leftarrow \Delta v'_{i',l} \times T_{\text{frag},1}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i). \qquad (48) $$

If we want to enhance the accuracy of the 2D approximation to better preserve the high precision in the DG scheme, we can interpolate the physical collision kernel, as illustrated in the right panel of Fig. 3. The error of the 2D cubic spline interpolation (green dashed line) of the discrete kernel array,

$$ \forall (i, j) \in [1, N]^2, K_{i,j} = \pi \sigma(x_i, x_j) \Delta v_{\text{Br}}(x_i, x_j), \qquad (49) $$

is only ~ 0.005 per cent. The 2D interpolation gives this high accuracy only in log-log space where the grid is regular. The inconvenience of using the interpolation for the coupling is that the approximation of the kernel has to be integrated to evaluate the terms in equations (A1)-(A3) after each hydrodynamic time-step. This would significantly reduce the global performance of the DG

scheme, which relies heavily on pre-computed integrals or efficient quadrature methods (see Section 5.3.1) to quickly generate the terms in equations (A1)-(A3).

In addition to the differential velocities of grains given by the hydrodynamic solver, sub-grid models for the differential velocities, such as Brownian motion, can be added in the DG scheme like

$$ \Delta v(x_i, x_j) = \sqrt{\Delta v_{\text{hydro}}(x_i, x_j)^2 + \Delta v_{\text{subgrid}}(x_i, x_j)^2}. \qquad (50) $$

If the sub-grid model is considered independent of time, such as the Brownian motion in equation (46), the integral terms for the DG scheme need only to be pre-computed once.

### 3.1.3 Integral of the flux

Now that we have a relation for the conservative flux, $F_{\text{frag,c}}[g](x, \tau)$, we define $F_{\text{frag,c}}$ the term with the integral of the conservative flux, in equations (38) and (41), which takes the form

$$ \begin{aligned} F_{\text{frag,c}}(j, k', \tau) = & -\frac{1}{2} \int_{I_j} \int_{x_{\text{min}}}^{x} \int_{x_{\text{min}}}^{x} \int_{x-z+x_{\text{min}}}^{x_{\text{max}}-z+x_{\text{min}}} x' b(x'; y, z) K(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \frac{g(y, \tau) g(z, \tau)}{yz} dy dz dx' dx \\ & -\frac{1}{2} \int_{I_j} \int_{x_{\text{min}}}^{x} \int_{x}^{x_{\text{max}}} \int_{x-z+x_{\text{min}}}^{x_{\text{max}}-z+x_{\text{min}}} x' b(x'; y, z) K(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \frac{g(y, \tau) g(z, \tau)}{yz} dy dz dx' dx \\ & + \int_{I_j} \int_{x_{\text{min}}}^{x} \int_{x-z+x_{\text{min}}}^{x_{\text{max}}-z+x_{\text{min}}} K(z, y) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times g(z, \tau) \frac{g(y, \tau)}{y} dy dz dx, \end{aligned} \qquad (51) $$
where $k' \in [0, k]$. For each term in the right-hand side of equation (51), the integrals over $z$ are split into two integrals: $z \in [x_{\min}, x_{j-1/2}]$ and $z \in [x_{j-1/2}, x]$ for the first term; $z \in [x, x_{j+1/2}]$ and $z \in [x_{j+1/2}, x_{\max}]$ for the second term; and $z \in [x_{\min}, x_{j-1/2}]$ and $z \in [x_{j-1/2}, x]$ for the third term. We apply the same method as in Section 3.1.2 to obtain the numerical integral of the flux

$$ \mathcal{F}_{\text{frag,c}}(j, k', \tau) = -\frac{1}{2} \sum_{l'=1}^{j-1} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \times \mathcal{T}_{\text{frag},1,A}(x_{\max}, x_{\min}, j, k', l', l, i', i) \\ -\frac{1}{2} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{j}^{i'}(\tau) g_{l}^{i}(\tau) \times \mathcal{T}_{\text{frag},1,B}(x_{\max}, x_{\min}, j, k', l, i', i) \\ -\frac{1}{2} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{j}^{i'}(\tau) g_{l}^{i}(\tau) \times \mathcal{T}_{\text{frag},2,A}(x_{\max}, x_{\min}, j, k', l, i', i) \\ -\frac{1}{2} \sum_{l'=j+1}^{N} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \times \mathcal{T}_{\text{frag},2,B}(x_{\max}, x_{\min}, j, k', l', l, i', i) \\ + \sum_{l'=1}^{j-1} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \times \mathcal{T}_{\text{coag},A}(x_{\max}, x_{\min}, j, k', l', l, i', i) \\ + \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{j}^{i'}(\tau) g_{l}^{i}(\tau) \times \mathcal{T}_{\text{coag},B}(x_{\max}, x_{\min}, j, k', l, i', i), \quad (52) $$

where the definition of the different terms and the method to compute them are detailed in Appendix B. Importantly, since the integral terms in $F_{\text{frag,c}}$ and $F_{\text{frag,c}}$ are evaluated analytically, the accuracy depends only on the order of polynomials to approximate $g$.

## 3.2 Scaling limiter

The DG scheme needs the use of a scaling limiter to preserve the positivity of the numerical solutions (Zhang & Shu 2010; Liu et al. 2019; Lombart et al. 2022). The scaling limiter is applied by a reconstruction step based on cell averaging. The reconstructed polynomials in each bin writes (Liu et al. 2019; Lombart et al. 2022)

$$ \left\{ \begin{array}{l} \forall j \in [1, N], \quad p_j(x, \tau) \equiv \psi_j(\tau) (g_j(x, \tau) - \bar{g}_j(\tau)) + \bar{g}_j(\tau) \\ \psi_j(\tau) \equiv \min \left\{ 1, \left| \frac{\bar{g}_j(\tau)}{m_j(\tau) - \bar{g}_j(\tau)} \right| \right\}, \end{array} \right. \quad (53) $$

where $m_j(\tau) \equiv \min_{x \in I_j} g_j(x, \tau)$ and $\bar{g}_j$ is the average of $g_j$ in bin $I_j$

$$ \bar{g}_j(\tau) = \frac{1}{h_j} \int_{I_j} g_j(x, \tau) dx = g_j^0(\tau). \quad (54) $$

Since $g_j$ is approximated by orthogonal polynomials, the scaling limiter coefficient $\psi_j$ is applied to all components of $p_j$ except the average value which writes $p_j^0(\tau) = g_j^0(\tau)$. In Lombart & Laibe (2021) and Lombart et al. (2022), there is a mistake in the application of the scaling limiter. We present here the correct implementation of the scaling limiter for the DG scheme by following the description in Guillet et al. (2019). For the initialization of the DG scheme, the $L^2$

projection of the initial mass distribution onto the Legendre basis can result in negative values for the numerical mass density. At that step, the use of the scaling limiter is required; therefore, the following replacement is applied on the components of $g_j$

$$ \forall j \in [1, N], \quad \forall i \ge 1, \quad g_j^i(\tau) \leftarrow \psi_j g_j^i(\tau). \quad (55) $$

Then, after each evolution of the component of $g_j$, meaning after each time-step in the SSPRK method (see Section 3.3), the replacement in equation (55) has to be applied to ensure the positivity of the numerical solution.

## 3.3 CFL criterion

The DG scheme is associated with the Strong Stability Preserving Runge-Kutta third-order method (SSPRK) to ensure that the high-order accuracy is preserved during time (Liu et al. 2019; Lombart & Laibe 2021; Lombart et al. 2022). The SSPRK time solver is stable under a suitable Courant-Friedrichs-Lewy (CFL) condition on the time-step. A precise estimation of the CFL condition has been investigated in several studies on coagulation and fragmentation processes (Filbet & Laurençot 2004; Liu et al. 2019; Lombart & Laibe 2021; Laibe & Lombart 2022; Lombart et al. 2022), but none so far dedicated to the evaluation of the CFL condition for the general non-linear fragmentation (probably due to the high degree of complexity in the flux given in equation 41). The method based on the Laplace transform used in Laibe & Lombart (2022) seems difficult to be applied to the general non-linear fragmentation. The method to determine the CFL criterion developed in Filbet & Laurençot (2004); Lombart et al. (2022) does not provide an analytical expression to evaluate the CFL criterion because of the non-linearity of the general non-linear fragmentation. Therefore, we propose to numerically evaluate the CFL condition. The DG scheme with SSPRK method for order 0 corresponds to the forward Euler discretization, i.e.

$$ \begin{aligned} g_j^{0,n+1} &= g_j^{0,n} + \frac{\Delta \tau}{\Delta x_j} [F_{\text{frag}}[g_j^n](x_{j+1/2}, \tau) - F_{\text{frag}}[g_j^n](x_{j-1/2}, \tau)], \quad (56) \\ \Delta \tau_{\text{CFL}} &\le \frac{g_j^{0,n} \Delta x_j}{|F_{\text{frag}}[g_j^n](x_{j+1/2}, \tau) - F_{\text{frag}}[g_j^n](x_{j-1/2}, \tau)|} \end{aligned} $$

for the nth time-step. The CFL condition is evaluated to ensure the positivity of the bin average at time-step $n+1$, i.e. $\bar{g}_j^{n+1} = g_j^{0,n+1} \ge 0$ (Liu et al. 2019; Lombart & Laibe 2021; Lombart et al. 2022). We obtain the condition

$$ \begin{aligned} \forall j \in [1, N], \\ g_j^{0,n+1} \ge 0 \\ \Rightarrow \Delta \tau_{\text{CFL}} \le \frac{g_j^{0,n} \Delta x_j}{|F_{\text{frag}}[g_j^n](x_{j+1/2}, \tau) - F_{\text{frag}}[g_j^n](x_{j-1/2}, \tau)|} \end{aligned} \quad (57) $$

In practice, we define

$$ \Delta \tau_{\text{CFL}} \equiv \min_j \left[ \frac{g_j^{0,n} \Delta x_j}{|F_{\text{frag}}[g_j^n](x_{j+1/2}, \tau) - F_{\text{frag}}[g_j^n](x_{j-1/2}, \tau)|} \right], \quad (58) $$

and the time solver is executed with the time-step $d\tau = C_{\text{CFL}}\Delta\tau_{\text{CFL}}$, where $C_{\text{CFL}} \in [0, 1]$ is the CFL coefficient to ensure stability, typically $C_{\text{CFL}} = 0.3$ for SSPRK order 3 (Gottlieb 2015).

# 4 NUMERICAL RESULTS

The DG scheme presented in Section 3 is tested against the analytical solution outlined in Section 2.7. Tests are performed with a limited
number of mass bins, i.e. $N = 20$, in order to reflect the constraints from 3D hydrodynamical simulations. We also performed simulations for a power-law mass distribution of fragments.

## 4.1 Evaluation of errors

Error measurements are performed to determine the experimental order of convergence (EOC) and the efficiency of the DG algorithm, similar to Liu et al. (2019) and Lombart et al. (2022). The continuous $L^1$ norm used to evaluate the errors is

$$ 
\begin{aligned}
e_{c,N}(\tau) &\equiv \sum_{j=1}^{N} \int_{I_j} |g_j(x, \tau) - g_{\text{exact}}(x, \tau)| \, dx, \\
&\approx \sum_{j=1}^{N} \frac{h_j}{2} \sum_{\alpha=1}^{R} \omega_{\alpha} |g_j(x_j^{\alpha}, \tau) - g_{\text{exact}}(x_j^{\alpha}, \tau)|,
\end{aligned} 
\quad (59) $$

where the integral is approximated by the Gauss-Legendre quadrature method. The terms $g_{\text{exact}}$ and $g_j$ are, respectively, the exact and the numerical solutions. $\omega_\alpha$ and $x_j^\alpha$ are the weights and the points in bin $I_j$ for the Gauss-Legendre quadrature method. We use $R = 16$ Gaussian points. The discrete $L^1$ error is evaluated at the geometric mean value, which is the mass of the representative grain in each bin. The $L^1$ norm is evaluated in logarithmic scale with the following change of variable:

$$ 
\begin{aligned}
\|f\|_1 &\equiv \int_{x_{\min}}^{x_{\max}} |f(x)| dx = \sum_{j=1}^{N} \int_{I_j} |f(x)| dx \\
&\approx \sum_{v \leftarrow \log(x)} \sum_{j=1}^{N} \int_{\log(x_{j-1/2})}^{\log(x_{j+1/2})} e^v |f(e^v)| dv,
\end{aligned} 
\quad (60) $$

where $f$ is an arbitrary function. Then, we apply the mid-point rule to evaluate the integral to obtain the discrete $L^1$ error,

$$ e_{d,N}(\tau) \equiv \sum_{j=1}^{N} \log \left( \frac{x_{j+1/2}}{x_{j-1/2}} \right) \hat{x}_j |g_j(\hat{x}_j, \tau) - g(\hat{x}_j, \tau)|, \quad (61) $$

where $\hat{x}_j = \sqrt{x_{j+1/2}x_{j-1/2}}$ is the geometric mean in bin $j$. The EOC is defined as

$$ \text{EOC} \equiv \frac{\ln \left( \frac{e_N(\tau)}{e_{2N}(\tau)} \right)}{\ln(2)}, \quad (62) $$

where $e_N$ is the error for $N$ bins and $e_{2N}$ for $2N$ bins. The error can be the continuous or discrete $L^1$ error. To avoid time-stepping errors, the errors are calculated after one time-step. The numerical total mass density of the system is the first moment of $g(x, \tau)$ and writes

$$ M_{1,N}(\tau) \equiv \int_{x_{\min}}^{x_{\max}} g(x, \tau) dx = \sum_{j=1}^{N} h_j g_j^0(\tau). \quad (63) $$

Therefore, the mass conservation is analysed with the absolute error of the total mass density given by

$$ e_{M_{1,N}}(\tau) \equiv \frac{|M_{1,N}(\tau) - M_1|}{M_1}, \quad (64) $$

where $M_1$ is the first moment of the exact solution $g_{\text{exact}}$, which is constant in time.

## 4.2 Implementation details

Simulations are performed for a mass range $x \in [x_{\min}, x_{\max}]$, with $x_{\min} = 10^{-6}$ and $x_{\max} = 10^3$, in order to follow the formation of small grains from the fragmentation of larger grains. Numerical solutions of the mass density are shown only for polynomial approximations of order $k \in \{0, 1, 2, 3\}$. Tests are performed with FORTRAN and the errors are calculated with PYTHON. The initial components $g_j^l$ are evaluated by the $L^2$ projection of the initial condition $g_0(x)$ on the Legendre polynomials basis in each bin (see equation 17 in Lombart & Laibe 2021)

$$ 
\left\{ 
\begin{array}{l}
g_0(x) \equiv x e^{-x}, \\
\forall j \in [1, N], \quad g_j^l(0) = \frac{2}{h_j d_i} \int_{I_j} g_0(x) \phi_i(\xi_j(x)) dx.
\end{array}
\right. 
\quad (65) $$

where $d_i$ is the normalization coefficient of the Legendre polynomial basis defined as

$$ d_i \equiv \frac{2}{2i + 1}. \quad (66) $$

The integral is evaluated by a Gauss-Legendre quadrature method with five points. To avoid numerical instabilities, a minimum physical threshold is set to $10^{-20}$. Therefore, during the simulations, any polynomials $g_j(x, \tau)$ with mean value below the threshold are changed into a constant polynomial with value $10^{-20}$. This step is applied after each update of the $g_j^l(\tau)$, i.e. for each sub time-step in the SSPRK solver (section 3.3 and equation 38 in Lombart & Laibe 2021). We choose the CFL coefficient $C_{\text{CFL}} = 0.3$ and the algorithm is run sequentially on the Apple M1 Max chip. We use the GFORTRAN v13.1.0 compiler.

## 4.3 Distribution of fragments from Feingold et al. (1988)

Numerical solutions for the constant kernel with the specific distribution of fragments, equation (29), are benchmarked against the analytical solution presented in Section 2.7. Tests are performed in quadruple precision, in order to maintain stability of the DG scheme. Simulations are performed from $\tau = 0$ to $\tau = 3 \times 10^{-3}$ with 100 constant time-steps $\Delta\tau = 3 \times 10^{-5}$. The coefficient $\gamma$ in equation (29) is set to $10^4$, meaning that the majority of the grains produced by fragmentation are of mass $10^{-4}$. The evolution of the numerical solutions with $N = 20$ bins and polynomials of order $k \in \{0, 1, 2, 3\}$ are shown in Fig. 4. Initially, the majority of the mass is represented by large particles of mass $x = 1$. Then, the fragmentation process occurs and a large number of fragments with mass $x = 10^{-4}$ are produced. At $\tau = 8.10^{-5}$, the number of fragments with mass $x = 10^{-4}$ is so large that the majority of the total mass starts to be represented by these fragments. We observe at that time a double-humped curve. The number of the large grains of mass $x = 1$ continue to decrease. At the final time, the analytical and numerical solutions converge to the mass distribution of fragments.

### 4.3.1 Positivity and mass conservation

Fig. 4 shows the numerical results of the mass density (linear scale for the first fourth rows and log scale for the last row) versus the mass $x$ in log scale. For all polynomial orders, the numerical solutions remain positive thanks to the combination of the CFL-limited SSPRK time solver (Section 3.3) and the scaling limiter (Section 3.2). Fig. 5(a) shows the absolute error $e_{M_{1,N}}$ from $\tau = 10^{-5}$ to $\tau = 3 \times 10^{-3}$. The mass is conserved for all orders $k$, as expected from the design of the DG scheme with the conservative flux.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Test for the constant kernel and the distribution of fragments given in equation (29). The first four rows show the numerical solution (solid lines) for $N = 20$ bins and the listed value of $k \in \{0, 1, 2, 3\}$ at time $\tau = 0$ (left column), $\tau = 9 \times 10^{-5}$ (middle column), and $\tau = 3 \times 10^{-3}$ (right column). In these rows, the mass density is in linear scale and the mass in logarithmic scale. The exact solution $g(x, \tau)$ is given by the blue dashed line. Vertical grey lines represent the boundaries of the bins. The last row shows the same numerical solutions in log-log scale. The accuracy improves with increasing values of $k$. Order 3 recovers the two peaks at $\tau = 3 \times 10^{-3}$ with an accuracy of order $\sim 0.1$ per cent for masses below the cusp, and $\sim 0.01$ per cent above.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Test for the constant kernel and the distribution of fragments given in equation (29) with $N = 20$ bins. Panel (a): Evolution of the absolute error $e_{M1,N}$ of the total mass. The mass is conserved for each order $k$. Panel (b): Evolution of the time-step $d\tau$. The time-step steadily increases with time as more fragments are produced until reaching the global time-step value $\Delta\tau = 3 \times 10^{-5}$. The plateau shows that from $\tau = 10^{-4}$, the CFL condition is greater than $3 \times 10^{-5}$. Panels (c),(d): Time evolution of the $L^1$ continuous and discrete errors, showing they remain bounded at large times. The variation of errors observed at $\tau = 10^{-3}$ is explained in Section 4.3.3. Panels (e),(f): The errors $e_{c,N}$ and $e_{d,N}$ are plotted versus the number of bins per decade. The experimental order of convergence is EOC $= k + 1$ for $e_{c,N}$ and for $e_{d,N}$, EOC $= k + 1$ for polynomials of odd orders and EOC $= k + 2$ for polynomials of even orders. An accuracy of order $\sim 0.1$ per cent is achieved with $\sim 4$ bins/decade or $\sim 1$ per cent with $\sim 2$ bins/decade for $k = 3$.</figcaption></figure>
### 4.3.2 Accuracy

Fig. 4 shows the accuracy of the numerical solutions improves with the order of the polynomials. The two plots in log-log scale at $\tau = 8 \times 10^{-5}$ and $\tau = 3 \times 10^{-3}$ show that the numerical diffusion is reduced (up to a factor of 10) in the decreasing parts (exponential decays) as the order of polynomials increases. The major part of the total mass is generally localized in the maximum of the curve. However, at $\tau = 8 \times 10^{-5}$, the curve exhibits two maxima, indicating that a significant proportion of the total mass is represented by both small and large grains. Near the peak at lower mass, numerical solutions with order $k = 3$ achieve absolute errors of order $\sim 0.1$ per cent while errors of order $\sim 1$ per cent are obtained with $k \in \{0, 1, 2\}$. Near the peak at higher mass, numerical solutions with order $k = 3$ achieve absolute errors of order $\sim 0.01$ per cent while errors of order $\sim 1$ per cent is obtained with $k = 0$.

### 4.3.3 Stability of the DG scheme

Figs 5(c) and (d) show the time evolution of the continuous $e_{c,N}$ and discrete $e_{d,N}L^1$ errors. With only $N = 20$ bins, both errors for each order remain bounded over the entire time interval. However, notice an unusual behaviour of the curves around $\tau \sim 10^{-3}$. Looking at the error analysis one time-step at a time, the increase in the error comes from the fact that the numerical solution near the lower peak is approximately polynomial in shape for $\tau < 10^{-3}$ but transitions to an exponential decay for $\tau > 10^{-3}$. The sudden introduction of a second exponential section in the solution approximated by polynomials leads to a short temporal increase in the error that lasts until the initial size distribution is sufficiently depleted.

### 4.3.4 Convergence of the DG scheme

The EOC is determined independently from the mass range by plotting the continuous and discrete $L^1$ errors as a function of the number of bins per decade $N_{\text{bin/dec}}$. Figs 5(e) and (f) show the numerical errors at time $\tau = 10^{-9}$ for several total bin numbers $N = 5, 10, 20, 40, 80$. In Fig. 5(e), the EOC for the continuous $L^1$ error is of order $k + 1$. In Fig. 5(f), the EOC for the discrete $L^1$ error is of order $k + 2$ for odd polynomials, and $k + 1$ for even polynomials. With $e_{d,N}$, an accuracy of order $\sim 0.1$ per cent is achieved with $\sim 4$ bins/decade with $k = 3$, $\sim 5$ bins/decade with $k = 2$, and more than 10 bins/decade for $k = 0, 1$. An accuracy of order $\sim 1$ per cent is achieved with $\sim 2$ bins/decade with $k = 3$, $\sim 3$ bins/decade with $k = 2$, and $\sim 8$ bins/decade with $k \in \{0, 1\}$.

### 4.3.5 Analysis of the CFL condition

The evolution of the time-step $d\tau$ under the CFL condition equation (57) is shown in Fig. 5(b). The time-step increases during time until a stable value $\sim 3 \times 10^{-5}$. The green points around $\tau \approx 10^{-4}$ are smaller than the CFL condition to obtain equally spaced global time-steps $\Delta\tau$. The plateau observed from $\tau = 5 \times 10^{-4}$ shows that the CFL condition is greater than $3 \times 10^{-5}$; therefore, the time solver runs at $d\tau = \Delta\tau$.

## 4.4 Power-law distribution of fragments

### 4.4.1 Power-law mass distribution

In astrophysics, the mass distribution of fragments is traditionally described by a power law:

$$n(m, t)dm \propto m^{\alpha}dm, \qquad (67)$$

where the determination of $\alpha$ remains an open problem. Several theoretical works obtained $\alpha = -11/6 \approx -1.8$ (Dohnanyi 1969; Williams & Wetherill 1994; Jones et al. 1996; Tanaka et al. 1996). Experimental studies of grain-grain collision found that $\alpha$ takes value from $-2$ to $-1$ (Blum & Münch 1993; Güttler et al. 2010; Deckers & Teiser 2014; Bukhari Syed et al. 2017). We present a numerical analysis to constrain the value of $\alpha$. Let us consider a mass distribution of fragments defined as

$$b(x; y, z) \equiv A(y, z)x^{\alpha}. \qquad (68)$$

The normalization coefficient $A(y, z)$ is obtained by applying local mass conservation to our physical mass domain

$$\int_{x_{\min}}^{y+z} xb(x; y, z)dx = y + z, \qquad (69)$$

giving

$$A(y, z) = \frac{(\alpha + 2)(y + z)}{(y + z)^{\alpha + 2} - x_{\min}^{\alpha + 2}} = \frac{(\alpha + 2)(y + z)}{x_{\min}^{\alpha + 2} \left[ \left( \frac{y+z}{x_{\min}} \right)^{\alpha + 2} - 1 \right]}, \qquad (70)$$

where $\alpha \in \mathbb{R}$, $(y+z, x_{\min}) \in \mathbb{R}_+^2$ and $y+z > x_{\min}$. The mass distribution of fragments is a positive function if $\alpha \in (-\infty, -2) \cup (-2, \infty)$. A singularity appears for $\alpha = -2$ for the normalization in equation (70). For that specific case, the coefficient $A$ has to be calculated from equation (68). The number of fragments produced after each collision is given by equation (16) with $0 \leftrightarrow x_{\min}$,

$$N_{\text{frag}}(y, z) = \frac{(\alpha + 2)(y + z) \left[ \left( \frac{y+z}{x_{\min}} \right)^{\alpha+1} - 1 \right]}{(\alpha + 1)x_{\min} \left[ \left( \frac{y+z}{x_{\min}} \right)^{\alpha+2} - 1 \right]}. \qquad (71)$$

The number of fragments is strictly positive for $\alpha \in (-\infty, -2) \cup (-2, -1) \cup (-1, \infty)$. Two singularities appear for $\alpha = -2$ and $\alpha = -1$. For the value $\alpha = -1$, equations (70)–(71) have to be calculated from equation (68). The physical condition on $N_{\text{frag}}$ is that at least two fragments are produced per collision, $N_{\text{frag}} \ge 2$. This inequality is difficult to solve in the general case with the variables $y, z, x_{\min}$, and $\alpha$. But we can give an estimation of $\alpha$ for given values of $x_{\min}$ and the ratio $(y+z)/x_{\min}$, i.e. the ratio of mass between the total mass of the colliding grains and the mass of the smallest grain considered. Conversely, by choosing $\alpha$, we can estimate $N_{\text{frag}}$ according to the ratio $(y+z)/x_{\min}$. The left panel in Fig. 6 shows the range $\alpha$ values as a function of mass ratio for which $N_{\text{frag}} \ge 2$. The physical condition $N_{\text{frag}} \ge 2$ is satisfied over a wide range of values for the ratio $(y+z)/x_{\min}$, when $\alpha \in (-\infty, -2) \cup (-2, -1) \cup (-1, 0)$. The important result is that the model can handle the physical value $\alpha = -11/6$ with the condition that $(y+z)/x_{\min} \ge 3.64$. A large number of fragments require a large mass ratio (see the right panel in Fig. 6).

### 4.4.2 DG scheme with Gauss quadrature

Applying the method used in Section 3.1.2 to evaluate the flux with the power-law distribution of fragments from equation (68),
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left: The range in values for $\alpha$, respecting the condition $N_{frag} \ge 2$, is plotted versus mass ratio $(y+z)/x_{min}$. We observe that the values $\alpha = -2$ and $\alpha = -1$ (dashed black line) correspond to the singularities from equations (70) and (71). Shaded regions indicate values of $\alpha$ that respect the condition $N_{frag} \ge 2$ for a given mass ratio. For $\alpha = -11/6$, the condition is respected only if the mass ratio is greater than 3.64. Right: For a given number of fragments $N_{frag} \in \mathbb{N}$, the mass ratio is calculated with $\alpha = -11/6$. A large number of fragments requires correspondingly high mass ratios (e.g. more than a 1000 for $N_{frag} = 100$ when $\alpha = -11/6$).</figcaption></figure>

the analytic integrals result in Gaussian hypergeometric functions ${}_2F_1$; however, the correct evaluation of this special function requires high accuracy and is very computationally inefficient. We therefore opt to evaluate the integrals numerically with a Gauss quadrature method. Since the limits in the integral on the variable $y$ in equations (42) and (51) depend on the variable $z$, the accurate evaluation of the integrals required a large number of Gauss points. For that reason, we use the Gauss-Kronrod quadrature method. To ensure the stability of the DG scheme, we require 15 Gauss points is 15 or 31 Gauss-Kronrod points. We present here an approach that differs from Liu et al. (2019) to evaluate the integral on $y$, for which the limits contain variables, i.e.

$$ \int_{x_{j-1/2}-z+x_{min}}^{x_{max}-z+x_{min}} g(y, \tau)dy, \qquad (72) $$

where $g$ is approximated by a piecewise polynomial function. The method used in Liu et al. (2019) is to find which bins contain the values $x_{j-1/2} - z + x_{min}$ and $x_{max} - z + x_{min}$ and then split the integral over the bins to approximate $g$ by the corresponding polynomials. Our approach is to sum the integral of $g_l$ over each bin $l$ with an automatic selection of the bins which intersect the range $[x_{j-1/2} - z + x_{min}, x_{max} - z + x_{min}]$. Equation (72) is approximated by

$$ \begin{aligned} & \int_{x_{j-1/2}-z+x_{min}}^{x_{max}-z+x_{min}} g(y, \tau)dy \\ & \approx \sum_{l=1}^{N} \left( \begin{aligned} & \mathbb{1}_{x_{l-1/2} < x_{max} - z + x_{min}} \mathbb{1}_{x_{j-1/2} - z + x_{min} < x_{l+1/2}} \\ & \times \min(x_{max} - z + x_{min}, x_{l+1/2}) \\ & \times \max(x_{j-1/2} - z + x_{min}, x_{l-1/2}) \end{aligned} \right) g_l(y, \tau)dy. \end{aligned} \qquad (73) $$

The property $[x_{l-1/2}, x_{l+1/2}] \cap [x_{j-1/2} - u + x_{min}, x_{max} - u + x_{min}]$ is ensured by the use of the operator $\mathbb{1}$. The integral on $y$ is evaluated

with the Gauss-Kronrod quadrature with the following formula:

$$ \int_{a}^{b} f(x)dx \approx \frac{b-a}{2} \sum_{\alpha=1}^{Q} \omega_{\alpha} f \left( \frac{b+a}{2} + \frac{b-a}{2} s_{\alpha} \right), \qquad (74) $$

where $Q$ is the number of Gauss-Kronrod points and $\omega_\alpha$ and $s_\alpha$ are the weights and node coefficients, respectively, of the Gauss-Kronrod quadrature method. Then, the integral on $z$ is evaluated with the Gauss-Kronrod method to obtain the first double integral for the flux in equation (42). The details on the evaluation of the flux $F_{frag,c}[g](x_{j-1/2}, \tau)$ and the integral of the flux $F_{frag,c}(j, k', \tau)$ are given in Appendix C.

### 4.4.3 Results

We present here the benchmark results of the DG scheme with a power-law mass distribution of fragments. We test our numerical algorithm using a multiplicative kernel with $\alpha = -11/6$ and $N = 20$. The numerical solutions for $k = 0, 1, 2, 3$ are benchmarked with a reference solution obtained by using $N = 160$ bins with $k = 3$. The simulations are performed from $\tau = 0$ to $\tau = 1$ with 100 time-steps $\Delta\tau$ of length $10^{-2}$. The simulations are performed in double precision. Fig. 7 shows the numerical results of the mass density function versus the mass in linear-log for the first fourth rows and log-log scale for the last row. Analyse of the performance, shown in Fig. 8, shows similar results to those found in Section 4.3. High accuracy of the numerical solutions is achieved with high order of polynomials. The experimental order of convergence is not shown since it is similar to the test in Fig. 5.

# 5 DISCUSSION

## 5.1 Conservative form

We have derived the first conservative form of the general non-linear fragmentation equation (equation 9) as a hyperbolic law equation (equation 22) with a fragmentation flux in mass space. The fragmentation flux, equation (25), includes a ‘coagulation’ flux
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Test for the multiplicative kernel and the power-law distribution of fragments in equation (70) with power-law index $\alpha = -11/6$. The first four rows show the numerical solution (solid lines) for $N = 20$ bins and the indicated value $k \in 0, 1, 2, 3$ at time $\tau = 0$ (left column), $\tau = 1.8 \times 10^{-5}$ (middle column) and $\tau = 1$ (right column). The mass density is in linear scale and the mass in logarithmic scale. The reference numerical solution $g(x, \tau)$ is given by the blue dashed line. Vertical grey lines represent the boundaries of the bins. The last row shows the same numerical solutions in log-log scale. The accuracy improves with increasing values of $k$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Test for the multiplicative kernel and the power-law mass distribution of fragments in equation (70) with $N = 20$ bins. The power-law coefficient is $\alpha = -11/6$. The results are similar to the constant kernel test in Fig. 5. The mass is conserved for each order $k$ in panel (a). Panel (b): evolution of the time-step $d\tau$. For each order $k$, the time-step globally increases with time as more fragments are produced. The numerical errors decrease with the order of polynomials $k$ in panels (c) and (d).</figcaption></figure>

that results from potential mass transfer events that take place when one grain breaks against another. Our new expression for the flux combines in a continuous way the gain and loss terms of equation (9). The conservative form is of particular interest since it allows the use of robust numerical schemes (e.g. finite-volume or DG schemes) that naturally conserve the mass to machine precision.

## 5.2 Performance

The high-order DG scheme (see Section 3) efficiently solves the general non-linear fragmentation equation with only $N = 20$ bins. High-order polynomial approximations can drastically reduce the numerical diffusion observed for low-order numerical schemes with the same bin resolution. While this improvement was modest (one order of magnitude) for the analytic test case, the accuracy increased by four orders of magnitude for the power-law test. Moreover, the improved accuracy is occurring near the peak of the mass density curve where the majority of the mass of the system resides. In the test illustrated in Fig. 4, we already achieve an accuracy of approximately 1 per cent with order $k = 2$ and $N = 20$ bins, an accuracy level comparable to those reported by 3D hydrodynamic codes (Teyssier 2002; Price et al. 2018). Thus, for the first time, we could treat dust

fragmentation in 3D hydrodynamic simulations together with gas and dust dynamics. The details of how to couple the algorithm with a 3D hydrodynamic code will be presented in a future work.

## 5.3 Limitations

### 5.3.1 DG scheme architecture

The DG method, with the analytical evaluation of the integrals presented in Section 3, requires the use of quadruple precision to handle the arithmetic of large numbers for polynomials of high order. Thus, in its current form, order $k = 3$ is the maximum limit of the algorithm to obtain a good balance between accuracy and computational efficiency. This issue mainly comes from the approximation of $g(y, \tau)$ over the entire mass range with the aid of Heaviside functions (see equation 43). To evaluate the integrals in equations (A1)–(A3) and equations (B1)–(B6), the difference of Heaviside functions is propagated through the integrals. Therefore, the algorithm needs to evaluate the subtraction of two integrals with potentially large values, as illustrated in Fig. 9. This method of approximation over the entire mass range explains why the accurate evaluation of the terms in equations (A1)–(A3) and
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Illustration of the method used to analytically evaluate the integral of $g(y, \tau)$ approximated on the entire mass range by $g_l(y, \tau)$ in equation (43). The integral over the range $[x_{j-1/2} - z + x_{\min}, x_{\max} - z + x_{\min}]$ is evaluated by summing the integral of $g_l(y, \tau)$ (purple rectangle for a given $l$) over all bins $l \in [1, N]$. The integral of $g_l(y, \tau)$ is obtained by the difference of two integrals (blue and orange area). This difference of integrals is propagated all along the derivation to obtain terms in equations (A1)–(A3) and (B1)–(B6). Therefore, high accuracy, such as quadruple precision, is required to correctly evaluate the terms for the flux and the integral of the flux.</figcaption></figure>

(B1)–(B6) requires quadruple precision, which is computationally inefficient.

One solution to this problem would be to rewrite the approximation of $g$ over the entire mass range as

$$ \forall y \in [x_{\min}, x_{\max}], \\ g(y, \tau) \approx \sum_{l=1}^{N} \sum_{i=0}^{k} g_l^i(\tau) \phi_i(\xi_l(y)) H(y - x_{l-1/2}) H(x_{l+1/2} - y), \quad (75) $$

which instead uses the product of two Heaviside functions. In this way, we avoid the difference of two large numbers. The details of the derivation to obtain the analytical integrals will be presented in a future work. However, as mentioned in Section 4.4.2, analytic evaluation of the integrals for the power-law mass distribution of fragments will probably still result in Gaussian hypergeometric functions. We will explore whether there is an efficient way to rewrite or simplify these solutions using MATHEMATICA.

The DG scheme with the Gauss quadrature (DGGQ) method in Section 4.4.2 is an interesting alternative to the DG scheme with analytical integrals. The DGGQ scheme has several advantages: (i) any collision kernel and mass distribution of fragments can be used, (ii) only double precision is needed, and (iii) it is relatively easy to implement. However, in order to reach high accuracy for the power-law mass distribution test in Section 4.4, we found we needed to use the Gauss-Kronrod method with at least 31 points. This number of points is determined experimentally. For instance, with 21 points and order 3 polynomials, the simulation does not reach the final time because the time-step under the CFL condition at some point tends to zero. This behaviour can be explained by the fact that the evaluation of the integrals is not sufficiently accurate with 21 points. A good balance has been found by using 31 points. This requires a significative amount of time to pre-compute the integrals used in the evaluation of $F_{frag,c}$ and $\mathcal{F}_{frag,c}$ in equations (42) and (51) – particularly for orders $k=2$ and $k=3$. In practice, to reduce the pre-computing time, the terms are generated with multithreading processes by using OPENMP. The required large number of Gauss

<table><caption>**Table 1.** Averaged elapsed wall time in seconds for a single time-step for each order of polynomials $k$.</caption><thead><tr><th>Order<br>of polynomials</th><th>Analytic test<br>time (s)</th><th>Power-law test<br>time (s)</th></tr></thead><tbody><tr><td>$k=0$</td><td>$10^{-3}$</td><td>$10^{-4}$</td></tr><tr><td>$k=1$</td><td>$10^{-2}$</td><td>$8 \times 10^{-4}$</td></tr><tr><td>$k=2$</td><td>$2 \times 10^{-2}$</td><td>$2 \times 10^{-3}$</td></tr><tr><td>$k=3$</td><td>$6 \times 10^{-2}$</td><td>$8 \times 10^{-3}$</td></tr></tbody></table>

points suggests that the Gauss quadrature is not well adapted to evaluating integrals with a variable in the limits. It is possible that other integration methods would better suit this application. For example, the double exponential formula or Tanh-Sinh quadrature (Takahasi & Mori 1974; Mori & Sugihara 2001; Muhammad & Mori 2005) may be a suitable candidate for the DG fragmentation scheme we present in this paper. Tanh-Sinh quadrature is an efficient quadrature method to evaluate an integral on a finite interval with exponential convergence. With this quadrature method, we are likely to reach better accuracy with less points than the Gauss-Kronrod quadrature method. This too will be explored in a future work.

### 5.3.2 Time execution

Table 1 gives the averaged elapsed wall-clock time for a single time-step for each order of polynomials $k$. We observe that the power-law mass distribution test executes faster than the analytic test since less terms are required to evaluate $F_{frag,c}$ and $\mathcal{F}_{frag,c}$ for each time-step (see Appendix C). A gain of a factor of 10 is reached with the DGGQ scheme for $k=3$. To couple the DG scheme to multifluid 3D hydrodynamic codes, the DG time solver has to be highly efficient. A first test to treat the general non-linear fragmentation equation in a 3D simulation of protoplanetary disc will be performed with the code PHANTOM (Price et al. 2018). A typical gas/dust simulation of disc uses 10 dust sizes (i.e. $N=10$) with $10^6$ SPH particles and 32 CPUs. In this configuration, one hydrodynamic time-step takes ~1 s. By coupling the DG scheme to PHANTOM, we aim to reach a one-to-one ratio in running time, meaning that the elapsed time for the DG scheme has to be ~1s. According to the CFL criterion in Section 3.3, the fragmentation process will probably impose a sub-cycling compared to the hydrodynamic time-step. Therefore, for the value of $d\tau = 8 \times 10^{-3}$ for $k=3$, the running time for one hydrodynamic time-step would be at least $10^6 \times 8 \times 10^{-3} = 8 \times 10^3$ s. There are several strategies that can be explored in order to reach this one-to-one ratio.

The first strategy is to use the matrix form of the general non-linear fragmentation equation in order to take advantage of the vast computational resources that have been developed for matrix operations. A similar approach was used by Sandu (2006) for the coagulation equation. The details of the matrix form is given in Appendix D. Thus, efficient matrix operation algorithms, included in the library BLAS, can be used to improve the time solver. Moreover, the matrix form in equation (D1) allows for a simple form of the Jacobian that can be paired with an implicit time solver (Sandu 2006), eliminating the need for subcycling altogether. Unfortunately, the implicit SSPRK order three method cannot be made unconditionally stable and has to follow a time-step restriction with an SSP coefficient $C \le 2$ (Gottlieb, Ketcheson & Shu 2009; Ketcheson, Macdonald & Gottlieb 2009; Gottlieb, Ketcheson & Shu 2011). Therefore, considering the time it takes to compute the Jacobian, this implicit
method does not provide a significant improvement over the explicit method. However, recently, Gottlieb et al. (2022) derived a high-order unconditionally SSP implicit Runge-Kutta method by considering restrictions on the second derivative of the variable g in the spatial operator L for the DG scheme (see the definition of L in Lombart & Laibe 2021). The linearized backward Euler scheme used in Sandu (2006) and the implicit SSP two derivatives Runge-Kutta order 3 (Gottlieb et al. 2022) will be applied for the DG scheme and tested in a future work in order to remove the need for any sub-cycling during the evolution of a simulation.

The second strategy is to use GPU parallelization to run all $10^6$ calls to the fragmentation solver for all SPH particles in ~ 1s. A future work will be concerned with the design of the DG scheme adapted for GPU programming language (e.g. CUDA and SYCL). Since the general non-linear fragmentation equation is a non-local partial differential equation, the main difficulty will be to precisely manage the memory access on the GPU device in order to obtain the best performance. The implicit solver will be implemented with the GPU version of the BLAS library. The GPU version of the DG scheme will allow the treatment of fragmentation in future exascale hydrodynamics codes, such as IDEFIX (Lesur et al. 2023b) and DYABLO (Aubert & Durocher 2021).

## 5.4 Physical mass distribution of fragments

Laboratory experiments and theoretical studies in astrophysics indicate that fragmentation results in small grains with a power-law mass distribution, accompanied by one or two additional remnants with non-power-law distributed masses (Güttler et al. 2010; Windmark et al. 2012; Blum 2018; Hirashita & Aoyama 2019; Lebreuilly et al. 2022; Hasegawa et al. 2023). To take into account the production of these outlying remnants, the mass distribution of fragments can be modified (Hirashita et al. 2021) as

$$b(x; y, z) \equiv Ax^{\alpha} + \delta(x - (y + z - m_{\text{frag}})), \quad (76)$$

where $\delta$ is the Dirac delta function, $m_{\text{frag}}$ is the total mass of fragments, and $y+z-m_{\text{frag}}$ is the mass of the remnant grain (assuming there is only one). In general, $m_{\text{frag}}$ is defined as a percentage of $y+z$, the mass of the two colliding grains. The normalization constant A is determined by the following local mass conservation equation:

$$\int_{x_{\min}}^{x_{\max}} x A x^{\alpha} dx = m_{\text{frag}} \Rightarrow A = \frac{(2 + \alpha)m_{\text{frag}}}{x_{\max}^{2+\alpha} - x_{\min}^{2+\alpha}}. \quad (77)$$

By giving the link between $x_{\min}$, $x_{\max}$, and $m_{\text{frag}}$, the number of fragments per collision can be obtained. For instance, if we take the values from Hirashita & Aoyama (2019),

$$\begin{gathered} x_{\min} = 10^{-6}x_{\max}, \quad x_{\max} = 0.02m_{\text{frag}}, \quad \alpha = -\frac{11}{6} \\ \Rightarrow N_{\text{frag}} \approx 10^6. \end{gathered} \quad (78)$$

So under these conditions, each collision produces approximately one million fragments and one remnant. To obtain two remnants, like the N-body simulations of collisions in Hasegawa et al. (2023), we can further modify equation (76) with another Dirac delta function,

$$b(x; y, z) \equiv Ax^{\alpha} + \delta(x - m_1) + \delta(x - (y + z - (m_1 + m_{\text{frag}}))), \quad (79)$$

where $m_1$ is the mass of the first remnant and $m_1 + m_{\text{frag}}$ is the mass of the second remnant. The application of the DG scheme with

the mass distributions of fragments equations (76) and (79) will be investigated in a future work.

## 5.5 Coagulation and fragmentation

The coagulation/fragmentation equation used in astrophysics is the combination of the Smoluchowski coagulation equation and the general non-linear fragmentation equation equation (9) and writes (Blum 2006; Banasiak et al. 2019; Barik & Giri 2020; Giri & Laurençot 2021a)

$$\begin{aligned} \frac{\partial f(x, \tau)}{\partial t} = &+ \frac{1}{2} \int_{0}^{x} (1 - P_{\text{frag}}(x - y, y, \Delta v)) \mathcal{K}(x - y, y) f(x - y, \tau) f(y, \tau) dy \\ &+ \frac{1}{2} \int_{0}^{\infty} \int_{0}^{\infty} \mathbb{1}_{y+z \ge x} b(x; y, z) P_{\text{frag}}(y, z, \Delta v) \\ & \times \mathcal{K}(y, z) f(y, \tau) f(z, \tau) dy dz \\ &- f(x, \tau) \int_{0}^{\infty} \mathcal{K}(x, y) f(y, \tau) dy, \end{aligned} \quad (80)$$

where $P_{\text{frag}}$ is the probability that two particles fragment. The first term in the right-hand side of equation (80) describes the formation of particles of mass x due to coagulation. The second term denotes the formation of particles of mass x due to the fragmentation of two colliding particles of masses y and z. The third term describes the loss of particles of mass x due to collisions leading to coagulation or fragmentation. Some mathematical papers proved the existence and uniqueness of mass-conserving solutions to equation (80) for a large class of collision kernels $\mathcal{K}$ but only for fragmentation without mass transfer (Barik & Giri 2020; Giri & Laurençot 2021a). Moreover, Feingold et al. (1988) derived a steady-state solution to the coagulation/fragmentation equation with the alternative fragmentation rate equation in equation (21) using a constant collision kernel, constant probability of fragmentation and the mass distribution of fragments in equation (29). In a future work, the DG scheme will be applied to the conservative form of the coagulation/fragmentation equation and benchmarked with the analytical steady-state solution given in Feingold et al. (1988).

## 5.6 Dust aggregates

The general non-linear fragmentation equation presented in Section 2 is a Smoluchowski-like equation, meaning that it respects the same assumptions as the Smoluchowski coagulation equation. In particular, grains are considered to be spheres of the same material and, consequently, so are the fragments. However, astrophysical grains are aggregates with a non-spherical shape (Suttner & Yorke 2001; Blum 2006; Okuzumi, Tanaka & Sakagami 2009; Blum 2018). To account for this, the evolution of the number density of grains would depend on two variables describing the properties of a grain, such as the mass x and the porosity $\Psi$. The evolution of this number density $f(x, \Psi, \tau)$ is governed by the bivariate coagulation/fragmentation equation (Kostoglou & Konstandopoulos 2001; Ormel, Spaans & Tielens 2007)

$$\frac{\partial f(x, \Psi_x, \tau)}{\partial \tau} = \left. \frac{\partial f(x, \Psi_x, \tau)}{\partial \tau} \right|_{\text{coag}} + \left. \frac{\partial f(x, \Psi_x, \tau)}{\partial \tau} \right|_{\text{frag}}, \quad (81)$$
with

$$ \left. \frac{\partial f(x, \Psi, \tau)}{\partial \tau} \right|_{\text{coag}} \equiv \frac{1}{2} \int_0^x \int_{\Psi_{\text{min}}}^{\Psi_{\text{max}}} \int_{\Psi_{\text{min}}} (1 - P_{\text{frag}}(x-y, \Psi_{x-y}; y, \Psi_y; \Delta v)) \times \mathcal{K}(x-y, \Psi_{x-y}; y, \Psi_y) f(y, \Psi_y, \tau) f(x-y, \Psi_{x-y}, \tau) \times \delta(\Psi_x - \Gamma(x-y, \Psi_{x-y}; y, \Psi_y)) d\Psi_{x-y} d\Psi_y dy - f(x, \Psi_x, \tau) \int_0^\infty \int_{\Psi_{\text{min}}}^{\Psi_{\text{max}}} \mathcal{K}(x, \Psi_x; y, \Psi_y) f(y, \Psi_y, \tau) d\Psi_y dy, \quad (82) $$

and

$$ \left. \frac{\partial f(x, \Psi_x, \tau)}{\partial \tau} \right|_{\text{frag}} \equiv \frac{1}{2} \int_0^\infty \int_0^\infty \int_{\Psi_{\text{min}}}^{\Psi_{\text{max}}} \int_{\Psi_{\text{min}}} \mathbb{1}_{y+z \ge x} b(x, \Psi_x|y, \Psi_y; z, \psi_z) \times \mathcal{K}(y, \Psi_y; z, \psi_z) f(y, \Psi_y, \tau) f(z, \psi_z, \tau) \times P_{\text{frag}}(y, \Psi_y; z, \psi_z, \Delta v) d\Psi_y d\psi_z dy dz - f(x, \Psi_x, \tau) \int_0^\infty \int_{\Psi_{\text{min}}}^{\Psi_{\text{max}}} \mathcal{K}(x, \Psi_x; y, \Psi_y) f(y, \Psi_y, \tau) d\Psi_y dy, \quad (83) $$

where $\Gamma(x, \Psi_x; y, \Psi_y)$ gives the porosity of the resulting aggregate from collision of two grains of mass and porosity $(x, \Psi_x)$ and $(y, \Psi_y)$ (Ormel et al. 2007). The expression of $\Gamma$ depends on the collision algorithm, such as particle-cluster aggregation (PCA) and cluster-cluster aggregation (CCA). The function $b(x, \Psi_x|y, \Psi_y; z, \psi_z)$ gives the distribution in mass and porosity of the fragments from two colliding grains. The mass of fragments generally follows a power law, but the distribution in porosity is poorly understood. Some recipes have been used to determine the porosity of the fragment, for instance in Hirashita et al. (2021). It may be possible to derive a conservative form of equation (81) by following Qamar & Warnecke (2007) and Das et al. (2023) and then to apply the high-order DG scheme to efficiently solve this 2D coagulation/fragmentation equation with a few number of bins in mass and in porosity.

## 5.7 Accounting for shape, porosity, and chemistry

The general non-linear fragmentation model only considers spherical grains with the same chemical composition. However, recent observations of protoplanetary discs and diffuse ISM require more complex grain models (Ysard 2019; Hensley & Draine 2023; Siebenmorgen 2023; Ysard et al. 2024, and references therein), since shape and composition strongly affect the optical properties of the dust. To match the complexity of these observations, we need dust models that can account for different aggregate shapes and/or chemical compositions. For the first time, the DG scheme applied to the conservative form of the coagulation and the fragmentation equations provide a robust numerical framework to handle the complexity of grain models with greater dimensionality (e.g. shape and porosity). For instance, following the evolution of the mass density distribution of dust aggregates requires the 2D coagulation and fragmentation equations with variables mass and porosity (Kostoglou & Konstandopoulos 2001; Okuzumi et al.

2009; Hirashita et al. 2021). In addition, the chemical composition can be treated by considering that each grain is composed of a fraction of different chemical species (Pilinis 1990; Jacobson et al. 1994; Sandu 2006). Then, the dust coagulation and fragmentation equation can be written in term of the mass density distribution for each species to obtain a system of coagulation and fragmentation equations which can be solved efficiently by the DG scheme. By combining chemical composition and aggregate shape, the system of multidimensional equations can be efficiently solved by the DG scheme since the mass distribution functions are approximated by polynomials with a few number of bins in each dimension, reducing the computational cost. The efficiency will be improved with the optimizations of the DG scheme mentioned in Section 5. In future studies, the multidimensional dust evolution model, accounting for the correct optical properties of grains, will be implemented in 3D hydrodynamics simulations to match observations, thanks to this high-order DG scheme.

## 6 CONCLUSION

The grain-grain collision outcomes leading to the formation of fragments are important for understanding population levels of small grains population in different astrophysical environments. Several physical processes, such as thermal balance and gas-dust dynamics, are strongly impacted by the evolution of the dust-size distribution, highlighting the need to accurately treat dust coagulation and fragmentation in 3D simulations.

We have presented the derivation of the conservative form of the general non-linear fragmentation equation utilizing a mass flux (see equations 22 and 25). The physical interpretation of the fragmentation equation is enriched by the formulation of the fragmentation flux, which contains a coagulation flux resulting from mass-transfer in sufficiently high-velocity fragmenting collisions that allow some grains to grow in mass (Blum 2018; Birnstiel 2024; Hasegawa et al. 2023). This conservative form enables the application of robust numerical schemes, such as the finite volume method or the DG scheme we presented in this work. The high-order DG scheme accurately solves the general non-linear fragmentation equation on a low-resolution mass grid of only ~ 20 bins, thus paving the way to address poly-disperse dust fragmentation in 3D hydrodynamic codes. The DG scheme meets all necessary requirements to be coupled to 3D codes: (i) a strictly positive mass density distribution ensured by the SSPRK third-order time solver combined with a slope limiter, (ii) conservation of the mass at machine precision, (iii) accuracy of ~ 0.1 – 1 per cent obtained by high-order discretization in mass and time space, and (iv) a fast algorithm facilitated by pre-computed analytic integrals or efficient numerical integration.

Accurately following the evolution of the dust-size distribution in 3D simulations due to fragmentation has long been an objective in astrophysics (Haworth et al. 2016). The DG scheme, applied on the conservative form of the fragmentation equation, provides a numerical framework to treat dust coagulation/fragmentation in 3D simulations and allow for dust models that consider more than just grain size (e.g. porosity and chemical composition).

## ACKNOWLEDGEMENTS

We thank the referee for providing very useful comments that helped us to improve the manuscript. The authors thank P.-K. Barik for discussion and contribution on the derivation of the conservative form. M. L. acknowledges support the funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research
and innovation programme (Grant agreement No. 101098309 - PEB-BLES). M. L. acknowledges support the funding from the National Science and Technology Council (NSTC 112-2636-M-003-001). MAH acknowledges support from the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) under Germany Excellence Strategy - EXC-2094-390783311 and partial funding by the (0:funding-sourc e 3:href="http://dx.doi.org/10.13039/501100001659")Deutsche Forschungsgemeinschaft/(0:funding-source) (DFG; German Research Foundation) – 325594231. We used MATHEMATICA (Wolfram Research, Inc. 2024).

DATA AVAILABILITY

Data are available on request.

REFERENCES

Andrews S. M., 2020, ARA&A, 58, 483

Aubert D., Durocher A., 2021, in Siebert A. et al., eds, SF2A-2021: Proceedings of the Annual meeting of the French Society of A&A, p. 473

Banasiak J., Lamb W., Laurençot P., 2019, in Burns J. et al., eds, Analytic Methods for Coagulation-Fragmentation Models, Volume I, Chapman and Hall/CRC, New-York

Barik P. K., Giri A. K., 2020, Zeitschrift für angewandte Mathematik und Physik, 71, 38

Birnstiel T., 2024, Annual Review of Astronomy and Astrophysics, 0066

Blum J., 2006, Adv. Phys., 55, 881

Blum J., 2018, Space Sci. Rev., 214, 52

Blum J., Münch M., 1993, Icarus, 106, 151

Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859

Brown P. S., 1986, J. Appl. Meteorol. Climatol., 25, 313

Bukhari Syed M., Blum J., Wahlberg Jansson K., Johansen A., 2017, ApJ, 834, 145

Cazaux S., Tielens A. G. G. M., 2004, ApJ, 604, 222

Chen L.-H., Hirashita H., Hou K.-C., Aoyama S., Shimizu I., Nagamine K., 2018, MNRAS, 474, 1545

Cheng Z., Redner S., 1990, J. Phys. Math. Gen., 23, 1233

Cockburn B., Shu C.-W., 1989, Math. Comput., 52, 411

Das A., Paul J., Heinrich S., Kumar J., 2023, Proc. R. Soc. A: Math. Phys. Eng. Sci., 479, 20220658

Deckers J., Teiser J., 2014, ApJ, 796, 99

Dohnanyi J. S., 1969, J. Geophys. Res., 74, 2531

Dražkowska J., Li S., Birnstiel T., Stammler S. M., Li H., 2019, ApJ, 885, 91

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

Ernst M. H., Pagonabarraga I., 2007, J. Phys. Math. Gen., 40, F331

Feingold G., Tzivion (Tzitzvashvili) S., Leviv Z., 1988, J. Atmos. Sci., 45, 3387

Filbet F., Laurençot P., 2004, SIAM J. Sci. Comp., 25, 2004

Garaud P., Meru F., Galvagni M., Olczak C., 2013, ApJ, 764, 146

Gillespie J. R., List R., 1978, Pure Appl. Geophys., 117, 599

Giri A. K., Laurençot P., 2021a, SIAM J. Math. Anal., 53, 4605

Giri A. K., Laurençot P., 2021b, J. Differ. Equ., 280, 690

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

Gottlieb S., 2015, in Kirby R. M., Berzins M., Hesthaven J. S., eds, Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2014. Lecture Notes in Computational Science and Engineering. Springer International Publishing, Cham, p. 17

Gottlieb S., Grant Z. J., Hu J., Shu R., 2022, SIAM J. Num. Anal., 60, 423

Gottlieb S., Ketcheson D. I., Shu C.-W., 2009, J. Sci. Comput., 38, 251

Gottlieb S., Ketcheson D., Shu C.-W., 2011, Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations, World Scientific, Singapore

Grabowski W. W., 2022, in Morbidelli R., ed., Rainfall. Elsevier, Netherlands, p. 49

Guillet T., Pakmor R., Springel V., Chandrashekar P., Klingenberg C., 2019, MNRAS, 485, 4209

Guillot T., Ida S., Ormel C. W., 2014, A&A, 572, A72

Güttler C., Blum J., Zsom A., Ormel C. W., Dullemond C. P., 2010, A&A, 513, A56

Hasegawa Y., Suzuki T. K., Tanaka H., Kobayashi H., Wada K., 2023, ApJ, 944, 38

Haworth T. J. et al., 2016, Publ. Astron. Soc. Aust., 33, e053

Hensley B. S., Draine B. T., 2023, ApJ, 948, 55

Hirashita H., Aoyama S., 2019, MNRAS, 482, 2555

Hirashita H., Hoang T., 2020, MNRAS, 494, 1058

Hirashita H., Il'in V. B., Pagani L., Lefèvre C., 2021, MNRAS, 502, 15

Hirashita H., Yan H., 2009, MNRAS, 394, 1061

Hoang T., 2019, ApJ, 876, 13

Hu Z., Srivastava R. C., 1995, J. Atmos. Sci., 52, 1761

Jacobson M. Z., 2011, J. Atmos. Sci., 68, 334

Jacobson M. Z., Turco R. P., Jensen E. J., Toon O. B., 1994, Atmos. Environ., 28, 1327

Jones A. P., Tielens A. G. G. M., Hollenbach D. J., 1996, ApJ, 469, 740

Ketcheson D. I., Macdonald C. B., Gottlieb S., 2009, Appl. Num. Math., 59, 373

Khain A. P., Pinsky M., 2018, Physical Processes in Clouds and Cloud Modeling. Cambridge Univ. Press, Cambridge

Kobayashi H., Tanaka H., 2010, Icarus, 206, 735

Kobayashi H., Tanaka H., 2021, ApJ, 922, 16

Kostoglou M., Karabelas A. J., 2000, J. Phys. Math. Gen., 33, 1221

Kostoglou M., Konstandopoulos A. G., 2001, J. Aerosol Sci., 32, 1399

Laibe G., Lombart M., 2022, MNRAS, 510, 5220

Lebreuilly U., Commerçon B., Laibe G. 2020, A&A, 641, A112

Lebreuilly U., Vallucci-Goy V., Guillet V., Lombart M., Marchand P., 2022, MNRAS, 518, 3326

Lesur G. et al., 2023a, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. p. 465

Lesur G. R. J., Baghdadi S., Wafflard-Fernandez G., Mauxion J., Robert C. M. T., Van den Bossche M., 2023b, A&A, 677, A9

List R., Gillespie J. R., 1976, J. Atmos. Sci., 33, 2007

Liu H., Gröpler R., Warnecke G., 2019, SIAM J. Sci. Comput., 41, B448

Lombart M., Hutchison M., Lee Y.-N., 2022, MNRAS, 517, 2012

Lombart M., Laibe G., 2021, MNRAS, 501, 4298

Low T. B., List R., 1982, J. Atmos. Sci., 39, 1607

McFarquhar G. M., 2004, J. Atmos. Sci., 61, 777

Mkee C. F., Ostriker E. C., 2007, ARA&A, 45, 565

Mori M., Sugihara M., 2001, J. Comp. Appl. Math., 127, 287

Muhammad M., Mori M., 2005, Japan J. Ind. Appl. Math., 22, 77

Okuzumi S., Tanaka H., Sakagami M.-a., 2009, ApJ, 707, 1247

Ormel C. W., Paszun D., Dominik C., Tielens A. G. G. M., 2009, A&A, 502, 845

Ormel C. W., Spaans M., Tielens A. G. G. M., 2007, A&A, 461, 215

Paszun D., Dominik C., 2009, A&A, 507, 1023

Pilinis C., 1990, Atmos. Environ., 24, 1923

Pinsky M., Khain A., Shapiro M., 2001, J. Atmos. Sci., 58, 742

Prat O. P., Barros A. P., 2007, J. Appl. Meteorol. Climatol., 46, 1480

Price D. J. et al., 2018, PASA, 35, e031

Pruppacher H., Klett J., 2010, Mysak L. A., Hamilton K., Microstructure of Atmospheric Clouds and Precipitation, Springer, Netherlands

Qamar S., Warnecke G., 2007, Chem. Eng. Sci., 62, 679

Robinson A.,Booth R. A.,Owen J. E.,2024, Monthly Notices of the Royal Astronomical Society, 1524, 0035

Safronov V. S., 1972, Evolution of the protoplanetary cloud and formation of the earth and planets by Safronov V.S. Translated from Russian. Jerusalem (Israel): Israel Program for Scientific Translations. Keter Publishing House, p. 212, Washington, D.C.

Sandu A., 2006, Aerosol Sci. Technol., 40, 261

Siebenmorgen R., 2023, A&A, 670, A115

Squire J., Hopkins P. F., 2020, MNRAS, 498, 1239

Stammler S. M., Birnstiel T., 2022, ApJ, 935, 35

Suttner G., Yorke H. W., 2001, ApJ, 551, 461

Takahasi H., Mori M., 1974, Publ. Res. Inst. Math. Sci., 9, 721

Tanaka H., Inaba S., Nakazawa K., 1996, Icarus, 123, 450

Testi L. et al., 2014, Protostars Planets VI, 339,University of Arizona Press,Arizona (US)
Teyssier R., 2002, A&A, 385, 337

Tu Y., Li Z.-Y., Lam K. H., 2022, MNRAS, 515, 4780

Vericel A., Gonzalez J.-F., Price D. J., Laibe G., Pinte C., 2021, MNRAS, 507, 2318

Visser R. G., Ormel C. W., 2016, A&A, 586, A66

Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T., 2009, ApJ, 702, 1490

Wang P. K., 2013, Physics and Dynamics of Clouds and Precipitation, Cambridge Univ. Press, Cambridge

Weidenschilling S. J., 1977, MNRAS, 180, 57

Williams D. R., Wetherill G. W., 1994, Icarus, 107, 117

Windmark F., Birnstiel T., G\d{a} 1/4 tler C., Blum J., Dullemond C. P., Henning T., 2012, A&A, 540, A73

Wolfram Research, Inc., 2024, Mathematica, Version 14.1, Wolfram Research, Inc., Champaign, Illinois, available: https://www.wolfram.com/mathematica

Yamasawa D., Habe A., Kozasa T., Nozawa T., Hirashita H., Umeda H., Nomoto K., 2011, ApJ, 735, 44

Youdin A. N., Goodman J., 2005, ApJ, 620, 459

Ysard N., Jones A.-P., Guillet V., Demyk K., Decleir M., Verstaete L., Choubani I., M.-A. Miville-Deschênes, Fanciullo L., 2024, Astronomy & Astrophysics, A34, 0004

Ysard N., 2019, Proc. Int. Astron. Union, 15, 53

Zhang X., Shu C.-W., 2010, J. Comput. Phys., 229, 3091

## APPENDIX A: TERMS TO EVALUATE THE FLUX

The terms $T_{frag,1}$, $T_{frag,2}$, and $T_{coag}$ in equation (44) take the form

$$ T_{frag,1}(x_{\max}, x_{\min}, j, l', l, i', i) \equiv \sum_{u=1}^{j-1} \int_{I_u} \int_{I_{l'}} \int_{x_{j-1/2}-z+x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \times \mathcal{K}(y, z) \frac{\phi_i(\xi_l(y))\phi_{i'}(\xi_{l'}(z))}{yz} \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx'}, \quad (A1) $$

$$ T_{frag,2}(x_{\max}, x_{\min}, j, l', l, i', i) \equiv \sum_{u=1}^{j-1} \int_{I_u} \int_{I_{l'}} \int_{x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \times \mathcal{K}(y, z) \frac{\phi_i(\xi_l(y))\phi_{i'}(\xi_{l'}(z))}{yz} \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx'}, \quad (A2) $$

$$ T_{coag}(x_{\max}, x_{\min}, j, l', l, i', i) \equiv \int_{I_{l'}}^{x_{\max}-z+x_{\min}} \int_{x_{j-1/2}-z+x_{\min}}^{x_{\max}-z+x_{\min}} \mathcal{K}(z, y) \phi_{i'}(\xi_{l'}(z)) \frac{\phi_i(\xi_l(y))}{y} \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydz}. \quad (A3) $$

## APPENDIX B: TERMS TO EVALUATE THE INTEGRAL OF THE FLUX

The terms in equation (52) are defined as

$$ T_{frag,1,A}(x_{\max}, x_{\min}, j, k', l', l, i', i) \equiv \int_{I_j}^{x} \int_{x_{\min}}^{x} \int_{I_{l'}}^{x_{\max}-z+x_{\min}} \int_{x-z+x_{\min}}^{x_{\max}-z+x_{\min}} x' b(x'; y, z) \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) $$

$$ \times \frac{\phi_i(\xi_l(y))\phi_{i'}(\xi_{l'}(z))}{yz} \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx'}x, \quad (B1) $$

$$ \begin{aligned} T_{frag,1,B}(x_{\max}, x_{\min}, j, k', l, i', i) \equiv & \int_{I_j}^{x} \int_{x_{\min}}^{x} \int_{x_{j-1/2}-z+x_{\min}}^{x_{\max}-z+x_{\min}} \int_{x-z+x_{\min}}^{x} x' b(x'; y, z) \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \frac{\phi_i(\xi_l(y))\phi_{i'}(\xi_{l'}(z))}{yz} \\ & \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx'}x, \end{aligned} \quad (B2) $$

$$ \begin{aligned} T_{frag,2,A}(x_{\max}, x_{\min}, j, k', l, i', i) \equiv & \int_{I_j}^{x} \int_{x_{\min}}^{x_{j+1/2}-z+x_{\min}} \int_{x}^{x_{\max}-z+x_{\min}} \int_{x_{\min}}^{x} x' b(x'; y, z) \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \frac{\phi_i(\xi_l(y))\phi_{i'}(\xi_{l'}(z))}{yz} \\ & \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx'}x, \end{aligned} \quad (B3) $$

$$ T_{frag,2,B}(x_{\max}, x_{\min}, j, k', l', l, i', i) \equiv \int_{I_j}^{x} \int_{x_{\min}}^{x} \int_{I_{l'}}^{x_{\max}-z+x_{\min}} \int_{x_{\min}}^{x} x' b(x'; y, z) \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \frac{\phi_i(\xi_l(y))\phi_{i'}(\xi_{l'}(z))}{yz} \\ & \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx'}x, \quad (B4) $$

$$ \begin{aligned} T_{coag,A}(x_{\max}, x_{\min}, j, k', l', l, i', i) \equiv & \int_{I_j}^{x} \int_{x_{\min}}^{x_{\max}-z+x_{\min}} \int_{x-z+x_{\min}}^{x_{\max}-z+x_{\min}} \mathcal{K}(z, y) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \phi_{i'}(\xi_{l'}(z)) \frac{\phi_i(\xi_l(y))}{y} \\ & \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx}, \end{aligned} \quad (B5) $$

$$ \begin{aligned} T_{coag,B}(x_{\max}, x_{\min}, j, k', l, i', i) \equiv & \int_{I_j}^{x} \int_{x_{j-1/2}-z+x_{\min}}^{x_{\max}-z+x_{\min}} \int_{x-z+x_{\min}}^{x} \mathcal{K}(z, y) \partial_x \phi_{k'}(\xi_j(x)) \\ & \times \phi_{i'}(\xi_j(z)) \frac{\phi_i(\xi_l(y))}{y} \\ & \times [H(y - x_{l-1/2}) - H(y - x_{l+1/2})] \mathrm{dydzdx}. \end{aligned} \quad (B6) $$

All the terms in equations (B1) to (B6) are analytically calculated with MATHEMATICA and translated into FORTRAN/C++. Then, $\mathcal{F}_{frag,c}$ is evaluated similarly to the numerical flux in equation (44). $T_{frag,1,A}$, $T_{frag,1,B}$, $T_{frag,2,A}$, $T_{frag,2,B}$, $T_{coag,A}$ and $T_{coag,B}$ are computed once at the beginning of the algorithm. $\mathcal{F}_{frag,c}$ is obtained by the sum of six terms in equation (52). The first term is obtained by computing the product of the subarray for index $(j, k')T_{frag,1,A}(x_{\max}, x_{\min}, j, k', l', l, i', i)$ with $g_{l'}^i(\tau)g_l^i(\tau)$ and summing over all elements. The same evaluation is applied to the other terms. The process is repeated for all $(j, k')$ to obtain $\mathcal{F}_{frag,c}$ on all the mass range for the DG scheme.
APPENDIX C: DG SCHEME WITH THE
GAUSS-KRONROD QUADRATURE

Here, we describe the implementation of the Gauss–Kronrod quadra-
ture method to evaluate the integrals in the term $F_{\text{frag,c}}[g](x_{j-1/2}, \tau)$
in equation (41) and the term $\mathcal{F}_{\text{frag,c}}(j, k', \tau)$ in equation
(51).

C1 Gauss-Kronrod quadrature for the flux

By using the evaluation of the integral on y in equation (73), the
numerical flux in equation (41) writes

$$
\begin{equation}
\begin{aligned}
F_{\text{frag,c}}[g](x_{j-1/2}, \tau) = & -\frac{1}{2} \sum_{l'=1}^{j-1} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times T_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) \\
& -\frac{1}{2} \sum_{l'=j}^{N} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times T_{\text{GK,frag},2}(x_{\text{max}}, x_{\text{min}}, l', l, i', i) \\
& + \sum_{l'=1}^{j-1} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) T_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i),
\end{aligned}
\tag{C1}
\end{equation}
$$

where GK stands for Gauss-Kronrod and the terms write

$$
\begin{align}
T_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) &\equiv \nonumber \\
& \quad \sum_{u=1}^{j-1} \int_{I_u} \int_{I_{l'}} \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - z + x_{\text{min}}} \mathbb{1}_{x_{j-1/2} - z + x_{\text{min}} < x_{l+1/2}} \nonumber \\
& \quad \times \min(x_{\text{max}} - z + x_{\text{min}}, x_{l+1/2}) \nonumber \\
& \quad \times \int_{\max(x_{j-1/2} - z + x_{\text{min}}, x_{l-1/2})} x' b(x'; y, z) \mathcal{K}(y, z) \nonumber \\
& \quad \times \frac{\phi_i(\xi_l(y)) \phi_{i'}(\xi_l'(z))}{yz} \mathrm{dy}z \mathrm{d}x', \tag{C2}
\end{align}
$$

$$
\begin{align}
T_{\text{GK,frag},2}(x_{\text{max}}, x_{\text{min}}, l', l, i', i) &\equiv \nonumber \\
& \quad \sum_{u=1}^{j-1} \int_{I_u} \int_{I_{l'}} \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - z + x_{\text{min}}} \nonumber \\
& \quad \times \min(x_{\text{max}} - z + x_{\text{min}}, x_{l+1/2}) \nonumber \\
& \quad \times \int_{\max(x_{\text{min}}, x_{l-1/2})} x' b(x'; y, z) \mathcal{K}(y, z) \nonumber \\
& \quad \times \frac{\phi_i(\xi_l(y)) \phi_{i'}(\xi_l'(z))}{yz} \mathrm{dy}z \mathrm{d}x', \tag{C3}
\end{align}
$$

$$
\begin{align}
T_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) &\equiv \nonumber \\
& \quad \int_{I_{l'}} \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - z + x_{\text{min}}} \mathbb{1}_{x_{j-1/2} - z + x_{\text{min}} < x_{l+1/2}} \nonumber \\
& \quad \times \min(x_{\text{max}} - z + x_{\text{min}}, x_{l+1/2}) \nonumber \\
& \quad \times \int_{\max(x_{j-1/2} - z + x_{\text{min}}, x_{l-1/2})} \mathcal{K}(y, z) \nonumber \\
& \quad \times \phi_{i'}(\xi_l'(z)) \frac{\phi_i(\xi_l(y))}{y} \mathrm{dy}z. \tag{C4}
\end{align}
$$

The Gauss-Kronrod quadrature is then applied to evaluate the
integrals. The first term writes

$$
\begin{equation}
\begin{split}
T_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) = & \\
& \sum_{u=1}^{j-1} \sum_{\gamma=1}^{Q} \sum_{\alpha=1}^{Q} \sum_{\beta=1}^{Q} \frac{h_u h_{l'} h_{f_1, j, l', l, \alpha}}{8} \\
& \times \mathbb{1}_{x_{j-1/2} < x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}} \mathbb{1}_{x_{j-1/2} - \hat{x}_{l'}^{\alpha} + x_{\text{min}} < x_{l+1/2}} \\
& \times \omega_{\gamma} \omega_{\alpha} \omega_{\beta} \hat{x}_{u}^{\gamma} b(\hat{x}_{u}^{\gamma}; \hat{x}_{f_1, j, l', l}^{\alpha, \beta}, \hat{x}_{l'}^{\alpha}) \mathcal{K}(\hat{x}_{f_1, j, l', l}^{\alpha, \beta}, \hat{x}_{l'}^{\alpha}) \\
& \times \frac{\phi_i\left(\frac{2}{h_l}(\hat{x}_{f_1, j, l', l}^{\alpha, \beta} - x_l)\right) \phi_{i'}(s_{\alpha})}{\hat{x}_{f_1, j, l', l}^{\alpha, \beta} \hat{x}_{l'}^{\alpha}},
\end{split}
\tag{C5}
\end{equation}
$$

where $Q$ is the number of Gauss–Kronrod points, $s$ and $\omega$ are the
node and the weight coefficients, respectively, and we define the
terms $\hat{x}_{l'}^{\alpha}$, $h_{f_1, j, l', l, \alpha}$, and $\hat{x}_{f_1, j, l', l}^{\alpha, \beta}$ as

$$
\begin{equation}
\begin{split}
\hat{x}_{l'}^{\alpha} &\equiv x_{l'} + \frac{h_{l'}}{2}s_{\alpha}, \\
h_{f_1, j, l', l, \alpha} &\equiv \min(x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}, x_{l+1/2}) \\
&\quad - \max(x_{j-1/2} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}, x_{l-1/2}), \\
\hat{x}_{f_1, j, l', l}^{\alpha, \beta} &\equiv \frac{1}{2} \left[ \min(x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}, x_{l+1/2}) \right. \\
&\qquad \left. + \max(x_{j-1/2} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}) \right] + \frac{h_{f_1, j, l', l, \alpha}}{2}s_{\beta}.
\end{split}
\tag{C6}
\end{equation}
$$

The second term $T_{\text{GK,frag},2}$ writes

$$
\begin{equation}
\begin{split}
T_{\text{GK,frag},2}(x_{\text{max}}, x_{\text{min}}, l', l, i', i) = & \\
& \sum_{u=1}^{j-1} \sum_{\gamma=1}^{Q} \sum_{\alpha=1}^{Q} \sum_{\beta=1}^{Q} \frac{h_u h_{l'} h_{f_2, l', l, \alpha}}{8} \\
& \times \mathbb{1}_{x_{j-1/2} < x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}} \\
& \times \omega_{\gamma} \omega_{\alpha} \omega_{\beta} \hat{x}_{u}^{\gamma} b(\hat{x}_{u}^{\gamma}; \hat{x}_{f_2, l', l}^{\alpha, \beta}, \hat{x}_{l'}^{\alpha}) \mathcal{K}(\hat{x}_{f_2, l', l}^{\alpha, \beta}, \hat{x}_{l'}^{\alpha}) \\
& \times \frac{\phi_i\left(\frac{2}{h_l}(\hat{x}_{f_2, l', l}^{\alpha, \beta} - x_l)\right) \phi_{i'}(s_\alpha)}{\hat{x}_{f_2, l', l}^{\alpha, \beta} \hat{x}_{l'}^{\alpha}},
\end{split}
\tag{C7}
\end{equation}
$$

where

$$
\begin{align}
h_{f_2, l', l, \alpha} &\equiv \min (x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}, x_{l+1/2}) - x_{l-1/2}, \nonumber \\
\hat{x}_{f_2, l', l}^{\alpha, \beta} &\equiv \frac{1}{2} \left[ \min (x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}, x_{l+1/2}) + x_{l-1/2} \right] \nonumber \\
&\quad + \frac{h_{f_2, l', l, \alpha}}{2} s_{\beta}. \tag{C8}
\end{align}
$$

The last term $T_{\text{GK,coag}}$ writes

$$
\begin{equation}
\begin{split}
T_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) = & \\
& \sum_{\alpha=1}^{Q} \sum_{\beta=1}^{Q} \frac{h_{l'} h_{c, j, l', l, \alpha}}{4} \\
& \times \mathbb{1}_{x_{j-1/2} < x_{\text{max}} - \hat{x}_{l'}^{\alpha} + x_{\text{min}}} \mathbb{1}_{x_{j-1/2} - \hat{x}_{l'}^{\alpha} + x_{\text{min}} < x_{l+1/2}} \\
& \times \omega_{\alpha} \omega_{\beta} \mathcal{K}(\hat{x}_{c, j, l', l}^{\alpha, \beta}, \hat{x}_{l'}^{\alpha}) \\
& \times \frac{\phi_i\left(\frac{2}{h_l}(\hat{x}_{c, j, l', l} - x_l)\right) \phi_{i'}(s_\alpha)}{\hat{x}_{c, j, l', l}^{\alpha, \beta}},
\end{split}
\tag{C9}
\end{equation}
$$

where $h_{c,j,l',l,\alpha} = h_{f_1,j,l',l,\alpha}$ and $\hat{x}_{c,j,l',l}^{\alpha,\beta} = \hat{x}_{f_1,j,l',l}^{\alpha,\beta}$.
C2 Gauss–Kronrod quadrature for the term integral of the flux

By using equation (73), the term $\mathcal{F}_{\text{frag,c}}$ in equation (51) writes

$$ 
\begin{aligned}
\mathcal{F}_{\text{frag,c}}(j, k', \tau) = & -\frac{1}{2} \sum_{l'=1}^{j} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times \mathcal{T}_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) \\
& -\frac{1}{2} \sum_{l'=j}^{N} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times \mathcal{T}_{\text{GK,frag},2}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) \\
& + \sum_{l'=1}^{j} \sum_{l=1}^{N} \sum_{i'=0}^{k} \sum_{i=0}^{k} g_{l'}^{i'}(\tau) g_{l}^{i}(\tau) \\
& \times \mathcal{T}_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i)
\end{aligned}
\quad (\text{C10})
$$

where the terms write

$$ 
\begin{aligned}
\mathcal{T}_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) &\equiv \\
& \int_{I_j} \sum_{u=1}^{j} \mathbb{1}_{x_{u-1/2} < x} \quad \int_{\max(x_{\text{min}}, x_{u-1/2})}^{\min(x, x_{u+1/2})} \quad \mathbb{1}_{x_{l'-1/2} < x} \\
& \times \int_{\max(x_{\text{min}}, x_{l'-1/2})}^{\min(x, x_{l'+1/2})} \quad \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - z + x_{\text{min}}} \\
& \times \mathbb{1}_{x-z + x_{\text{min}} < x_{l+1/2}} \\
& \times \int_{\min(x_{\text{max}} - z + x_{\text{min}}, x_{l+1/2})}^{\min(x_{\text{max}} - z + x_{\text{min}}, x_{l'-1/2})} x' b(x'; y, z) \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\
& \times \frac{\phi_i(\xi_l(y)) \phi_{i'}(\xi_l'(z))}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \mathrm{d}x,
\end{aligned}
\quad (\text{C11})
$$

$$ 
\begin{aligned}
\mathcal{T}_{\text{GK,frag},2}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) &\equiv \\
& \int_{I_j} \sum_{u=1}^{j} \mathbb{1}_{x_{u-1/2} < x} \quad \int_{\max(x_{\text{min}}, x_{u-1/2})}^{\min(x, x_{u+1/2})} \quad \mathbb{1}_{x < x_{l'+1/2}} \\
& \times \int_{\max(x, x_{l'-1/2})}^{\min(x_{\text{max}}, x_{l'+1/2})} \quad \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - z + x_{\text{min}}} \\
& \times \int_{\max(x_{\text{min}}, x_{l-1/2})}^{\min(x_{\text{max}} - z + x_{\text{min}}, x_{l+1/2})} x' b(x'; y, z) \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\
& \times \frac{\phi_i(\xi_l(y)) \phi_{i'}(\xi_l'(z))}{yz} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \mathrm{d}x,
\end{aligned}
\quad (\text{C12})
$$

$$ 
\begin{aligned}
\mathcal{T}_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) &\equiv \\
& \int_{I_j} \mathbb{1}_{x_{l'-1/2} < x} \quad \int_{\max(x_{\text{min}}, x_{l'-1/2})}^{\min(x, x_{l'+1/2})} \quad \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - z + x_{\text{min}}} \\
& \times \mathbb{1}_{x-z + x_{\text{min}} < x_{l+1/2}} \\
& \times \min(x_{\text{max}} - z + x_{\text{min}}, x_{l+1/2}) \\
& \times \int_{\max(x-z+x_{\text{min}}, x_{l-1/2})}^{\min(x-z+x_{\text{min}}, x_{l'-1/2})} \mathcal{K}(y, z) \partial_x \phi_{k'}(\xi_j(x)) \\
& \times \frac{\phi_i(\xi_l(y)) \phi_{i'}(\xi_l'(z))}{y} \mathrm{d}y \mathrm{d}z \mathrm{d}x' \mathrm{d}x,
\end{aligned}
\quad (\text{C13})
$$

We then apply the Gauss–Kronrod quadrature method to evaluate the integrals. The first term writes

$$ 
\begin{aligned}
\mathcal{T}_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) = & \sum_{\lambda=1}^{Q} \sum_{u=1}^{j} \sum_{\gamma=1}^{Q} \sum_{\alpha=1}^{Q} \sum_{\beta=1}^{Q} \frac{h_j h_{u,\lambda} h_{l',\lambda} h_{f_1, j, l', l, \lambda, \alpha}}{16} \\
& \times \mathbb{1}_{x_{u-1/2} < \xi_j^{\lambda}} \mathbb{1}_{x_{l'-1/2} < \xi_j^{\lambda}} \\
& \times \mathbb{1}_{x_{l-1/2} < x_{\text{max}} - \xi_j^{\alpha} + x_{\text{min}}} \mathbb{1}_{\xi_j^{\lambda} + x_{\text{min}} < x_{l+1/2}} \\
& \times \omega_{\lambda} \omega_{\gamma} \omega_{\alpha} \omega_{\beta} \hat{x}_{u}^{\gamma} b(\hat{x}_{u}^{\gamma}, \hat{x}_{f_1, l', l}^{\lambda, \alpha, \beta}, \hat{x}_{l'}^{\alpha}) \mathcal{K}(\hat{x}_{f_1, l', l}^{\lambda, \alpha, \beta}, \hat{x}_{l'}^{\alpha}) \\
& \times \frac{\phi_i\left(\frac{2}{h_l}(\hat{x}_{f_1, l', l}^{\alpha, \beta} - x_l)\right) \phi_{i'}(s_{\alpha})}{\hat{x}_{f_1, l', l}^{\alpha, \beta} \hat{x}_{l'}^{\alpha}} \\
& \times \partial_x \phi_{k'}(\xi_j(x)) \Big|_{x=\xi_j^{\lambda}},
\end{aligned}
\quad (\text{C14})
$$

where we define the terms

$$ 
\begin{aligned}
\hat{x}_j^\lambda &\equiv x_j + \frac{h_j}{2} s_\lambda, \\
\hat{x}_u^\gamma &\equiv x_u + \frac{h_u}{2} s_\gamma, \\
\hat{x}_{l'}^\alpha &\equiv x_{l'} + \frac{h_{l'}}{2} s_\alpha, \\
h_{u,\lambda} &\equiv \min(\hat{x}_j^\lambda, x_{u+1/2}) - x_{u-1/2}, \\
h_{l',\lambda} &\equiv \min(\hat{x}_j^\lambda, x_{l'+1/2}) - x_{l'-1/2}, \\
h_{f_1, j, l', l, \lambda, \alpha} &\equiv \min(x_{\text{max}} - \hat{x}_{l'}^\alpha + x_{\text{min}}, x_{l+1/2}) \\
&\quad - \max(\hat{x}_j^\lambda - \hat{x}_{l'}^\alpha + x_{\text{min}}, x_{l-1/2}), \\
\hat{x}_{f_1, l', l}^{\lambda, \alpha, \beta} &\equiv \frac{1}{2} [\min(x_{\text{max}} - \hat{x}_{l'}^\alpha + x_{\text{min}}, x_{l+1/2}) \\
&\quad + \max(\hat{x}_j^\lambda - \hat{x}_{l'}^\alpha + x_{\text{min}}, x_{l-1/2})] \\
&\quad + \frac{h_{f_1, j, l', l, \lambda, \alpha}}{2} s_\beta.
\end{aligned}
\quad (\text{C15})
$$
The second term $\mathcal{T}_{\text{GK,frag},2}$ writes

$$ \begin{aligned} \mathcal{T}_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) = & \sum_{\lambda=1}^{Q} \sum_{u=1}^{j} \sum_{\gamma=1}^{Q} \sum_{\alpha=1}^{Q} \sum_{\beta=1}^{Q} \frac{h_j h_{u,\lambda} h_{j',\lambda} h_{f_2,j',l',\alpha}}{16} \\ & \times \mathbb{1}_{x_{u-1/2} < \hat{x}_j^\lambda} \mathbb{1}_{\hat{x}_j^\lambda < x_{j'+1/2}} \mathbb{1}_{x_{j'-1/2} < x_{\text{max}} - \hat{x}_j^\alpha} + x_{\text{min}} \\ & \times \omega_{\lambda} \omega_{\gamma} \omega_{\alpha} \omega_{\beta} \hat{x}_u^\gamma b \left( \hat{x}_u^\gamma; \hat{x}_{f_1,j',l',}^{\alpha,\beta}, \hat{x}_{j',l'}^\alpha \right) \mathcal{K} \left( \hat{x}_{f_2,j',l',}^{\alpha,\beta}, \hat{x}_{j',l'}^\alpha \right) \\ & \times \frac{\phi_i \left( \frac{2}{h_l} \left( \hat{x}_{f_2,l',l}^{\alpha,\beta} - x_l \right) \right) \phi_{j'}(s_\alpha)}{\hat{x}_{f_2,l',l}^{\alpha,\beta} \hat{x}_{j',l'}^\alpha} \\ & \times \partial_x \phi_{k'}(\xi_j(x)) \Big|_{x=\hat{x}_j^\lambda} \end{aligned} \quad (C16) $$

where we define the terms

$$ \begin{aligned} h_{l',\lambda} &\equiv \min(x_{\text{max}}, x_{l'+1/2}) - \min(\hat{x}_j^\lambda, x_{l'-1/2}), \\ h_{f_2,j',l',\alpha} &\equiv \min(x_{\text{max}} - \hat{x}_l^\alpha + x_{\text{min}}, x_{l+1/2}) \\ &\quad - \max(x_{\text{min}}, x_{l-1/2}), \\ \hat{x}_{f_1,j',l',}^{\alpha,\beta} &\equiv \frac{1}{2} \left[ \min(x_{\text{max}} - \hat{x}_l^\alpha + x_{\text{min}}, x_{l+1/2}) \right. \\ &\quad \left. + \max(x_{\text{min}}, x_{l-1/2}) \right] + \frac{h_{f_1,j',l',\lambda,\alpha}}{2} s_\beta. \end{aligned} \quad (C17) $$

The last term $\mathcal{T}_{\text{GK,coag}}$ writes

$$ \begin{aligned} \mathcal{T}_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) = & \sum_{\lambda=1}^{Q} \sum_{\alpha=1}^{Q} \sum_{\beta=1}^{Q} \frac{h_j h_{l',\lambda} h_{c,j',l',\lambda,\alpha}}{8} \\ & \times \mathbb{1}_{x_{l'-1/2} < \hat{x}_j^\lambda} \mathbb{1}_{x_{l'-1/2} < x_{\text{max}} - \hat{x}_j^\alpha} + x_{\text{min}} \mathbb{1}_{\hat{x}_j^\lambda - \hat{x}_j^\alpha + x_{\text{min}} < x_{l+1/2}} \\ & \times \omega_{\lambda} \omega_{\alpha} \omega_{\beta} \mathcal{K} \left( \hat{x}_{c,l',l}^{\lambda,\alpha,\beta}, \hat{x}_{j',l'}^{\alpha} \right) \\ & \times \frac{\phi_i \left( \frac{2}{h_l} \left( \hat{x}_{c,l',l}^{\lambda,\alpha,\beta} - x_l \right) \right) \phi_{j'}(s_\alpha)}{\hat{x}_{c,l',l}^{\lambda,\alpha,\beta}} \\ & \times \partial_x \phi_{k'}(\xi_j(x)) \Big|_{x=\hat{x}_j^\lambda}, \end{aligned} \quad (C18) $$

where $h_{c,j',l',\lambda,\alpha} = h_{f_1,j',l',\lambda,\alpha}$ and $\hat{x}_{c,l',l}^{\lambda,\alpha,\beta}$.

## APPENDIX D: MATRIX FORM

Originally used by Sandu (2006) for the coagulation equation, the matrix form of equation (38) writes

$$ \mathbf{A} \frac{\mathrm{d}g(\tau)}{\mathrm{d}\tau} = \begin{bmatrix} \mathbf{g}^T(\tau) \cdot \mathbf{B}^1 \cdot \mathbf{g}(\tau) \\ \vdots \\ \mathbf{g}^T(\tau) \cdot \mathbf{B}^{N \times (k+1)} \cdot \mathbf{g}(\tau) \end{bmatrix}, \quad (D1) $$

where $\mathbf{g} = [g_m]_{0 \le m \le N \times (k+1)}$ with the index $m = i + j \times (k+1)$ for which $i \in \llbracket 0, k \rrbracket$ and $j \in \llbracket 1, N \rrbracket$. The matrix **A** and the 3-tensor **B** write

$$ \begin{aligned} \mathbf{A} &= \operatorname{diag}[a_m]_{0 \le m \le N \times (k+1)}, \\ \mathbf{B} &= [b_{m,n}^s]_{0 \le m \le N \times (k+1), 0 \le n \le N \times (k+1), 0 \le s \le N \times (k+1)}, \end{aligned} \quad (D2) $$

where $\forall m \in \llbracket 1, N \times (k+1) \rrbracket$, $i = m/N$, $j = m \pmod{N}$, $a_m = \frac{2}{h_{c,j'}^{l'+1}}$ and $b_{m,n}^s$ gathers all the integrals for the flux and the integral of the flux terms. The 3-tensor **B** is a sparse tensor with elements

$$ \begin{aligned} \forall (m, n, s) \in \llbracket 0, N \times (k+1) \rrbracket^3, \\ b_{m,n}^s = \frac{1}{2} \left( \mathbb{1}_{l' \le j} \mathcal{T}_{\text{GK,frag},1}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) \right) \\ + \frac{1}{2} \left( \mathbb{1}_{j \le l'} \mathcal{T}_{\text{GK,frag},2}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) \right) \\ - \left( \mathbb{1}_{l' \le j} \mathcal{T}_{\text{GK,coag}}(x_{\text{max}}, x_{\text{min}}, j, k', l', l, i', i) \right. \\ \left. - \phi_k(1) \mathbb{1}_{j \le N-1} \right) \\ \times \left( \left[ \mathbb{1}_{l' \le j} T_{\text{frag},1}(x_{\text{max}}, x_{\text{min}}, j+1, l', l, i', i) \right. \right. \\ \left. \left. + \mathbb{1}_{j+1 \le l'} T_{\text{frag},2}(x_{\text{max}}, x_{\text{min}}, j+1, l', l, i', i) \right] \right. \\ \left. - \mathbb{1}_{l' \le j} T_{\text{coag}}(x_{\text{max}}, x_{\text{min}}, j+1, l', l, i', i) \right), \\ + \phi_k(-1) \\ \times \left( \left[ \mathbb{1}_{l' \le j-1} T_{\text{frag},1}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) \right. \right. \\ \left. \left. + \mathbb{1}_{j \le l'} T_{\text{frag},2}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) \right] \right. \\ \left. - \mathbb{1}_{l' \le j-1} T_{\text{coag}}(x_{\text{max}}, x_{\text{min}}, j, l', l, i', i) \right), \\ k' = s/N, \quad j = s \pmod{N}, \\ i' = m/N, \quad l' = m \pmod{N}, \\ i = n/N, \quad l = n \pmod{N}. \end{aligned} \quad (D3) $$

The presence of the operator $\mathbb{1}$ in $b_{m,n}^s$ implies that the components of the 3-tensor **B** are sparse matrices.

This paper has been typeset from a TeX/LA TeX file prepared by the author.