# Introducing cuDisc: a 2D code for protoplanetary disc structure and
evolution calculations

Alfie Robinson,1 Richard A. Booth2 and James E. Owen1
1Astrophysics Group, Imperial College London, Prince Consort Road, London SW7
2AZ, UK
2School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a new 2D axisymmetric code, cuDisc, for studying protoplanetary
discs, focusing on the self-consistent calculation of dust dynamics, grain
size distribution and disc temperature. Self-consistently studying these
physical processes is essential for many disc problems, such as structure
formation and dust removal, given that the processes heavily depend on one
another. To follow the evolution over substantial fractions of the disc
lifetime, cuDisc uses the CUDA language and libraries to speed up the code
through GPU acceleration. cuDisc employs a second-order finite-volume Godonuv
solver for dust dynamics, solves the Smoluchowski equation for dust growth and
calculates radiative transfer using a multi-frequency hybrid ray-tracing/flux-
limited-diffusion method. We benchmark our code against current state-of-the-
art codes. Through studying steady-state problems, we find that including 2D
structure reveals that when collisions are important, the dust vertical
structure appears to reach a diffusion-settling-coagulation equilibrium that
can differ substantially from standard models that ignore coagulation. For low
fragmentation velocities, we find an enhancement of intermediate-sized dust
grains at heights of $\sim 1$ gas scale height due to the variation in
collision rates with height, and for large fragmentation velocities, we find
an enhancement of small grains around the disc mid-plane due to collisional
“sweeping” of small grains by large grains. These results could be important
for the analysis of disc SEDs or scattered light images, given these
observables are sensitive to the vertical grain distribution.

###### keywords:

protoplanetary discs – methods: numerical – stars: pre-main sequence

††pubyear: 2015††pagerange: Introducing cuDisc: a 2D code for protoplanetary
disc structure and evolution calculations–A

##  1 Introduction

The past few decades have seen the study of young planetary systems and their
formation environments become a rapidly evolving field with major interest
within the astrophysics community. This is largely due to an unprecedented
wealth of observations made possible by recent observatories such as ALMA
(Wootten & Thompson, 2009) and the continual discovery of diverse exoplanetary
systems (e.g Mayor et al., 2011; Batalha et al., 2013; Winn & Fabrycky, 2015;
Madhusudhan, 2019; Zhu & Dong, 2021). Protoplanetary discs, the discs of gas
and dust that form around young stars, are the birthplace of such planetary
systems and have, therefore, been subject to extensive theoretical interest.
Various questions relating to their nature have arisen in recent years that
remain at least partly unanswered by current theoretical models. Examples of
such problems include: the nature of the mechanisms behind “sub-structure”
formation in protoplanetary discs, as observations have shown these objects to
exhibit diverse features from axisymmetric rings and gaps to non-axisymmetric
arcs (Andrews, 2020; Bae et al., 2022); the connection between the spatial
distribution and evolution of chemical species in discs to the eventual
compositions of planetary cores and atmospheres (Öberg et al., 2011; Booth et
al., 2017; Madhusudhan, 2019; Eistrup, 2023); and mechanisms for the dispersal
of protoplanetary discs after their observed lifetimes of $\sim$ a few Myr
(Ercolano & Pascucci, 2017; Owen & Kollmeier, 2019).

Exploring these problems requires sophisticated numerical modelling, given the
plethora of physical processes that govern the structure and evolution of
protoplanetary discs. Our understanding of discs has primarily been advanced
through the use of state-of-the-art codes for studying the dynamics and
thermodynamics of the gas and dust that comprise the disc material. 2D and 3D
simulations have typically been used to study discs on short, dynamical time-
scales, given their computational cost, whilst 1D models have often been used
to study discs on longer, secular time-scales. Work done using these models
has hugely advanced our understanding of protoplanetary discs; however, it has
become evident that for certain problems, the interplay of each of the facets
of disc physics - dynamics, thermodynamics and the dust size distribution -
must be studied self-consistently over secular time-scales. The 1D code DustPy
(Stammler & Birnstiel, 2022) is the current state-of-the-art for studying
problems of this nature; however, it cannot be used if the problem depends on
the intricacies of the disc vertical structure. Examples of such problems
include temperature instabilities (Watanabe & Lin, 2008; Wu & Lithwick, 2021;
Fuksman & Klahr, 2022) where 2D temperature solvers have been used but dust
dynamics and growth neglected, snow line instabilities (Owen, 2020) where 1D
temperature and dynamics solvers have been used, and problems relating to disc
dispersal such as the removal of dust from discs via radiation pressure (Owen
& Kollmeier, 2019; Krumholz et al., 2020), which has been studied with both 1D
and 2D solvers but without dust growth/fragmentation, or entrainment in
photoevaporative or magnetic winds (Franz et al., 2020; Booth & Clarke, 2021;
Hutchison & Clarke, 2021; Rodenkirch & Dullemond, 2022), again with 1D and 2D
approaches but with dust grain growth neglected.

This paper presents a new code, cuDisc, that includes 2D multi-fluid dust
dynamics coupled to both 2D temperature evolution and 1D secular gas
evolution. Dust and gas can be evolved, whilst simultaneously evolving the
dust grain size distribution and the resulting temperature structure. This
self-consistent calculation means fewer assumptions about the system’s state
for a particular scenario must be made. The code uses a 2D grid in the
poloidal plane, meaning that assumptions about vertical structure in the dust
do not have to be made as the structure can be calculated self-consistently.
In order to allow evolution calculations to be run for significant fractions
of the typical disc lifetime ($\sim$ Myr) in computationally feasible time-
scales, cuDisc utilises GPU acceleration via the employment of the CUDA
language and libraries. In the rest of this paper, Sections 2, 3, 4 and 5
outline the numerical methods used by cuDisc whilst Section 6 outlines the
code structure and Section 7 shows some example science calculations made
using cuDisc where we study grain growth in a steady-state transition disc.

##  2 Grid Structure

cuDisc is a 2D code in the poloidal plane, with axisymmetry adopted about the
disc’s rotation axis. The grid is structured in a fashion that makes certain
features of disc physics easier to calculate; cell interfaces are defined
along lines of constant polar angle above the mid-plane, $\theta$, and
constant cylindrical radius, $R$, as shown in Fig. 1. This means that the
basis vectors of the grid coordinate system are spherical radius,
$\boldsymbol{\hat{r}}$, and cylindrical height, $\boldsymbol{\hat{Z}}$,
although vector quantities (such as velocity) are dealt with in the standard
cylindrical components. This grid structure allows for ray-tracing from the
central star for calculating quantities such as optical depth whilst also
making calculations that require vertical integration easy, such as computing
hydrostatic equilibrium. Cell spacing can be arbitrary in principle, but the
standard implementation is logarithmic in $R$ with a power law in $\theta$.
The individual cell structure is shown in Fig. 2. Physical quantities are
stored at cell centres, and edge coordinates are required for the
reconstruction of the cell-centred quantities to the cell interfaces for the
advection routine (see later sections). Cell volumes, $V_{i,j}$, and interface
areas $A^{R}_{i,j}$ & $A^{Z}_{i,j}$ are given by

$V_{i,j}=\frac{1}{3}d(R^{3})_{i}d(\tan\theta)_{j},$ (1)
$A^{R}_{i,j}=(R^{e}_{i})^{2}d(\tan\theta)_{j},$ (2)
$A^{Z}_{i,j}=\frac{1}{2}\frac{d(R^{2})_{i}}{\cos\theta^{e}_{j}}.$ (3)

These geometric factors are written in such a way as to avoid divergence as
$R\rightarrow 0$. Explicitly, $d(R^{n})_{i}=(R^{e}_{i+1})^{n}-(R^{e}_{i})^{n}$
and $d(\tan\theta)_{j}=\tan\theta^{e}_{j+1}-\tan\theta^{e}_{j}$. The $Z$
coordinates are then given by $Z^{e}_{i,j}=R^{c}_{i}\tan\theta^{e}_{i,j}$ and
$Z^{c}_{i,j}=R^{c}_{i}\tan\theta^{c}_{i,j}$. Due to the non-orthogonality of
the grid coordinate bases, when calculating fluxes through interfaces, the
dot-product of velocity with the interface normals must be calculated to only
account for the velocity component orthogonal to the interface.

At each edge of the active cell domain, ghost cells are employed in order to
set the boundary conditions. The conditions along each of the four boundaries
can be set independently as either: outflow, where ghost cell quantities are
set by the values in the first adjacent active cell, zero, where ghost cell
quantities are set to floor values, or closed, where ghost cell quantities are
set by the values in the adjacent active cell but with zero flux over the
boundary. The default setup assumes vertical symmetry in the disc about the
mid-plane and, therefore, uses a minimum $\theta$ of 0 with a closed boundary
condition and outflow boundary conditions for the three other boundaries.

Figure 1: The computational grid used in cuDisc. Cell boundaries exist along
lines of constant radius and constant angle above the mid-plane, $\theta$.
Figure 2: The cell structure and indexing implemented in cuDisc. The indices
in the $R$ and $Z$ directions are $i$ and $j$, respectively.

##  3 Dynamics Solver

In cuDisc, the dust species are evolved by treating each grain size as a
pressureless fluid and solving their associated advection-diffusion equations.
We employ a second-order finite-volume Godonuv scheme for solving the set of
equations (Stone & Gardiner, 2009). We write our equations in terms of vector
fields of the conserved quantities, $\boldsymbol{Q}$, their associated fluxes,
$\boldsymbol{F}(\boldsymbol{Q})$, and any source terms, $\boldsymbol{S}$. For
our system, these fields are given by

$\boldsymbol{Q}_{i}=\begin{pmatrix}\rho_{i}\\\ \rho_{i}v_{R,i}\\\
\rho_{i}v_{\phi,i}R\\\ \rho_{i}v_{Z,i}\end{pmatrix},$ (4)

where $\rho_{i}$ and $\boldsymbol{v}_{i}=(v_{R,i},v_{\phi,i},v_{Z,i})$ are the
volume density and velocity of dust species $i$ respectively,

$\boldsymbol{F}_{i}=\begin{pmatrix}\rho_{i}\boldsymbol{v}_{i}+\boldsymbol{F}_{\text{diff},i}\\\
v_{R,i}(\rho_{i}\boldsymbol{v}_{i}+\boldsymbol{F}_{\text{diff},i})\\\
v_{\phi,i}R(\rho_{i}\boldsymbol{v}_{i}+\boldsymbol{F}_{\text{diff},i})\\\
v_{Z,i}(\rho_{i}\boldsymbol{v}_{i}+\boldsymbol{F}_{\text{diff},i})\end{pmatrix},$
(5)

where the two terms for each conserved quantity are the fluxes generated by
advection, $\rho_{i}\boldsymbol{v}_{i}$, and diffusion,
$\boldsymbol{F}_{\text{diff},i}$,

$\boldsymbol{S}_{i}=\begin{pmatrix}0\\\
\dfrac{\rho_{i}v_{\phi,i}^{2}}{R}-\rho_{i}\Omega^{2}R+f_{\text{drag},i,R}+f_{\text{ext},i,R}\vspace{5pt}\\\
f_{\text{drag},i,\phi}\vspace{5pt}\\\
-\rho_{i}\Omega^{2}Z+f_{\text{drag},i,Z}+f_{\text{ext},i,Z}\end{pmatrix},$ (6)

where $\Omega$ is the Keplerian angular velocity,
$\sqrt{GM_{*}/(R^{2}+Z^{2})^{3/2}}$, and $\boldsymbol{f}_{\text{drag}}$ and
$\boldsymbol{f}_{\text{ext}}$ are source terms that arise due to dust-gas drag
and any external forces (e.g. radiation pressure). We formulate diffusion in
the momentum equations as the diffusive flux acting to diffuse the advective
quantities. We do not currently consider the other terms discussed in recent
works, such as Huang & Bai (2022)black; these being the advection of the
diffusive flux and the time-dependent diffusive flux. Diffusion related terms
arise out of modelling turbulence by writing the density and velocity
components as a short-term average plus a short-term fluctuation, then
averaging the resulting equations over the short-term, keeping any terms that
include correlations between fluctuations (Reynolds averaging, see e.g. Cuzzi
et al. 1993). Note that we solve an azimuthal equation even though all $\phi$
derivatives are zero by axisymmetry, as the advection of angular momentum
throughout the $(R,Z)$ plane could be important in some problems; this is
sometimes referred to as 2.5D. Also, note that we write this equation in
angular momentum conserving form as this removes the Coriolis force term,
which has been shown to lead to loss of angular momentum conservation (see
e.g. Kley, 1998).

The quantities, $\boldsymbol{Q}$, are updated by solving the advection-
diffusion equations given by

$\frac{\partial\boldsymbol{Q}}{\partial
t}+\nabla\cdot\boldsymbol{F}=\boldsymbol{S}.$ (7)

We solve this set of equations in two stages via operator splitting: a
transport step where the advection-diffusion equations are solved as
homogeneous hyperbolic equations and a source step where the quantities are
updated through the source terms. For the transport step, the equations are
integrated over volume to find the flux-conservative form. Discretising this
gives the equation used to evolve a quantity $Q$ in cell $i,j$ to time $n+1$,

$Q^{n+1}_{i,j}=Q^{n}_{i,j}-\frac{\Delta
t}{V_{i,j}}(FA)^{n+\frac{1}{2}}_{i,j},$ (8)

where $(FA)_{i,j}$ is the total area-weighted-flux exiting through the
interfaces between cell $i,j$ and adjacent cells, given by

$(FA)_{i,j}=F^{R}_{i+1,j}A^{R}_{i+1,j}-F^{R}_{i,j}A^{R}_{i,j}\\\
+F^{Z}_{i,j+1}A^{Z}_{i,j+1}-F^{Z}_{i,j}A^{Z}_{i,j}.$ (9)

The time index of $n+\frac{1}{2}$ indicates that the fluxes should be
calculated at the half time-step to represent the average flux over the full
time-step $\Delta t$. Following Stone & Gardiner (2009), we calculate these
half-time fluxes in the following manner:

  1. 1.

The primitive quantities, $(\rho_{i},\boldsymbol{v}_{i})$, are reconstructed
at both sides of the cell interfaces from the conserved quantities at cell
centres using a first-order donor cell method and a Riemann solver is used to
calculate first-order advective fluxes through each interface at time $n$.

  2. 2.

The diffusive fluxes at each interface at time $n$ are calculated from the
cell-centred conserved quantities and added to the advective fluxes to get the
total area-weighted flux exiting each cell at time $n$, $(FA)^{n}_{i,j}$.

  3. 3.

The conserved quantities are advanced by a half time-step using these fluxes:

$Q^{n+\frac{1}{2}}_{i,j}=Q^{n}_{i,j}-\frac{1}{2}\frac{\Delta
t}{V_{i,j}}(FA)^{n}_{i,j},$ (10)

  4. 4.

These half-updated quantities are given a half time-step source update
(described in detail in subsection 3.1).

  5. 5.

Primitive quantities at time $n+\frac{1}{2}$ are now reconstructed from these
half-updated quantities using a second-order piecewise linear method with a
modified van Leer limiter detailed by Mignone (2014). These half-time
primitive quantities are used to calculate second-order advective fluxes
through each interface at time $n+\frac{1}{2}$.

  6. 6.

The diffusive fluxes at each interface at time $n+\frac{1}{2}$ are calculated
from the half-updated cell-centred quantities and added to the second-order
advective fluxes to get the total area-weighted flux exiting each cell at time
$n+\frac{1}{2}$, $(FA)^{n+\frac{1}{2}}_{i,j}$.

In steps (i) and (v), where the advective fluxes are calculated, the Riemann
problem has to be solved at the cell interfaces. For pressureless fluids,
collisions between fluids with different velocities lead to “delta shocks”,
which travel at the Roe-averaged velocity (Roe, 1981),

$\hat{v}=\frac{\sqrt{\rho_{l}}\boldsymbol{v}_{l}\cdot\boldsymbol{\hat{n}}+\sqrt{\rho_{r}}\boldsymbol{v}_{r}\cdot\boldsymbol{\hat{n}}}{\sqrt{\rho_{l}}+\sqrt{\rho_{r}}},$
(11)

where $r$ and $l$ denote the values of the reconstructed quantities to the
left and right of the interface in question. The reconstructed velocities are
dotted with the interface normal, $\boldsymbol{\hat{n}},$ to find the velocity
orthogonal to each interface; this is important given that we use a non-
orthogonal grid. Numerically, we capture this process according to Leveque
(2004) by taking the upwind flux depending on the sign of the Roe-averaged
velocity:

$F_{\text{adv}}=\begin{cases}Q_{l}\boldsymbol{v}_{l}\cdot\boldsymbol{\hat{n}}&\text{if
}\hat{v}>0\\\ Q_{r}\boldsymbol{v}_{r}\cdot\boldsymbol{\hat{n}}&\text{if
}\hat{v}<0\\\
\frac{1}{2}(Q_{l}\boldsymbol{v}_{l}\cdot\boldsymbol{\hat{n}}+Q_{r}\boldsymbol{v}_{r}\cdot\boldsymbol{\hat{n}})&\text{if
}\hat{v}=0\\\ 0&\text{if
}\boldsymbol{v}_{l}\cdot\boldsymbol{\hat{n}}<0<\boldsymbol{v}_{r}\cdot\boldsymbol{\hat{n}},\end{cases}$
(12)

The diffusive fluxes are then added to these advective fluxes to find the
total flux.

###  3.1 Dust source terms and diffusion

The source terms given by Eqn. 6 include curvature terms, gravitational terms,
drag terms and any other external forces. Drag is caused by dust-gas
interactions and is given by

$\boldsymbol{f}_{\text{drag}}=-\frac{\rho}{t_{s}}(\boldsymbol{v}-\boldsymbol{v}_{g}),$
(13)

where $\boldsymbol{v}_{g}$ is the gas velocity and $t_{s}$ is the stopping
time that conveys how well-coupled dust grains are to the gas flow; short
stopping times imply well-coupled dust grains that are entrained in the gas
flow. For the problems that we intend to study using cuDisc, the largest dust
grains have radii $\sim$ cm size, meaning that we are in the regime of Epstein
drag. This gives the stopping time as

$t_{s}=\frac{\rho_{m}s}{\rho_{g}v_{\text{th}}},$ (14)

where $\rho_{m}$ is the dust grain internal density, $s$ is the dust grain
radius, $\rho_{g}$ is the gas mass density and $v_{\text{th}}$ is the thermal
velocity of the gas, usually given by $(8/\pi)^{1/2}c_{s}$, where $c_{s}$ is
the isothermal gas sound speed.

Including source term updates, the full time-step is given by:

  1. 1.

$\boldsymbol{Q}^{*}_{i,j}=\boldsymbol{Q}^{n}_{i,j}-\dfrac{1}{2}\dfrac{\Delta
t}{V_{i,j}}\boldsymbol{(FA)}^{n}_{i,j}$

  2. 2.

$\boldsymbol{Q}^{**}_{i,j}=\boldsymbol{Q}^{*}_{i,j}+\dfrac{1}{2}\Delta
t\boldsymbol{S}_{\text{exp}}(\boldsymbol{P}^{n}_{i,j})$

  3. 3.

$\boldsymbol{Q}^{**}_{i,j}\rightarrow\boldsymbol{P}^{**}_{i,j}$

  4. 4.

$\boldsymbol{P}^{n+\frac{1}{2}}_{i,j}=\left(\rho^{**}_{i,j},\;\boldsymbol{v}^{**}_{i,j}-\dfrac{\Delta
t/2}{\Delta
t/2+t_{s}}(\boldsymbol{v}^{**}_{i,j}-\boldsymbol{v}_{g,i,j})\right)$

  5. 5.

$\boldsymbol{Q}^{{\dagger}}_{i,j}=\boldsymbol{Q}^{n}_{i,j}-\dfrac{\Delta
t}{V_{i,j}}\boldsymbol{(FA)}^{n+\frac{1}{2}}_{i,j}$

  6. 6.

$\boldsymbol{Q}^{{\dagger}{\dagger}}_{i,j}=\boldsymbol{Q}^{{\dagger}}_{i,j}+\Delta
t\boldsymbol{S}_{\text{exp}}(\boldsymbol{P}^{n+\frac{1}{2}}_{i,j})$

  7. 7.

$\boldsymbol{Q}^{{\dagger}{\dagger}}_{i,j}\rightarrow\boldsymbol{P}^{{\dagger}{\dagger}}_{i,j}$

  8. 8.

$\boldsymbol{P}^{n+1}_{i,j}=\left(\rho^{{\dagger}{\dagger}}_{i,j},\;\boldsymbol{v}^{{\dagger}{\dagger}}_{i,j}-\dfrac{\Delta
t}{\Delta
t+t_{s}}(\boldsymbol{v}^{{\dagger}{\dagger}}_{i,j}-\boldsymbol{v}_{g,i,j})\right)$

where $\boldsymbol{Q}$ & $\boldsymbol{P}$ represent conserved quantities and
primitive quantities respectively, $\boldsymbol{S}_{\text{exp}}$ is the vector
of source terms that are solved explicitly (curvature, gravity and external
forces) and quantities with an asterisk ($*$) or a dagger (${\dagger}$)
signify intermediate states. Steps (iii) and (vii) signify conversions of
conserved quantities to primitive quantities. Steps (iv) and (viii) are the
drag updates to the primitive velocities; these are solved implicitly, as if
they were solved explicitly, then short stopping times would prohibitively
limit the length of an explicit time-step. As the implicit update only
includes the drag terms, the Jacobian for the implicit system of equations is
linear and diagonal, meaning the solution has a simple, exact analytic form.
Boundary conditions are set before step (i) and again after step (iv).

The diffusive mass flux discussed in the previous section arises due to the
diffusion of dust grains via turbulent gas motions. This diffusive velocity is
calculated in a way analogous to molecular diffusion Clarke & Pringle (1988);
Takeuchi & Lin (2002), where the diffusive flux is given by

$\boldsymbol{F}_{\text{diff},i}=-\frac{\rho_{g}\nu}{\text{Sc}_{\textcolor{black}{i}}}\nabla\left(\frac{\rho_{i}}{\rho_{g}}\right),$
(15)

where $\nu$ is the kinematic viscosity of the disc, and Scblacki is the
Schmidt number, a dimensionless number that represents the strength of the
dust-gas coupling black for a given grain. blackIn the calculations presented
in this paper, Sci is set to unity for all dust species, however the code
allows any choice for the form of Sci, which may vary arbitrarily with both
position and dust properties. Given that we formulate the diffusion in the
momentum equations as the diffusive flux acting to diffuse the advective
quantities, the diffusive momentum flux at an interface is given by
$\boldsymbol{F}_{\text{diff},i}$ multiplied by the left or right advective
velocity depending on the sign of the diffusive flux itself,

$F_{\text{diff,mom}}=\begin{cases}F_{\text{diff}}\boldsymbol{v}_{l}\cdot\boldsymbol{\hat{n}}&\text{if
}F_{\text{diff}}>0\\\
F_{\text{diff}}\boldsymbol{v}_{r}\cdot\boldsymbol{\hat{n}}&\text{otherwise.}\\\
\end{cases}$ (16)

These diffusive fluxes are then added to the advective fluxes defined earlier
to give the total flux.

The time interval $\Delta t$ for each time-step is determined from the
Courant-Friedrichs-Lewy (CFL; Courant et al. 1928) condition:

$\Delta t_{\text{CFL}}=\text{min}\left(\Delta
t^{\text{adv}}_{\text{CFL}},\Delta t^{\text{diff}}_{\text{CFL}}\right),$ (17)

where $\Delta t^{\text{adv}}_{\text{CFL}}$ & $\Delta
t^{\text{diff}}_{\text{CFL}}$ are the CFL time-steps for advection and
diffusion respectively, given by

$\Delta
t^{\text{adv}}_{\text{CFL}}=C^{\text{adv}}\text{min}\left(|dR_{i}/v^{s}_{R,i,j}|,|dZ_{i,j}/v^{s}_{Z,i,j}|\right),$
(18)

where the minimum is over all cells $i,j$ and species $s$, and

$\Delta
t^{\text{diff}}_{\text{CFL}}=C^{\text{diff}}\text{min}\left(|(dR_{i})^{2}/(\nu_{i,j}/\text{Sc}_{s})|,|(dZ_{i,j})^{2}/(\nu_{i,j}/\text{Sc}_{s})|\right).$
(19)

Both time-steps have associated Courant numbers, $C^{\text{adv}}$ &
$C^{\text{diff}}$, which are safety factors set by default to 0.4 and 0.2,
respectively.

###  3.2 Diffusion

For both the dust dynamics and the radiative diffusion problem presented in
Section 4, we need to calculate diffusive fluxes. In general, a diffusive flux
can be written as

$\boldsymbol{F}=-D\nabla u,$ (20)

where $D$ is the diffusion constant and $u$ is the quantity being diffused.

Calculating the diffusive flux requires care in cuDisc given the non-
orthogonal coordinate system in use. This is because using the standard
difference formulae would give gradients of $u$ that are not perpendicular to
the interface in question and therefore using them directly would provide
incorrect fluxes across the interfaces. To deal with this we need a method
that calculates the 2D gradient of $u$ that we can dot-product with the
interface normal to find the flux across the interface. To this end, we use
the second-order accurate, conservative method detailed by Wu et al. (2012)
for diffusion problems on arbitrary mesh structures. Here, we present this
method specifically for our mesh.

To compute the gradients of $u$ within a cell, Wu et al. introduce
interpolation points on the cell interfaces. The location of these
interpolation points and accompanying values of $u$ are determined through the
conditions (i) $u$ must be continuous everywhere and (ii) the flux across the
interface must be continuous, i.e. $D\boldsymbol{\hat{n}}\cdot\nabla u$ must
be continuous at the interface, $\boldsymbol{\hat{n}}$ being the interface
normal vector. In this paper, we describe how these interpolation points are
found; for readers interested in why these criteria lead to the interpolation
points, we refer to the appendix of Wu et al. (2012). The flux at an interface
is built up from one-sided fluxes on each side of an interface. Take a radial
interface $\sigma$ between cells $k$ and $l$ (see Fig. 3). The net flux across
the face from the perspective of cell $k$ is given by

Figure 3: The cell stencil used for calculating diffusive fluxes over
interfaces. The cells are shown as orthogonal for illustrative ease. The total
flux $\tilde{F}_{\sigma}$ through the interface $\sigma$ that connects cells
$k$ and $l$ is calculated using the one-sided fluxes, $F_{k\sigma}$ and
$F_{l\sigma}$. The same can be said for the total flux $\tilde{F}_{\lambda}$
for interface $\lambda$ but with the one-sided fluxes $F_{a\lambda}$ and
$F_{l\lambda}$. The cell colours indicate which cells are used to construct
the one-sided fluxes: cells with yellow for $F_{k\sigma}$, cyan for
$F_{l\sigma}$, green for $F_{k\lambda}$, and red for $F_{l\lambda}$.

$\tilde{F}_{\sigma}=w_{k\sigma}F_{k\sigma}-w_{l\sigma}F_{l\sigma},$ (21)

where $w_{k\sigma}$ and $w_{l\sigma}$ are weights assigned to the one-sided
flux contributions from the cells on each side of the interface, $F_{k\sigma}$
and $F_{l\sigma}$, where these fluxes are directed out of cells $k/l$ across
the interface $\sigma$. The weights are defined as

$w_{k\sigma}=\frac{\mu_{k\sigma}}{\mu_{k\sigma}+\mu_{l\sigma}},\hskip
10.0ptw_{k\sigma}+w_{l\sigma}=1,$ (22)

where $\mu_{k\sigma}=\frac{d_{k\sigma}}{D_{k}}$, $d_{k\sigma}$ being the
shortest (perpendicular) distance between the cell centre and the plane of the
interface $\sigma$ and $D_{k}$ the diffusion constant associated with cell
$k$.

To include all flux orthogonal to the interface $\sigma$, the one-sided flux
$F_{k\sigma}$ is calculated by decomposing the interface normal vector
$\boldsymbol{n}_{\sigma}$ into two vectors, $\boldsymbol{v}_{kl}$ and
$\boldsymbol{v}_{km}$, that point from cell centre $k$ to the interpolation
points on the cell interface in question and the clockwise-adjacent interface
that joins cell $k$ to cell $m$, as shown in Fig. 4. To find the interpolation
point $\boldsymbol{v}_{kl}$ we first find two vectors that connect cell centre
$k$ to the points on the interface that are closest to cell centres $k$ and
$l$, shown as $\boldsymbol{dv}_{kl,1}$ and $\boldsymbol{dv}_{kl,2}$ in Fig. 4.
$\boldsymbol{v}_{kl}$ is then constructed using a weighted sum of
$\boldsymbol{dv}_{kl,1}$ and $\boldsymbol{dv}_{kl,2}$,

$\boldsymbol{v}_{kl}=w_{k\sigma}\boldsymbol{dv}_{kl,1}+w_{l\sigma}\boldsymbol{dv}_{kl,2}.$
(23)

The vector $\boldsymbol{v}_{km}$ is constructed similarly but with vectors and
weights associated with interface $\tau$ that connects cells $k$ and $m$. The
interface normal is then decomposed into these two vectors by solving:

$c_{kl}\boldsymbol{v}_{kl}+c_{km}\boldsymbol{v}_{km}=\boldsymbol{n}_{\sigma}$
(24)

for the coefficients $c_{kl}$ and $c_{km}$. These coefficients allow us to
calculate the gradients across the interface without explicitly calculating
the values of $u$ at the interpolation points. The one-sided flux $F_{k}$ is
then calculated via the weighted sum of these gradients:

$F_{k\sigma}=-D_{k}\left[w_{l\sigma}c_{kl}(u_{l}-u_{k})+w_{m\tau}c_{km}(u_{m}-u_{k})\right].$
(25)

To calculate the net flux $\tilde{F}_{\sigma}$ from Eqn. 21, $F_{l\sigma}$ is
constructed similarly to $F_{k\sigma}$ but from the perspective of cell $l$.
For polar/height interfaces such as interface $\lambda$ in Fig. 3, the process
is the same except that the stencil for constructing one-sided fluxes uses the
interface in question and the anti-clockwise adjacent interface as opposed to
clockwise.

Figure 4: The vector decomposition used for the cell interface normals when
calculating diffusive fluxes. Green vectors show those that point from cell
centre $k$ to the points on the interface $\sigma$ that are closest to cell
centres $k$ and $l$. The red vector $\boldsymbol{v}_{kl}$ is then constructed
via a weighted sum of the two green vectors. $\boldsymbol{v}_{km}$ is
constructed from vectors for the interface between cells $k$ and $m$ but are
not shown for simplicity.

###  3.3 Gas

In cuDisc, the 1D viscous evolution equation for gas surface density is
solved, and the 2D gas distribution is calculated through hydrostatic
equilibrium. Currently, there is no feedback from the dust to the gas;
however, it may be added in the future. The viscous evolution equation is
given by (Pringle, 1981)

$\frac{\partial\Sigma_{g}}{\partial t}=\frac{3}{R}\frac{\partial}{\partial
R}\left(R^{1/2}\frac{\partial}{\partial
R}(R^{1/2}\nu\Sigma_{g})\right)+S_{g},$ (26)

where $S_{g}$ represents any source terms (e.g. mass loss via photoevaporative
winds). The equation is solved in two stages via operator splitting; in the
first stage, the first term on the right-hand side is solved as a diffusion
problem using an explicit forward-time-centered-spatial scheme, and in the
second stage, the source terms are solved using explicit Eulerian updates.

The equations we use to calculate the gas velocities are found from viscous
theory (see e.g Urpin, 1984; Balbus & Papaloizou, 1999), written in the
coordinate grid basis of cuDisc,
$(\boldsymbol{\hat{r}},\boldsymbol{\hat{Z}})$. These velocities are needed
explicitly for their influence on dust dynamics - they are not used for
evolving the gas. Under the assumption that $v_{Z}\ll v_{R},v_{\phi}$
(justified for thin discs as $v_{Z}\sim(H/R)v_{R}$ where $H$ is the gas scale
height, see Urpin 1984), we find

$v_{\phi}=\bigg{[}\bigg{.}\frac{GM_{*}\cos^{2}{\theta}}{r}+\frac{r}{\rho_{g}}\left(\frac{\partial
P}{\partial r}-\sin{\theta}\frac{\partial P}{\partial
Z}\right)\bigg{.}\bigg{]}^{1/2},$ (27)

$v_{R}=r\cos{\theta}\left[\frac{\partial(rv_{\phi})}{\partial
r}-r\sin{\theta}\frac{\partial v_{\phi}}{\partial Z}\right]^{-1}\times\\\
\left\\{\frac{1}{\rho_{g}}[\nabla\cdot\boldsymbol{T}]_{\phi}-v_{Z,\text{cyl}}\frac{\partial
v_{\phi}}{\partial Z}\right\\},$ (28)

where the derivatives are performed along the grid directions of spherical $r$
and cylindrical $Z$. $v_{Z,\text{cyl}}$ is the physical velocity in the
$Z$-direction. blackWe do not include a form for the vertical velocity due to
viscosity here as it depends on assumptions made about the turbulence model
(see e.g. Philippov & Rafikov, 2017); it can also be governed by other
processes such as disc winds. For these reasons, by default (and for the tests
in this paper) the vertical gas velocity is assumed to be 0. Different models
of vertical gas motions can be included depending on the problem in question.
$\theta$ is the angle above the disc mid-plane, $P$ is the gas pressure, given
by the ideal gas law, and $\boldsymbol{T}$ is the viscous stress tensor. To
calculate $[\nabla\cdot\boldsymbol{T}]_{\phi}$, the $r\phi$ and $Z\phi$
components of $\boldsymbol{T}$ are required. Assuming a Navier-Stokes-like
viscosity, these have the forms:

$T^{r\phi}=\frac{\rho_{g}\nu}{\cos^{2}{\theta}}\left[r\frac{\partial}{\partial
r}\left(\frac{v_{\phi}}{r}\right)-\sin{\theta}\frac{\partial
v_{\phi}}{\partial Z}\right],$ (29)
$T^{Z\phi}=\frac{\rho_{g}\nu}{\cos^{2}{\theta}}\left[\frac{\partial
v_{\phi}}{\partial Z}-r\sin{\theta}\frac{\partial}{\partial
r}\left(\frac{v_{\phi}}{r}\right)\right],$ (30)

and the tensor divergence in the $\phi$ direction is explicitly written as

$[\nabla\cdot\boldsymbol{T}]_{\phi}=\frac{1}{r^{3}}\frac{\partial(r^{3}T^{r\phi})}{\partial
r}+\frac{\partial T^{Z\phi}}{\partial Z}.$ (31)

These derivatives are calculated using finite differencing on the grid
variables to calculate cell-centred values for the gas velocities.

###  3.4 Kinematic viscosity

The kinematic viscosity $\nu$ is important for governing the gas’ evolution
and dust diffusion. It is often parameterised following Shakura & Sunyaev
(1973) as $\alpha c_{s}H$ where $H$ is the gas scale height, and $\alpha$ is a
constant that is used to control the strength of the effective viscosity
generated by turbulent processes. With a constant $\alpha$, this assumes that
the viscosity varies with temperature changes in the disc. However, $\alpha$
need not be assumed constant; it can depend on other quantities such as
ionisation fraction or magnetic field strength (see e.g Jankovic et al.,
2021). A different way to parameterise the viscosity that removes this
variation with the temperature profile calculated by cuDisc is
$\nu=\nu_{0}(R/R_{0})$, where the linear dependence on radius comes from
assuming a mid-plane temperature profile proportional to $R^{-1/2}$. Both
methods of setting the viscosity are available in cuDisc. By default, $v_{0}$
is calculated assuming a mid-plane temperature of 100 K at 1 AU, but the user
can edit this.

###  3.5 Hydrostatic equilibrium

The 2D gas density structure is calculated by solving the hydrostatic
equilibrium equation,

$\frac{dP}{dZ}=-\frac{GM_{*}}{(R^{2}+Z^{2})^{\frac{3}{2}}}\rho_{g}Z.$ (32)

First, we rewrite the equation using the ideal gas law to obtain

$d\log P=\frac{GM_{*}}{c_{s}^{2}}d\left(\frac{1}{r}\right),$ (33)

where we have utilised $ZdZ/r^{3}=dr/r^{2}=-d(1/r)$, $r$ being the spherical
radius, and $c_{s}$ is the isothermal sound speed. We then take the
exponential of Eqn. 33 and in each cell $i,j$ calculate

$\frac{P_{i,j}}{P_{i,j-1}}=\exp\left[GM_{*}\left<\frac{1}{c_{s}^{2}}\right>_{i,j}\left(\frac{1}{r_{i,j}}-\frac{1}{r_{i,j-1}}\right)\right],$
(34)

where $\left<\frac{1}{c_{s}^{2}}\right>_{i,j}$ is the average reciprocal of
the sound speed squared,

$\left<\frac{1}{c_{s}^{2}}\right>_{i,j}=\frac{1}{2}\left[\frac{1}{c_{s,i,j}^{2}}+\frac{1}{c_{s,i,j-1}^{2}}\right].$
(35)

The pressure in each cell relative to the cell at $j=0$ is then found by
multiplicatively scanning over $Z$,

$\frac{P_{i,j}}{P_{i,j=0}}=\prod_{k=1}^{j}\frac{P_{i,k}}{P_{i,k-1}}.$ (36)

This method enforces positivity of the calculated pressure. The pressure is
then converted to density and normalised via the gas surface density,

$\rho_{g,i,j}=\frac{\Sigma_{g,i}(RdR)_{i}}{\sum_{k}\omega_{i,k}V_{i,k}}\omega_{i,j},$
(37)

where $\omega_{i,j}$ is the unnormalised density and $(RdR)_{i}$ is calculated
as
$(RdR)_{i}=\frac{1}{2}d(R^{2})_{i}=\frac{1}{2}[(R^{e}_{i+1})^{2}-(R^{e}_{i})^{2}]$.

###  3.6 Dynamics Tests

To demonstrate that the code is second-order accurate, we now show test
problems run using 2D Gaussian pulses with the analytic form

$\rho(x,y,t)=\frac{A}{t}\exp\left(-\frac{(x-x_{p}(t))^{2}+(y-y_{p}(t))^{2}}{4Dt}\right),$
(38)

where $x_{p}(t)=x_{0}+v_{x}(t-t_{0})$ with an equivalent expression for
$y_{p}$. For these tests, we used the following parameter choices: $A=D=1$,
$x_{0}=30$, $y_{0}=0$, $t_{0}=0.1$, $v_{x}=5$ and $v_{y}=2$. These pulses were
advected in both $R$ and $Z$ directions whilst undergoing diffusion from
$t=t_{0}$ to $t=1$. Cells were logarithmically-spaced in $R$ between 10 and 50
and linearly-spaced in $\theta$ between $-\pi/6$ and $\pi/6$. Cartesian areas
and volumes were used, and the boundaries were set to outflow conditions. An
example run can be seen in Fig. 5. The L2 error norm,
black$\sqrt{\sum_{ij}(\rho_{ij}-\rho_{\text{an},ij})^{2}/(N\times N)}$, was
calculated using the analytic solution, $\rho_{\text{an}}$ compared to the
computed solution, $\rho$. Fig. 6 shows how this error decreases with number
of cells $N$ on a 2D $N\times N$ grid. As shown, the L2 error is proportional
to $1/N^{2}$, demonstrating that our method is second-order accurate.

Figure 5: Advection-diffusion of a 2D gaussian pulse using the dynamics
routine in cuDisc. The left and right panels show the initial and final
states, respectively. This simulation had a resolution of 512 $\times$ 512
cells. Figure 6: The L2 error norm for advection-diffusion of a 2D gaussian
pulse using the dynamics routine in cuDisc. Slopes of $1/N$ and $1/N^{2}$,
where the 2D grid has size $N\times N$, are shown for reference.

To test our code on disc-related problems, we compared simulated steady-state
results to analytic results given by Takeuchi & Lin (2002). The gas was
initialised using their analytic density and velocity profiles. Three dust
species (10 $\mu$m, 100 $\mu$m and 1 mm) were initialised as being well-mixed
with the gas (i.e. the same density and velocity profiles) with a standard
dust-to-gas ratio of 0.01 applied to the density. In order to study vertical
settling, radial fluxes were set to 0 for comparison with Takeuchi & Lin’s
analytic profiles, and the code was run for 100,000 years to ensure a steady-
state in the dust, blackas the settling timescale for the smallest grains in
this set-up was 90,000 years. black300 cells were used in the vertical
dimension. Fig. 7 shows the computed equilibrium dust-to-gas profiles of the
dust species compared to the analytic solutions, showing strong agreement.
blackA small discrepancy is visible for the smallest grains at large $Z$; this
is due to our inclusion of the advection of momentum and because we do not
make a thin disc approximation in the gravitational potential (Takeuchi & Lin
invoke $Z<<R$ to simplify the problem). With radial velocities turned on, we
also match the analytic profiles for radial velocity as a function of height,
as shown in Fig. 8.

Figure 7: Steady-state dust-to-gas ratios as a function of height above the
disc mid-plane (in units of the gas scale height, $h_{g}$) at a radius of 10
AU for three dust species. The analytic solutions given by Takeuchi & Lin
(2002) are overplotted for comparison. Figure 8: Steady-state radial
velocities as a function of height above the disc mid-plane (in units of the
gas scale height, $h_{g}$) at a radius of 10 AU for three dust species. The
analytic solutions given by Takeuchi & Lin (2002) are overplotted for
comparison.

##  4 Temperature solver

Radiative transfer is a complex problem, and various numerical methods can be
employed whose results differ depending on the weights attributed to accuracy
versus computation time. For certain problems, calculated temperatures may
need to be highly accurate and need not evolve (e.g. comparing line fluxes to
observations), allowing long computation times to be a viable option. Other
problems may not require such high accuracy in the calculation of temperature
but instead may require fast calculation to assess the thermal evolution of a
system.

Since the aim of cuDisc is to study the evolution of protoplanetary discs over
long timescales ($\sim 10^{5}$ yr or longer), we require a fast temperature
solver. For this reason, our method uses a hybrid of ray-tracing and multi-
band flux-limited diffusion as opposed to exact techniques such as Monte
Carlo. We benchmark our solver against RADMC-3D (Dullemond et al., 2012), a
Monte Carlo code that requires longer computation times for more accurate
results.

###  4.1 Hybrid ray-tracing + multi-band flux-limited diffusion

In cuDisc, we use a hybrid radiative transfer method that employs ray-tracing
from the star for heating due to stellar radiation and multi-band flux-limited
diffusion (FLD) for treating re-emitted/scattered radiation. FLD is a model
developed by Levermore & Pomraning (1981) that aims to solve the radiative
transfer problem by treating the transport of radiative energy as a diffusion
problem. The diffusive flux is limited to ensure that radiation does not
propagate faster than the speed of light. The limiter is chosen to tend
towards the correct physical forms in limits of high and low optical depth.
The following scheme used in cuDisc is based on similar schemes detailed by
Kuiper et al. (2010); Commerçon et al. (2011); Bitsch et al. (2013). In
conservative form, the equations governing the internal energy density,
$\epsilon$, and radiative energy density in a given wavelength band $b$,
$E^{b}_{R}$, are given by

$\frac{\partial\epsilon}{\partial
t}=-\sum_{b}\left(\sum_{s}\rho_{s}\kappa^{b}_{P,s}\right)\left(f^{b}(T)B(T)-cE_{R}^{b}\right)+S_{\text{heat}},$
(39) $\frac{\partial E_{R}^{b}}{\partial
t}+\nabla\cdot\boldsymbol{F}^{b}=\left(\sum_{s}\rho_{s}\kappa_{P,s}^{b}\right)\left(f^{b}(T)B(T)-cE_{R}^{b}\right)+S_{\text{sca}}^{b},$
(40)

where $\rho_{s}$ is the mass density of species $s$, $\kappa_{P,s}^{b}$ is the
Planck mean opacity in wavelength band $b$ for species $s$, $B(T)=4\sigma
T^{4}$ accounts for radiative cooling, $S_{\text{heat}}$ accounts for heating
via stellar irradiation and viscous processes and $S_{\text{sca}}^{b}$ is the
scattered stellar radiation. As Eqn. 40 is for a particular wavelength band, a
Planck factor $f^{b}$ is included that accounts for the fraction of total
radiated thermal energy that is emitted in the particular band in question. We
do not include advective terms in our formulation. $\boldsymbol{F}^{b}$ is the
flux of radiative energy in a given band, and this term is calculated in flux-
limited diffusion as

$\boldsymbol{F}^{b}=-\frac{c\lambda^{b}}{\rho\kappa_{R}^{b}}\nabla E_{R}^{b},$
(41)

where $\kappa_{R}^{b}$ is the Rosseland mean opacity in the given band and
$\lambda^{b}$ the flux-limiter. In our multi-band scheme, the Planck opacities
in Eqns. 39 & 40 are calculated using absorption opacities, whilst the
Rosseland opacity is calculated using the total opacity, absorption plus
scattering. The choice of flux-limiter used for cuDisc is that given by Kley
(1989),

$\lambda^{b}=\begin{cases}\dfrac{2}{3+\sqrt{9+10R_{b}^{2}}},&\text{for
}R_{b}<2\vspace{0.2cm}\\\ \dfrac{10}{9+10R_{b}+\sqrt{81+180R_{b}}},&\text{for
}R_{b}>2\end{cases}$ (42)

where

$R_{b}=\frac{1}{\rho\kappa_{R}^{b}}\frac{|\nabla E_{R}^{b}|}{E^{b}_{R}}.$ (43)

Taking the limits of high and low optical depth ($\rho\kappa_{R}\gg 1$ &
$\rho\kappa_{R}\ll 1$) the flux tends to the diffusion limit,
${\boldsymbol{F}=-(c/3\rho\kappa_{R})\nabla E_{R}}$, and the free-streaming
limit, ${\boldsymbol{F}=-cE_{R}(\nabla E_{R}/|\nabla E_{R}|)}$, respectively.

The heating term $S_{\text{heat}}$ is split into stellar heating and viscous
heating. The viscous heating is calculated explicitly as

$S_{\text{visc}}=\frac{9}{4}\alpha\rho_{g}c_{s}^{2}\Omega,$ (44)

where $\alpha$ parameterises the effective turbulent viscosity of the disc
(Shakura & Sunyaev, 1973), $c_{s}$ is the isothermal sound speed, and $\Omega$
is the Keplerian angular velocity. Stellar heating in each cell is calculated
via ray tracing from the star using

$S_{*}=\sum_{B}\frac{L_{B}}{4\pi r^{2}\Delta
r}e^{-\tau_{B}}(1-e^{-\rho\kappa_{B}\Delta r})(1-a_{B}),$ (45)

where the sum is over the wavelength bands used for stellar radiation, $L_{B}$
is the stellar luminosity in each band, and $\Delta r$ is the (spherical)
radial length of the cell. By default, the number of bands used for stellar
heating ($\sim 100)$ is larger than the number used for the FLD routine
detailed above ($\sim 20$) as the heating calculation is computationally
cheap. For each wavelength band, the $e^{-\tau_{B}}$ term accounts for the
flux that has been attenuated (removed through extinction; i.e. absorption
plus scattering opacity) by the material along the radial path between the
star and the cell of interest, whilst $(1-e^{-\rho\kappa_{B}\Delta r})$
accounts for the amount of remaining flux that is attenuated in the cell.
$\kappa_{B}$ is the extinction of stellar radiation at the centre of each
wavelength band. The term $(1-a_{B})$ includes the albedo, $a_{B}$, to take
into account the effect of scattering. The albedo is calculated as the ratio
of scattering opacity to total extinction at each wavelength. The diffusive
flux generated by scattering at each wavelength is therefore given by

$S_{\text{sca},B}=\frac{L_{B}}{4\pi r^{2}\Delta
r}e^{-\tau_{B}}(1-e^{-\rho\kappa_{B}\Delta r})a_{B}.$ (46)

This is then binned into the wavelength bins used for the FLD routine and
included in Eqn. 40.

Figure 9: Comparison of mid-plane and surface temperatures calculated by
cuDisc and RADMC-3D for the Pinte et al. (2009) benchmark adapted to have a
grain distribution. The $x$-axis is plotted with units of cylindrical radius
minus the inner disc radius (0.1 AU). The left and right plots show models
with total dust masses of $1\times 10^{-2}$ $M_{\oplus}$ and 10 $M_{\oplus}$,
respectively.

The internal energy is given by ${\epsilon=\rho c_{V}T}$, where $\rho$ is the
total density of all species and $c_{V}$ is the specific heat at constant
volume, meaning that Eqn. 39 can be written in terms of $B$,

$\frac{\rho c_{V}}{16\sigma T^{3}}\frac{\partial B(T)}{\partial
t}=-\rho\kappa_{P}(B(T)-cE_{R})+S.$ (47)

This allows the system of equations to be set up and solved implicitly for the
quantities $B$ and $E_{R}$. This method ensures the correct equilibrium is
reached for large $\Delta t$ \- this being the standard case for our
simulations. We also note that the Commerçon et al. (2011) linearization
writes $B(T)$ in terms of $T$, instead of $T$ in terms of $B(T)$ as we do. The
discretised forms of Eqns. 47 & 40 for cell $i,j$ at time-step $n$ are

$\frac{\rho^{n}_{i,j}c_{V}}{16\sigma(T^{n}_{i,j})^{3}\Delta
t}B^{n+1}_{i,j}V_{i,j}\\\ \hskip
28.45274pt+\sum_{b}\left(\sum_{s}\rho^{n}_{s,i,j}\kappa_{P,s,i,j}^{b}\right)\left(f^{b}_{i,j}B^{n+1}_{i,j}-cE^{b,n+1}_{R,i,j}\right)V_{i,j}\\\
=\frac{\rho^{n}_{i,j}c_{V}}{16\sigma(T^{n}_{i,j})^{3}\Delta
t}B^{n}_{i,j}V_{i,j}+S^{n}_{\text{heat},i,j}V_{i,j},$ (48)

and

$\frac{E^{b,n+1}_{R,i,j}}{\Delta
t}V_{i,j}-\left(\sum_{s}\rho^{n}_{s,i,j}\kappa_{P,s,i,j}\right)\left(f^{b}_{i,j}B^{n+1}_{i,j}-cE^{b,n+1}_{R,i,j}\right)V_{i,j}+\\\
\sum^{i+1}_{k=i-1}\sum^{j+1}_{l=j-1}\tilde{D}^{b}_{kl}E^{b,n+1}_{R,k,l}=\frac{E^{b,n}_{R,i,j}}{\Delta
t}V_{i,j}+S^{b,n}_{\text{sca},i,j}V_{i,j}.$ (49)

$\tilde{D}^{b}_{kl}$ is the diffusion matrix, constructed by summing all terms
that apply to cell $k,l$ when finding $\tilde{F}_{\sigma}A_{\sigma}$ (Eqn. 21)
for each of the four interfaces around cell $i,j$. The construction of the
elements of $\tilde{D}_{kl}$ can be found in more detail in Appendix A. We
time-lag the diffusion constant in Eqn. 41 and the Planck factor, $f^{b}$
(i.e. they are evaluated at $t^{n}$ instead of $t^{n+1}$). The system of
equations defined by Eqns. 48 & 49 are solved using incomplete-LU
preconditioning and the Bi-Conjugate Gradient Stabilized (BiCGStab) solver
described in Naumov (2011). The convergence criterion is set by comparing the
maximum fractional residual between the left and right-hand sides of each
equation at each iteration to a user-controlled relative tolerance, set by
default to $10^{-4}$.

###  4.2 Opacities

Dust absorption and scattering opacities can be set in various ways in cuDisc.
Simple, functional forms that approximate the opacity dependence on grain size
and wavelength may be used, or opacity tables generated using packages such as
the DSHARP opacity library (Birnstiel et al., 2018) can be read in. If tables
are used, cuDisc can interpolate the opacity data to the user-defined grain
sizes and wavelength bins using piecewise-cubic Hermite interpolation (Fritsch
& Carlson, 1980) in both grain size and wavelength space. In general, more
wavelength bins are used for the cheaper stellar heating calculation in Eqn.
45 than are used in the FLD scheme. For conversion between these two
wavelength grids, the opacities are binned by calculating the Planck mean
opacity over the range of wavelengths corresponding to each coarse bin - these
partial means result in a better estimate of the temperature in optically thin
regions when a small number of bins are used.

###  4.3 Radiative transfer tests

To test our temperature solver, we ran modified versions of standard benchmark
tests used in the literature for comparing temperature solvers: the Pascucci
et al. (2004) and Pinte et al. (2009) tests. Fig. 9 shows the comparison of
mid-plane and surface temperatures calculated using cuDisc and the Monte Carlo
radiative transfer code RADMC-3D for the Pinte test, modified to include a
distribution of grain sizes. The density of each grain species $i$ was set by

$\rho_{i}(R,Z)={F_{i}}\frac{\Sigma_{\text{tot}}(R)}{\sqrt{2\pi}h_{i}(R)}\exp\left(-\frac{Z^{2}}{2h_{i}^{2}(R)}\right),$
(50)

where $F_{i}$ is the fraction of the total surface density
$\Sigma_{\text{tot}}$ distributed to each grain size and $h_{i}$ the scale
height for each grain. The grain size distribution was set according to the
MRN (Mathis et al., 1977) distribution with a maximum grain size of 0.5 cm,
and the total surface density was set according to Pinte et al. (2009) as

$\Sigma_{\text{tot}}=\Sigma_{0}(R/\text{AU})^{-1.625},$ (51)

where $\Sigma_{0}$ was set by the total desired disc mass. The scale height
for each grain size was set as the approximate height reached once grains are
settled (Dubrulle et al., 1995),

$h_{i}=h_{g}\sqrt{\frac{1}{1+\text{St}_{i}/\alpha}},$ (52)

where $h_{g}$ is the gas scale height, Sti is the grain Stokes number, set
simply in this problem as St${}_{i}=0.05a_{i}(R/\text{AU})$ where $a_{i}$ is
the grain radius, and $\alpha$ is the gas turbulence parameter, set to
$1\times 10^{-3}$. The gas scale height was set according to Pinte et al.
(2009) as $h_{g}=10(R/100\text{ AU})^{1.125}$ AU. The opacities of the dust
grains were set using the Birnstiel et al. (2018) tool for calculating
opacities, together with the dielectric constants given by Draine (2003).
Scattering was treated as being isotropic as cuDisc is currently unable to
treat anisotropic scattering. Fig. 9 shows the results for two different total
disc masses ($1\times 10^{-2}$ & 10 $M_{\oplus}$). These masses corresponded
to mid-plane optical depths from the star to the observer at 0.81 $\mu$m of 87
and $8.7\times 10^{4}$, respectively. These values are just over an order of
magnitude smaller than the optical depths for the Pinte test at the same disc
masses, as our use of a grain distribution lowers the overall opacity of the
disc due to the presence of large grains. The inner and outer bounds of the
grids were set to (0.1, 400) AU in $R$ and (0, $\pi/4$) rad in $\theta$. 150
equally-spaced cells were used in $\theta$, whilst the number of cells in $R$
was adapted depending on the disc mass. For both disc masses, 200
logarithmically-spaced cells were used between 0.15 & 400 AU, whilst to
minimise the optical depth in cells very close to the inner edge, we follow
the method detailed by Ramsey & Dullemond (2015) and use 40 & 96
logarithmically-spaced cells between 0.1 & 0.15 AU for the lower & higher disc
mass respectively. 100 logarithmically-spaced wavelengths between 0.1 & 3000
$\mu$m were used for stellar heating and binned down to 20 bands for the FLD
calculations. Apart from the binned wavelength bands (which are not required),
these same grid structures were used in RADMC-3D for consistency. Resolution
tests were performed, which confirmed that convergence was reached for the
grid resolutions quoted here.

We find that our results in the optically thin surface layers are accurate to
within a few percent for the entirety of the disc apart from the very inner
region between 0.1 and 0.2 AU. At the mid-plane, the accuracy varies but
remains within $\sim$20 % for most of the disc - a level of accuracy high
enough for our problems.

##  5 Dust growth and fragmentation

Dust coagulation is performed using the method outlined in Brauer et al.
(2008) but without the vertical integration that they use to convert the
problem to 1D along the disc mid-plane. In this method, the Smoluchowski
coagulation equation (Smoluchowski, 1916) becomes

$\dot{\rho}_{k}=\sum_{i=0}^{N-1}\sum_{j=0}^{i}\dot{\rho}_{\text{gain},ijk}-\sum_{i=0}^{N-1}\dot{\rho}_{\text{loss},ik}$
(53)

where dots represent time derivatives, $N$ is the total number of dust
species, $\rho_{k}$ is the density of the $k$-th dust species, and
$\dot{\rho}_{\text{gain},ijk}$ and $\dot{\rho}_{\text{loss},ik}$ are the gains
in density due to collision products of other grains (through coagulation and
fragmentation) and losses in density due to collisions between the $k$-th dust
species and all others respectively. The loss is given by

$\dot{\rho}_{\text{loss},ik}=m_{k}n_{i}n_{k}K_{ik}(p_{{\rm coag},ij}+p_{{\rm
frag},ij})$ (54)

where $m_{k}$ is the mass of the $k$-th grain, $K_{ij}$ is referred to as the
kernel, and $p_{{\rm coag},ij}$ and $p_{{\rm frag},ij}$ are the probabilities
of coagulation and fragmentation, which follow

$p_{\text{coag},ij}+p_{\text{frag},ij}=1$ (55)

in the absence of bouncing (as assumed here). The form of this kernel is taken
from Birnstiel et al. (2010) and is given by

$K_{ij}=\sigma_{ij}\Delta v_{ij},$ (56)

where $\sigma_{ij}$ is the collision cross-section of the grains and $\Delta
v_{ij}$ is the relative velocity of the grains, given by

$\Delta v_{ij}^{2}=v_{\text{BM},ij}^{2}+v_{\text{turb},ij}^{2}+(\Delta
v_{\text{lam},ij})^{2},$ (57)

where $v_{\text{BM},ij}$, $v_{\text{turb},ij}$, and $\Delta v_{\text{lam},ij}$
are the relative velocities due to Brownian motion, gas turbulence, and
laminar dust motion respectively. Here $\Delta v_{\text{lam},ij}$ is the root-
mean-square relative velocity between the dust particles, calculated using the
velocities determined by the dynamical solver. Brownian motion most strongly
affects small grains and is given by

$v_{\text{BM},ij}=\sqrt{\frac{8k_{B}T(m_{i}+m_{j})}{\pi m_{i}m_{j}}},$ (58)

where $k_{B}$ is the Boltzmann constant and $T$ the disc temperature. The form
of $v_{\text{turb},ij}$ is taken from Ormel & Cuzzi (2007) and is given by

$v_{\text{turb},ij}=\sqrt{\alpha}c_{s}V_{\text{rel},ij},$ (59)

where $V_{\text{rel},ij}$ is the relative velocity between grains scaled to
the eddy velocity $\sqrt{\alpha}c_{s}$,

$V_{\text{rel},ij}=\sqrt{\Delta V_{\text{I}}^{2}+\Delta V_{\text{II}}^{2}},$
(60)

where $\Delta V_{\text{I}}^{2}$ & $\Delta V_{\text{II}}^{2}$ are the relative
velocities induced by “slow” class I and “fast” class II eddies, respectively
given by Eqns. 17 & 18 in Ormel & Cuzzi (2007). For computational speed, we
use a quadratic fit to find the Stokes number of the boundary between class I
& II eddies for a given pair of grains, St∗ (Eqn. 21d in Ormel & Cuzzi
(2007)). This fit is accurate to within 2%.

In a collision, grains can coagulate to form more massive grains or fragment
to produce a remnant grain and a set of fragments that are distributed over
lower-mass grains. Combining all of these collision products gives us the
density gain for a dust species as

$\dot{\rho}_{\text{gain},ijk}=n_{i}n_{j}K_{ij}\left[p_{\text{coag},ij}(m_{i}+m_{j})C_{ijk}\right.\\\
\left.+p_{\text{frag},ij}(m_{\text{rem},ij}R_{ijk}+m_{\text{frag},ij}F_{ijk})\right],$
(61)

where $C_{ijk}$, $R_{ijk}$ and $F_{ijk}$ are the fractions of the mass
associated with coagulation products, remnant products and fragment products
that are distributed into bin $k$, and $m_{\text{rem},ij}$ and
$m_{\text{frag},ij}$ are the masses of remnants and fragments produced by the
collision. Note to avoid double counting we replace
$K_{ij}\rightarrow\frac{1}{2}K_{ij}$ when $i=j$.

The particle velocities are assumed to follow a Maxwell-Boltzmann distribution
around $\Delta v_{ij}$ to calculate the fragmentation probability. This is a
simplified version of the calculations done by Garaud et al. (2013), who
separated the stochastic (turbulent and Brownian) motions from the laminar
motions; we are therefore assuming that stochastic motions are larger than any
laminar relative motions. This is also the approach taken by Stammler &
Birnstiel (2022). In this way, the fragmentation probability is given by

$p_{\text{frag},ij}=\left[\frac{3}{2}\left(\frac{v_{\text{frag}}}{\Delta
v_{ij}}\right)^{2}+1\right]\exp\left[-\frac{3}{2}\left(\frac{v_{\text{frag}}}{\Delta
v_{ij}}\right)^{2}\right].$ (62)

The coefficients $C_{ijk}$, $R_{ijk}$ and $F_{ijk}$ add the collision products
into the appropriate mass bins. Since our standard choice of the mass grid is
logarithmic, these products do not map directly to other mass bins and must be
split over the two bins on either side of the product mass. Defining the
integers $\ell(m)$ and $u(m)$ as: $m_{\ell(m)}$ being the largest mass such
that $m_{\ell(m)}<m$ and $m_{u(m)}$ being the smallest mass such that
$m_{u(m)}>m$, we may write

$C_{ijk}=\begin{cases}\varepsilon(m_{i}+m_{j})&k=\ell(m_{i}+m_{j})\\\
1-\varepsilon(m_{i}+m_{j})&k=u(m_{i}+m_{j})\\\ 0&\text{otherwise,}\end{cases}$
(63)

and

$R_{ijk}=\begin{cases}\varepsilon(m_{\text{rem},ij})&k=\ell(m_{\text{rem},ij})\\\
1-\varepsilon(m_{\text{rem},ij})&k=u(m_{\text{rem},ij})\\\
0&\text{otherwise,}\end{cases}$ (64)

where

$\varepsilon(m)=\frac{m_{u(m)}-m}{m_{u(m)}-m_{\ell(m)}}.$ (65)

The distribution of fragments, $F_{ijk}$, is calculated following Rafikov et
al. (2020). First, we assume that the number density distribution of fragments
is given by a power law,

$n(m)dm\propto m^{-\eta}dm$ (66)

where $n(m)dm$ is the number of particles per unit volume in the mass range
$[m,m+dm]$ and $\eta$ is the fragmentation parameter that controls how
fragments are distributed amongst the smaller mass grains; by default, we take
the standard value of 11/6 found in the steady-state solutions of Dohnanyi
(1969); Tanaka et al. (1996). Converting to mass density and integrating over
the bins to find the fraction of mass distributed into each bin gives

$N_{k}(m_{\rm
max})=\frac{\left(m^{e}_{k+1}\right)^{2-\eta}-\left(m^{e}_{k}\right)^{2-\eta}}{m_{\rm
max}^{2-\eta}-\left(m^{e}_{0}\right)^{2-\eta}},$ (67)

where $m_{\rm max}$ is the mass of the largest fragment and where
$m^{e}_{k+1}$ and $m^{e}_{k}$ are the upper and lower edges of the
$k^{\text{th}}$ mass bin.

We then place each $m_{\text{frag},ij}$ into the bin $l$ by finding the
smallest $l$ such that $m^{e}_{l+1}>m_{\text{frag},ij}$ and setting $m_{{\rm
max},ij}=m^{e}_{l+1}$. From this we arrive at $F_{ijk}=N_{k}(m_{{\rm
max},ij})$. However, $F_{ijk}$ is not used directly. Instead, we first
determine the total amount of fragments falling into each $l$-bin before
distributing them according to $N_{k}(m^{e}_{l+1})$. This results in a much
more efficient algorithm, as outlined in Rafikov et al. (2020).

The proportion of mass ending up as fragments is defined as a multiple of the
mass of the smaller of the two colliding grains. Denoting the more massive
particle as the target with mass, $m_{\text{tar}}$ and the smaller as the
impactor with mass, $m_{\text{im}}$, we define the $\chi_{\text{im}}$ such
that
$m_{\text{frag},ij}=\text{min}\left(\chi_{\text{im}}m_{\text{im}},m_{\text{tar}}\right)$.
The remnant mass is then given by

$m_{\text{rem}}=\begin{cases}m_{\text{tar}}-\chi_{\text{im}}m_{\text{im}}&\text{for
}\chi_{\text{im}}m_{\text{im}}<m_{\text{tar}}\\\
0&\text{otherwise},\end{cases}$ (68)

where $\chi_{\text{im}}$ is a factor controlling the amount of material the
impactor can remove relative to its mass, usually taken to be unity. The
fragment mass is then given by the total mass of the colliding grains minus
the remnant mass.

To integrate Eqn. 53, cuDisc employs the Bogacki-Shampine 3(2) embedded Runge-
Kutta method with error estimation (Bogacki & Shampine, 1989) that adaptively
calculates the time-step required to meet user-specified relative and absolute
tolerances. As default, these are a 1% relative tolerance and an absolute
tolerance of $1\times 10^{-10}$ multiplied by the sum over all grain densities
in a cell.

###  5.1 Coagulation tests

To test our coagulation routine, we produced direct comparisons to the 1D disc
evolution code DustPy using a vertically integrated kernel (Birnstiel et al.,
2010), which can be seen in Fig. 10. The vertically integrated kernel is used
to study the growth and fragmentation of dust surface densities instead of
volume densities and is constructed by dividing the coagulation kernel (Eqn.
56) by $\sqrt{2\pi(h_{i}^{2}+h_{j}^{2})}$, where $h_{i}$ & $h_{j}$ are the
scale heights of the dust grains, assuming a Gaussian vertical density
profile. The relative vertical velocity is set according to the difference in
terminal settling velocities of the dust grains in question, taking the dust
grain scale height as its representative height $Z$ above the mid-plane. Each
simulation was run with dynamics turned off, at a single radius in the disc,
with the same initial conditions. The simulations were run until a steady-
state was reached in the dust grain size distribution. To compare the two
codes, we plot the mass-grid independent densities for each grain, given by

$\sigma_{d}(m)=m\frac{\partial\Sigma_{d}}{\partial m}.$ (69)

200 grain species were used for both simulations, as this is comparable to the
typical amount used for science calculations. The two codes show strong
agreement. The only slight difference at $\sim 7\times 10^{-4}$ cm is due to
different methods used by the two codes to fit the relative turbulent
velocities described in Ormel & Cuzzi (2007). The difference occurs because
the approximation used in DustPy is not continuous in this region.

We also ran tests using a constant kernel, the results of which were
consistent with the analytic solutions.

Figure 10: Comparison of steady-state dust grain size distributions calculated
by cuDisc and DustPy. The simulation parameters used were $R=20$ AU, $T=36$ K,
$\mu=2.4$, $\alpha=10^{-3}$ and $v_{\text{frag}}=100$ cm s-1.

##  6 Code structure

The overall time evolution is controlled by the CFL condition on the dust
dynamics. Gas, temperature and coagulation updates are then performed when
required depending on the evolution time-scales associated with each physical
process. A schematic of the code structure can be seen in Fig. 11. By default,
gas surface density updates and subsequent hydrostatic equilibrium
calculations are also performed at each CFL time-step due to their low
computational cost. For temperature calculations, an update is performed on
the first time-step of the simulation, and the percentage change with respect
to the initial state is calculated. How often the temperature solver should be
employed is then controlled by choosing a percentage tolerance for temperature
updates that allows the code to calculate the approximate time period that can
be allowed to elapse before the next update is required; by default, this is
set to 0.1%. However, this tolerance should be chosen depending on the
expected time-scales of temperature evolution for each specific problem.
Coagulation updates are managed differently, given that the coagulation
routine is explicit and, as such, performs its own sub-steps within each
global time-step (i.e. sub-cycling). To control how often these updates are
performed, the final sub-time-step within the coagulation routine is extracted
and used to compare to global time-steps. By default, the coagulation routine
is triggered when one sub-time-step taken from the last coagulation step has
elapsed in global simulation time; however, this can be changed by the user.

Figure 11: The schematic for cuDisc. Dust & gas dynamics are updated every
time-step whilst temperature and coagulation updates are performed when
certain conditions are met. These conditions can be controlled by the user in
order to align with the timescales relevant to the physical problem at hand.

##  7 Grain growth in a steady-state transition disc

To show a simple science test case and compare it to other works, we ran
simulations of a disc with an inner hole evolving towards a steady-state,
representative of a “transition” disc structure (e.g. Owen, 2016). To compare
to DustPy, one simulation was run with a vertically isothermal temperature
profile set to the mid-plane temperature adopted by DustPy, whilst one had the
2D temperature solver switched on. We will refer to these simulations as
isothermal and non-isothermal from here on. The gas surface density was set to
have a mass within 100 AU of 0.017 $M_{\odot}$ with a power law radial
dependence and a sharp exponential cut off at 5 AU; i.e.
$\Sigma_{g}=\Sigma_{0}R_{\text{AU}}^{-1}\exp{[-(5/R_{\text{AU}})^{10}]}$. For
these tests, the gas surface density was not evolved; however, the vertical
gas density profile was updated for the non-isothermal test to maintain
hydrostatic equilibrium. Radial gas velocities were set to 0, whilst azimuthal
velocities were updated with the temperature to maintain force balance. 150
logarithmically-spaced cells were used in $R$ between 3 & 20 AU, and a total
of 150 cells were used in $\theta$ between 0 & $\pi/4$, with 92 linearly-
spaced cells between 0 & $\pi/9$ with double the resolution of the remaining
linearly-spaced 58 cells. This sub-division was used to enhance the resolution
at the mid-plane. The viscosity, $\nu$, was set according to Shakura & Sunyaev
(1973) via $\alpha c_{s}H$. The parameters chosen for each of the simulations
are given in Table 1. The stellar radiation was assumed to be a blackbody at
the effective temperature given in Table 1. Note that simulations were run for
fragmentation velocities of 100 and 1000 cm s-1. For the 100 cm s-1 runs, 100
dust species with grain sizes logarithmically spaced between 0.1 $\mu$m and
0.5 cm were initialised with a size distribution according to the MRN
distribution (Mathis et al., 1977) and spatially distributed as being well-
mixed with the gas at a local dust-to-gas ratio of 0.01. For the 1000 cm s-1
runs, the mass grid was adjusted to 135 logarithmically spaced grain sizes
between 0.1 $\mu$m and 20 cm to account for the larger grains present.
black100 & 135 mass bins were used respectively to give $m_{n+1}/m_{n}=1.38$
for both sets of simulations. This choice maintains $>7$ bins per mass decade,
a requirement for accuracy in the coagulation routine (Ohtsuki et al., 1990).
For these calculations, changes in dust composition throughout the disc due to
sublimation were neglected. In this section we discuss the 100 cm s-1 runs,
whilst the 1000 cm s-1 runs are discussed in Section 7.2.

Fig. 12 shows the 2D dust density profiles of three grain sizes after 1 Myr of
evolution, at which time the simulations had reached a steady state, for both
the isothermal disc and non-isothermal disc, whilst Fig. 13 shows the
temperature profile of the non-isothermal disc after 1 Myr. Dust settling is
apparent from the stratification of the different grains. The non-isothermal
disc also exhibits an increase in the proportion of large grains due to the
cool interior that allows for the coagulation of grains up to larger sizes.
This is because the strength of the gas turbulent velocity is proportional to
the sound speed when assuming a Shakura & Sunyaev (1973) viscosity, and
reduced turbulent velocities allow particles to grow larger before reaching
the fragmentation limit. The temperature structure of the non-isothermal disc
also exhibits a super-heated surface layer with temperatures greater than the
blackbody equilibrium temperature; this is due to the small grains that occupy
the upper regions of the disc whose opacities mean they are good absorbers of
stellar optical photons, but bad emitters of their own infra-red photons (see
e.g. Chiang & Goldreich, 1997). Above the super-heated layer, the dust
densities drop to the floor value and the temperature is set by the opacities
given to the gas. By default, these opacities are set as the dust-to-gas ratio
floor value multiplied by the opacities of the smallest dust grains, as this
leads to a smooth temperature transition in the upper regions of the disc. In
reality, the temperature in the gaseous atmosphere above the dust disc is
controlled by other processes (see e.g. Woitke et al., 2009). In principle,
these processes can be included in the framework if necessary.

Parameter | Value
---|---
$\alpha$ | $10^{-3}$
$\mu$ | 2.4
$v_{\text{frag}}$ | 100, 1000 cm s-1
$M_{*}$ | 1 $M_{\odot}$
$R_{*}$ | 1.7 $R_{\odot}$
$T_{*,\text{eff}}$ | 4500 K
$M_{\text{disc}}$ | 0.017 $M_{\odot}$
black$\Sigma_{0}$ |  black250 g cm-2
black$\rho_{m}$ |  black1.6 g cm-3
Dust opacities | DSHARP mix
| (Birnstiel et al., 2018)
Table 1: Parameters used for the transition disc simulations detailed in
Section 7. Figure 12: Dust profiles for three grain sizes after 1 Myr of
evolution for a vertically isothermal disc (left) and a disc with the 2D
temperature solver switched on (right). These runs were set with
$v_{\text{frag}}=100$ cm s-1. Figure 13: Temperature profile after 1 Myr of
evolution for a disc with the 2D temperature solver switched on and
$v_{\text{frag}}=100$ cm s-1.  Figure 14: Vertically-integrated mass-grid
independent dust densities at a radius of 6 AU for each $v_{\text{frag}}=100$
cm s-1 simulation.

Figs. 14 & 15 compare the vertically integrated dust grain size distributions
for the different simulations at a radius of 6 AU and for the entire disc,
respectively. Vertically-integrated mass-grid independent densities were
calculated from the volume densities computed by cuDisc through

$\sigma_{d}(m)=\int_{-\infty}^{\infty}m\frac{\partial\rho_{d}}{\partial m}dz.$
(70)

Fig. 15 shows that the spatial distribution of dust is very similar for all
runs; the maximum grain size is set by fragmentation throughout the entire
disc due to the fairly low choice of 100 cm s-1 as the fragmentation velocity.
The non-isothermal run is slightly more peaked at 6 AU. This is because the
equilibrium spatial distribution is set by the balance between diffusion and
radial drift, and, as previously mentioned, the non-isothermal disc has lower
turbulence and, therefore, less diffusion of the dust out of the dust trap.
Fig. 14 shows that the simulations follow very similar distributions up to
$\sim 30\mu$m, with the non-isothermal run having an overall higher density
across all grain sizes due to the increased surface density at the dust trap.
The DustPy simulation exhibits a much more prominent trough at 30 - 40 $\mu$m
than the cuDisc simulations, and the upper cut-off in the size distribution
differ across the simulations; the DustPy and non-isothermal runs show similar
cut-offs whilst the isothermal run is lower.

Figure 15: Dust grain size distributions as a function of disc radius for each
$v_{\text{frag}}=100$ cm s-1 simulation. The fragmentation limit (given by
Birnstiel et al. 2012) is over-plotted on each distribution for reference.
Figure 16: Dust grain size distributions at 1, 2 and 3 gas scale heights above
the mid-plane for each $v_{\text{frag}}=100$ cm s-1 simulation. The heights in
$Z$ of the plotted gas scale heights in each plot are different as they are
calculated from the mid-plane temperatures of the cuDisc runs. For the
isothermal simulation, this is the same as DustPy, but for the non-isothermal
simulation, the gas scale heights are lower due to the cooler interior of the
non-isothermal disc.

###  7.1 Vertical structure comparison

To compare the vertical structure, we converted the dust surface densities
from DustPy into 2D volume densities by calculating the diffusion-settling
equilibrium in the vertical axis for each dust grain (see Takeuchi & Lin
2002), as might be done when taking a DustPy output and calculating a
radiative transfer simulation. Fig. 16 shows the grain size distributions at
1, 2 and 3 gas scale heights, where the gas scale height is calculated from
the mid-plane temperatures. This means the scale heights are the same for the
isothermal and DustPy runs, but lower for the non-isothermal run given the
cooler interior temperature. It can now be seen that the upper cut-off grain
size for the isothermal and DustPy runs are very similar at the mid-plane but
differ at increasing height. This arises because the fragmentation limit sets
the maximum grain size and the fragmentation limit is proportional to gas
density. The density decreases as we move away from the mid-plane in the 2D
cuDisc runs, leading to larger Stokes numbers and, therefore, larger turbulent
velocities. This lowers the fragmentation threshold on grain size. After
vertical integration, this leads to the lower cut-off grain size we see in
Fig. 14. For the non-isothermal simulation, we see a similar result but with
the cut-off at larger sizes due to the cooler interior. After vertical
integration, this leads to the larger maximum grain size seen in Fig. 14.
There are also fewer particles at high altitudes for the same reason: the
cooler temperature lowers the scale height of the disc, bringing the dust
closer to the mid-plane.

The trough in the grain size distribution arises at the particle size where
turbulent motions start to become a strong source of relative velocities
between particles. The location of the trough is dependent on the Reynolds
number, which represents the strength of turbulent motions. In cuDisc, the
Reynolds number varies as a function of height because it depends on the gas
density and sound speed; however, in DustPy, it is assumed to be equal to the
mid-plane value everywhere. This causes the trough location to vary as a
function of height in cuDisc, leading to a smearing out of the feature in the
vertically integrated density.

Figure 17: The vertical structure of three different grain sizes at a radius
of 6 AU for the two $v_{\text{frag}}=100$ cm s-1 cuDisc simulations over-
plotted on the calculated DustPy profiles. The gas scale heights are
calculated from the mid-plane temperatures of the respective cuDisc
simulations; as the non-isothermal simulation has a cooler mid-plane, the gas
scale height is lower in the right-hand panel.

Fig. 17 shows the density of three different grain sizes, $\sim 0.5$, $\sim
50$ and $\sim 150$ $\mu$m, as functions of height. The seemingly lower scale
height of the 150 $\mu$m grains in the cuDisc runs when compared to DustPy is
due to the decrease in the maximum grain size allowed by fragmentation with
height in the cuDisc runs discussed earlier. This decrease in maximum
allowable grain size can be seen by noting that 150 $\mu$m is around the peak
of the grain size distribution at the mid-plane in Fig. 16 but that at 1 gas
scale height, the peak has moved to $\sim$100 $\mu$m. This effect is less
noticeable for the small grains (i.e. 0.5 $\mu$m in Fig. 17) as the small
grain distribution is less affected by the change in the upper cut-off grain
size. The non-isothermal vertical profiles are more condensed than DustPy,
again due to the cooler interior temperatures leading to lower scale heights.

The vertical profiles found using cuDisc also show increases in the amount of
intermediate size (10-100 $\mu$m) grains away from the mid-plane, with some
grain sizes having a higher density at around a gas scale height than at the
mid-plane. To investigate this, we need to look at how the collision rates
vary with height. Fig. 18 shows the collision rate, Stokes numbers and
relative turbulent velocity of grains with sizes 19 and 24 $\mu$m as functions
of height for the non-isothermal run. Re is the Reynolds number that describes
the strength of turbulent viscosity over molecular viscosity, i.e.
Re$=\nu_{t}/\nu_{\text{mol}}$, whilst St∗ is the Stokes number for the
boundary between slow and fast eddies (class I and class II in the terminology
used by Ormel & Cuzzi 2007) for a particular particle. Slow eddies have turn-
over times longer than particle stopping times and induce large-scale
systematic motion in the grains, whilst fast eddies have turn-over times
smaller than the particle stopping times and, therefore, induce stochastic
motions in the particles (Völk et al., 1980). Slow eddies do not drive large
relative velocities between similarly sized grains, whereas fast eddies do
(Ormel & Cuzzi, 2007). Close to the mid-plane, the Stokes numbers of the
grains are smaller than the Stokes number associated with the smallest
turbulent eddies (Re-1/2), meaning the particles only experience slow eddies;
therefore, the relative velocities are low. As we move to higher regions of
the disc, the particle Stokes numbers are in the intermediate regime, and
particles feel the impact of fast eddies. These stochastic motions lead to
larger relative velocities between similarly sized particles. The densities of
the intermediate-sized grains are enhanced at these heights due to the
increase in collision rates in combination with the decrease in maximum grain
size with height discussed previously. Higher up in the disc, settling causes
a decrease in dust density that leads to a sharp decrease in collisions
regardless of large turbulent velocities, counteracting the enhancement seen
at around one gas scale height.

Figure 18: Collision rates, Stokes numbers and relative turbulent velocities
of two different grain sizes, 19 and 24 $\mu$m, as a function of height. St∗
is the boundary between class I and class II eddies for a particular particle.
Its minimum value is equal to the Stokes number of the smallest eddies,
Re-1/2, whilst its maximum value is equal to the Stokes number associated with
the largest eddies (one orbit), 1.

###  7.2 Particle sweeping

Figure 19: Dust grain size distributions as a function of disc radius for each
$v_{\text{frag}}=1000$ cm s-1 simulation. The fragmentation and drift limits
(given by Birnstiel et al. 2012) are over-plotted for reference in solid white
and dashed white respectively. Figure 20: The vertical structure of three
different grain sizes for the $v_{\text{frag}}=1000$ cm s-1 simulations. The
top and bottom rows show the profiles within the dust trap ($\sim$5.7 AU), and
away from the dust trap ($\sim 10$ AU) respectively. The density of the 0.5 cm
grain is at the floor value at 10 AU in all of the simulations due to being
above the drift limit. The gas scale heights are calculated from the mid-plane
temperatures of the respective cuDisc simulations.

These results differ from those seen by Krijt & Ciesla (2016), who found that
small grains become trapped in the disc mid-plane due to colliding with, and
sticking to, larger grains before vertical mixing can distribute them higher
in the disc. This “sweeping-up” of small grains by large grains decreases the
abundance of small grains at high altitudes in regions where particle
collision rates at the mid-plane are high compared to the rate of diffusion
due to turbulent processes. A reason why we do not see this effect may be that
the largest grains in our simulations are $\sim 150$ $\mu$m - these grains are
less strongly settled than the cm-sized grains reached in the Krijt & Ciesla
(2016) models, lessening the effect that leads to trapping in the mid-plane.
To investigate this, we ran the same set of simulations with a higher
fragmentation velocity of 1000 cm s-1. Fig. 19 shows the vertically-integrated
density as a function of radius and grain size for this set of simulations,
now with the addition of the approximate drift limit. Both the cuDisc and the
DustPy simulations show very similar dust distributions, with the majority of
the disc exhibiting drift-limited growth and the largest grains in the dust
trap reaching $\sim$ a few cm. Turning to the vertical structure, Fig. 20
shows the density profiles of three different dust species from each
simulation at the peak of the dust trap ($\sim$5.7 AU) on the top row and away
from the trap ($\sim 10$ AU) on the bottom row. The cuDisc dust profiles in
the trap clearly exhibit the trapping of smaller grains around the mid-plane,
as found by Krijt & Ciesla (2016). However, the effect is not visible at 10
AU, where the isothermal profiles closely follow the Takeuchi & Lin (2002)
analytic profiles calculated from the DustPy densities, and the non-isothermal
profiles differ at a few scale heights due to the hot surface layers that
increase the turbulent velocities, lofting the smaller grains to higher
heights. The lack of sweeping at these radii is to be expected if the cause of
the effect is the presence of large cm-sized grains, because here the the
maximum grain size is lower due to the disc being in the drift-limited regime.

Krijt & Ciesla (2016) found that the effect of sweeping should be important in
regions where the collisional timescale of grains is less than the timescale
associated with diffusion to higher regions of the disc. In terms of disc
quantities, this criterion can be written as
$\alpha/(\Sigma_{d}/\Sigma_{g})<1$. Fig. 21 shows the radial regions of the
discs in both the low and high fragmentation velocity simulations where this
criterion is satisfied. In the high fragmentation velocity disc, the condition
is only met around the peak of the dust trap - this concurs with our findings,
as shown in Fig. 20. For the low fragmentation velocity case, the condition is
satisfied for the entirety of the disc; but we find no evidence of trapping.
We suspect this is due to the difference in scale height between the largest
and smallest grains in the simulations as discussed above; for the low
fragmentation velocity simulation, the difference between the scale heights of
the 0.5 and 150 micron grains is a factor of $\sim 2$, whilst, for the high
fragmentation velocity simulation, the difference in scale height between the
0.5 micron and 1 cm grains is a factor of $\sim 80$. In the low fragmentation
velocity case, this means that any sweeping by the largest grains does not
manifest itself as an appreciable change in the vertical density distribution
of the smaller grains, as the sweeping grains and swept-up grains have similar
scale heights. The effect is noticeable, however, in the high fragmentation
velocity case, where the largest grains are much more settled than the
smallest grains.

Figure 21: The regions of the disc within which the Krijt & Ciesla (2016)
criterion for mid-plane trapping of small grains is satisfied, shown as
hatched regions, for both the low and high fragmentation velocity cuDisc
simulations.

###  7.3 Diffusion-settling-coagulation equilibrium

These results indicate that the diffusion-settling equilibrium profiles
(Takeuchi & Lin, 2002) for the dust vertical structure do not fully describe
the dust population in regions where collisions are important. In these
regions, we suspect that a diffusion-settling-coagulation equilibrium is
established; the form of which varies depending on the fragmentation velocity.
For a low fragmentation velocity of 100 cm s-1, we find enhancements of
intermediate-sized dust grains at $\sim$ one gas scale height, whilst for a
high fragmentation velocity of 1000 cm s-1, we find that large grains sweeping
up smaller grains in the dust trap leads to the enhancement of small grains at
the disc mid-plane. These findings may have implications for the analysis of
observed disc spectral energy distributions (SEDs) as one must assume a disc
structure to estimate the emission layers of different-sized dust grains (see
e.g. D’Alessio et al., 2006). Scattered light images could also be affected,
as mid-plane densities calculated from the observed small grain distribution
at large disc heights are dependent on the assumed vertical structure of said
grains. However, our results also demonstrate that if one only cares about
overall, general, disc evolution and not specific problems that require a
resolved vertical dimension (e.g. winds or temperature instabilities caused by
vertical structure), the 1D results from DustPy are a good approximation to
the full 2D problem.

##  8 Summary

cuDisc is a new protoplanetary disc code that aims to allow long time-scale
calculations of discs with self-consistently calculated dynamics,
thermodynamics and dust grain growth and fragmentation. Modelling these
physical processes alongside one another allows us to answer important
problems in disc physics, such as structure formation due to instabilities and
dust removal. With the use of GPU acceleration, cuDisc enables simulations to
run for large fractions of the disc lifetime, removing the need for running
simplified secular models. We have shown that 2D structure affects the dust
spatial and grain size distributions even for simple systems, such as dust
trapped in the pressure bump of a transition disc. This may be important when
analysing SEDs and scattered light images, given the need to assume some model
of grain vertical structure. We also find that for studying overall disc
evolution, 1D models can quite accurately concur with 2D models. More features
can and will be added to the code as we continue to develop it, such as more
detailed dust microphysics and ice-vapour chemistry, making many other
problems in disc physics able to be investigated through the use of cuDisc.

## Acknowledgements

blackWe would like to thank the referee for their detailed and helpful
comments. RAB and JEO are supported by Royal Society University Research
Fellowships. This work has received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 853022, PEVAP) and a Royal Society Enhancement
Award. Some of this work was performed using the Cambridge Service for Data
Driven Discovery (CSD3), part of which is operated by the University of
Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital
funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC
operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.
For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC-BY) licence to any Author Accepted Manuscript version arising.

## Data Availability

The GitHub repository for cuDisc can be found at
https://github.com/cuDisc/cuDisc/.

## References

  * Andrews (2020) Andrews S. M., 2020, Annual Review of Astronomy and Astrophysics, 58, 483
  * Bae et al. (2022) Bae J., Isella A., Zhu Z., Martin R., Okuzumi S., Suriano S., 2022, arXiv e-prints, p. arXiv:2210.13314
  * Balbus & Papaloizou (1999) Balbus S. A., Papaloizou J. C. B., 1999, ApJ, 521, 650
  * Batalha et al. (2013) Batalha N. M., et al., 2013, ApJS, 204, 24
  * Birnstiel et al. (2010) Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79
  * Birnstiel et al. (2012) Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148
  * Birnstiel et al. (2018) Birnstiel T., et al., 2018, ApJ, 869, L45
  * Bitsch et al. (2013) Bitsch B., Crida A., Morbidelli A., Kley W., Dobbs-Dixon I., 2013, A&A, 549, A124
  * Bogacki & Shampine (1989) Bogacki P., Shampine L., 1989, Applied Mathematics Letters, 2, 321
  * Booth & Clarke (2021) Booth R. A., Clarke C. J., 2021, MNRAS, 502, 1569
  * Booth et al. (2017) Booth R. A., Clarke C. J., Madhusudhan N., Ilee J. D., 2017, MNRAS, 469, 3994
  * Brauer et al. (2008) Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859
  * Chiang & Goldreich (1997) Chiang E. I., Goldreich P., 1997, ApJ, 490, 368
  * Clarke & Pringle (1988) Clarke C. J., Pringle J. E., 1988, MNRAS, 235, 365
  * Commerçon et al. (2011) Commerçon B., Teyssier R., Audit E., Hennebelle P., Chabrier G., 2011, A&A, 529, A35
  * Courant et al. (1928) Courant R., Friedrichs K., Lewy H., 1928, Mathematische Annalen, 100, 32
  * Cuzzi et al. (1993) Cuzzi J. N., Dobrovolskis A. R., Champney J. M., 1993, Icarus, 106, 102
  * D’Alessio et al. (2006) D’Alessio P., Calvet N., Hartmann L., Franco-Hernández R., Servín H., 2006, ApJ, 638, 314
  * Dohnanyi (1969) Dohnanyi J. S., 1969, J. Geophys. Res., 74, 2531
  * Draine (2003) Draine B. T., 2003, ApJ, 598, 1026
  * Dubrulle et al. (1995) Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237
  * Dullemond et al. (2012) Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, record ascl:1202.015
  * Eistrup (2023) Eistrup C., 2023, ACS Earth and Space Chemistry, 7, 260
  * Ercolano & Pascucci (2017) Ercolano B., Pascucci I., 2017, Royal Society Open Science, 4, 170114
  * Franz et al. (2020) Franz R., Picogna G., Ercolano B., Birnstiel T., 2020, A&A, 635, A53
  * Fritsch & Carlson (1980) Fritsch F. N., Carlson R. E., 1980, SIAM Journal on Numerical Analysis, 17, 238
  * Fuksman & Klahr (2022) Fuksman J. D. M., Klahr H., 2022, ApJ, 936, 16
  * Garaud et al. (2013) Garaud P., Meru F., Galvagni M., Olczak C., 2013, ApJ, 764, 146
  * Huang & Bai (2022) Huang P., Bai X.-N., 2022, ApJS, 262, 11
  * Hutchison & Clarke (2021) Hutchison M. A., Clarke C. J., 2021, MNRAS, 501, 1127
  * Jankovic et al. (2021) Jankovic M. R., Owen J. E., Mohanty S., Tan J. C., 2021, MNRAS, 504, 280
  * Kley (1989) Kley W., 1989, A&A, 208, 98
  * Kley (1998) Kley W., 1998, A&A, 338, L37
  * Krijt & Ciesla (2016) Krijt S., Ciesla F. J., 2016, ApJ, 822, 111
  * Krumholz et al. (2020) Krumholz M. R., Ireland M. J., Kratter K. M., 2020, MNRAS, 498, 3023
  * Kuiper et al. (2010) Kuiper R., Klahr H., Dullemond C., Kley W., Henning T., 2010, A&A, 511, A81
  * Leveque (2004) Leveque R. J., 2004, Journal of Hyperbolic Differential Equations, 01, 315
  * Levermore & Pomraning (1981) Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321
  * Madhusudhan (2019) Madhusudhan N., 2019, Annual Review of Astronomy and Astrophysics, 57, 617
  * Mathis et al. (1977) Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
  * Mayor et al. (2011) Mayor M., et al., 2011, arXiv e-prints, p. arXiv:1109.2497
  * Mignone (2014) Mignone A., 2014, Journal of Computational Physics, 270, 784
  * Naumov (2011) Naumov M., 2011, NVIDIA White Paper
  * Öberg et al. (2011) Öberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJ, 743, L16
  * Ohtsuki et al. (1990) Ohtsuki K., Nakagawa Y., Nakazawa K., 1990, Icarus, 83, 205
  * Ormel & Cuzzi (2007) Ormel C. W., Cuzzi J. N., 2007, A&A, 466, 413
  * Owen (2016) Owen J. E., 2016, Publ. Astron. Soc. Australia, 33, e005
  * Owen (2020) Owen J. E., 2020, MNRAS, 495, 3160
  * Owen & Kollmeier (2019) Owen J. E., Kollmeier J. A., 2019, MNRAS, 487, 3702
  * Pascucci et al. (2004) Pascucci I., Wolf S., Steinacker J., Dullemond C. P., Henning T., Niccolini G., Woitke P., Lopez B., 2004, A&A, 417, 793
  * Philippov & Rafikov (2017) Philippov A. A., Rafikov R. R., 2017, ApJ, 837, 101
  * Pinte et al. (2009) Pinte C., Harries T. J., Min M., Watson A. M., Dullemond C. P., Woitke P., Ménard F., Durán-Rojas M. C., 2009, A&A, 498, 967
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Rafikov et al. (2020) Rafikov R. R., Silsbee K., Booth R. A., 2020, ApJS, 247, 65
  * Ramsey & Dullemond (2015) Ramsey J. P., Dullemond C. P., 2015, A&A, 574, A81
  * Rodenkirch & Dullemond (2022) Rodenkirch P. J., Dullemond C. P., 2022, A&A, 659, A42
  * Roe (1981) Roe P., 1981, Journal of Computational Physics, 43, 357
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Smoluchowski (1916) Smoluchowski M. V., 1916, Zeitschrift fur Physik, 17, 557
  * Stammler & Birnstiel (2022) Stammler S. M., Birnstiel T., 2022, ApJ, 935, 35
  * Stone & Gardiner (2009) Stone J. M., Gardiner T., 2009, New Astronomy, 14, 139
  * Takeuchi & Lin (2002) Takeuchi T., Lin D. N. C., 2002, ApJ, 581, 1344
  * Tanaka et al. (1996) Tanaka H., Inaba S., Nakazawa K., 1996, Icarus, 123, 450
  * Urpin (1984) Urpin V. A., 1984, Soviet Ast., 28, 50
  * Völk et al. (1980) Völk H. J., Jones F. C., Morfill G. E., Roeser S., 1980, A&A, 85, 316
  * Watanabe & Lin (2008) Watanabe S.-i., Lin D. N. C., 2008, ApJ, 672, 1183
  * Winn & Fabrycky (2015) Winn J. N., Fabrycky D. C., 2015, Annual Review of Astronomy and Astrophysics, 53, 409
  * Woitke et al. (2009) Woitke P., Kamp I., Thi W. F., 2009, A&A, 501, 383
  * Wootten & Thompson (2009) Wootten A., Thompson A. R., 2009, IEEE Proceedings, 97, 1463
  * Wu & Lithwick (2021) Wu Y., Lithwick Y., 2021, ApJ, 923, 123
  * Wu et al. (2012) Wu J., Gao Z., Dai Z., 2012, Journal of Computational Physics, 231, 7152
  * Zhu & Dong (2021) Zhu W., Dong S., 2021, Annual Review of Astronomy and Astrophysics, 59, 291

##  Appendix A Diffusion Matrix

The full form of $\tilde{D}_{kl}$ is found by writing out $\tilde{F}_{\sigma}$
(Eqn. 21) for each interface $\sigma$ about cell $i,j$ and summing all terms
that apply to cell $k,l$ in the 3$\times$3 stencil of cells centred on cell
$i,j$. For this exercise, we will change cell references from the form $i,j$
to the single letter form $k$, for ease of reading. As an example, referring
to Fig. 3, when calculating the fluxes over the interfaces $\sigma$ and
$\lambda$, $\tilde{F}_{\sigma}$ and $\tilde{F}_{\lambda}$ respectively,
according to Eqn. 25 we find

$\tilde{F}_{\sigma}=-D_{k}w_{k\sigma}\left[w_{l\sigma}c_{kl}(u_{l}-u_{k})+w_{m\tau}c_{km}(u_{m}-u_{k})\right]+\\\
D_{l}w_{l\sigma}\left[w_{k\sigma}c_{lk}(u_{k}-u_{l})+w_{p\nu}c_{lp}(u_{p}-u_{l})\right],$
(71)
$\tilde{F}_{\lambda}=-D_{a}w_{a\lambda}\left[w_{l\lambda}c_{al}(u_{l}-u_{a})+w_{m\beta}c_{am}(u_{m}-u_{k})\right]+\\\
D_{l}w_{l\lambda}\left[w_{a\lambda}c_{l}a(u_{l}-u_{a})+w_{b\epsilon}c_{lb}(u_{b}-u_{l})\right],$
(72)

where $\nu$, $\beta$ and $\epsilon$ are the interfaces between cells $l$ &
$p$, $a$ & $m$ and $l$ & $b$ respectively.

Taking for example cell $m$, the value in the diffusion matrix $\tilde{D}_{m}$
is therefore given by the sum of all terms that are multiplied by $u_{m}$,
multiplied by their respective interface areas A, i.e.

$\tilde{D}_{m}=-D_{k}w_{k\sigma}w_{m\tau}c_{km}A_{\sigma}-D_{a}w_{a\lambda}w_{m\beta}c_{am}A_{\lambda}.$
(73)