# Developing a non-Newtonian fluid model for dust, for application to astrophysical flows

Elliot M. Lynch, Guillaume Laibe

► To cite this version:

Elliot M. Lynch, Guillaume Laibe. Developing a non-Newtonian fluid model for dust, for application to astrophysical flows. Journal of Fluid Mechanics, 2024, 1001, pp.A37. 10.1017/jfm.2024.1088. hal-04879729

HAL Id: hal-04879729

https://hal.science/hal-04879729v1

Submitted on 10 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

# Developing a Non-Newtonian Fluid Model for Dust, for Application to Astrophysical Flows

Elliot M. Lynch ¹†, and Guillaume Laibe,¹

¹Univ Lyon, Univ Lyon 1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France

(Received xx; revised xx; accepted xx)

In the astrophysics community it is common practice to model collisionless dust, entrained in a gas flow, as a pressureless fluid. However a pressureless fluid is fundamentally different from a collisionless fluid - the latter of which generically possess a non-zero anisotropic pressure or stress tensor. In this paper we derive a fluid model for collisionless dust, entrained in a turbulent gas, starting from the equations describing the motion of individual dust grains. We adopt a covariant formulation of our model to allow for the geometry and coordinate systems prevalent in astrophysics, and provide a closure valid for the accretion disc context. We show that the continuum mechanics properties of a dust fluid corresponds to a higher-dimensional anisotropic Maxwell fluid, after the extra dimensions are averaged out, with a dynamically important rheological stress tensor. This higher-dimensional treatment has the advantage of keeping the dust velocity and velocity of the fluid seen, and their respective moments, on the same footing. This results in a simplification of the constitutive relation describing the evolution of the dust Rheological stress.

**Key words:** Particle/fluid flows – Multiphase and Particle-laden Flows, Turbulent mixing – Mixing, Gas dynamics – Compressible Flows

**MSC Codes** (Optional) Please enter your MSC Codes here

## 1. Introduction

The dynamics of dust in turbulent flows is important to a wide array of astrophysical, geophysical and engineering applications. In the case of astrophysical applications, dusty astrophysical fluids often combine a high Mach number with subsonic turbulence which feeds off of a Rayleigh stable shear flow. The dust number density is typically much lower than that of the gas, such that dust-dust collisions are infrequent. However, dust particles are typically too numerous to be kept track of individually. As such there is a need to be able to model the dynamics of weakly collisional/collisionless dust in turbulent gasses effectively.

The most physically accurate method of evolving dust grains in fluids is an N-body

<footnote>† Email address for correspondence: elliot.lynch@ens-lyon.fr</footnote>
approach where each solid particle is evolved independently (although this approach can
still exhibit spurious trapping behaviour, Commerçon et al. 2023). However, this approach is
typically prohibitively expensive for practical computations in the astrophysics setting, due to
the large range of lengthscales and number of dust particles involved, except on the smallest
of scales. Two common methods are used to make modelling dust dynamics computationally
tractable. One is to significantly reduce the number of dust grains compared with reality,
or to treat N-body particles as a dust aggregate; for instance the dust module in Athena,
(Bai & Stone 2010; Zhu et al. 2014), and PLUTO, (Mignone et al. 2019), and superparticle
implementations by Youdin & Johansen (2007); Balsara et al. (2009) and Yang & Johansen
(2016). This is commonly used when there is no back reaction or interaction between dust
grains as the number of particles required to achieve convergence will be much lower. In
accretion disc simulations, making use of such methods, it is common to employ of
order 10 particles per cell (Laibe & Price 2012a), which is not sufficient to adequately sample
the particle velocity distribution (Peirano et al. 2006). On smaller scales, in particular for the
small/incompressible shearing box (Latter & Papaloizou 2017), adequate particle resolution
may be possible with current computational resources and would provide an excellent check
on models capable of simulating the global disc scale. The second method is to treat the
dust as a continuous fluid (Barrière-Fouchet et al. 2005; Laibe & Price 2012a,b, 2014; Lin
& Youdin 2017; Lin 2019; Bi et al. 2021). In this paper we shall derive such a fluid model,
starting from a stochastic differential equation for the motion of individual grains entrained
in a turbulent gas flow.

The most common model of a dust fluid (in the astrophysics community) is to model it as a pressureless fluid coupled to the gas via the drag terms (as has been done in Barrière-Fouchet et al. 2005; Laibe & Price 2012a,b, 2014; Lin & Youdin 2017; Lin 2019; Bi et al. 2021). The justification for treating the dust as a pressureless fluid is that when the dust number density is much lower than that of the gas, dust-dust collisions are unimportant to the dust dynamics (although could be important for fragmentation/coagulation) which is dominated by gravity and the dust-gas interaction. As dust collisions are unimportant the dust, according to the literature, can be treated as pressureless. Unfortunately this argument for pressureless dust is flawed due to a misunderstanding about the micro-physical origin of pressure in a fluid.

The issue with this argument is that it conflates fluid pressure with collisionality. However
fluid pressure is not a measure of fluid collisionality, but instead is a measure of the mean
squared (density weighted) velocity dispersion of the particles. Crucially a collisionless fluid
can have a non-zero velocity dispersion, and will thus have a non-zero pressure tensor. In fact
weakly collisional/collisionless fluids often have large anisotropic pressure tensors and the
hydrodynamical description of the fluid breaks down, not because fluid properties such as
pressure and density are not defined, but because of the difficulty in truncating the moment
expansion, used to derive hydrodynamics from kinetic theory, at finite order (Chapman &
Cowling 1990; Grad 1948, 1949; Bobylev 1982, 2018). Collisions in a fluid are not the source
of pressure - instead the effect of collisions is to ensure that the moment expansion truncates
by damping higher order moments, along with isotropising the fluid pressure tensor (e.g.
Levermore 1996, see also Boltzman's H-theorem). In conclusion, while there is a strong
argument that dust in astrophysical fluids (and many geophysical fluids) can be approximated
as being collisionless, we cannot conclude, a priori, that the dust pressure is negligible. In
addition to this pressure from the particle motion, in turbulent gas-dust mixture there is an
additional dust Reynolds stress from the turbulent motion.

Stochastic differential equations (SDEs) have been used to model turbulent motion in
fluids (e.g. Thomson 1987; Pope 1987; Sawford 1991; Minier et al. 2004). Various authors
have extended such stochastic models of turbulent fluids to describe the motion of dust
grains entrained in the flow (e.g. Dubrulle et al. 1995; Minier 2001; Carballido et al. 2006;
Youdin & Lithwick 2007; Minier et al. 2014; Minier 2015; Ormel & Liu 2018; Laibe et al. 2020; Booth & Clarke 2021). Dubrulle et al. (1995), Carballido et al. (2006), Fromang & Papaloizou (2006), Ormel & Liu (2018), Laibe et al. (2020) and Booth & Clarke (2021) used their models to calculate the steady state vertical structure of a dust layer in an astrophysical disc. Youdin & Lithwick (2007) calculated the dust velocity correlations in a rotating shear flow and, importantly for our work, calculated a dust fluid model by preforming a moment expansion of the Fokker-Planck equation associated with the stochastic dust motion.

In this paper we develop a dust fluid model starting from a system of stochastic differential equations (SDEs) describing the motion of a single dust grain in a turbulent gas. To do this we shall preform a moment expansion of the Fokker-Planck equation associated with the SDEs, similar to that preformed Youdin & Lithwick (2007) but without the restrictive assumption that the correlation time is the shortest timescale in the problem, and adopt a closure valid for the accretion disc context. This approach differs from the more commonly adopted method of Reynolds averaging the pressureless two-fluid model and including a closure relation motivated by the interaction of dust grains with individual turbulent eddies (e.g. adopted by Binkert 2023). Our approach makes use of a novel 6-dimensional formulation, which keeps the dust velocity and velocity of the fluid seen, along with their moments, on the same footing. In this formulation the dust Kinetic tensor, Reynolds stress for the fluid seen and dust-gas cross correlation tensor combine into a single 6-dimensional stress tensor, which is advected by the flow. We adopt a covariant formulation of the dust fluid equation so that the model can be adapted to non-Cartesian coordinates often adopted in astrophysics problems. This will also allow for the adoption of orbital coordinates systems (e.g. Ogilvie & Latter 2013b; Ogilvie & Barker 2014), which will facilitate the study of distorted (elliptical or warped) dust discs. Finally, we explore the physical properties of our dust fluid model and consider the behaviour of the dust stress tensor in a rotating shear flow. Studying the behaviour of the dust fluid in rotating shear flows allows us to connect our model to problems in astrophysical and experimental fluid dynamics (accretion discs and dusty Taylor-Couette flows respectively). This may provide a basis to experimentally test the model in the lab.

In Section 3 we consider a SDE for motion of a single dust grain in a turbulent gas disc. In Section 4 we derive the dust fluid equations by performing a moment expansion of the Fokker-Planck equation associated with the SDE introduced in Section 3 and discuss our closure scheme. Section 5-7 describes the physics of the model. Section 5 discusses the dust fluid physics and highlights key properties of the model. Section 6 considers the hyperbolic structure, and wave modes, of the dust fluid equations. Sections 7 looks at the behaviour of the dust rheological† stress tensor in a rotating shear flows. In Section 8 we suggest possible refinements that could be made to the model. We present our conclusions in Section 9 and further mathematical derivations are given in the Appendices.

## 2. Overview of Astrophysical Flows

In this section we shall briefly outline the key properties of the astrophysical fluids, which are the primary motivation for developing this model, for the benefit of non-astrophysicists. The primary flow of interest are protoplanetary discs and other dusty accretion discs, with an additional interest in dusty quasi-spherical flows present in star formation and dusty planetary atmospheres. Focusing on accretion discs - these are disc like structures of gas and solid matter in approximately Keplerian rotation about 1 (or more) central object which dominate the gravitational field. The gas in such a system has the following properties

<footnote>† whenever we speak of the dust rheology or rheological stress we are referring to the rheology of the dust fluid and not the, entirely separate, rheology of the individual solid dust grains.</footnote>
* The flow in the inertial frame, stationary with respect to the centre of mass of the system, is highly hypersonic. However, in the fluid frame it principally behaves like a subsonic shear flow in a rapidly rotating frame.

* The geometry of the flow naturally lends itself to using cylindrical or spherical coordinates, both for simplifying analytical treatment and for improved angular momentum conservation, diffusivity and speed of numerical schemes.

* The discs are Rayleigh stable, however they can exhibit subsonic hydrodynamical or magnetohydrodynamical turbulence. Magnetohydrodynamical turbulence in discs (due to the magneto-rotational instability Balbus & Hawley 1991; Hawley & Balbus 1991; Hawley et al. 1995) is much stronger than hydrodynamical turbulence (e.g. Vertical Shear Instability Nelson et al. (2013); Lin & Youdin (2015); Flock et al. (2017); Svanberg et al. (2022) or parametric instability Papaloizou (2005a,b); Ogilvie & Latter (2013a); Barker & Ogilvie (2014)). However, discs that are cool enough for the presence of dust are typically too cool to be well ionised, which tends to suppress the action of the magnetic fields. Thus turbulence in such discs is expected to be hydrodynamical and very subsonic.

* The disc is stratified with a pressure scale height $H \sim R/M$, where R is the cylindrical distance from the central object and M is the Mach number. This vertical confinement gives the disc a shallow-water like character and is also important for setting the maximum size of turbulent eddies. The rapid rotation means the eddies (inertial waves) are predominantly vertical with vertical extend $\sim$ the scale height.

* Characteristic timescales are the orbital period of $\sim 1$ day-$10^3$ years (depending on the position in the disc). Characteristic lengthscale are the scale height $H \le 0.1R$ and cylindrical radius $R \sim 0.1-100$ Au (astronomical units) $\sim 10^7-10^{10}$ km.

* Molecular viscosity is typically sufficient low that it can be neglected (Although the Kolmogorov scale is $\sim 10$m Armitage (2020)).

The typical properties of dust in protoplanetary discs and prestellar cores are

* The dust is polydispersed with size $\sim \mu$m - 10 cm and forms a near continuous distribution in size space, however we will only consider the monodispersed case in this paper. For computational reasons most simulations of dusty accretion discs are monodispersed at present. The monodispersed case is also of observational interest as observations tend to be sensitive to a narrow range in size space dependant on the observational wavelength.

* Dust to gas mass ratio is typically $\ge 0.01$, with the vast majority of the mass in the largest grains (Testi et al. 2014).

* Total number of grains $\ge 1$ mm is $\sim 10^{32}$. The dust number density is $n \sim 10^{-9}cm^{-3}$ this corresponds to $\sim 10^{27}$ particles per cubic scale height (Testi et al. 2014; Lesur et al. 2022).

* The mean free path for dust-dust collision is $\sim 10^5$ km, with the collision timescale being typically much longer than the stopping time.

### 3. Stochastic differential equation for dust particle motion in a dust disc.

Consider a dust grain entrained in a gas flow, in the Epstein regime, where the gas velocity at the dust grain position is denoted $v^g$. The position $x$ and velocity $v$ for a dust particle, subject to force per unit mass, $f$, and gas drag, are given by the following set of differential equations:

$$dx_i = v_i dt, \qquad (3.1)$$

$$dv_i = f_i dt - \frac{1}{t_s} (v_i - v_i^g) dt, \qquad (3.2)$$
where $t_s$ is the stopping time for the dust particle under consideration. Typically we take the force per unit mass to be due to gravity with $f_i = -\nabla_i \phi$, where $\phi$ is the gravitational potential. Here $x_i, v_i, v_i^g$ and $f_i$ are the covariant components of the vectors $\mathbf{x}, \mathbf{v}, \mathbf{v}^g$ and $\mathbf{f}$ respectively. These are related to the contravariant components, $x^i, v^i, v_i^g$ and $f^i$ via the metric tensor $\gamma_{ij}$, where $x_i = \gamma_{ij}x^j, v_i = \gamma_{ij}v^j, v_i^g = \gamma_{ij}v_i^g$ and $f_i = \gamma_{ij}f^j$ and we have adopted the Einstein summation convention such that pairs of matching covariant, contravariant indices are implicitly summed over (see e.g. Hobson et al. 2006, for details).

The stopping time, in the Epstein regime, for a spherical dust grain of size $s$ and grain density $\rho_{\text{grain}}$ in a gas of density $\rho_g$ is

$$t_s = \frac{\rho_{\text{grain}} s}{\rho_g c_s} \sqrt{\frac{\pi \gamma}{8}}, \qquad (3.3)$$

where $\gamma$ is the adiabatic index of the gas and $c_s$ is the gas sound speed (Epstein 1924; Baines et al. 1965; Whipple 1972). The relative importance of gas drag is dictated by a comparison between the stopping time and some characteristic timescale of the fluid flow, $t_f$. This is encapsulated by the Stokes number $St = t_s/t_f$ which is a dimensionless number which controls how strongly the gas and dust are coupled. In rotating shear flows, with angular velocity $\Omega$, it is typical to take $t_f = \Omega^{-1}$ (although in some applications it can be useful to instead set $t_f$ to be the timescale associated with the fluid shear).

A commonly used model for the stochastic gas velocity, subject to homogeneous turbulence, is to model it as a Ornstein-Uhlenbeck process,

$$dv_i^g = -\frac{1}{t_c}v_i^g dt + \sqrt{\frac{2\alpha}{t_c}} c_s dW_i, \qquad (3.4)$$

where $t_c$ is the correlation time (or “eddy turnover” time) of the turbulence, $c_s$ is the gas sound speed, $\alpha$ is a dimensionless measure of the strength of the fluid turbulence and $W_i$ is a Wiener process. This model of turbulence regards the turbulent flow as a member of a statistical ensemble of similar flows (Thomson 1987), with each “draw” following a fluid element in a single realisation of the flow.

As with the stopping time it is useful to introduce a dimensionless correlation time $\tau_c = t_c/t_f$. Some authors define the Stokes number to be $St = t_s/t_c$, however this only really makes sense in homogeneous turbulence applications where $t_c$ is the only fluid timescale.

For more complex fluid flows, in the infinite Reynolds number limit, we can model turbulence as undergoing an Ornstein-Uhlenbeck walk about the mean flow. In this model the gas velocity evolves according to

$$dv_i^g = f_i^g dt - \frac{1}{t_c}(v_i^g - u_i^g)dt + \sqrt{\frac{2\alpha}{t_c}} c_s dW_i, \qquad (3.5)$$

where $f_i^g$ is the force per unit mass on the gas and $u_i^g = \mathbb{E}_g(v_i^g)$ is the mean gas velocity at the dust location. This mean gas velocity needs to be solved for separately, for which we use the Equations A 10-A 12 in Appendix A. In the absence of back reaction the force per unit mass on the gas is due to gravity and pressure gradients with $f_i^g = -\nabla_i \phi - \rho_g^{-1} \nabla_i p_g$ where $p_g$ is the gas pressure and $\rho_g$ is the gas density. With this choice of $f_i^g$, Equation 3.5 amounts to modelling the pressure fluctuation and dissipation terms as being responsible for the Ornstein-Uhlenbeck terms present above (Pope 2000). $f_i^g, \alpha, t_c, c_s$ and $\mathbf{u}^g$ are all functions of space and, in general, time. For instance, in accretion discs, $t_c$ is typically proportional to the orbital period and is thus an increasing function of cylindrical radius. Likewise, the
soundspeed and $\alpha$ vary (typically slowly) throughout the disc, although $\alpha$ is often assumed to be constant. All these quantities must be evaluated at the dust particle position. In principle one may be able to incorporate the effects of back reaction into $f_i^g$, we shall give a brief discussion of this possibility in Section 8.

Combining the model for the gas and dust, we arrive at a system of stochastic differential equations describing the motion of a dust grain in a turbulent gas,

$$dx_i = v_i dt, \qquad (3.6)$$

$$dv_i = f_i dt - \frac{1}{t_s}(v_i - v_i^g)dt, \qquad (3.7)$$

$$D_d v_i^g = f_i^g dt - \frac{1}{t_c}(v_i^g - u_i^g)dt + \sqrt{\frac{2\alpha}{t_c}} c_s dW_i. \qquad (3.8)$$

Now one can regards each ‘draw’ as selecting, and following, a single dust grain entrained with the turbulent flow. The gas fluid elements do not, in general, follow the dust grains, so we must correct for the fact we are taking a sample of the gas along the trajectory of the dust. Following Minier et al. (2004, 2014) we take the operator $D_d$ to be

$$D_d v_i^g = dv_i^g - (u^k - u_g^k) \nabla_k u_i^g dt. \qquad (3.9)$$

This can be thought of as a separate ‘advection’ step which, on average, corrects for the difference in the gas and dust trajectories. One can more compactly write these equations in terms of the dynamics of a particle in 6-dimensions, subject to drag, stochastic forcing and force per unit mass $F_\alpha$ (which contains contributions from the force on the dust and gas $f_i$, $f_i^g$, along with the shift correction, $(u^k - u_g^k)\nabla_k u_i^g$),

$$dX_\alpha = V_\alpha dt, \\ dV_\alpha = F_\alpha dt - C_{\alpha\beta} (V^\beta - U_g^\beta) dt + \sigma_{\alpha\beta} dW^\beta, \qquad (3.10)$$

where we have adopted the convention that Greek indices are over the 6-dimensional space and Latin indices are taken over the 3-dimensional space. These 6-dimensional indices are raised and lowered with a 6-dimensional metric tensor $g_{\alpha\beta}$, constructed from $\gamma_{ij}$, which will be properly defined in the next section. We have introduced $U_g^\beta$, the mean gas velocity “seen” by the dust; the 6×6 drag tensor $C_{\alpha\beta}$, which incorporate both the gas-dust drag on the stopping time along with the return of the stochastic gas velocity towards the mean on the turbulent correlation time which in the 6D picture acts like a “drag” between the gas components of the velocity and the mean gas flow. We have also introduced $\sigma_{\alpha\beta}$ which controls the strength of the stochastic forcing in each component of the momentum equation - i.e. it’s the 6D form of the last term in equation 3.8. In addition to simplifying the subsequent derivations, Equation 3.10 allows us to derive the fluid model for more general drag and turbulence models without increasing the complexity. For instance the subsequent derivations works equally well for anisotropic stochastic driving.

One can also include anisotropic correlation times as seen in some two-phase turbulence models (e.g. Minier et al. 2004, 2014), based on the analysis of Csanady (1963), which attempts to incorporate the effects of spatial correlation on the fluid seen by the dust particles. We have chosen not to include this correction as the proposed form of the correction in the literature (as described in Minier et al. 2004, 2014) predicts that rapidly drifting particle in rotating shear flows experience the same turbulence as that in homogeneous-isotropic turbulence. This likely arises due to the Csanady correction neglecting the anisotropy in the
correlation length induced by the shear. It is possible that the two-step stochastic model (as discussed in Minier & Henry 2023) will better account for the effects of spatial correlations and improve the modelling of dusty anisotropic turbulence in the future.

### 3.1. Geometry of the 6-dimensional space

The 3 additional dimensions in the 6-D system are a set of dummy gas degrees of freedom corresponding to the gas displacement. These should not be thought of as the gas position vector as the gas is coincident with the dust. These additional dimensions are, in a sense, non-physical and, in order that the 6D system agree with the 3D system, the 6D system must possess translational invariance along these dummy directions. The coordinate basis of the gas displacement are independent of the basis of the dust position vector. However it is useful to choose the basis of the gas displacement dimensions such that it reflects the underlying (physical) 3D coordinate system.

To construct this coordinate system we first consider the coordinates of the underlying 3D system with metric tensor $\gamma_{ij}$ and associated Christoffel symbols $\Gamma_{ij}^k$. Introducing basis vectors for the 6-D system, $\{\hat{e}_\alpha\}$, and the notation $\alpha_d \in \{1, 2, 3\}$ and $\alpha_g \in \{4, 5, 6\}$ such that $\hat{e}_{\alpha_d}$ give the basis vectors of the dust position vector and $\hat{e}_{\alpha_g}$ gives the basis vectors of the gas displacement vector. Additionally it is useful to introduce the bijection $*: \{1...6\} \to \{1...6\}$, which interchanges the 'dummy gas' and position indices with $1, 2, 3 \mapsto 4, 5, 6$ and $4, 5, 6 \mapsto 1, 2, 3$.

Throughout this work we shall make use of symmetrising/anti-symmetrising operations on the tensor indices with $E_{(\alpha_1 \cdots \alpha_n)}$ and $E_{[\alpha_1 \cdots \alpha_n]}$, for some tensor $\mathbf{E}$, being symmetrisation and anti-symmetrisation of the indices in brackets, where $E_{(\alpha\beta)} = \frac{1}{2}(E_{\alpha\beta} + E_{\beta\alpha})$ and $E_{[\alpha\beta]} = \frac{1}{2}(E_{\alpha\beta} - E_{\beta\alpha})$. The operation $*$ does not commute with symmetrisation/anti-symmetrisation operations, but instead follows the obvious order of operations such that

$$E_{(\alpha, \beta^*)} = \frac{1}{2}(E_{\alpha\beta^*} + E_{\beta^*\alpha}), \qquad (3.11)$$

$$E_{(\alpha, \beta)^*} = \frac{1}{2}(E_{\alpha\beta^*} + E_{\beta\alpha^*}), \qquad (3.12)$$

with equivalent expressions for antisymmetrisation.

The physical solutions must be independent of the gas displacement, we can therefore integrate out the dummy gas dimensions. Introducing an integral over the dummy gas directions,

$$\bar{\tau} := \int \cdot J_g d^3 x_g, \qquad (3.13)$$

where $J_g$ is the Jacobian determinant of the dummy gas coordinates. Thus for $\mathbf{E}$, an arbitrary tensorial quantity, we have

$$\overline{\nabla}_{\alpha_g} \mathbf{E} = 0. \qquad (3.14)$$

For Cartesian gas displacement coordinates this integrating out of the non-physical space is straight-forward. Unfortunately if the coordinate system describing the dust position is non-Cartesian then we need to rotate the “dummy” components of vectors so that they reflect the underlying 3D coordinate system (e.g. when calculating the gas drag). It is instead useful to setup the geometry of our 6D space so that the rotation happens automatically. To do this we first introduce the metric tensor of the 6D coordinate system
$$ g_{\alpha\beta} = \begin{cases} \gamma_{\alpha\beta}, & \alpha, \beta \in \{1, 2, 3\}, \\ \gamma_{\alpha^*\beta^*}, & \alpha, \beta \in \{4, 5, 6\}, \\ 0, & \text{otherwise.} \end{cases} \qquad (3.15) $$

We also introduce a metric connection $\overline{\nabla}_{\alpha}$ that is responsible for rotating the dummy gas coordinate system. We require that this connection satisfy the following properties:

(i) $\overline{\nabla}_{\alpha}$ is a metric connection, so that $\overline{\nabla}_{\alpha}g_{\beta\gamma} = \overline{\nabla}_{\alpha}g^{\beta\gamma} = 0$.

(ii) Translational invariance with respect to the gas displacement such that $\overline{\nabla}_{\alpha_g} \mathbf{E}(\mathbf{x}_d) = 0$ for tensoral quantity $\mathbf{E}$.

(iii) Alignment of the dummy gas coordinates with the position coordinates. For vectors $\mathbf{A}$, $\mathbf{B}$ and $\tilde{\mathbf{B}}$ with $B^{\alpha_d} = 0$ and $\tilde{B}^{\alpha} = B^{\alpha^*}$ then we require $(\mathbf{A} \cdot \overline{\nabla}\mathbf{B})^{\beta} = (\mathbf{A} \cdot \overline{\nabla}\tilde{\mathbf{B}})^{\beta^*}$.

Property (ii) ensures that $\overline{\nabla}_{\alpha_g} \mathbf{E}(\mathbf{x}_d) = \overline{\nabla}_{\alpha_g} \mathbf{E}$, where $\nabla_i$ is the covariant derivative, and allows us to carry out the integral over the dummy gas directions by replacing covariant derivatives with $\overline{\nabla}_i$. Property (iii) is required to ensure that the geometric terms in Lagrangian time derivatives act the same on the dust and gas components of the 6D vectors. This can be seen considering $\mathbf{A} = \mathbf{U}$ and considering the action of the Lagrangian time derivative, $\mathbf{D} = \partial_t + \mathbf{U} \cdot \overline{\nabla}$, on the vectors $\mathbf{B}$ and $\tilde{\mathbf{B}}$. As $B^{\alpha} = \tilde{B}^{\alpha^*}$ one requires $(\mathbf{D}\mathbf{B})^{\alpha} = (\mathbf{D}\tilde{\mathbf{B}})^{\alpha^*}$ this requires condition (iii) as $\mathbf{U}$ is arbitrary and $(\partial_t \mathbf{B})^{\alpha} = (\partial_t \tilde{\mathbf{B}})^{\alpha^*}$.

The connection which satisfies these properties, given the metric tensor (3.15), acts on the basis vectors $\{\hat{e}_{\alpha}\}$ as follows,

$$ \overline{\nabla}_{\alpha_d} \hat{e}_{\beta_d} = \mathcal{T}_{\alpha_d \beta_d}^{\gamma_d} \hat{e}_{\gamma_d}, \qquad \overline{\nabla}_{\alpha_d} \hat{e}_{\beta_g} = \mathcal{T}_{\alpha_d \beta_g^*}^{\gamma_g} \hat{e}_{\gamma_g}, \qquad (3.16) $$

with $\overline{\nabla}_{\alpha_g} \hat{e}_{\beta} = 0$. As $\mathcal{T}_{ij}^k$ are the Christoffel symbols components for the 3D coordinate system associated with the metric $\gamma_{ij}$ it is straight-forward to show that this connection satisfies property 1. Property 2 follows from $\overline{\nabla}_{\alpha_g} \hat{e}_{\beta} = 0$. Finally for property 3,

$$ (\mathbf{A} \cdot \overline{\nabla} \mathbf{B})^{\beta} = A^{\alpha} \overline{\nabla}_{\alpha} B^{\beta} = A^{\alpha} \partial_{\alpha} B^{\beta} + A^{\alpha} \mathcal{T}_{\alpha\gamma}^{\beta} B^{\gamma}, \qquad (3.17) $$

$$ \begin{aligned} (\mathbf{A} \cdot \overline{\nabla} \tilde{\mathbf{B}})^{\beta} &= A^{\alpha} \overline{\nabla}_{\alpha} \tilde{B}^{\beta} = A^{\alpha} \partial_{\alpha} \tilde{B}^{\beta} + A^{\alpha} \mathcal{T}_{\alpha\gamma}^{\beta*} \tilde{B}^{\gamma} \\ &= A^{\alpha} \partial_{\alpha} B^{\beta*} + A^{\alpha} \mathcal{T}_{\alpha\gamma*}^{\beta*} B^{\gamma*} = (\mathbf{A} \cdot \overline{\nabla} \mathbf{B})^{\beta*}. \end{aligned} \qquad (3.18) $$

While this connection has the advantage of keeping the gas coordinate system aligned and avoids the necessity of including rotation matrices in the equation of motion it does have one major drawback in that it is not torsion free (since it isn't the Levi-Civita connection). This torsion arises when $\overline{\nabla}_{\alpha_d} \hat{e}_{\beta_g} \neq 0$ as $\overline{\nabla}_{\alpha_t} \hat{e}_{\beta_d} = 0$, by construction, and is associated with the rotation of the dummy gas coordinate system. The torsion tensor, $S_{\alpha\beta}^{\gamma}$, is given by

$$ S_{\alpha\beta}^{\gamma} \hat{e}_{\gamma} = \overline{\nabla}_{\alpha} \hat{e}_{\beta} - \overline{\nabla}_{\beta} \hat{e}_{\alpha}, \qquad (3.19) $$

making use of the properties of the connection the torsion tensor components are

$$ S_{\alpha_d \beta_g}^{\gamma_g} = -S_{\beta_g \alpha_d}^{\gamma_g} = \mathcal{T}_{\alpha_d \beta_g^*}^{\gamma_g}, \qquad (3.20) $$

with all other components zero.
Finally, after specialising to the oriented 6-D geometry one can write $U_g^\alpha$ in terms of the mean gas velocity in the gas frame, $u_g^i$,

$$U_g^\alpha = \begin{cases} u_g^\alpha, & \alpha \in \{1, 2, 3\}, \\ u_g^{\alpha*}, & \alpha \in \{4, 5, 6\}, \end{cases} \qquad (3.21)$$

while the drag and diffusion tensors can be written in terms of the metric tensor. The 6-D force per unit mass is

$$F_{\alpha} = \begin{cases} f_{\alpha}, & \beta \in \{1, 2, 3\}, \\ f_{\alpha*}^{g} + (U^{\beta} - U_{g}^{\beta}) \nabla_{\beta} U_{\alpha}^{g}, & \alpha \in \{4, 5, 6\}, \end{cases} \qquad (3.22)$$

while the drag tensor is

$$C_{\alpha\beta} = \begin{cases} \frac{1}{t_s} g_{\alpha\beta}, & \alpha, \beta \in \{1, 2, 3\}, \\ -\frac{1}{t_s} g_{\alpha\beta*}, & \alpha \in \{1, 2, 3\}, \quad \beta \in \{4, 5, 6\}, \\ 0, & \alpha \in \{4, 5, 6\}, \quad \beta \in \{1, 2, 3\}, \\ \frac{1}{t_c} g_{\alpha\beta}, & \alpha, \beta \in \{4, 5, 6\}, \end{cases} \qquad (3.23)$$

while the diffusion tensor is

$$D_{\alpha\beta} = \begin{cases} \frac{\alpha c_s^2}{t_c} g_{\alpha\beta}, & \alpha, \beta \in \{4, 5, 6\}, \\ 0, & \text{otherwise.} \end{cases} \qquad (3.24)$$

This diffusion tensor is applicable to isotropic diffusivity. More generally one can include an anisotropic diffusivity by introducing an $\alpha$-tensor $a_{\alpha\beta}$, in which case the diffusion tensor will be

$$D_{\alpha\beta} = \begin{cases} \frac{c_s^2}{t_c} a_{\alpha\beta}, & \alpha, \beta \in \{4, 5, 6\}, \\ 0, & \text{otherwise.} \end{cases} \qquad (3.25)$$

If one were to instead use the more usual Levi-Civita connection the above expressions would be considerably more complex as they would need to include the rotation of the dummy gas directions.

# 4. Derivation of the dust-fluid model.

## 4.1. Derivation of the Fokker-Planck equation

In order to derive the dust fluid model we must first obtain the Fokker-Planck equation associate with Equation 3.10, and then perform a moment expansion to derive the fluid model. To do this consider an arbitrary ($C^2$) function of the dust particle position, velocity and stochastic gas displacement, $A = A(X, V)$. By use of Ito's lemma this evolves according to

$$dA = \frac{\partial A}{\partial X_{\alpha}} dX_{\alpha} + \frac{\partial A}{\partial V_{\alpha}} dV_{\alpha} + \frac{1}{2} \frac{\partial^2 A}{\partial V_{\alpha} \partial V_{\beta}} \langle dV_{\alpha}, dV_{\beta} \rangle, \qquad (4.1)$$

where the angle bracket $\langle \cdot, \cdot \rangle$ denotes the covariance. Where the covariance of a Wiener process $dW_\alpha$ is given by
$$ \langle dW^\alpha, dW^\beta \rangle = g^{\alpha\beta} dt. \quad (4.2) $$

This leads to the following for the covariance of the velocity,

$$ \langle dV_\alpha, dV_\beta \rangle = \sigma_{\alpha\mu}\sigma_{\beta\nu}\langle dW^\mu dW^\nu \rangle = 2D_{\alpha\beta}dt, \quad (4.3) $$

where we have introduced the diffusion tensor, $D_{\alpha\beta} = \frac{1}{2}g^{\mu\nu}\sigma_{\alpha\mu}\sigma_{\beta\nu}$. Substituting Equation 3.10 into Equation 4.1, we arrive at

$$ \begin{aligned} dA = \frac{\partial A}{\partial X_\alpha}V_\alpha dt + \frac{\partial A}{\partial V_\alpha} \left[ F_\alpha - C_{\alpha\beta} (U^\beta - U_g^\beta) \right] dt \\ + \frac{\partial A}{\partial V_\alpha} \sigma_{\alpha\beta} dW^\beta + D_{\alpha\beta} \frac{\partial^2 A}{\partial V_\alpha \partial V_\beta} dt. \end{aligned} \quad (4.4) $$

The expectation of A is given by,

$$ \mathbb{E}[A] = \int p^L(\mathbf{X}, \mathbf{V}, t, \mathbf{X}_0, \mathbf{V}_0, t_0) A(\mathbf{X}, \mathbf{V}) d^6\mathbf{X} d^6\mathbf{V}, \quad (4.5) $$

where $p^L(\mathbf{X}, \mathbf{V}, t, \mathbf{X}_0, \mathbf{V}_0, t_0)$ is the probability for the system to arrive at state $(\mathbf{X}, \mathbf{V}, t)$ from an initial state $(\mathbf{X}_0, \mathbf{V}_0, t_0)$. $\mathbb{E}[dA]$ is given by,

$$ \mathbb{E}[dA] = \int dp^L(\mathbf{X}, \mathbf{V}, t, \mathbf{X}_0, \mathbf{V}_0, t_0) A(\mathbf{X}, \mathbf{V}) d^6\mathbf{X} d^6\mathbf{V}. \quad (4.6) $$

Substituting Equation 4.4 into the above and, after appropriate integration by parts (assuming appropriate regularity conditions for $p$, namely that $p$ and $\frac{\partial p}{\partial V_\alpha}$ vanish as $V^\beta \to \infty$) we arrive at

$$ \begin{aligned} \int A \left\{ dp^L + \frac{\partial}{\partial X_\alpha} (p^L V_\alpha) dt + \frac{\partial}{\partial V_\alpha} \left[ \left(F_\alpha - C_{\alpha\beta} (U^\beta - U_g^\beta)\right) p^L \right] dt \right. \\ \left. - D_{\alpha\beta} \frac{\partial^2 p^L}{\partial V_\alpha \partial V_\beta} dt \right\} d^6 \mathbf{X} d^6 \mathbf{V}, \end{aligned} \quad (4.7) $$

provided that $\int_{\partial} p^L A \mathbf{V} \cdot dS d^6 \mathbf{V} = 0$, where $\int_{\partial} dS$ denotes an integral over the spatial boundaries. i.e. the expected net flux of A through the domain boundaries is zero.

As A is arbitrary (baring being $C^2$, and the boundary conditions) we arrive at the Fokker-Planck equation for $p$,

$$ \frac{\partial p^L}{\partial t} + \frac{\partial}{\partial X^\alpha} (p^L V^\alpha) + \frac{\partial}{\partial V_\alpha} \left[ \left( F_\alpha - C_{\alpha\beta} (V^\beta - U_g^\beta) \right) p^L \right] = D_{\alpha\beta} \frac{\partial^2 p^L}{\partial V_\alpha \partial V_\beta}. \quad (4.8) $$

This equation gives an evolutionary equation for the Lagrangian transition PDF, describing the probability of finding a particle at $\mathbf{X}, \mathbf{V}$ at time $t$ conditional on it being located at $\mathbf{X}_0, \mathbf{V}_0$ at time $t_0$. The fluid model will consist of a set of Eulerian fields located at a given position in Space and must be obtained from the Eulerian Mass Density Function (MDF), $p(\mathbf{X}, \mathbf{V}, t)$, which is the expected mass density of particles at $\mathbf{X}, \mathbf{V}$ at time $t$ (Pope 1985, 2000; Minier & Peirano 2001). This will contain contributions from particles with different initial condition $(\mathbf{X}_0, \mathbf{V}_0)$, arriving from differing trajectories. This can be obtained from the Eulerian mass
density function at $t_0$, $p(\mathbf{X}_0, \mathbf{V}_0, t_0)$, by using the transition PDF and integrating over the initial positions and velocities (Pope 1985, 2000; Minier & Peirano 2001),

$$p(\mathbf{X}, \mathbf{V}, t) = \int p^L(\mathbf{X}, \mathbf{V}, t, \mathbf{X}_0, \mathbf{V}_0, t_0) p(\mathbf{X}_0, \mathbf{V}_0, t_0) d^6\mathbf{X}_0 d^6\mathbf{V}. \quad (4.9)$$

We can obtain the Fokker-Planck equation for $p$ by multiplying Equation 4.8 by $p(\mathbf{X}_0, \mathbf{V}_0, t_0)$ and integrating over the initial position and velocities. This leaves the form of the Fokker-Planck Equation unchanged and we obtain

$$\frac{\partial p}{\partial t} + \frac{\partial}{\partial X^{\alpha}}(pV^{\alpha}) + \frac{\partial}{\partial V_{\alpha}}\left[(F_{\alpha} - C_{\alpha\beta}(V^{\beta} - U_{g}^{\beta}))p\right] = D_{\alpha\beta}\frac{\partial^2 p}{\partial V_{\alpha} \partial V_{\beta}}. \quad (4.10)$$

## 4.2. Moment expansion of the Fokker-Planck equation

Fluid dynamical models can be derived from the Fokker-Planck equation via a moment expansion, in a similar manner to that done in kinetic theory. In performing this moment expansion we wish to arrive at a set of PDEs in space and time from the initial PDE in $(t, \mathbf{X}, \mathbf{V})$. This means we need to compute a moment expansion in $\mathbf{V}$. A similar procedure was carried out by Youdin & Lithwick (2007). Defining the velocity moments of $p$ as follows,

$$\rho_6 := \int p d^6 \mathbf{V}, \quad (4.11)$$

$$\rho_6 U_{\alpha} := \int p V_{\alpha} d^6 \mathbf{V}, \quad (4.12)$$

$$\Pi_{\alpha\beta} := \int p(V_{\alpha} - U_{\alpha})(V_{\beta} - U_{\beta})d^6\mathbf{V}, \quad (4.13)$$

$$\Pi_{\alpha_1 \dots \alpha_k} := \int p(V_{\alpha_1} - U_{\alpha_1}) \dots (V_{\alpha_k} - U_{\alpha_k}) d^6 \mathbf{V}. \quad (4.14)$$

Note that this moment expansion is in the 6-dimensional space, so that $\rho_6$ is the 6-dimensional mass density and $\Pi_{\alpha\beta}$ is the 6-dimensional rheological stress tensor† and we have chosen a normalisation such that

$$\int p d^6 \mathbf{V} d^3 \mathbf{x}_g = \rho_d, \quad (4.15)$$

where $\rho_d$ is the dust density (i.e. the density of the dust phase, this is equal to the grain density, $\rho_{\text{grain}}$, multiplied by the dust volume fraction). We have opted to normalise with respect to the dust mass density rather than the dust number density so that $\int \Pi_{\alpha\beta} d^3 \mathbf{x}_g$ has the same units as the gas pressure.

Taking the zeroth velocity moment of Equation 4.10 we arrive at the (6-D) dust continuity equation,

$$\dot{\rho}_6 + \nabla_{\alpha} [\rho_6 U^{\alpha}] = 0. \quad (4.16)$$

<footnote>† We have chosen to call the second velocity moment the rheological stress tensor rather than the dust pressure tensor as it contains contributions from both the dust pressure (particle velocity dispersion) and the dust Reynolds stress. These two stresses are indistinguishable due to the way we have formulated the averaging. This can run into issues when dust-dust collisions are included as the dust collisional velocity is principally sensitive to the particle (rather than turbulent) velocity dispersion (Fox 2014; Capecelatro et al. 2016b).</footnote>
The first V moment of Equation 4.10 leads to the (6-D) dust momentum equation,

$$ \frac{\partial}{\partial t} [\rho_6 U_\alpha] + \nabla^\beta [\Pi_{\alpha\beta} + \rho_6 U_\alpha U_\beta] - \rho_6 F_\alpha + \rho_6 C_{\alpha\beta} (U^\beta - U_g^\beta) = 0. \quad (4.17) $$

Taking the second V moment yields a constitutive relation for the (6-D) dust stress tensor,

$$ \frac{\partial}{\partial t} [\Pi_{\alpha\beta} + \rho_6 U_\alpha U_\beta] + \nabla^\gamma [\Pi_{\alpha\beta\gamma} + 3U_{(\alpha}\Pi_{\beta\gamma)} + \rho_6 U_\alpha U_\beta U_\gamma] \\ - 2\rho_6 U_{(\alpha} [F_\beta) - C_{\beta)\gamma} (U^\gamma - U_g^\gamma)] + 2\Pi^\gamma_{(\alpha} C_{\beta)\gamma} = 2\rho_6 D_{(\alpha\beta)}. \quad (4.18) $$

Here we have made use of the notation for the symmetrisation of the tensor indices. As we shall make extensive use of this notation we give explicit expressions for the symmetrised terms in the above equation as an illustrative example, $U_{(\alpha}\Pi_{\beta\gamma)} = \frac{1}{3}(U_\alpha\Pi_{\beta\gamma} + U_\beta\Pi_{\alpha\gamma} + U_\gamma\Pi_{\alpha\beta})$ and $2U_{(\alpha} [F_\beta) - C_{\beta)\gamma} (U^\gamma - U_g^\gamma)] = U_\alpha [F_\beta - C_{\beta)\gamma} (U^\gamma - U_g^\gamma)] + U_\beta [F_\alpha - C_{\alpha\gamma} (U^\gamma - U_g^\gamma)]$.

Higher velocity moments can be computed in a similar manner. Making use of the expressions for the velocity moments of the terms of the Fokker-Planck equation given in Appendix C, we can take the k-th velocity moment of the Fokker-Planck equation to obtain

$$ \frac{\partial \Pi_{\alpha_1 \cdots \alpha_k}}{\partial t} + k \Pi_{(\alpha_1 \cdots \alpha_{k-1}} \left[ D U_{\alpha_k} + \nabla_{\alpha_k} \phi + C_{\alpha_k}^{\gamma} (U_{\gamma} - U_{\gamma}^{g}) \right] + \nabla^{\sigma} \left[ \Pi_{\alpha_1 \cdots \alpha_k} \sigma + U_{\sigma} \Pi_{\alpha_1 \cdots \alpha_k} \right] \\ + k \Pi_{\sigma(\alpha_1 \cdots \alpha_{k-1}} \left[ \nabla^{\sigma} U_{\alpha_k} + C_{\alpha_k}^{\sigma} \right] = k(k-1) \Pi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_k} \quad (4.19) $$

Making use of the dust momentum equation this simplifies to

$$ (D + \nabla_{\sigma} U^{\sigma}) \Pi_{\alpha_1 \cdots \alpha_k} + k \Pi_{\sigma(\alpha_1 \cdots \alpha_{k-1}} \nabla^{\sigma} U_{\alpha_k)} + \nabla^{\sigma} \Pi_{\alpha_1 \cdots \alpha_k} \sigma \\ = -k \left[ \Pi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} C_{\alpha_k)} \sigma - \frac{1}{\rho} \Pi_{(\alpha_1 \cdots \alpha_{k-1}} \nabla^{\sigma} \Pi_{\alpha_k)} \sigma - (k-1) \Pi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_k)} \right]. \quad (4.20) $$

Taking $k=2$ in the above equation we recover the constitutive relation for $\Pi_{\alpha\beta}$ (to obtain this we note that $\Pi_\alpha = 0$ by the definition of $U^\alpha$).

It is useful to define various tensor advection operators $\mathcal{D}$, $\mathcal{D}_1$ and $\mathcal{D}_2$. When acting on the k-th velocity moment these are given by

$$ \mathcal{D} \Pi_{\alpha_1 \cdots \alpha_k} = D \Pi_{\alpha_1 \cdots \alpha_k} + k \Pi_{\gamma(\alpha_1 \cdots \alpha_{k-1}} \nabla_{\alpha_k)} U^{\gamma} + \Pi_{\alpha_1 \cdots \alpha_k} \nabla_{\gamma} U^{\gamma}, \quad (4.21) $$

$$ \mathcal{D}_1 \Pi_{\alpha_1 \cdots \alpha_k} = D \Pi_{\alpha_1 \cdots \alpha_k} + \Pi_{\alpha_1 \cdots \alpha_k} \nabla_{\gamma} U^{\gamma}, \quad (4.22) $$

$$ \mathcal{D}_2 \Pi_{\alpha_1 \cdots \alpha_k} = D \Pi_{\alpha_1 \cdots \alpha_k} + k \Pi_{\gamma(\alpha_1 \cdots \alpha_{k-1}} \nabla^{\gamma} U_{\alpha_k)} + \Pi_{\alpha_1 \cdots \alpha_k} \nabla_{\gamma} U^{\gamma}. \quad (4.23) $$

The first of these is closely related to the convective Maxwell derivative, with $\mathcal{D}\Pi_{\alpha_1 \cdots \alpha_k=0}$ implying that the tensorial quantity $\rho_6^{-1}\Pi_{\alpha_1 \cdots \alpha_k=0}$ (i.e. the k-th velocity correlation) is passively advective by the flow. The other operators $\mathcal{D}_1$ and $\mathcal{D}_2$ are defined for convenience. This highlights one advantage of the 6-d formalisation as couplings between the dust kinetic tensor ($T_{\alpha_d\beta_d}$), cross correlation tensor ($T_{\alpha_d\beta_g}$) and fluid seen Reynolds stress ($R_{\alpha_g\beta_g}$), are shown to arise from the advection of the dust Rheological stress by the 6D flow.

Rearranging the continuity, momentum and constitutive equations, and making use of the operator $\mathcal{D}_2$, we obtain
$$D\rho_6 = -\rho_6\nabla_\alpha u^\alpha, \qquad (4.24)$$

$$\rho_6 DU_\alpha = \rho_6 F_\alpha - \nabla^\beta \Pi_{\alpha\beta} - \rho_6 C_{\alpha\beta} (U^\beta - U_g^\beta), \qquad (4.25)$$

$$\mathcal{D}_2 \Pi_{\alpha\beta} = -\nabla^\gamma \Pi_{\alpha\beta\gamma} - 2 \left( \Pi^\gamma_{(\alpha} C_{\beta)\gamma} - \rho_6 D_{(\alpha\beta)} \right). \qquad (4.26)$$

As the right hand side of Equation 4.26 is symmetrised this ensures that $\Pi_{\alpha\beta}$ remains symmetric for symmetric initial conditions. Using a similar argument to that advanced in Ogilvie (2003); Lynch & Ogilvie (2021), $\Pi_{\alpha\beta}$ is positive semi-definite for positive semi-definite initial conditions (See Appendix B.1 for a details). The evolutionary equation for the k-th velocity moment simplifies to

$$\mathcal{D}_2 \Pi_{\alpha_1 \dots \alpha_k} = - \nabla^\gamma \Pi_{\gamma \alpha_1 \dots \alpha_k} - k \left[ \Pi^\sigma_{(\alpha_1 \dots \alpha_{k-1}} C_{\alpha_k)\sigma} - \frac{1}{\rho} \Pi_{(\alpha_1 \dots \alpha_{k-1}} \nabla^\sigma \Pi_{\alpha_k)\sigma} - (k-1) \Pi_{(\alpha_1 \dots \alpha_{k-2}} D_{\alpha_{k-1})\sigma} \right] \qquad (4.27)$$

Alternatively one can write the constitutive equation in terms of the operator $\mathcal{D}$ and obtain the following alternative form of Equation 4.26,

$$\mathcal{D} \Pi_{\alpha\beta} = -2 \left( \Pi^\gamma_{(\alpha} A_{\beta)\gamma} - \rho_6 D_{(\alpha\beta)} \right), \qquad (4.28)$$

where we have defined

$$A_{\alpha\beta} = C_{\alpha\beta} - \omega^\gamma \varepsilon_{\alpha\beta\gamma}, \qquad (4.29)$$

where $\omega^\gamma$ is the dust-fluid vorticity and

$$\omega^\gamma \varepsilon_{\gamma\alpha\beta} = \nabla_\alpha U_\beta - \nabla_\beta U_\alpha. \qquad (4.30)$$

The evolutionary equation for the k-th velocity moment can be similarly rewritten. In the full 6-D model, with the Levi-Civita connection, Equation 4.28 is the more useful form of the constitutive relation as it is independent of the Christoffel symbol components (by symmetry) and it is more connected to the underlying physics of the rheological stress tensor where the operator $\mathcal{D}$ is responsible for passively advecting the pressure tensor and the drag, vorticity and turbulent “heating” on the right hand side of Equation 4.28 act like sources/sinks for the stress tensor. Unfortunately in the presence of torsion, the constitutive equation based on Equation 4.28 ends up more complicated to manipulate than that based on Equation 4.26 owing to the addition of terms involving the torsion tensor. As such we shall stick to Equation 4.26 for the constitutive relation from this point onwards.

Finally, for the purposes of developing the closure scheme for the moment expansion, it is useful to express the evolutionary equation for the k-th velocity moment in terms of the operator $\mathcal{D}_1$,

$$\mathcal{D}_1 \Pi_{\alpha_1 \dots \alpha_k} = - \nabla^\sigma \Pi_{\alpha_1 \dots \alpha_k \sigma} - k \left[ \Pi^\sigma_{(\alpha_1 \dots \alpha_{k-1}} B_{\alpha_k)\sigma} \right. \\ \left. - \frac{1}{\rho} \Pi_{(\alpha_1 \dots \alpha_{k-1}} \nabla^\sigma \Pi_{\alpha_k)\sigma} - (k-1) \Pi_{(\alpha_1 \dots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_k)} \right], \qquad (4.31)$$

where we have introduced $B_{\alpha\beta} = C_{\alpha\beta} + \nabla_\beta U_\alpha$.
### 4.3. Closure scheme.

As is usual for a moment expansion we now have an infinite tower of velocity moments that is not useful for practical computations and must now consider a closure scheme. In this section we shall show that when the fluid is thermally stable, and the turbulent velocity small relative the fluid velocity, the third velocity moment typically decays until it is asymptotically small relative to the stress tensor and we can therefore drop the $\nabla^\gamma \Pi_{\alpha\beta\gamma}$ in the constitutive relation and close the moment expansion at the 2nd velocity moment.

#### 4.3.1. Well coupled ordering scheme.

Previous authors have noted that when the dust is well coupled to the gas ($\text{St} \ll 1$) it can be approximated with a fluid description. We can consider such a 'well coupled' ordering scheme by introducing a small parameter $\epsilon > 0$, which can be regarded as a characteristic Stokes number such that $\text{St} = O(\epsilon)$. We consider units such that $U^\alpha = O(1)$, $D_1 = O(1)$ and sufficiently weak turbulence heating such that $D_{\alpha\beta} = O(\epsilon^2)$. In our units the spatial gradients are limited such that $\nabla^\sigma = O(\epsilon^{-1})$ (In that the magnitude of the spatial gradients cannot significantly exceed $\epsilon^{-1}$, they can, however, be $\ll \epsilon^{-1}$).

Introducing rescaled velocity moment, $\tilde{\Pi}_{\alpha_1 \cdots \alpha_k}$, such that $\Pi_{\alpha_1 \cdots \alpha_k} = \epsilon^{\delta_k} \tilde{\Pi}_{\alpha_1 \cdots \alpha_k}$, and stretched/rescaled variable $\tilde{X} = X/\epsilon$, such that $\nabla = \epsilon^{-1} \tilde{\nabla}$, then we arrive at a rescaled equation for the k-th velocity moment

$$ \epsilon^{\delta_k} \mathcal{D}_1 \Pi_{\alpha_1 \cdots \alpha_k} = -\epsilon^{\delta_{k+1}-1} \nabla^\sigma \Pi_{\alpha_1 \cdots \alpha_k} \sigma - k \left[ \epsilon^{\delta_k-1} \Pi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1})} B_{\alpha_k} \sigma - \frac{\epsilon^{\delta_{k-1}+\delta_2-1}}{\rho} \Pi_{(\alpha_1 \cdots \alpha_{k-1})} \nabla^\sigma \Pi_{\alpha_k} \sigma - (k-1) \epsilon^{\delta_{k-2}+2} \Pi_{(\alpha_1 \cdots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k} \right]. \quad (4.32) $$

Proposing $\delta_k = 3\text{ceil}(k/2)$, we can rearrange the above to obtain, for even $k$,

$$ \epsilon \mathcal{D}_1 \Pi_{\alpha_1 \cdots \alpha_k} + k \left[ \Pi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1})} B_{\alpha_k} \sigma - (k-1) \Pi_{(\alpha_1 \cdots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k} \right] \\ = \epsilon^3 \left[ -\nabla^{\sigma} \Pi_{\alpha_1 \cdots \alpha_k} \sigma + \frac{k}{\rho} \Pi_{(\alpha_1 \cdots \alpha_{k-1})} \nabla^{\sigma} \Pi_{\alpha_k} \sigma \right], \quad (4.33) $$

For $k=2$ the left hand side corresponds to the constitutive model with $\Pi_{\alpha\beta\gamma} = 0$. For odd $k$ we instead have

$$ \epsilon \mathcal{D}_1 \Pi_{\alpha_1 \cdots \alpha_k} = -\nabla^\sigma \Pi_{\alpha_1 \cdots \alpha_k} \sigma - k \left[ \Pi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1})} B_{\alpha_k} \sigma - \frac{1}{\rho} \Pi_{(\alpha_1 \cdots \alpha_{k-1})} \nabla^\sigma \Pi_{\alpha_k} \sigma - (k-1) \Pi_{(\alpha_1 \cdots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k} \right]. \quad (4.34) $$

Thus we find that the correction to the evolutionary equation for the second velocity moment $\Pi_{\alpha\beta}$ is suppressed by a factor of $\epsilon^3$, relative to the leading order terms. Crucially this strong suppression means that Stokes numbers slightly less than one may still be well approximated by our fluid model, provided that we retain the $O(\epsilon)$ advection term ($\mathcal{D}_1 \Pi_{\alpha\beta}$) which will no longer be negligible.

According to Equation 4.34 the evolution of the third velocity moment will depend on
gradients of the 4-th velocity moment at leading order. Thus, we gain no advantages if we were to truncate the expansion at the third velocity moment over truncating at the second.

### 4.3.2. Near Maxwellian ordering scheme

We now wish to consider a situation where the dust distribution function is initially close to a Maxwellian velocity distribution and determine under what circumstances the departure from a Maxwellian velocity distribution remains small. Consider an asymmetric Maxwellian velocity distribution,

$$f = \frac{|A|^{1/2}}{(2\pi)^{n/2}} \exp\left(-\frac{1}{2}Q^{\alpha\beta}(V_{\alpha} - U_{\alpha})(V_{\beta} - U_{\beta})\right), \quad (4.35)$$

with second velocity moment

$$W_{\alpha\beta} = \int (V_{\alpha} - U_{\alpha})(V_{\beta} - U_{\beta})fd^{n}V. \quad (4.36)$$

This is related to $Q^{\alpha\beta}$ through $Q^{\alpha\sigma}W_{\sigma\beta} = \delta_{\beta}^{\alpha}$. More generally we define the k-th velocity moment for the Maxwellian velocity distribution as

$$W_{\alpha_1 \dots \alpha_k} = \int (V_{\alpha_1} - U_{\alpha_1}) \dots (V_{\alpha_k} - U_{\alpha_k}) f d^n V. \quad (4.37)$$

For odd $k$, $W_{\alpha_1 \dots \alpha_k} = 0$. Using standard results for Maxwellian distributions (e.g. Withers 1985) we obtain the following relationship between the k-th and $(k-2)$-th velocity moment;

$$W_{\alpha_1 \dots \alpha_k} = (k-1)W_{(\alpha_1 \dots \alpha_{k-2})} W_{\alpha_{k-1} \alpha_k}. \quad (4.38)$$

By symmetry of the velocity moments we also have $W_{\alpha_1 \dots \alpha_k} = W_{(\alpha_1 \dots \alpha_k)} = (k-1)W_{(\alpha_1 \dots \alpha_{k-2})} W_{\alpha_{k-1} \alpha_k}$.

Starting from the assumption that the 2nd velocity moment evolves according to

$$DW_{\alpha_1 \alpha_2} = -2 \left[ W^{\sigma}_{(\alpha_1} B_{\alpha_2)} \sigma - D_{\alpha_1 \alpha_2} \right], \quad (4.39)$$

we wish to show that the k-th velocity moment evolves according to

$$DW_{\alpha_1 \dots \alpha_k} = -k \left[ W^{\sigma}_{(\alpha_1 \dots \alpha_{k-1})} B_{\alpha_k)} \sigma - (k-1)W_{(\alpha_1 \dots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k)} \right] \quad (4.40)$$

Assuming this is the case for the $(k-2)$-th velocity moment then we can substitute Equation 4.38 into the above equation to obtain

$$ \begin{aligned} DW_{\alpha_1 \dots \alpha_k} &= (k-1)W_{(\alpha_1 \alpha_2)} DW_{\alpha_3 \dots \alpha_k} + (k-1)W_{(\alpha_1 \dots \alpha_{k-2})} DW_{\alpha_{k-1} \alpha_k} \\ &= -(k-1)(k-2)W_{(\alpha_1 \alpha_2)} \left[ W^{\sigma}_{\alpha_3 \dots \alpha_{k-1}} B_{\alpha_k)} \sigma - (k-3)W_{\alpha_3 \dots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_k)} \right] \\ &\quad - 2(k-1)W_{(\alpha_1 \dots \alpha_{k-2})} \left[ W^{\sigma}_{\alpha_{k-1} \alpha_k)} \sigma - D_{\alpha_{k-1} \alpha_k)} \right] \\ &= -(k-1) \left[ (k-2)W_{\sigma(\alpha_1 \dots \alpha_{k-3})} B^{\sigma}_{\alpha_{k-2}} W_{\alpha_{k-1} \alpha_k)} + 2W_{(\alpha_1 \dots \alpha_{k-2})} W^{\sigma}_{\alpha_{k-1} \alpha_k)} B_{\alpha_k)} \sigma \right] \\ &\quad + (k-1)(k-3) \left[ (k-2)W_{(\alpha_1 \dots \alpha_{k-4})} D_{\alpha_{k-3} \alpha_{k-2}} W_{\alpha_{k-1} \alpha_k)} + 2W_{(\alpha_1 \dots \alpha_{k-4})} W_{\alpha_{k-3} \alpha_{k-2}} D_{\alpha_{k-1} \alpha_k)} \right] \\ &= -k \left[ W^{\sigma}_{(\alpha_1 \dots \alpha_{k-1})} B_{\alpha_k)} \sigma - (k-1)W_{(\alpha_1 \dots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k)} \right]. \end{aligned} \quad (4.41) $$
Thus we see that if the (k-2)-th velocity moment evolves according to Equation 4.40, and the 2nd velocity moment evolves according to Equation 4.39, then the k-th velocity moment also evolves according to Equation 4.40. Starting with the 4-th velocity moment we see that, given Equations 4.39 and 4.41, it evolves according to Equation 4.40. We can thus proceed by induction to arbitrary k, and conclude that $W_{\alpha_1 \dots \alpha_k}$ evolve according to Equation 4.40.

Consider a dust fluid which varies on some short lengthscale $L_{\text{dust}}$ embedded with a gas that varies on a long lengthscale $L_{\text{gas}}$. This introduces a separation of scales for which we introduce $\xi$ for coordinates describing variation on the short dust lengthscale and $\mathbf{x}$ describing variation on the gas lengthscale. Naturally, the properties of the gas depend only on $\mathbf{x}$ (and time). We propose a nearly Maxwellian dust velocity distribution with the following asymptotic scheme

$$ \Pi_{\alpha_1 \dots \alpha_k} = \epsilon^k \rho(\xi, \mathbf{x}) W_{\alpha_1 \dots \alpha_k}(\mathbf{x}) + \epsilon^{k+1} \Sigma_{\alpha_1 \dots \alpha_k}(\xi, \mathbf{x}), \quad (4.42) $$

$$ U = U_0(\mathbf{x}) + \epsilon^k u_0(\xi, \mathbf{x}), \quad (4.43) $$

$$ \nabla = \epsilon^{-1} \frac{\partial}{\partial \xi} + \frac{\partial}{\partial \mathbf{x}}, \quad (4.44) $$

where $\epsilon$ is treated as a book-keeping parameter. Strictly speaking one should also expand the density, however the $O(\epsilon)$ terms due to the effects of the Non-Maxwellian velocity perturbation can be absorbed into the definition of $\Sigma_{\alpha_1 \dots \alpha_k}$. While we can often treat $\kappa = 2$ (i.e. the part of the mean velocity that varies on the dust lengthscale is $O(\epsilon^2)$) we shall assume $\kappa = 1$ throughout as this will allow for a wider range of dust flows.

Substituting Equations 4.42-4.44 into Equation 4.31, and making use of $D_1\rho = 0$, the evolutionary equation for the k-th velocity moment becomes:

$$ \begin{aligned} D_1 \Pi_{\alpha_1 \dots \alpha_k} &= \epsilon^k \rho D W_{\alpha_1 \dots \alpha_k} + \epsilon^{k+1} u_0^\alpha \frac{\partial}{\partial x^\alpha} W_{\alpha_1 \dots \alpha_k} + \epsilon^{k+1} D_1 \Sigma_{\alpha_1 \dots \alpha_k} \\ &= -\epsilon^k \left( \frac{\partial}{\partial \xi^\sigma} + \epsilon \frac{\partial}{\partial x^\sigma} \right) (\rho W_{\alpha_1 \dots \alpha_k \sigma}) - \epsilon^{k+1} \left( \frac{\partial}{\partial \xi^\sigma} + \epsilon \frac{\partial}{\partial x^\sigma} \right) \Sigma_{\alpha_1 \dots \alpha_k \sigma} \\ &\quad - k \epsilon^k \left[ \rho W_{(\alpha_1 \dots \alpha_{k-1})}^\sigma + \epsilon \Sigma_{(\alpha_1 \dots \alpha_{k-1})}^\sigma \right] B_{\alpha_k} \sigma \\ &\quad - (\rho W_{(\alpha_1 \dots \alpha_{k-1})} + \epsilon \Sigma_{(\alpha_1 \dots \alpha_{k-1})}) \left( \frac{\partial}{\partial \xi^\sigma} + \epsilon \frac{\partial}{\partial x^\sigma} \right) (\rho W_{\alpha_k} \sigma + \epsilon \Sigma_{\alpha_k} \sigma) \\ &\quad - (k-1) (\rho W_{(\alpha_1 \dots \alpha_{k-2})} + \epsilon \Sigma_{(\alpha_1 \dots \alpha_{k-2})}) D_{\alpha_{k-1} \alpha_k} \end{aligned} \quad (4.45) $$

where, here, $D = \partial_t + U_0^\alpha \nabla_\alpha$ is the Lagrangian time derivative with respect to the leading order flow described by $U_0$.

Making use of Equation 4.40 for the evolution of $W_{\alpha_1 \dots \alpha_k}$, along with the recurrence relation for $W_{\alpha_1 \dots \alpha_k}$ (Equation 4.38) and rearranging we obtain and equation for the evolution of the non-Maxwellian part of the velocity moment,
$$
\begin{align}
& \mathcal{D}_1 \Sigma_{\alpha_1 \cdots \alpha_k} + u_0^\alpha \frac{\partial}{\partial x^\alpha} W_{\alpha_1 \cdots \alpha_k} + k \rho W_\sigma (\alpha_1 \frac{\partial}{\partial x^\sigma} W_{\alpha_2 \cdots \alpha_k}) + \frac{\partial}{\partial \xi^\sigma} \Sigma_{\alpha_1 \cdots \alpha_k \sigma} \nonumber \\
& \qquad + k \left[ \Sigma_{(\alpha_1 \cdots \alpha_{k-1}) \alpha_k}^\sigma B_{\alpha_k} \sigma - \rho W_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial \xi^\sigma} \Sigma_{\alpha_k} \sigma \right. \nonumber \\
& \qquad \left. - \Sigma_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial \xi^\sigma} \rho W_{\alpha_k} \sigma - (k-1) \Sigma_{(\alpha_1 \cdots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k} \right] \tag{4.46} \\
& = \epsilon \left[ k \rho W_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial x^\sigma} \Sigma_{\alpha_k} \sigma + k \Sigma_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial x^\sigma} (\rho W_{\alpha_k} \sigma + \epsilon \Sigma_{\alpha_k} \sigma) \right. \nonumber \\
& \qquad \left. + k \Sigma_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial x^\sigma} \Sigma_{\alpha_k} \sigma - \frac{\partial}{\partial x^\sigma} \Sigma_{\alpha_1 \cdots \alpha_k \sigma} \right], \nonumber
\end{align}
$$

where the terms on the right hand side are all sub-leading. Dropping these subleading terms we obtain

$$
\begin{equation}
\begin{split}
0 = \mathcal{D}_1 \Sigma_{\alpha_1 \cdots \alpha_k} &+ u_0^\alpha \frac{\partial}{\partial x^\alpha} W_{\alpha_1 \cdots \alpha_k} + k \rho W_\sigma (\alpha_1 \frac{\partial}{\partial x^\sigma} W_{\alpha_2 \cdots \alpha_k}) + \frac{\partial}{\partial \xi^\sigma} \Sigma_{\alpha_1 \cdots \alpha_k \sigma} \\
& + k \left[ \Sigma_{(\alpha_1 \cdots \alpha_{k-1}) \alpha_k}^\sigma B_{\alpha_k} \sigma - \rho W_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial \xi^\sigma} \Sigma_{\alpha_k} \sigma \right. \\
& \qquad \left. - \Sigma_{(\alpha_1 \cdots \alpha_{k-1})} W_{\alpha_k} \sigma \frac{\partial}{\partial \xi^\sigma} \rho - (k-1) \Sigma_{(\alpha_1 \cdots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_k} \right].
\end{split}
\tag{4.47}
\end{equation}
$$

This confirms that the asymptotic ordering scheme (Equations 4.42-4.44) is self consistent and the non-Maxwellian terms are suppressed by a factor of $\epsilon \sim L_{\text{dust}}/L_{\text{gas}}$ relative to the Maxwellian terms. However, for the purposes of the equation of motion the pressure gradients are the more important quantity. For the nearly Maxwellian velocity distribution considered here the stress gradients are

$$
\begin{align}
\nabla_{\beta} \Pi^{\beta \alpha} &= \epsilon \left( \frac{\partial}{\partial \xi^{\beta}} + \frac{\partial}{\partial x^{\beta}} \right) \rho W^{\beta \alpha} + \epsilon^2 \left( \frac{\partial}{\partial \xi^{\beta}} + \frac{\partial}{\partial x^{\beta}} \right) \Sigma_{\beta \alpha} \tag{4.48} \\
&= \epsilon W^{\beta \alpha} \frac{\partial}{\partial \xi^{\beta}} \rho + O(\epsilon^2). \nonumber
\end{align}
$$

Thus the effects of the non-Maxwellian terms are $O(\epsilon^2)$, and are thus small relative to the acceleration and gravity, which are taken to be $O(1)$, when the dust layer is dynamically cool.

### 4.3.3. Are the ordering schemes attractors?

We have two separate situations where we can truncate the moment expansion by neglecting the third (and higher) velocity moment(s). The first is when $\text{St} \lesssim 1$, meaning that the dust is tightly coupled to the gas and the higher order velocity moments are suppressed by interaction with the gas. The second is for dynamically cool dust layers where $L_{\text{dust}} \ll L_{\text{gas}}$ (with the lengthscale typically being the dust and gas scale heights), where the non-Maxwellian velocity moments are suppressed by the confinement of the dust. This latter scenario is of interest for
dust with $St > 1$ in gas flows which are not strongly stirred, in the presence of vertical gravity, as these would be expected to settle into a hydrostatically supported dust layer which is much thinner than a hydrostatically supported gas flow. Of course the existence of a consistent asymptotic scaling does not guarantee that the fluid regime is an attractor. While a complete exploration of when this state becomes an attractor, and thus allow for a fluid treatment of the dust, is beyond the scope of this work, in this section we shall present an argument showing that velocity moments which start far from this asymptotic scaling are expected to damp towards this scaling, subject to the dust fluid being thermally stable.

Consider a situation where the either the well coupled or near-Maxwellian ordering scheme holds. We wish to explore what happens where some perturbation increases the k-th velocity moment sufficiently such that it breaks the ordering scheme. If the k-th velocity moment is large, while all other velocity moments keep the same ordering as in the fluid ordering schemes, then the only terms that are important in the evolutionary equation for k-th velocity moment are those involving $\Pi_{\alpha_1 \cdots \alpha_k}$. Thus evolution of the k-th velocity moment is approximately described by

$$ \mathcal{D}_1 \Pi_{\alpha_1 \cdots \alpha_k} = -k \Pi_{(\alpha_1 \cdots \alpha_{k-1})}^{\sigma} B_{\alpha_k} \sigma. \quad (4.49) $$

Defining $W_{\alpha_1 \cdots \alpha_k} = \rho^{-1} \Pi_{\alpha_1 \cdots \alpha_k}$ this simplifies to

$$ DW_{\alpha_1 \cdots \alpha_k} = -k W_{(\alpha_1 \cdots \alpha_{k-1})}^{\sigma} B_{\alpha_k} \sigma. \quad (4.50) $$

We wish to show that $W_{\alpha_1 \cdots \alpha_k}$ decays subject to certain constraints on $B_{\alpha\beta}$. To do this we make use of the adjoint problem,

$$ DY^{\alpha} = B_{\beta}^{\alpha} Y^{\beta}, \quad (4.51) $$

in order to relate the evolutionary equation for $W_{\alpha_1 \cdots \alpha_k}$ for arbitrary $k$ to that with $k=2$. This allows us to relate the behaviour of Equation 4.50 to properties of the constitutive equation, in particular the thermal stability of the flow.

Equations 4.50 and 4.51 are related by an invariant scalar $\chi = W_{\alpha_1 \cdots \alpha_k} Y^{\alpha_1} \cdots Y^{\alpha_k}$, with

$$ \begin{aligned} D\chi &= Y^{\alpha_{n+1}} \cdots Y^{\alpha_k} DQ_{\alpha_{n+1} \cdots \alpha_k} + (k-n)Q_{\alpha_{n+1} \cdots \alpha_k} Y^{\alpha_{n+1}} \cdots Y^{\alpha_{k-1}} DY^{\alpha_k} \\ &= -(k-n)Y^{\alpha_{n+1}} \cdots Y^{\alpha_k} Q_{\sigma\alpha_{n+1} \cdots \alpha_{k-1}} M_{\alpha_k}^{\sigma} + (k-n)Q_{\sigma\alpha_{n+1} \cdots \alpha_{k-1}} Y^{\alpha_{n+1}} \cdots Y^{\alpha_{k-1}} M_{\beta}^{\sigma} Y^{\beta} \\ &= 0. \end{aligned} \quad (4.52) $$

Consider now $k=2$ and define the associated scalar, $\zeta = Q_{\alpha\beta}Y^{\alpha}Y^{\beta}$, we also assume $Q_{\alpha\beta}$ is positive definite at $t=t_0$. Without loss of generality we can take $Y^{\alpha} = y_0^{\alpha}$ at $t=t_0$, where $|y_0|=1$, such that

$$ \zeta = (Q_{\alpha\beta} Y^{\alpha} Y^{\beta})|_{t=t_0} = Q_{\alpha\beta} y_0^{\alpha} y_0^{\beta} > 0. \quad (4.53) $$

At $t=t_1 > t_0$ we write $Y^{\alpha} = \mathcal{Y}y^{\alpha}$, with $|\mathcal{Y}|=1$, such that

$$ \zeta = (Q_{\alpha\beta} Y^{\alpha} Y^{\beta})|_{t=t_1} = \mathcal{Y}^2 Q_{\alpha\beta} y^{\alpha} y^{\beta}. \quad (4.54) $$

As $Q_{\alpha\beta}$ is positive semi-definite, for all $t$, $Q_{\alpha\beta}y^{\alpha}y^{\beta} \ge 0$. Using the fact that $\zeta$ is constant, we obtain
$$ \mathcal{Y}^2 = \frac{Q_{\alpha\beta} y_0^{\alpha} y_0^{\beta}}{Q_{\alpha\beta} y^{\alpha} y^{\beta}}, \qquad (4.55) $$

In order for the fluid to be thermally stable $Q_{\alpha\beta}$ must ultimately decay towards zero. If this were not the case then there would exist components of $\Pi_{\alpha\beta}$ where heating by the disc turbulence is not balanced by cooling from the $\Pi_{(\alpha}^{\sigma} B_{\beta)}\sigma$ term, and would thus experience thermal runaway. Thus for $\delta > 0$, there exists a $t = t_{\text{cool}} > t_0$ such that the components of $Q_{\alpha\beta}$ at $t = t_{\text{cool}}$ satisfy $|Q_{\alpha\beta}| < \delta$. It should be noted that certain components of $Q_{\alpha\beta}$ can experience transient growth (e.g. due to the shearing out of the initial conditions), but must ultimately decline in order to ensure thermal stability. As $\delta$ is arbitrary we can choose $\delta$ small enough such that

$$ Q_{\alpha\beta} y^{\alpha} y^{\beta} \le \sum_{\alpha,\beta} |Q_{\alpha\beta}||y^{\alpha}||y^{\beta}| < \delta \sum_{\alpha,\beta} |y^{\alpha}||y^{\beta}| < Q_{\alpha\beta} y_0^{\alpha} y_0^{\beta}, \qquad (4.56) $$

for $t > t_{\text{cool}}$. From this we can conclude that $\mathcal{Y} > 1$, for $t > t_{\text{cool}}$. By choosing $t$ large enough we can make $\mathcal{Y}$ arbitrarily large. Typically one expects $t_{\text{cool}}$ to be of order the cooling/settling time in the fluid as this decay is linked to the dynamical cooling of the dust fluid.

Now consider the scalar $\chi$ associated with $Q_{\alpha_{n+1}\cdots\alpha_k}$. As $\chi$ is constant we have,

$$ \begin{aligned} Q_{\alpha_{n+1} \cdots \alpha_k} \Big|_{t=t_0} y_0^{\alpha_{n+1}} \cdots y_0^{\alpha_k} &= (Q_{\alpha_{n+1} \cdots \alpha_k} Y^{\alpha_{n+1}} \cdots Y^{\alpha_k})_{t=t_0} \\ &= (Q_{\alpha_{n+1} \cdots \alpha_k} Y^{\alpha_{n+1}} \cdots Y^{\alpha_k})_{t=t_1} \qquad (4.57) \\ &= \mathcal{Y}^k Q_{\alpha_{n+1} \cdots \alpha_k} \Big|_{t=t_1} y^{\alpha_{n+1}} \cdots y^{\alpha_k}. \end{aligned} $$

Rearranging this we obtain,

$$ Q_{\alpha_{n+1} \cdots \alpha_k} \Big|_{t=t_1} y^{\alpha_{n+1}} \cdots y^{\alpha_k} = \mathcal{Y}^{-k} Q_{\alpha_{n+1} \cdots \alpha_k} \Big|_{t=t_0} y_0^{\alpha_{n+1}} \cdots y_0^{\alpha_k} \le \mathcal{Y}^{-k} \sum_{\alpha_{n+1} \cdots \alpha_k} \left| Q_{\alpha_{n+1} \cdots \alpha_k} \Big|_{t=t_0} \right|. \qquad (4.58) $$

Again, by choosing $t_1$ large enough we can take $\mathcal{Y}$ to be as large as we like, this means that $Q_{\alpha_{n+1} \cdots \alpha_k} \Big|_{t=t_1} y^{\alpha_{n+1}} \cdots y^{\alpha_k}$ can be made arbitrarily small. As we can do this for any unit vector $y^\alpha$, and $Q_{\alpha_{n+1} \cdots \alpha_k}$ is symmetric, we conclude that the components of $Q_{\alpha_{n+1} \cdots \alpha_k}$ will become arbitrarily small at late times

Thus we can conclude, from the above argument, that thermal stability of the fluid flow is a necessary, and sufficient condition for $W_{\alpha_1 \cdots \alpha_k}$ to decay. This means that the k-th velocity moment, $\Pi_{\alpha_1 \cdots \alpha_k}$, decays when its evolution can be well approximated by Equation 4.49 and the fluid is thermally stable. This implies that the dust fluid ordering schemes are stable to (nonlinear) perturbations to the higher order velocity moments, which should damp until they are compatible with the fluid ordering scheme derived above on approximately the cooling/settling time of the dust fluid.

The above argument shows that thermal stability is a necessary condition for the fluid dust description to remain valid. It is not, however, a sufficient condition as the argument only applies to (nonlinear) perturbations to the k-th velocity moment in isolation. Thus, in principle, there could exist perturbations to the multiple orders of velocity moment simultaneously that can be sustained and will not damp towards the fluid ordering scheme. For now we shall work under the assumption that thermal stability is sufficient to ensure the damping of higher-order velocity moments, however the exploration of the stability of
the fluid description against more general perturbations should be explored if the dust fluid model finds widespread use.

### 4.4. Obtaining the dust fluid equations

As a result of the asymptotic argument presented above, we can take $\Pi_{\alpha\beta\gamma} = 0$ and only consider the first 3 moments of the Fokker-Planck equation. This yields a continuity, momentum and constitutive relation for a 6-dimensional dust fluid. This dust fluid has a high degree of symmetry as physical properties must be independent of the gas displacement $\{\mathbf{x}_g\}^\dagger$. One can, therefore, integrate out these redundant degrees of freedom.

In integrating out the dummy gas degrees of freedom, we replace the connections $\nabla_\alpha$ with the connections $\bar{\nabla}_\alpha$ which ensure that the components of vectors associated with the integrated out directions remain correctly aligned. Our normalisation means that $\rho_d = \bar{\rho}_6$, we also introduce the rheological stress tensor $T_{\alpha\beta} = \bar{\Pi}_{\alpha\beta}$ and we can always choose the size of the dummy gas dimensions such that $U^\alpha = \bar{U}^\alpha$. With these choices the dust fluid equations are

$$ \bar{D}\rho_d = -\rho_d\bar{\nabla}_\alpha U^\alpha, \quad (4.59) $$

$$ \rho_d \bar{D} U_\alpha = \rho_d F_\alpha - \bar{\nabla}^\beta T_{\alpha\beta} - \rho_d C_{\alpha\beta} (U^\beta - U_g^\beta), \quad (4.60) $$

$$ \bar{D}_2 T_{\alpha\beta} = -2 \left( T^\gamma_{(\alpha} C_{\beta)\gamma} - \rho_d D_{(\alpha\beta)} \right), \quad (4.61) $$

where

$$ \bar{D} = \partial_t + U^\alpha \bar{\nabla}_\alpha = \begin{pmatrix} \partial_t + u_d^i \nabla_i & \mathbf{0} \\ \mathbf{0} & \partial_t + u_d^i \nabla_i \end{pmatrix}, \quad (4.62) $$

which corresponds to the usual (3-dimensional) Lagrangian time derivative with respect to the mean dust flow applied to the dust and dummy gas components of the (6-D) tensor independently. $C_{\alpha\beta}$ and $D_{\alpha\beta}$ are given by Equations 3.23 and 3.24 (for isotropic stochastic driving). Finally the operator $\bar{D}_2$, when acting on $T_{\alpha\beta}$, is given by

$$ \bar{D}_2 T_{\alpha\beta} = \bar{D} T_{\alpha\beta} + 2 T_{\gamma d(\alpha} \bar{\nabla}^{\gamma d} U_{\beta)} + T_{\alpha\beta} \bar{\nabla}_{\gamma d} U^{\gamma d}. \quad (4.63) $$

Finally to highlight the effects of the torsion we consider its contribution to the dust fluid vorticity,

$$ \omega^\gamma \epsilon_{\gamma\alpha\beta} = 2\bar{\partial}_{[\alpha} U_{\beta]} + S_{\alpha\beta}^{\gamma g} U_{\gamma g}, \quad (4.64) $$

meaning

$$ \omega^\gamma \epsilon_{\gamma\alpha g\beta d} = -\bar{\partial}_{\beta d} U_{\alpha g} - \mathcal{T}_{\beta d}^{\gamma g} \alpha_g^* U_{\gamma g}. \quad (4.65) $$

This additional contribution to the vorticity is associated with the rotation of the gas displacement vectors.

For Cartesian coordinates $(x, y, z, x_g, y_g, z_g)$, in Euclidean space, with $C_{\alpha\beta}$ and $D_{\alpha\beta}$ are given by Equations 3.23 and 3.24, then Equations 4.59-4.61 are explicitly

<footnote>† For the Stochastic model considered here. If the fluid equations were to be derived based on a two-step Stochastic model, as outlined in Minier & Henry (2023), then the dummy gas variable would influence the fluid model, which may allow the effects of spatial correlations to be included.</footnote>
$$ (\partial_t + U^{\alpha_d} \partial_{\alpha_d}) \rho_d = -\rho_d \partial_{\alpha_d} U^{\alpha_d}, \qquad (4.66) $$

$$ \rho_d (\partial_t + U^{\alpha_d} \partial_{\alpha_d}) U_{\alpha_d} = -\rho_d \partial_{\alpha_d} \phi - \partial_{\beta_d} T_{\alpha_d}^{\beta_d} - \frac{\rho_d}{t_s} (U_{\alpha_d} - U_{\alpha_g}), \qquad (4.67) $$

$$ \begin{aligned} \rho_d (\partial_t + U^{\alpha_d} \partial_{\alpha_d}) U_{\alpha_g} = & -\rho_d \partial_{\alpha_g^*} \phi - f_d \partial_{\alpha_g^*} p_g + \rho_d (U^{\beta_d} - U_g^{\beta_d}) \partial_{\beta_d} U_{\alpha_d}^g \\ & - \partial_{\beta_d} T_{\alpha_g}^{\beta_d} - \frac{\rho_d}{t_c} (U_{\alpha_g} - U_{\alpha_g}^g), \end{aligned} \qquad (4.68) $$

$$ \begin{aligned} (\partial_t + U^{\alpha_d} \partial_{\alpha_d}) T_{\alpha_d \beta_d} + T_{\alpha_d}^{\gamma_d} \partial_{\gamma_d} U_{\beta_d} + T_{\beta_d}^{\gamma_d} \partial_{\gamma_d} U_{\alpha_d} + T_{\alpha_d \beta_d} \partial_{\gamma_d} U^{\gamma_d} \\ = & -\frac{2}{t_s} T_{\alpha_d \beta_d} + \frac{1}{t_s} T_{\alpha_d \beta_d^*} + \frac{1}{t_s} T_{\alpha_d^* \beta_d}, \end{aligned} \qquad (4.69) $$

$$ \begin{aligned} (\partial_t + U^{\alpha_d} \partial_{\alpha_d}) T_{\alpha_d \beta_g} + T_{\alpha_d}^{\gamma_d} \partial_{\gamma_d} U_{\beta_g} + T_{\beta_g}^{\gamma_d} \partial_{\gamma_d} U_{\alpha_d} + T_{\alpha_d \beta_g} \partial_{\gamma_d} U^{\gamma_d} \\ = & -\left(\frac{1}{t_c} + \frac{1}{t_s}\right) T_{\alpha_d \beta_g} + \frac{1}{t_s} T_{\alpha_d^* \beta_g}, \end{aligned} \qquad (4.70) $$

$$ \begin{aligned} (\partial_t + U^{\alpha_d} \partial_{\alpha_d}) T_{\alpha_g \beta_g} + T_{\alpha_g}^{\gamma_d} \partial_{\gamma_d} U_{\beta_g} + T_{\beta_g}^{\gamma_d} \partial_{\gamma_d} U_{\alpha_g} + T_{\alpha_g \beta_g} \partial_{\gamma_d} U^{\gamma_d} \\ = & -\frac{2}{t_c} \left( T_{\alpha_g \beta_g} - \alpha c_s^2 \rho_d \delta_{\alpha_g \beta_g} \right). \end{aligned} \qquad (4.71) $$

In the 3-dimensional picture we have $u_d^i = U^i$, $u_s^i = U^{i*}$ are the (3D) dust velocity and fluid seen respectively and $p_{ij} = T_{ij}$, $\tau_{ij} = T_{ij^*}$, $\sigma_{ij} = T_{i^*j^*}$ are the dust kinetic tensor, dust-gas correlation tensor and Reynolds stress of the fluid seen respectively. Equations 4.66-4.71 are equivalent to

$$ (\partial_t + u_d^i \partial_i) \rho_d = -\rho_d \partial_i u_d^i, \qquad (4.72) $$

$$ \rho_d (\partial_t + u_d^j \partial_j) u_i^d = -\rho_d \partial_i \phi - \partial_j p_i^j - \frac{\rho_d}{t_s} (u_i^d - u_i^s), \qquad (4.73) $$

$$ \rho_d (\partial_t + u_d^j \partial_j) u_i^s = -\rho_d \partial_i \phi - f_d \partial_i p_g + \rho_d (u_d^j - u_g^j) \partial_j u_i^g - \partial_j \tau_i^j - \frac{\rho_d}{t_c} (u_i^s - u_i^g), \qquad (4.74) $$

$$ (\partial_t + u_d^k \partial_k) p_{ij} + p_i^k \partial_k u_j^d + p_j^k \partial_k u_i^d + p_{ij} \partial_k u_d^k = -\frac{2}{t_s} p_{ij} + \frac{1}{t_s} \tau_{ij} + \frac{1}{t_s} \tau_{ji}, \qquad (4.75) $$

$$ (\partial_t + u_d^k \partial_k) \tau_{ij} + p_i^k \partial_k u_j^s + \tau_j^k \partial_k u_i^d + \tau_{ij} \partial_k u_d^k = -\left(\frac{1}{t_c} + \frac{1}{t_s}\right) \tau_{ij} + \frac{1}{t_s} \sigma_{ij}, \qquad (4.76) $$

$$ (\partial_t + u_d^k \partial_k) \sigma_{ij} + \tau_i^k \partial_k u_j^s + \tau_j^k \partial_k u_i^s + \sigma_{ij} \partial_k u_d^k = -\frac{2}{t_c} \left( \sigma_{ij} - \alpha c_s^2 \rho_d \delta_{ij} \right). \qquad (4.77) $$

# 5. Properties of the dust fluid model

We will now describe some key features of our dust fluid model.

## 5.1. The mean gas velocity experienced by the dust is different to that experienced by the gas

Unlike the pressureless, non-turbulent models the dust experiences a different mean gas velocity to the gas. Part of this is due to the “crossing trajectory effect” (e.g. see Minier 2001; Minier et al. 2004, 2014) where the mean gas velocity “seen” by the dust is that following the Lagrangian trajectory traced by the dust, rather than that traced by the fluid particles. In addition to this, the dust experience a subsample of the gas velocity field rather than the gas
velocity field itself. This distinction is vital for producing dust dispersion by the turbulence. If the dust experienced the same gas velocity distribution as the gas then a local dust density maxima of perfectly coupled dust would not spread in homogeneous gas turbulence. The dust to gas density ratio (in the 3D picture), in such a setup, evolves according to

$$ \partial_t (\rho_d / \rho_g) = -\nabla_i (\rho_d u_d^i / \rho_g) \quad (5.1) $$

$$ = -\nabla_i (\rho_d u_s^i / \rho_g) \quad (5.2) $$

Thus, in order that the gas turbulence disperse the dust, we require the velocity of the fluid seen $u_s^i \neq u_g^i = 0$. The subsampling of the gas velocity distribution means the larger number of dust grains at the centre of the overdensity experience more “draws” from the gas velocity distribution and thus experience a greater gas dispersion (this would equally be true for “marked” gas fluid elements). This means the dust experiences a mean gas flow directed away from the maxima due to the resulting gradient in the cross pressure.

## 5.2. Anisotropic Dust Rheological Stress Tensor

The most important feature of the dust fluid model is the fluid stress is not zero, and can be dynamically important. In fact one expects dust settling/drift to concentrate dust until dust stress gradients become dynamically important. Also present is a form of “cross-pressure”, arising from correlations between the dust and gas motion, which modifies the mean-gas velocity experienced by the dust.

This rheological stress tensor is anisotropic in the presence of strong shear or rotation. In general, the gas turbulence heats the dust and isotropises the dust stress tensor on timescales longer than the correlation time. However, in strong shear flows the velocity dispersion induced in the dust by the turbulence is sheared out resulting in an anisotropic stress tensor (just as happens for the gas Reynolds tensor). The flow vorticity also provides an additional anisotropic heating term in the dust. In Section 7 we explore this effect further by considering the dust stress tensor in a rotating shear flow.

As we shall show in the next section, the presence of a nonzero elastic stress means the dust fluid supports waves, specifically seismic waves.

## 5.3. Viscoelasticity

The dust fluid exhibits viscoelastic behaviour (See Appendix B.2) This behaviour is easiest to see when $t_s \sim t_c = O(\text{De})$, where $\text{De} = t_r/t_f$ is the Deborah number of the dust fluid, which is the ratio of the characteristic relaxation time $t_r \sim t_s \sim t_c$ to the characteristic fluid timescale $t_f$. When $\text{De} \gg 1$ the dust stress tensor evolves according to

$$ \overline{D}_2 T_{\alpha\beta} = \overline{D} T_{\alpha\beta} - 2T^{\gamma}_{(\alpha} \varepsilon_{\beta)\gamma\sigma} \omega^{\sigma} = 0. \quad (5.3) $$

This corresponds to an elastic stress with a vortical heating (or “gyroscopic motion” Gavrilyuk & Gouin 2012) term, and evolves in an identical manner to a Reynolds stress in the absence of source terms. When $\text{De} \ll 1$ the stress tensor is approximately

$$ T_{\alpha\beta} = p_d \left( 1 + \frac{t_s}{t_c} \Theta^g_{\alpha\beta} \right) g_{\alpha\beta} + \frac{1}{2} p_x (g_{\alpha\beta^*} + g_{\alpha^*\beta}) - 2\mu^{\mu\nu} \overline{\nabla}_{\mu} U_{\nu} + O(\text{De}^2). \quad (5.4) $$

At leading order this consists of an isotropic, isothermal, effective, dust pressure with sound speed $\sqrt{\frac{\alpha}{1+t_s/t_c}} c_s$, a cross pressure $p_x = p_d$, from the dust-gas velocity correlations, and an additional pressure like contribution to the dummy gas components of the Rheological
stress. The next terms in the expansion are an anisotropic viscous stress characterised by the viscosity tensor $\mu_{\alpha\beta}^{\mu\nu}$; including a “cross” viscosity, which likely encapsulates the the decorrelation of the dust and gas velocities in the presence of shear. Explicit expressions for $\mu_{\alpha\beta}^{\mu\nu}$ are given in Appendix B.2. For weak gas turbulence $\alpha \ll 1$ the the viscous terms, for small dust grains, are typically negligible and the dust primarily behaves like an inviscid isothermal gas with a lower temperature than than the gas phase. The difference between $U_{\alpha g}^g$ and $U_{\alpha d}$ (mean gas velocities experienced by the gas and dust respectively) as a result of the cross pressure term allows for dust diffusion to occur in this limit.

The local expression, Equation 5.4, arises due to the fact that in the De $\ll 1$ limit $\mathbf{v} - \mathbf{v}^g$ and $\mathbf{v}^g - \mathbf{u}^g$, in the original stochastic-differential equations, are “fast variables” with no memory of the previous fluid state (Minier 2016). For $\tau_c \ll 1$, but $\text{St} \sim 1$, only $\mathbf{v}^g - \mathbf{u}^g$ is a fast variable and we have a local closure for $T_{\alpha d\beta g}$ and $P_{\alpha g\beta g}$ but not $P_{\alpha d\beta d}$, which then has a fluid memory of order the stopping time. In the large Deborah number limit the fast terms in Equation 3.8 are negligible, resulting in a fully nonlocal behaviour for $T_{\alpha\beta}$ (Equation 5.3).

### 5.4. Eddy-Knudsen number effect

While it might be expected that small dust grains should closely follow the gas with the dust velocity correlations being set by the gas velocity correlations. This turns out to only be the case when the dust sees the turbulence as a continuum. This is explored further in Appendix B.3. Whether the dust see the turbulence as a continuum or is sensitive to individual eddies is determined by a form of “eddy Knudsen number”:

$$ \text{Kn}_e = \frac{\lambda}{L} = \frac{t_c \Delta U^*}{L}, \qquad (5.5) $$

where $\lambda = t_c \Delta U^*$ is the mean free path of a dust grain in the turbulent flow representing the lengthscale a dust grain is transported by a single eddy, $\Delta U^*$ is a characteristic velocity difference between the dust and gas and $L$ is a characteristic lengthscale of variations in the fluid flow.

When $\text{Kn}_e \ll 1$ a dust grain interacts with many turbulent eddies over the lengthscale on which the dust fluid varies, meaning the dust experiences the turbulence as a continuum of stochastic perturbations. When $\text{Kn}_e \gtrsim 1$ the dust is instead strongly affected by individual eddies (in a similar manner to how weakly collisional gases can be strongly perturbed by individual collisions). Thus, in this regime, the dust is sensitive to individual eddies. In the short stopping time limit the equation for the dust stress simplify to (see Appendix B.3)

$$ \begin{aligned} t_c \tilde{D} \left( \frac{T_{\alpha g \beta_g}}{\rho_d} \right) &+ 2 \frac{t_c}{\rho_d} T_{\gamma_d}(\alpha_g) \nabla^k u_{\beta_g}^g + 2 \left( \frac{T_{\alpha_g \beta_g}}{\rho_d} - \alpha c_s^2 g \alpha_g \beta_g \right) \\ & = - \text{Kn}_e \left[ \frac{\Delta U^\gamma}{\Delta U^*} L \bar{\nabla}_\gamma \left( \frac{T_{\alpha_g \beta_g}}{\rho_d} \right) + 2 \frac{T_\gamma(\alpha_g)}{\rho_d} L \bar{\nabla}^\gamma \frac{\Delta U_{\beta_g}}{\Delta U^*} \right]. \end{aligned} \qquad (5.6) $$

When $\text{Kn}_e \to 0$ this matches the equation governing the evolution of the gas velocity correlations, meaning the dust velocity correlations are indeed set by those of the gas. This is no longer the case when $\text{Kn}_e \sim 1$ and the dust velocity correlations can depart strongly from those of the gas, even when the mean velocity of the gas and dust remain tightly coupled.

## 6. Hyperbolic structure and linear waves

In this section we will rearrange the equations into hyperbolic form, which is useful for some types of numerical solver and for understand the wave modes in the system. We wish to find
a state vector **W**, matrices **A**<sub>*i*</sub> and source vector such that the dust-fluid equations take the form

$$ \frac{\partial \mathbf{W}}{\partial t} + \mathbf{A}^{\alpha} \overline{\nabla}_{\alpha} \mathbf{W} = \mathbf{S}, \qquad (6.1) $$

To start we rearrange the equations so that all the source/sink terms are on the right hand side,

$$ \dot{\rho}_d + U^\alpha \overline{\nabla}_\alpha \rho_d + \rho_d \overline{\nabla}_\alpha U^\alpha = 0, \qquad (6.2) $$

$$ \dot{T}_{\alpha\beta} + U^{\gamma} \overline{\nabla}_{\gamma} T_{\alpha\beta} + 2T_{\gamma(\alpha} \overline{\nabla}^{\gamma} U_{\beta)} + T_{\alpha\beta} \overline{\nabla}_{\gamma} U^{\gamma} = -2 \left( T^{\gamma}_{(\alpha} C_{\beta)\gamma} - \rho_d D_{\alpha\beta} \right), \qquad (6.3) $$

$$ \dot{U}_{\alpha} + U^{\beta} \overline{\nabla}_{\beta} U_{\alpha} + \frac{1}{\rho_d} \overline{\nabla}^{\beta} T_{\beta\alpha} = F_{\alpha} - C_{\alpha\beta} (U^{\beta} - U_g^{\beta}), \qquad (6.4) $$

where we have exchanged the momentum and constitutive relation as it will make $\mathbf{A}^{\alpha}$ easier to diagonalise. The state vector for this system is

$$ \mathbf{W} = \begin{pmatrix} \rho_d \\ \mathbf{T} \\ \mathbf{U} \end{pmatrix}, \qquad (6.5) $$

The source vector is

$$ \mathbf{S} = \begin{pmatrix} 0 \\ -\mathbf{C}\mathbf{T} - \mathbf{T}\mathbf{C}^{\mathrm{T}} + 2\rho_d \mathbf{D} \\ \mathbf{F} - \mathbf{C}(\mathbf{U} - \mathbf{U}_g) \end{pmatrix}. \qquad (6.6) $$

The matrices $\mathbf{A}^{\alpha}$ are given by

$$ \mathbf{A}^{\alpha} = \begin{pmatrix} U^{\alpha} & \mathbf{0} & \rho_d \hat{\mathbf{e}}^{\alpha} \\ 0 & U^{\alpha} \mathbf{I} & \mathbf{M}^{\alpha} \\ 0 & \frac{1}{\rho_d} \hat{\mathbf{e}}^{\alpha} & U^{\alpha} \mathbf{I} \end{pmatrix}, \qquad (6.7) $$

where **I** denote the identity matrix and

$$ (\mathbf{M}^{\alpha})_{\sigma\beta\gamma} = T_{\sigma\beta} \delta^{\alpha}_{\gamma} + 2T^{\alpha}_{(\sigma} \delta^{\gamma)}_{\beta)}, \qquad (6.8) $$

such that

$$ (\mathbf{M}^{\alpha} \overline{\nabla}_{\alpha} \mathbf{U})_{\sigma\beta} = (\mathbf{M}^{\alpha})_{\sigma\beta\gamma} \overline{\nabla}_{\alpha} U^{\gamma}, \qquad (6.9) $$

For the system to be hyperbolic we must show that all the eigenvalues of $\mathbf{A}^{\alpha} \hat{n}_{\alpha}$ are real for unit vector $\hat{n}_{\alpha}$, and the eigenvectors span the 28 dimensional state-space. Without loss of generality we can orient our coordinate system such that $\hat{n} = \hat{e}^1$ to point along the positive x-direction. Physically we must remember that the dummy gas and position dimensions are distinct, however we do not need to consider the case where **n** has non-zero components in the 'dummy gas' directions as we require physical quantities to be independent of $\mathbf{x}_g$.

It is useful to separate out the velocity into the velocity along the x direction (along the direction of propagation), $U_1$, and the velocity in the other directions, $U_{\alpha}$ (where, for the rest of this section, we take the indices $\alpha, \beta$ and $\gamma$ to run over 2, ..., 6.). We similarly separate out the stress tensor into compression along the x direction, $T_{11}$, shear components in the x direction $T_{1\alpha}$ and the components of the stress in other directions, $T_{\alpha\beta}$. We shall split the
momentum and constitutive relation in a similar manner. This results in the following state vector

$$ \mathbf{W} = \begin{pmatrix} \rho \\ T_{11} \\ T_{1\alpha} \\ T_{\alpha\beta} \\ U_1 \\ U_{\alpha} \end{pmatrix}, \qquad (6.10) $$

The eigenvalues, $v$, for $A^\alpha n_\alpha$ can be derived from the determinant of the following matrix

$$ \mathbf{A}^1 - v \mathbf{I} = \begin{pmatrix} u^x - v & 0 & \mathbf{0} & \mathbf{0} & \rho_d & \mathbf{0} \\ 0 & u^x - v & \mathbf{0} & \mathbf{0} & 3T_{11} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & (u^x - v)\mathbf{I} & \mathbf{0} & 2T_{1\alpha} & T_{11}\mathbf{I} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & (u^x - v)\mathbf{I} & T_{\alpha\beta} & 2T_{1(\alpha}\hat{\mathbf{e}}_{\beta)}^\mathbf{T} \\ 0 & \frac{1}{\rho_d} & \mathbf{0} & \mathbf{0} & (u^x - v) & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \frac{1}{\rho_d}\mathbf{I} & \mathbf{0} & \mathbf{0} & (u^x - v)\mathbf{I} \end{pmatrix}, \quad (6.11) $$

which has a characteristic equation

$$ (v - u^x)^{16} \left( (v - u^x)^2 - 3 \frac{T_{11}}{\rho_d} \right) \left( (v - u^x)^2 - \frac{T_{11}}{\rho_d} \right)^5 = 0. \qquad (6.12) $$

This results in 16 nonpropagating (in the fluid frame) wavemodes, with wavespeed $v = u^x$. These consist of the entropy wave with eigenvector $\begin{pmatrix} 1 \\ 0.27 \end{pmatrix}$ and 15 'stress' waves with eigenvectors $\begin{pmatrix} 0.7 \\ 0.6 \end{pmatrix}$, where we have introduced the notation $0_n = 0$, with $n$ denoting the number of zeros in the column.

Two of the propagating waves can be identified as P-waves, with wavespeed $v = u^x \pm \sqrt{\frac{3T_{11}}{\rho_d}}$. The P-waves are analogous to sound-waves, but with an anisotropic soundspeed, with seismic wavespeed anisotropy being a well known phenomena in geophysics (Thomsen 1986). These wavemodes have eigenvectors,

$$ \begin{pmatrix} \pm \rho_d \\ \pm 3T_{11} \\ \pm 3T_{1\alpha} \\ \pm \left( T_{\alpha\beta} + \frac{2}{T_{11}} T_{1(\alpha} T_{\beta)1} \right) \\ -\sqrt{\frac{3T_{11}}{\rho_d}} \\ -\sqrt{\frac{3}{\rho_d T_{11}}} T_{1\alpha} \end{pmatrix}, \qquad (6.13) $$

Finally there are 10 propagating waves which can be identified as S-waves, with wavespeed $v = u^x \pm \sqrt{\frac{T_{11}}{\rho_d}}$. As is typical for elastic media, the S-waves have slower wavespeeds than the P-waves. These wavemodes have eigenvectors,
$$ \begin{pmatrix} 0 \\ 0 \\ \pm T_{11} \hat{\mathbf{e}}_{\gamma} \\ \pm 2T_{1(\alpha)} \delta_{\beta}^{\gamma} \\ 0 \\ -\sqrt{\frac{T_{11}}{\rho_d}} \hat{\mathbf{e}}_{\gamma} \end{pmatrix}. \qquad (6.14) $$

These eigenvectors span the 28-dimensional state-space of the dust fluid model.

Upon decomposing the velocity and pressure tensor, the source vector is given by

$$ S = \begin{pmatrix} 0 \\ -2T_1^1 C_{11} - 2T_1^\gamma C_{1\gamma} \\ -2T_{(1}\alpha) C_{\alpha)1} - 2T_{(1}\gamma}^\gamma C_{\alpha)\gamma} \\ -2T_{(\alpha}^1 C_{\beta)1} - 2T_{(\alpha}^\gamma C_{\beta)\gamma} + 2\rho_d D_{\alpha\beta} \\ F_1 - C_{11}(U^1 - U_g^1) - C_{1\gamma}(U^\gamma - U_g^\gamma) \\ F_\alpha - C_{\alpha 1}(U^1 - U_g^1) - C_{\alpha\gamma}(U^\gamma - U_g^\gamma) \end{pmatrix}, \qquad (6.15) $$

Thus we see that there are no source terms for the entropy wave. Turbulent diffusion (D) and the drag dependant coupling between pressure tensor components are sources/sinks of the stress waves. Finally the force per unit mass F and drag terms are sources/sinks of the P and S-waves. In practice whether the wavemodes can propagate in the dust fluid will depend on these source/sink terms as strong damping (such as by drag) may cause the waves to be evanescent in certain regions of parameter space.

While the aforementioned wavemodes represent all the waves present in the bulk. The dust fluid can support additional wavemodes when it occupies a thin layer, or other gravitationally confined structure. In such a situation the disc posses dust breathing modes associated with periodic oscillations of the dust scale-height. These are analogous to the surface waves in seismology.

# 7. Rheological Stress in a Rotating Shear Flow

## 7.1. Steady state

In order to better understand the behaviour of the rheology, we consider the specific example of a steady rotating shear flow in the kinematic limit (i.e. we impose a rotation profile in the dust and gas and neglect the modification to the fluid flow from the resulting stress gradients). Rotating shear flows are of particular interest in astrophysics as they are important for understanding accretion discs. They are also a common setup in experimental fluid dynamics (e.g. Taylor-Couette flows). To study this problem we adopt (6-D) cylindrical polar coordinates (R, φ, z, R_g, φ_g, z_g) with metric tensor components

$$ g_{RR} = g_{R_g R_g} = g_{zz} = g_{z_g z_g} = 1, \quad g_{\phi\phi} = g_{\phi_g \phi_g} = R^2. \qquad (7.1) $$

and connection coefficients,

$$ \Gamma_{\phi\phi}^{R} = \Gamma_{\phi\phi}^{R_g} = -R, \qquad (7.2) $$

$$ \Gamma_{\phi R}^{\phi} = \Gamma_{R \phi}^{\phi} = \Gamma_{\phi R_g}^{\phi_g} = \Gamma_{R \phi_g}^{\phi_g} = 1/R, \qquad (7.3) $$
with all other components zero. The fluid flow consists of a rotating shear flow where both the dust and gas rotate on cylinders with angular velocity $\Omega = \Omega(R)$. This leads to the 6-D mean velocity of the dust fluid of

$$U^{\gamma} = \Omega(R)(\hat{e}_{\phi}^{\alpha} + \hat{e}_{\phi g}^{\alpha}). \qquad (7.4)$$

We additionally assume that the fluid is vertically homogeneous. By specifying that the dust mean velocity should exactly follow that of the gas we are implicitly taking the zero Eddy-Knudsen number limit. Thus the stress for small Stokes dusts is entirely specified by the velocity correlations in the gas.

With this geometry, and imposed velocity, the operator $\overline{D}_2$ (when acting on $T_{\alpha\beta}$) is

$$\overline{D}_2 T_{\alpha\beta} = \overline{D} T_{\alpha\beta} + 2T^{R}_{(\alpha}(\hat{e}_{\beta}^{\phi} + \hat{e}_{\beta}^{\phi_g}) \partial_R (R^2\Omega) - 2\Omega\Gamma^{\gamma}_{\phi(\alpha} T_{\beta)\gamma} - 2R^2\Omega(\Gamma^{\phi}_{\gamma(\alpha} + \Gamma^{\phi_g}_{\gamma(\alpha}})T^{\gamma}_{\beta)} \qquad (7.5)$$

We are interested in the steady state solution to the stress tensor with $\overline{DT}_{\alpha\beta} = 0$. We thus have the following for the constitutive relation in the steady rotating shear flow,

$$4R(\Omega-A)T^{R}_{(\alpha}(\hat{e}_{\beta}^{\phi}+\hat{e}_{\beta}^{\phi_g})-2\Omega\Gamma^{\gamma}_{\phi(\alpha}T_{\beta)\gamma}-2R^2\Omega(\Gamma^{\phi}_{\gamma(\alpha}+\Gamma^{\phi_g}_{\gamma(\alpha}})T^{\gamma}_{\beta)}=-2(T^{\gamma}_{(\alpha}C_{\beta)\gamma}-\rho_d D_{\alpha\beta}), \qquad (7.6)$$

where we have introduced Oort's first constant $A = -(R/2)\Omega_R$, which is a measure of the fluid shear rate. The Rayleigh stability criterion corresponds to $A/\Omega < 1$.

Explicitly, the “dust-dust” components of Equation 7.6, which can be thought of as the equations governing the behaviour of the 3D dust stress, are

$$-\frac{4\Omega}{R}T_{R\phi} = -\frac{2}{t_s}(T_{RR} - T_{RR_g}), \qquad (7.7)$$

$$-\frac{2\Omega}{R}T_{\phi\phi} + 2R(\Omega - A)T_{RR} = -\frac{2}{t_s}T_{R\phi} + \frac{1}{t_s}(T_{R\phi_g} + T_{R_g\phi}), \qquad (7.8)$$

$$4R(\Omega - A)T_{R\phi} = -\frac{2}{t_s}(T_{\phi\phi} - T_{\phi\phi_g}). \qquad (7.9)$$

The “dust-gas” components Equation 7.6, which govern the behaviour of the “cross-stress” - from the cross correlation between the dust and gas velocities, are

$$-\frac{2\Omega}{R}T_{R_g\phi} - \frac{\Omega}{R}(T_{R\phi} + T_{R\phi_g}) = -\left(\frac{1}{t_s} + \frac{1}{t_c}\right)T_{RR_g} + \frac{1}{t_s}T_{R_g R_g}, \qquad (7.10)$$

$$-\frac{2\Omega}{R}T_{\phi\phi_g} - \Omega R(T_{RR} - T_{RR_g}) + 2R(\Omega - A)T_{RR} = -\left(\frac{1}{t_s} + \frac{1}{t_c}\right)T_{R\phi_g} + \frac{1}{t_s}T_{R_g\phi_g}, \qquad (7.11)$$

$$-\frac{\Omega}{R}(T_{\phi\phi} + T_{\phi\phi_g}) + 2R(\Omega - A)T_{RR_g} = -\left(\frac{1}{t_s} + \frac{1}{t_c}\right)T_{R_g\phi} + \frac{1}{t_s}T_{R_g\phi_g}, \qquad (7.12)$$

$$2R(\Omega - A)(T_{R\phi} + T_{R\phi_g}) + R\Omega(T_{R_g\phi} - T_{R\phi}) = -\left(\frac{1}{t_s} + \frac{1}{t_c}\right)T_{\phi\phi_g} + \frac{1}{t_s}T_{\phi_g\phi_g}. \qquad (7.13)$$

Finally the “gas-gas” components of Equation 7.6, which describe the behaviour the (3D) gas Reynolds stress along the dust trajectory, are
<figure><img src="image_2.png" /><figcaption>Figure 1: Horizontal stress tensor components in a rotating shear flow, as a function of St and A/Ω, with different dimensionless correlation times</figcaption></figure>

$$ -\frac{2\Omega}{R}(T_{R_g}\phi_g + T_{R_g}\phi) = -\frac{2}{t_c}(T_{R_g R_g} - \alpha c_s^2 \rho_d), \quad (7.14) $$

$$ -\frac{\Omega}{R}(T_{\phi\phi_g} + T_{\phi_g\phi_g}) - \Omega R(T_{RR_g} - T_{R_g R_g}) + 2R(\Omega - A)T_{RR_g} = -\frac{2}{t_c}T_{R_g\phi_g}, \quad (7.15) $$

$$ 2R\Omega T_{R_g\phi_g} + 2R(\Omega - 2A)T_{R\phi_g} = -\frac{2}{t_c}(T_{\phi_g\phi_g} - \alpha c_s^2 \rho_d R^2). \quad (7.16) $$

It is straightforward, if rather laborious, to invert the above equations. However, the resulting expressions are somewhat cumbersome and not particularly informative. We shall instead use a symbolic algebra package to obtain expressions for the pressure tensor components that can be used in numerical computations (we provide a Mathematica script to do this in the supplementary materials). We can then numerically explore the behaviour of the dust pressure tensor.

Figure 1 shows how the horizontal stress tensor components change with A/Ω and St for different dimensionless correlation times. Figure 2 shows the locations in the parameter space where stress tensor components pass through zero - indicating that the rotating shear flow no longer posses a steady solution. Together these shows the general behaviour of the dust stress tensor. This shows that the dust stress tends to get weaker at larger Stokes and tends to isotropy in the absence of shear A/Ω → 0. There are singularities/zeros of the pressure
<figure><img src="image_3.png" /><figcaption>Figure 2: Locations in the parameter space of zeros of the stress tensor components, indicating a lack of steady solutions. Black Solid: $\tau_c = 10^{-3}$, Black Dashed: $\tau_c = 10^{-2}$, Black Dotted: $\tau_c = 10^{-2}$, Blue Solid: $\tau_c = 1$, Blue Dashed: $\tau_c = 10^3$. The region of the parameter space above these lines contains no (physical) steady solutions.</figcaption></figure>

tensor at large $A/\Omega$ and St associated with the breakdown of the fluid dust description. For $\tau_c \le 1$ these asymptote to the Rayleigh stability criterion for large St, for small Stokes dust drag/cooling helps to regularise the behaviour of the pressure tensor allow for steady solutions at higher $A/\Omega$. For large $\tau_c$ small Stokes numbers are less effective at regularising the behaviour of the stress tensor and we see a zero of the stress tensor at $A/\Omega \sim 2$ for small Stokes. This occurs because of a breakdown of the turbulence model for large $\tau_c$ and $A/\Omega$ (See Appendix A.2 and A.3).

Figure 3 shows the stress tensor components for different Stokes numbers in a rotating shear flow with $\tau_c = 1$. The left plot shows the Rayleigh stable Keplerian shear flow with $A = \frac{3}{4}\Omega$ (The Reynolds stress of the gas associated with this flow is shown in the left hand plot of Figure 7 in Appendix A.3). For small Stokes numbers the tight coupling to the gas means the dust stress is set by the velocity correlations in the gas. As the Stokes number increases there is a competition between the isotropising effect of the turbulence and the shearing out of the radial component of the stress tensor, leading to an increasingly anisotropic flow.

The right hand plot shows a Rayleigh unstable shear flow with $A = 1.1\Omega$ (The Reynolds stress of the gas is shown in the right hand plot of Figure 7). Again at low Stokes the dust stress is set by the gas velocity correlations. The stress tensor components diverge as the Stokes number increases, and one approaches the breaks down in the fluid description, before vanishing indicating a lack of accessible steady flow solution.

Figure 4 shows how the speed of the P-wave varies with Stokes number and direction. S-waves velocities are the same as those of the P-waves but rescaled by $1/\sqrt{3}$. As the dust rheological stress becomes increasingly anisotropic the P-waves and S-waves propagate more radially than azimuthally. In the right hand plot of the Rayleigh unstable flow the P (and S) waves cannot propagate at large enough Stokes numbers - further increasing the Stokes number leads to a breakdown of the fluid description.

## 7.2. Accretion flow solutions

Section 7.1 gives the effects of the leading order velocity field on the pressure tensor in a rotating shear flow and neglects the presence of an accretion flow driven by the $R_{r\phi}$ Reynolds stress and the effects of gas pressure gradients. To study this effect we implement a 1D hydro-solver to solve the dust-fluid equations in aligned-cylindrical coordinates. We consider an axisymmetric dust flow around a Keplerian gravitational potential, neglecting the effects of vertical gravity.

For the gas properties we solve Equations A 10-A 12 perturbatively, with the leading order terms being that due to gravity and circular-Keplerian motion. We then consider the 1st order correction to the gas velocity due to the gas pressure and turbulence. This has the
<figure><img src="image_4.png" /></figure>

Figure 3: Horizontal stress tensor components for a fluid with a dimensionless correlation time, $\tau_c = 1$. Left: A Rayleigh stable Keplerian shear flow where $\Omega \propto R^{-3/2}$ ($A = \frac{3}{4}\Omega$); Right: Rayleigh unstable shear flow with $A = 1.1\Omega$. The gas Reynolds stress for these two cases is shown in Figure 7 in Appendix A.3. For the Keplerian shear flow the rheological stress tensor becomes increasingly anisotropic as the Stokes number increases. In the Rayleigh unstable case there is a maximum Stokes number above which the fluid dust description breaks down as the constitutive equation becomes thermally unstable.

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

Figure 4: Illustration of the directional dependence of the P-wave velocity at different Stokes numbers for the disc considered in Figure 3. The distance from the origin is proportional to the wavespeed of the P-wave, while the angle to the x-axis is the direction of propagation. The axes are aligned such that the x-axis is directed in the radial direction, while the y-axis points in the direction of the fluid motion. On this figure the sound wave in the gas would be a circle of unit radius. For the Keplerian shear flow the P-waves slow down and preferentially propagate radially as the Stokes number increases. For the Rayleigh unstable flow the P-wave velocity is highly anisotropic, up to the breakdown in the fluid description where there is no longer a steady background on which the P-wave can propagate.

effect of driving a slow, radial accretion flow and making the gas rotation sub-Keplerian in the presence of a negative pressure gradient. We consider constant $\alpha$ and sound speed throughout.

As a starting point we consider the case of a constant gas surface density, and neglect the slow change in the gas density due to accretion. This is not self-consistent as the timescale on which the gas density evolves is typically expected to be comparable to the dust sound crossing time of interest. We also consider a more realistic example, where we solve for the steady-state, turbulent gas profile, resulting in a gas surface density of $\rho_g \propto R^{-3/2}$. Further details on the gas-flow considered are given in Appendix A.4.

We solve the dust fluid equations in aligned-cylindrical coordinates using an implementation
<figure><img src="image_6.png" /><figcaption>Figure 5: Density plots for the rotating shearflow, (left) constant density flow, (right) steady state accretion flow. Both are for Keplerian shear flows, solved with the HLL solver with $\tau_c = 0.1$. The dashed line indicates the initial density profile, while the solid black line shows the final density profile. Grey lines show the evolution of the density profile, spaced every 100 inner orbits (left), and every 200 inner orbits (right).</figcaption></figure>

tation of the HLL (Harten et al. 1983) and Roe (Roe 1981) approximate-Reimann solvers, and constant reconstruction. We use an operator-split Van-Leer integrator (Van Leer 2006; Eleuterio 1999), and use an RK(2)3 integrator to integrate the source terms.

We take units such that the radius of the inner boundary is 1 and $GM = 1$, resulting in the circular Keplerian frequency on the inner boundary also being unity. The domain spans $R \in [1, 5]$. In these units we consider a gas disc with constant sound speed $c=0.2$, and turbulence with $\alpha = 0.02$ and $\tau_c = 0.1$ or $\tau_c = 0.01$. This is adopted for computational convenience (principally difficulties with ensuring positivity of the stress tensor and so that the dust sound crossing time is not too long), and is not reflective of realistic disc turbulence (particularly for dust hosting discs). We consider dust with stopping time $t_s = 0.01$ at the reference orbit $R = 1$ and reference gas density $\rho_g = 1$. For the constant gas density case we start with a constant dust density $\rho_d = 1$, while for the steady state we start with a step profile,

$$ \rho_d = 0.1 + 0.9[\tanh(2R + 1) + 1]. \qquad (7.17) $$

In both cases we take the initial dust velocity to be equal to the gas velocity and an isotropic dust stress with $T_{\alpha\beta} = \alpha c_s^2 g_{\alpha\beta}$. Notably this initial dust stress does not correspond to the anisotropic stress expected in a steady rotating shear flow (as shown in Section 7.1). We adopt zero-gradient boundary conditions with a diode inner boundary for $U^R$, i.e. we set $U^R(1) = 0$ whenever it is positive thereby only allowing an outflow on the inner boundary. We add wavekilling zones to our simulations, applying a large artificial viscosity near the boundary, which decreases to zero within a distance of 0.2 of the boundary. This is done to stop grid scale oscillations being excited by the boundary, particularly when using the Roe solver. For the constant gas density case we integrate for 1000 inner orbits, while for the steady state accretion flow we integrate for 3000 inner orbits. In both cases this does not reach a steady state due to the long relaxation time in the outer disc.

Figure 5 shows the density of both cases at different times, integrated with the HLL solver for $\tau_c = 0.1$. In the steady state case the initial dust step in the outer disc has drifted in due to the drag from the sub-Keplerian gas, before the dust density start to relax towards the steady state for the induced drift velocity Figure 6 shows the dust stress tensor for the final output of both these simulations, compared against the steady state solutions of Section 7.1. The
<figure><img src="image_7.png" /><figcaption>Figure 6: Diagonal dust-stress components for the flows considered in Figure 5. Full colour lines are the simulation, while greyed out/transparent lines are the solutions from Section 7.1. These are closer in the inner disc where the velocity is closer to Keplerian motion. The outer disc in the constant gas density case is still far from the steady state profile.</figcaption></figure>

Roe solver is not able to reach as large a correlation time as the HLL solver (which is more diffusive), however we carry out the same simulations with a correlation time of $\tau_c = 0.01$ where we find agreement between the two solvers.

The simulations exhibit several of the expected feature of the dust fluid in an accretion disc. In the steady accretion flow, Figure 5 shows the initial step profile drifting inwards due to the action of gas drag with the sub-Keplerian gas flow. This is well established behaviour for dust in accretion flows and occurs even in the absent of dust-pressure. For the constant gas density case the dust is dragged inwards by the accretion flow and builds up on the boundary. This appears to be due to the adopted boundary conditions being partially permeable to the dust, with the zero-gradient in the radial velocity resulting in a lower outflow velocity than expected for a continuation of the accretion flow. The presence of this partial obstruction is then communicated into the disc by diffusion of the dust due to the gas turbulence. Preliminary tests simulating the dust fluid in accretion flows with a gas-pressure maximum suggests, as expected, the gas pressure maximum is no longer a dust trap. This differs from the behaviour found in pressureless dust models, with the dust now able to reach the inner boundary given enough time. This effect will have major implications for the transport of solids in protoplanetary discs and needs to be explored more thoroughly in the future.

The numerical implementation of the dust-fluid solver, presented here, is far from a practical implementation and requires numerous improvements to be useful. For most purposes the HLL solver is too diffusive and one should either further develop the Roe solver to be more stable at large correlation times or implement a HLLC (Toro et al. 1994) type solver. Currently the solver struggles to maintain positivity of the stress tensor, particularly $T_{\phi\phi}$, at larger correlation times, with the correlation time reachable by the HLL solver being $\tau_c = 0.1$, rather than the $\tau_c = 1-5$ expected of realistic disc turbulence. One possible reason for this is using the total second velocity moment, $T_{\alpha\beta} + \rho_d U_\alpha U_\beta$ as a conservative variable. This leads to errors due to the large difference in scale between the orbital motion, $\rho_d U_\phi U_\phi$ and the dust stress tensor component $T_{\phi\phi}$. On obvious improvement would be to implement a FARGO (Masset 2000) like advection scheme and subtract off the orbital motion, to reduce the error on $T_{\phi\phi}$. Is it also possible that a finite difference scheme will perform better than
the finite-volume Reimann solvers employed here when the orbital motion is treated in this manor (Benítez-Llambay & Masset 2016). Finally much work needs to be done to determine appropriate boundary conditions, as the current zero-gradient boundaries excite grid scale oscillations that are currently dealt with by adding viscosity close to the boundary. This probably indicates that such a choice of boundary are ill posed in the dusty, rotating shear flow setting.

## 8. Discussion

In this paper we have chosen not to include the back reaction of the dust on the gas, despite it's importance even in dust poor flows. One can include back reaction in an ad-hoc manner by adding a dust drag term to $f_i^g$,

$$f_i^g = -\nabla_i \phi - \frac{1}{\rho_g} \nabla_i \rho_g - \frac{1}{t_s} f_d (u_i^g - u_i^d), \qquad (8.1)$$

where $u_i^d$ is the mean dust velocity and $f_d = \rho_d/\rho_g$ is the dust to gas density ratio. Here the back reaction on the gas is only between the mean velocities and does not account for the effect of the dust on the gas turbulence, through corresponding source terms in the evolutionary equation for the gas Reynolds stress. A more self consistent approach would be to allow back reaction with the stochastic gas and dust velocities, which leads to the obvious improvement on Equation 8.1 by replacing the mean flow back reaction term with $-\frac{1}{t_s} f_d (v_i^g - v_i^d)$ (e.g. as done in Minier & Peirano 2001; Minier et al. 2004; Minier 2015). However, $v^g$ is the velocity of a fluid element (seen), while $v^d$ is the velocity of an individual dust particle, and one expects there to be multiple dust particles within a given fluid element - it is thus not clear whether this modification is self consistent (see also the discussion in Minier & Peirano 2001).

As above, in addition to adding in the dust drag, dust loading can affect the fluid turbulence. An alternative way to include this effect is make $t_c$ and $\alpha$ dependant on collective dust properties, the most important effect likely being a dependence on the dust to gas ratio $f_d = \rho_d/\rho_g$. It would, thus, be useful to have a more rigorous treatment of back reaction, this would also be important for ensuring total energy in the gas+dust system is conserved (particularly to ensure the turbulence doesn't act as an infinite source of energy).

A rigorous treatment of energy conservation allows for energy to be exchanged between the gas turbulence, gas thermal energy and dust mechanical potential energy, along with the mean flow of both phases. This is particularly important in the compressible setting as the pressure is dynamical, and affected by the gas temperature, rather than being a Lagrange multiplier enforcing incompressibility (the coupling has been considered in the incompressible setting e.g. by Fox 2014). This coupling naturally leads to the damping of the gas turbulence due to dust loading as energy is transformed from the turbulent fluctuations into heating the dust and is then transformed to the gas thermal energy due to gas drag leading to turbulent gas motions being converted to heat on $\sim t_s f_d^{-1}$. Finally this more complete modelling of the energy exchange between the three stores of energy may allow for more complex phenomena like intermittence and predator-prey dynamics which are known to be important in many turbulence processes (e.g. Diamond et al. 1994).

Our model considers a situation where dust-dust collisions are rare due to the low number density of the dust relative to the gas. As the dust density increases, dust-dust collisions can become important, particularly for small grains. Inclusion of dust-dust collisions would allow contact between the dust fluid theory and existing work on dust pressure in weakly collisionless dust systems (Goldreich & Tremaine 1978; Borderies et al. 1985; Araki &
Tremaine 1986; Latter & Ogilvie 2008; Larue et al. 2023). The inclusion of dust collisions in a SDE model for particle laden flows, and associate moment closure, has been studied in Innocenti et al. (2019, 2021) and Capecelatro et al. (2016a,b). This includes a separation of the dust Reynolds stress from the pressure tensor, which is important when dust-dust collisions are included as the dust collision velocity is principally sensitive to the particle velocity dispersion rather than the turbulent velocity dispersion (Fox 2014; Capecelatro et al. 2016b). Gas kinetic effects can also be important for smaller dust grains in regions of low gas density, where finite Knudsen number effects become important. Here dust gas collisions are infrequent enough, and impart sufficient momentum on the dust grain, that are an additional source of stochasticity on top of the gas turbulence which will act to heat the dust. Both dust-dust collisions and kinetic gas effects are likely important in systems with very large dust to gas ratios - particularly when the gas is produced by sublimating/colliding dust.

Finally more sophisticated models of gas turbulence's (e.g. Sawford 1991; Pope 2002) include two timescales (correlation time and Kolmogorov time) and maybe used to derive finite Reynolds number effects (formally our model is for turbulence with an infinite Reynolds number). As Reynolds numbers in astrophysical (and geophysical) gasses are very large this effect is likely only important in a limited region of parameter space - but maybe needed to obtain the correct collisional velocity for small grains (e.g. to reproduce the results of Ormel & Cuzzi 2007), or for the experimental verification of the model.

## 9. Conclusion

In this paper we have derived a fluid model for collisionless dust in a turbulent gas, starting from a system of stochastic differential equations describing the motion of a single dust grain. To allow for the coordinate systems and geometries common to astrophysics, we have adopted a covariant form for our dust-fluid model. We show that the continuum mechanics properties of dust in a turbulent gas corresponds to a 6-dimensional anisotropic Maxwell fluid with a dynamically important rheological stress tensor. The 6-dimensional formulation keeps the dust and fluid seen velocities, and their respective moments, on the same footing. The coupling between the dust kinetic tensor, dust-gas cross pressure and fluid seen Reynolds stress are obtained from the advection of the 6-dimensional dust stress tensor by the fluid flow.

In summary our conclusions are:

* We have developed a dust fluid model, using a closure valid in the accretion disc context, and demonstrated that the self-consistency of the moment truncation used to obtain the fluid description is closely related to the thermals stability of the fluid.

* Collisionless dust has a non-zero anisotropic rheological stress which can be dynamically important, such as in dusty atmospheres where the dust is in hydrostatic equilibrium between the dust stress, vertical gravity and the gas Reynolds stress.

* The dust can support seismic (P and S) waves

* Whether velocity correlations of small dust grains are set by the gas velocity correlations is determined by a form of Eddy-Knudsen number. Which can lead to small dust grains not being well mixed with the gas.

We have suggested several potential extension to our model some of which we intend to pursue in future work.

**Acknowledgements.** The authors wish to thank Richard Booth and Michaël Bourgoin for invaluable help with the literature along with Henrik Latter, Tobias Heinamann, Geoffroy Lesur and Francesco Lovascio for discussions which greatly clarified the physical interpretation of our model. We also thank Cathie Clarke, Andrew Sellek and the group of the CRAL for useful discussions on this work. We thank all three referees for a thorough review, bringing several suggestions that put the manuscript into its present form.
**Funding.** The authors would like to thank the European Research Council (ERC). This research was supported by the ERC through the CoG project PODCAST No 864965. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 823823.

**Declaration of interests.** The authors report no conflict of interest.

**Data availability statement.** The data underlying this article will be shared on reasonable request to the corresponding author.

# Appendix A. Model for the gas

## A.1. Formulation

In our model for the turbulent gas an individual fluid element evolves according to the following set of stochastic differential equations,

$$d\mathbf{x} = \mathbf{v}dt, \qquad (A1)$$

$$d\mathbf{v} = \mathbf{f}_g dt - \frac{1}{t_c}(\mathbf{v} - \mathbf{u})dt + \sqrt{\frac{2\alpha}{t_c}}c_s d\mathbf{W}, \qquad (A2)$$

where $(\mathbf{x}, \mathbf{v})$ are the position and velocity of the fluid elements, $\mathbf{F}$ is the force per unit mass on the gas in the absence of turbulence and the turbulence results in an Ornstein-Uhlenbeck walk around the mean fluid flow with correlation time $t_c$. $\mathbf{W}$ is a Wiener process, with $c_s$ the gas sound speed and $\alpha$ a dimensionless measure of the strength of the turbulence. The fluid flow is a member of a statistical ensemble of similar flows with each fluid element following a single realisation of the flow (Pope 1985; Thomson 1987).

The Fokker-Planck equation associated with these equations can be derived in a similar way to that of the dust,

$$\frac{\partial p}{\partial t} + \frac{\partial}{\partial x^i} [v^i p] + \frac{\partial}{\partial v^i} \left[ p f_g^i - \frac{1}{t_c} (v^i - u^i) p \right] = \frac{\alpha c_s^2}{t_c} \frac{\partial^2 p}{\partial v^2}. \qquad (A3)$$

Taking the 0th, 1st and 2nd velocity moments of this equation we arrive at

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x^i} [u_i \rho] = 0, \qquad (A4)$$

$$\frac{\partial}{\partial t} [\rho u_i] + \frac{\partial}{\partial x^j} [R_{ij} + \rho u_i u_j] - \rho f_i^g = 0, \qquad (A5)$$

$$\frac{\partial}{\partial t} [R_{ij} + \rho u_i u_j] + \frac{\partial}{\partial x^k} [R_{ijk} + 3u_{(i}R_{jk)} + \rho u_i u_j u_k] - 2\rho u_{(i}f_{j)}^g + \frac{2}{t_c} R_{ij} = 2 \frac{\alpha c_s^2}{t_c} \rho g_{ij}. \qquad (A6)$$

These can be rearranged to obtain

$$D\rho = -\rho \nabla_i u^i, \qquad (A7)$$

$$\rho D u_i = -\nabla^j R_{ij} + \rho f_i^g, \qquad (A8)$$

$$(D + \nabla_k u^k) R_{ij} + 2R_{k(i} \nabla^k u_{j)} = -\nabla^k R_{ijk} - \frac{2}{t_c} (R_{ij} - \alpha c_s^2 \rho g_{ij}), \qquad (A9)$$

where, in this Appendix, $D = \partial_t + u^i \nabla_i$ is the Lagrangian derivative with respect to the
gas flow. Our closure scheme for this model assumes $R_{ijk} = 0$. We shall show, in the next section, this can be justified on the basis of a near-Maxwellian ordering scheme for the velocity moments, similar to the dust. Finally using a similar argument to that presented in Appendix B.1, for the dust, we can show that (when $R_{ijk} = 0$) $R_{ij}$ is positive semi-definite for positive semidefinite initial conditions. Thus the equations to be solved for the gas phase are

$$D\rho = -\rho\nabla_i u^i, \qquad (A 10)$$

$$\rho Du_i = -\nabla^j R_{ij} + \rho f_i^g, \qquad (A 11)$$

$$(D + \nabla_k u^k) R_{ij} + 2R_{k(i} \nabla^k u_{j)} = -\frac{2}{t_c} (R_{ij} - \alpha c_s^2 \rho g_{ij}). \qquad (A 12)$$

Equivalently, one can perform a Reynolds decomposition of this flow resulting in the equivalent set of equations, which are closer to the formulation of Thomson (1987),

$$d\mathbf{x} = (\mathbf{v}_t + \mathbf{u})dt, \qquad (A 13)$$

$$d\mathbf{v}_t = -\frac{1}{t_c}(\mathbf{v}_t - \mathbf{v}_{hs})dt + \sqrt{\frac{2\alpha}{t_c}}c_s d\mathbf{W}, \qquad (A 14)$$

where the total gas velocity is equal to the sum of the mean velocity **u** and turbulent velocity **v**$_t$, $\mathbf{v} = \mathbf{u} + \mathbf{v}_t$. The mean velocity obeys the usual Reynolds averaged equation

$$\rho D \mathbf{u} = \rho \mathbf{f}_g - \nabla \cdot \mathbf{R}, \qquad (A 15)$$

and we have introduced $\mathbf{v}_{hs} = \frac{t_c}{\rho} \nabla \cdot \mathbf{R} - t_c \mathbf{v}_t \cdot \nabla \mathbf{u}$. Minier et al. (2014) has shown that these two formulations are equivalent.

## A.2. Justification of the closure scheme

In this section we will determining a closure for the gas phase in our model. This will exploit the separation of scales between the hypersonic background motion and the highly subsonic turbulence and show that there exist a well defined asymptotic scaling in which the departure from an anisotropic Maxwellian velocity distribution is small. This is similar to the near-Maxwellian ordering scheme for the dust fluid considered in Section 4.3.2. As we did in the dust fluid one can obtain an evolutionary equation for the k-th velocity moment of the gas turbulence,

$$\mathcal{D}_1 R_{i_1 \dots i_k} + k R_j (i_1 \dots i_{k-1} \nabla^j u_{i_k}) = - \nabla^j R_{i_1 \dots i_k j} + \frac{k}{\rho} R_{(i_1 \dots i_{k-1} \nabla^j R_{i_k}) j} \\ - \frac{k}{t_c} (R_{i_1 \dots i_k} - (k-1) \alpha c_s^2 R_{(i_1 \dots i_{k-2} g_{i_{k-1} i_k})}), \qquad (A 16)$$

where we have introduced the differential operator $\mathcal{D}_1 = D + \nabla_i u^i$.

Consider a high Mach-number gas flow with subsonic turbulence and introduce two (potentially) small parameters $\delta$, which is of order $1/M$ with $M$ the Mach number, and $\alpha < 1$ which is a measure of the strength of the turbulence. We introduce two lengthscales a long lengthscale $L = O(1)$ (with long lengthscale variable $\mathbf{x}$), and ‘short’ lengthscale $l = O(\delta)$ (with short lengthscale variable $\xi$). To leading order the gas has a Maxwellian velocity distribution where the mean velocity is $O(1)$ and the gas sound speed is $O(\delta)$. This
ordering scheme is subtly different to the near-Maxwellian ordering scheme of the dust fluid as we generally have $L \le l \le L_{\text{dust}}$. The small parameter $\epsilon$ in the dust fluid problem is $O(\alpha^{1/2}\delta)$ in the gas ordering scheme. Our ordering scheme will be valid provided that the turbulent velocities are small compared with the typical fluid velocities and correspond to $\epsilon \ll 1$. This can either be due to the flow having a high Mach-number ($\delta \ll 1$), as is typical in astrophysics, or when the turbulence is weak ($\alpha \ll 1$).

At leading order we consider a gas with a Maxwellian velocity distribution which varies on the long lengthscale $L$, but having a gas density that can vary on the short lengthscale $l$. At higher order the distribution function has a non-Maxwellian velocity component, which is allowed to vary on the short lengthscale. The nearly Maxwellian asymptotic scheme is

$$R_{i_1 \cdots i_k} = \alpha^{\text{ceil}(k/2)} \delta^k \rho(\xi, \mathbf{x}) W_{i_1 \cdots i_k}(\mathbf{x}) + \alpha^{\text{ceil}((k+1)/2)} \delta^{k+1} \Sigma_{i_1 \cdots i_k}(\xi, \mathbf{x}), \quad (A 17)$$

$$\mathbf{u} = \mathbf{u}_0(\mathbf{x}) + \alpha^2 \delta^2 \mathbf{u}_1(\xi, \mathbf{x}), \quad (A 18)$$

$$\nabla = \delta^{-1} \frac{\partial}{\partial \xi} + \frac{\partial}{\partial \mathbf{x}}, \quad (A 19)$$

where $W_{i_1 \cdots i_k}$ are the Maxwellian velocity correlations and have the same properties as to there dust counterparts and evolve according to

$$D W_{i_1 \cdots i_k} = -k \left[ W^{j}_{(i_1 \cdots i_{k-1})} B_{i_k,j} - (k-1) W_{(i_1 \cdots i_{k-2})} D_{i_{k-1} i_k} \right]. \quad (A 20)$$

As we did with the dust, we can absorb perturbations to the gas density, from the non-Maxwellian terms, into the definition of $\Sigma_{i_1 \cdots i_k}$.

Substituting the ordering scheme (Equations A 17-A 19) into Equation A 16 we arrive at

$$\begin{aligned} \mathcal{D}_1 R_{i_1 \cdots i_k} &= \alpha^{\text{ceil}(k/2)} \delta^k \rho D W_{i_1 \cdots i_k} + \alpha^{\text{ceil}(k/2)+2} \delta^{k+2} u_1^j \frac{\partial}{\partial x^j} W_{i_1 \cdots i_k} + \alpha^{\text{ceil}((k+1)/2)} \delta^{k+1} \mathcal{D}_1 \Sigma_{i_1 \cdots i_k} \\ &= -\alpha^{\text{ceil}((k+1)/2)} \delta^k \left( \frac{\partial}{\partial \xi^j} + \delta \frac{\partial}{\partial x^j} \right) (\rho W_{i_1 \cdots i_k j}) - \alpha^{\text{ceil}(k/2)+1} \delta^{k+1} \left( \frac{\partial}{\partial \xi^j} + \delta \frac{\partial}{\partial x^j} \right) \Sigma_{i_1 \cdots i_k j} \\ &\quad - k \delta^k \left[ \alpha^{\text{ceil}((k+1)/2)} \rho W_{i_1 \cdots i_k} + \alpha^{\text{ceil}(k/2)+1} \delta \Sigma_{i_1 \cdots i_k}^\sigma \right. \\ &\quad \left. - \left( \alpha^{\text{ceil}((k+1)/2)} W_{(i_1 \cdots i_{k-1})} + \alpha^{\text{ceil}(k/2)+1} \delta \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \right) \left( \frac{\partial}{\partial \xi^j} + \delta \frac{\partial}{\partial x^j} \right) (\rho W_{i_k})_j + \alpha \delta \Sigma_{i_k})_j \right. \\ &\quad \left. - (k-1) \left( \alpha^{\text{ceil}(k/2)} \rho W_{(i_1 \cdots i_{k-2})} + \alpha^{\text{ceil}((k+1)/2)} \delta \Sigma_{(i_1 \cdots i_{k-2})} \right) g_{i_{k-1} i_k} \right], \end{aligned} \quad (A 21)$$

where, here, $D = \partial_t + u_0^i \nabla_i$ and we have made use of $\mathcal{D}_1 \rho_g = (D + \nabla_i u_0^i) \rho_g = 0$.

Making use of Equation A 20 for the evolution of $W_{i_1 \cdots i_k}$, along with the recurrence relation for $W_{i_1 \cdots i_k}$ and rearranging we obtain and equation for the evolution of the non-Maxwellian part of the turbulent velocity moment,
$$
\begin{align*}
& \alpha^{\text{ceil}((k+1)/2)} \mathcal{D}_1 \Sigma_{i_1 \cdots i_k} + \alpha^{\text{ceil}(k/2)+1} \frac{\partial}{\partial \xi^j} \Sigma_{i_1 \cdots i_k} j + k \left[ \begin{aligned}
& \alpha^{\text{ceil}(k/2)+1} \Sigma_{i_1 \cdots i_k}^{\sigma} \\
& - \alpha^{\text{ceil}((k+1)/2)+1} W_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial \xi^j} \Sigma_{i_k} j + \alpha^{\text{ceil}((k+1)/2)} \rho W_{j(i_k)} \frac{\partial}{\partial x^j} W_{i_1 \cdots i_{k-1}} \\
& - \alpha^{\text{ceil}(k/2)+1} \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial \xi^j} (\rho W_{i_k}) j - (k-1) \alpha^{\text{ceil}((k+1)/2)} \Sigma_{(i_1 \cdots i_{k-2})} g_{i_{k-1} i_k}
\end{aligned} \right] \\
& = \delta \left[ k \alpha^{\text{ceil}((k+1)/2)+1} W_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial x^j} \Sigma_{i_k} j + k \alpha^{\text{ceil}(k/2)+1} \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial x^j} (\rho W_{i_k}) j \right. \\
& \qquad \left. + k \alpha^{\text{ceil}(k/2)+2} \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \left( \frac{\partial}{\partial \xi^j} + \delta \frac{\partial}{\partial x^j} \right) \Sigma_{i_k} j \right. \\
& \qquad \left. + \alpha^{\text{ceil}(k/2)+2} u_1^j \frac{\partial}{\partial x^j} W_{i_1 \cdots i_k} + \alpha^{\text{ceil}(k/2)+1} \frac{\partial}{\partial x^\sigma} \Sigma_{\alpha_1 \cdots \alpha_k \sigma} \right]. \tag{A 22}
\end{align*}
$$

Keeping only leading order terms in $\delta$ then, for even $k$, we have

$$
\begin{align*}
& \mathcal{D}_1 \Sigma_{i_1 \cdots i_k} + \frac{\partial}{\partial \xi^j} \Sigma_{i_1 \cdots i_k} j + k \left[ \Sigma_{i_1 \cdots i_k}^{\sigma} - \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial \xi^j} (\rho W_{i_k}) j - (k-1) \Sigma_{(i_1 \cdots i_{k-2})} g_{i_{k-1} i_k} \right] \\
& - \delta \left[ k \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial x^j} (\rho W_{i_k}) j + \frac{\partial}{\partial x^\sigma} \Sigma_{\alpha_1 \cdots \alpha_k \sigma} \right] \\
& = \delta \alpha \left[ k W_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial x^j} \Sigma_{i_k} j + k \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \left( \frac{\partial}{\partial \xi^j} + \delta \frac{\partial}{\partial x^j} \right) \Sigma_{i_k} j + u_1^j \frac{\partial}{\partial x^j} W_{i_1 \cdots i_k} \right], \tag{A 23}
\end{align*}
$$

while for odd $k$, we have

$$
\begin{align*}
& \mathcal{D}_1 \Sigma_{i_1 \cdots i_k} + \alpha \frac{\partial}{\partial \xi^j} \Sigma_{i_1 \cdots i_k} j + k \left[ \begin{aligned}
& \alpha \Sigma_{i_1 \cdots i_k}^{\sigma} - \alpha W_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial \xi^j} \Sigma_{i_k} j + \rho W_{j(i_k)} \frac{\partial}{\partial x^j} W_{i_1 \cdots i_{k-1}} \\
& - \alpha \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial \xi^j} (\rho W_{i_k}) j - (k-1) \Sigma_{(i_1 \cdots i_{k-2})} g_{i_{k-1} i_k}
\end{aligned} \right] \\
& = \delta \alpha \left[ k \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \frac{\partial}{\partial x^j} (\rho W_{i_k}) j + k \alpha \rho^{-1} \Sigma_{(i_1 \cdots i_{k-1})} \left( \frac{\partial}{\partial \xi^j} + \delta \frac{\partial}{\partial x^j} \right) \Sigma_{i_k} j + \frac{\partial}{\partial x^\sigma} \Sigma_{\alpha_1 \cdots \alpha_k \sigma} \right]. \tag{A 24}
\end{align*}
$$

In both cases the right-hand side can be dropped at leading order if either $\alpha$ or $\delta$ are small. This confirms that the asymptotic ordering scheme (Equations A 17-A 19) is self-consistent and the non-Maxwellian terms are suppressed by a factor of $\sim \alpha/M$ relative to the Maxwellian terms. For the purposes of the effect on the gas equation of motion one must consider the effect on the Reynolds stress gradients. For the nearly Maxwellian velocity distribution considered the gradients of the Reynolds Stress are
$$ \nabla_j R^{ij} = \alpha \delta \left( \frac{\partial}{\partial \xi^j} + \frac{\partial}{\partial x^j} \right) \rho W^{ij} + \alpha^2 \delta^2 \left( \frac{\partial}{\partial \xi^j} + \frac{\partial}{\partial x^j} \right) \Sigma_{ij}. \quad (A 25) $$

Thus the effects of the non-Maxwellian terms are $O(\alpha^2\delta^2)$ and are thus small relative to the acceleration and gravity, which are taken to be $O(1)$, or the pressure gradients which are $O(\delta)$

As with the dust fluid, the existence of a consistent asymptotic scaling does not guarantee that the it is an attractor. We shall not repeat the argument here, but one can follow a similar line of reasoning to Section 4.3.3 to demonstrate that turbulent velocity moments which start far from the asymptotic scaling are expected to damp towards the scaling, subject to the equation governing the evolution of the Reynolds stress having a stable equilibrium. As with the dust fluid, this is not sufficient to completely show that the near-Maxwellian ordering is an attractor as it does not account for the possibility of (nonlinear) perturbations to multiple velocity moments mutually supporting each other against decay. As with the dust fluid a more complete analysis of when the near-Maxwellian ordering scheme acts as an attractor must be left for future work.

### A.3. Reynolds stress in a rotating shear flow

In this section we shall derive the steady state Reynolds stress in a rotating shear flow. This will aid our discussion of the behaviour of the dust rheological stress in a rotating shear flow in Section 7 along with illustrating some key properties of our turbulence model. The steady state Reynolds stress in the gas satisfies the following equation,

$$ \nabla_k (u^k R_{ij}) + 2R_{k(i} \nabla^k u_{j)} = -\frac{2}{t_c} (R_{ij} - \alpha c_s^2 \rho_g g_{ij}). \quad (A 26) $$

Adopting cylindrical polar coordinates $(R, \phi, z)$, with the usual expressions for $g_{ij}$ and $\Gamma_{ij}^k$. We consider a rotating shear flow with velocity, $u^k = \Omega(R)\hat{e}_\phi$. Substituting this into Equation A 26, and neglecting gradients in $R_{ij}$ at leading order, we have

$$ -2\Omega\Gamma_{\phi(i}^s R_{j)s} + 4R(\Omega - A)R^R_{(i}\hat{e}_{j)\phi} - 2R^k_{(i}\Gamma_{j)k}^\phi R^2\Omega = -\frac{2}{t_c}(R_{ij} - \alpha c_s^2 \rho_g g_{ij}), \quad (A 27) $$

where we have in introduced Oort's constant $A = -\frac{R}{2}\Omega_R$. Explicitly the components of the Reynolds stress equations are

$$ -\frac{4\Omega}{R}R_{R\phi} = -\frac{2}{t_c}(R_{RR} - \alpha c_s^2 \rho_g), \quad (A 28) $$

$$ -\frac{2\Omega}{R}R_{\phi\phi} + 2R(\Omega - A)R_{RR} = -\frac{2}{t_c}R_{R\phi}, \quad (A 29) $$

$$ 4R(\Omega - A)R_{R\phi} = -\frac{2}{t_c}(R_{\phi\phi} - \alpha c_s^2 \rho_g R^2). \quad (A 30) $$

We can rearrange these to obtain the following expressions for the Reynolds stress components in the rotating shear flow,
<figure><img src="image_8.png" /><figcaption>Figure 7: Reynolds stress components (in units of $\alpha c_s^2$) for two different rotating shear flows. The left is for the, Rayleigh stable, circular Keplerian rotational profile where $\Omega \propto R^{-3/2}$. Here the Reynolds stress becomes increasingly anisotropic as $\tau_c$ increases, although the cross term $R_{R\phi}$ is maximum near $\tau_c = 1$. The right is for a marginally Rayleigh unstable flow where $A = 1.1\Omega$. Here we see that the Reynolds stress diverges as $\tau_c \to \sqrt{5/2}$, where the moment expansion used to derive the turbulent stress model breaks down.</figcaption></figure>

$$R_{RR} = \alpha c_s^2 \rho_g \frac{1 + 2\tau_c^2 \left(2 - \frac{A}{\Omega}\right)}{1 + 4\tau_c^2 \left(1 - \frac{A}{\Omega}\right)}, \quad (A 31)$$

$$R_{R\phi} = \alpha c_s^2 R \left(\frac{A}{\Omega}\right) \frac{\tau_c}{1 + 4\tau_c^2 \left(1 - \frac{A}{\Omega}\right)} \rho_g, \quad (A 32)$$

$$R_{\phi\phi} = \alpha c_s^2 \rho_g R^2 \frac{1 + 2\tau_c^2 \left(1 - \frac{A}{\Omega}\right) \left(2 - \frac{A}{\Omega}\right)}{1 + 4\tau_c^2 \left(1 - \frac{A}{\Omega}\right)}. \quad (A 33)$$

This solution breaks down in the presence of strong shear, when $\frac{A}{\Omega} \ge 1 + \frac{1}{4\tau_c^2}$, as either $R_{RR}$ or $R_{\phi\phi}$ will be negative. This means there is no stable equilibrium for the Reynolds stress. As discussed in the previous section, this will also result in a breakdown of the moment expansion used to derive the equation governing the evolution of $R_{ij}$ as higher order moments can grow to become important. In Figure 7 we show the Reynolds stress components for the Rayleigh stable Keplerian shear flow and a Rayleigh unstable flow with $A = 1.1\Omega$.

### A.4. Accretion flow solutions - gas

In this section we derive the background gas solutions used in the numerical modelling of Section 7.2. Consider a Keplerian shear flow where the, circular, Keplerian motion is take to be order 1 and $u^R$, $R$, $p$ and partial time derivative $\partial_t$ are $O(\epsilon)$, for some small parameter $\epsilon$. We will neglect vertical gravity throughout, such that the Keplerian potential $\Phi = -1/R$ is a function of the Cylindrical radius, $R$, only (where we have adopted units in which $GM = 1$). We consider velocity, in cylindrical coordinates, of the form
$$u^R = \epsilon \tilde{u}^R, \qquad (A 35)$$

$$u^\phi = \Omega_K + \epsilon \tilde{u}^\phi, \qquad (A 36)$$

$$u^z = 0. \qquad (A 37)$$

The gas equations for an axisymmetric, vertically invariant flow, neglecting vertical gravity are

$$\dot{\rho} + R^{-1} \partial_R (R \rho_g \tilde{u}^R) = 0, \qquad (A 38)$$

$$2\rho R \Omega_K \tilde{u}^\phi = \partial_R p + (\nabla \cdot \mathbf{R})_R, \qquad (A 39)$$

$$\frac{1}{2} \rho \tilde{u}^R R \Omega_K = -(\nabla \cdot \mathbf{R})_\phi. \qquad (A 40)$$

The vertical momentum equation is trivially solved by $u^z = 0$. The only time derivative remaining is that in the continuity equation, responsible for the slow accretion of the gas. At this order the Reynolds stress is given by Equation A 31-A 33 and $R_{zz} = \alpha c_s^2 \rho_g$. For the Keplerian shear flow we take $A/\Omega = 3/4$ and $\alpha, c_s$ and $\tau_c$ to be constants. Substituting these into the dust stress gradients we have

$$(\nabla \cdot \mathbf{R})_R = \frac{1}{2} \alpha \frac{2 + 5\tau_c^2}{1 + \tau_c^2} \partial_R p + \frac{15}{8} \frac{p}{R} \frac{\alpha \tau_c^2}{1 + \tau_c^2}, \qquad (A 41)$$

$$(\nabla \cdot \mathbf{R})_\phi = \frac{3}{4} \frac{\alpha \tau_c}{1 + \tau_c^2} R^{-1} \partial_R [p R^2], \qquad (A 42)$$

noting that at this order $R_{Rz} = R_{\phi z} = 0$. This results in the following corrections to the gas velocity,

$$\tilde{u}^\phi = \frac{1}{2\rho R \Omega_K} \left[ \left( 1 + \frac{1}{2}\alpha \frac{2 + 5\tau_c^2}{1 + \tau_c^2} \right) \partial_R p + \frac{15}{8} \frac{p}{R} \frac{\alpha \tau_c^2}{1 + \tau_c^2} \right], \qquad (A 43)$$

$$\tilde{u}^R = -\frac{3}{2\rho R^2 \Omega_K} \frac{\alpha \tau_c}{1 + \tau_c^2} \partial_R [p R^2]. \qquad (A 44)$$

We consider two scenarios. One where we neglect the gas continuity equation and consider a fixed dust density. This approximation can be reasonable if the dust drift timescale or the characteristic lengthscale of the dust fluid is short, however we adopt it here mostly for illustrative purposes. The second scenario is to solve for a steady accretion flow with $\dot{\rho} = 0$ this leads to a gas pressure profile of

$$p = -\frac{4}{3} \mathcal{F} \frac{1 + \tau_c^2}{\alpha \tau_c} \Omega_K + C R^{-2}, \qquad (A 45)$$

where $\mathcal{F}$ and $C$ are constants. For our dust-fluid simulations in Section 7.2 we take $\rho \propto R^{-3/2}$, which is compatible with the steady state pressure profile above.
# Appendix B. Properties of the Dust Fluid Model

## B.1. Realisability

A necessary property of the constitutive relation is that the stress tensor be realisable from a second velocity moment of some distribution function. A similar property must hold for the dust Reynolds stress, the proof of which proceeds the same as the proof for the dust rheological stress tensor - to avoid repeating our selves we shall only cover the latter. As $\Pi_{\alpha\beta} = \int p(V_{\alpha} - U_{\alpha})(V_{\beta} - U_{\beta})d^6V$, $\Pi_{\alpha\beta}$ must be positive semi-definite. Thus for all positive semidefinite initial conditions $\Pi_{\alpha\beta}(0)$ our constitutive model most conserve the positive semi-definite character of $\Pi_{\alpha\beta}$. This is similar to the requirements for constitutive models of the MRI (Ogilvie 2003; Lynch & Ogilvie 2021).

Following (Lynch & Ogilvie 2021) we introduce the quadratic form $Q = \Pi_{\alpha\beta}Y^{\alpha}Y^{\beta}$, if the stress tensor is positive semi-definite then $Q \ge 0$ for all vectors $Y^{\alpha}$ at all points in the fluid. We will show by contradiction that an initially positive semi-definite $\Pi_{\alpha\beta}$ cannot evolve into one that is not positive semi-definite. Suppose, to the contrary, that some point in the flow $Q < 0$ for some vector $X^{\alpha}$ at some time after the initial state. Let us consider a smooth, evolving vector field $Y^{\alpha}$ that matches the vector $X^{\alpha}$ at the given point and time. The corresponding quadratic form $Q$ is then a scalar field that evolves according to

$$ \begin{aligned} \mathcal{D}Q &= Y^{\alpha}Y^{\beta}\mathcal{D}_2\Pi_{\alpha\beta} + \Pi_{\alpha\beta}\mathcal{D}_2(Y^{\alpha}Y^{\beta}) \\ &= -2Y^{\alpha}Y^{\beta}(\Pi^{\gamma}_{\alpha}\bar{A}_{\beta\gamma} - \rho_d D_{\alpha\beta}) + \Pi_{\alpha\beta}\mathcal{D}_2Y^{\alpha}Y^{\beta}, \end{aligned} \quad (B 1) $$

where, when operating on $Y^{\alpha}$, the differential operator $\mathcal{D}_2$ is

$$ \mathcal{D}_2Y^{\alpha} = DY^{\alpha} - Y^{\gamma}\nabla^{\alpha}U_{\gamma} - \frac{1}{2}Y^{\alpha}\nabla_{\gamma}U^{\gamma}. \quad (B 2) $$

By assumption, $Q$ is initially positive and evolves continuously to a negative value at the given later time. Therefore $Q$ must pass through zero at some intermediate time, which we denote by $t = 0$ without loss of generality. We can also assume, without loss of generality, that the vector field evolves according to $\mathcal{D}_2Y^{\alpha} = 0$, which means that it is advected by the flow. The equation for $Q$ then becomes

$$ DQ = -2Y^{\alpha}Y^{\beta}(\Pi^{\gamma}_{\alpha}C_{\beta\gamma} - \rho_d D_{\alpha\beta}) , \quad (B 3) $$

At $t = 0$, $Q = 0$ and, as $\Pi_{\alpha\beta}$ is positive semi-definite, one can show that $Y^{\alpha}\Pi_{\gamma\alpha} = 0$ so that the time derivative of $Q$ is given by

$$ DQ|_{t=0} = 2Y^{\alpha}Y^{\beta}\rho_d D_{\alpha\beta} \ge 0, \quad (B 4) $$

provided that $D_{\alpha\beta}$ is also positive semidefinite (which is guaranteed as $D_{\alpha\beta} = \frac{1}{2}g^{\mu\nu}\sigma_{\alpha\mu}\sigma_{\beta\nu}$). This contradicts the assumption that $Q$ passes through zero from positive to negative at $t = 0$. We thus conclude that $\Pi_{\alpha\beta}$ remains positive semi-definite provided it is initially.

## B.2. Viscoelasticity

In this appendix we shall explore the viscoelastic behaviour of the dust rheological stress. It's easiest to see the viscoelastic behaviour of the model when $t_s \sim t_c$. Introducing a characteristic relaxation time of the dust fluid $t_r \sim t_s \sim t_c$, and a characteristic fluid timescale $t_f$. We introduce the Deborah number $De = t_r/t_f$, with $C_{\alpha\beta}$ and $D_{\alpha\beta}$ being $O(De^{-1})$ while $\rho_d$ and $\bar{D}_2$ are $O(1)$. We can rewrite the constitutive relation as
$$ \bar{D}_2 T_{\alpha\beta} = -\frac{2}{\mathrm{De}} \left( T^\gamma_{(\alpha} C_{\beta)\gamma} - \rho_d D_{\alpha\beta} \right), \quad (B 5) $$

where we now treat De as a bookkeeping parameter to keep track of terms in an expansion in Deborah number. The constitutive model now has a similar form to classic viscoelastic models (e.g. the Oldroyd-B model, Oldroyd (1950)). The elastic limit can be recovered by taking $De \to \infty$ leading to,

$$ \bar{D}_2 T_{\alpha\beta} = \bar{D} T_{\alpha\beta} - 2 T^\gamma_{(\alpha} \varepsilon_{\beta)\gamma\sigma} \omega^\sigma = 0. \quad (B 6) $$

This corresponds to an elastic flow, with a source term from the flow vorticity. It is equivalent to the evolution of the Reynolds stress in the absence of source terms (e.g. see Gavrilyuk & Gouin 2012).

If we instead take the short Deborah number limit we can develop a series solution to Equation B 5 (similar to Lynch & Ogilvie 2021), which takes the form

$$ T_{\alpha\beta} = \sum_{n=0}^{\infty} \mathrm{De}^n T_{\alpha\beta}^{(n)}, \quad (B 7) $$

where

$$ T^{(0)\gamma}_{(\alpha} C_{\beta)\gamma} = \rho_d D_{\alpha\beta} \quad (B 8) $$

and

$$ T^{(n)\gamma}_{(\alpha} C_{\beta)\gamma} = -\frac{1}{2} \bar{D}_2 T^{(n-1)}_{\alpha\beta}, \quad n > 0. \quad (B 9) $$

Both these equations take the form $T^{(n)\gamma}_{(\alpha} C_{\beta)\gamma} = Q^{(n)}_{\alpha\beta}$, which can be inverted to obtain

$$ T^{(n)}_{\alpha_d \beta_d} = t_s Q^{(n)}_{\alpha_d \beta_d} + \frac{t_s t_c}{t_s + t_c} \left( Q^{(n)}_{(\alpha_d \beta_d^*)} + Q^{(n)}_{(\alpha_d^* \beta_d)} + \frac{t_c}{t_s} Q^{(n)}_{\alpha_d^* \beta_d^*} \right) \quad (B 10) $$

$$ T^{(n)}_{\alpha_d \beta_g} = \frac{t_s t_c}{t_s + t_c} \left( 2Q^{(n)}_{(\alpha_d \beta_g)} + \frac{t_c}{t_s} Q^{(n)}_{\alpha_d^* \beta_g} \right), \quad (B 11) $$

$$ T^{(n)}_{\alpha_g \beta_g} = t_c Q^{(n)}_{\alpha_g \beta_g}. \quad (B 12) $$

Substituting $Q^{(0)}_{\alpha\beta} = D_{\alpha\beta}$, and making use of the properties of $D_{\alpha\beta}$ we obtain

$$ T^{(0)}_{\alpha_d \beta_d} = \frac{t_c^2}{t_s + t_c} D_{\alpha_d^* \beta_d^*} = \frac{t_c}{t_s + t_c} \alpha c_s^2 \rho_d g_{\alpha_d \beta_d} \quad (B 13) $$

$$ T^{(0)}_{\alpha_d \beta_g} = \frac{t_c^2}{t_s + t_c} D_{\alpha_d^* \beta_g} = \frac{t_c}{t_s + t_c} \alpha c_s^2 \rho_d g_{\alpha_d^* \beta_g}, \quad (B 14) $$

$$ T^{(0)}_{\alpha_g \beta_g} = t_c D_{\alpha_t \beta_t} = \alpha c_s^2 \rho_d g_{\alpha_g \beta_g}. \quad (B 15) $$

To calculate $T^{(1)}_{\alpha\beta}$ we first need to calculate $\bar{D}_2 T^{(0)}_{\alpha\beta}$. For simplicity we shall assume that $\alpha$, $c_s$, $t_s$, $t_c$ are constant, and that the metric tensor is time independent then we obtain
$$ (\overline{\mathcal{D}_2} \mathbf{T}^{(0)})_{\alpha_d \beta_d} = \frac{2t_c}{t_s + t_c} \alpha c_s^2 \rho_d \bar{\nabla}_{(\alpha_d} U_{\beta_d)}, \quad (B 16) $$

$$ (\overline{\mathcal{D}_2} \mathbf{T}^{(0)})_{\alpha_d \beta_g} = \frac{t_c}{t_s + t_c} \alpha c_s^2 \rho_d (\bar{\nabla}_{\alpha_d} U_{\beta_g} + \bar{\nabla}_{\beta_g} U_{\alpha_d}) \quad (B 17) $$

$$ (\overline{\mathcal{D}_2} \mathbf{T}^{(0)})_{\alpha_g \beta_g} = \frac{t_c}{t_s + t_c} \alpha c_s^2 \rho_d (\bar{\nabla}_{\alpha_g} U_{\beta_g} + \bar{\nabla}_{\beta_g} U_{\alpha_g}). \quad (B 18) $$

Substituting this into Equations B 10-B 12 we obtain

$$ T_{\alpha_d \beta_d}^{(1)} = -\frac{1}{2} \frac{t_c^2}{t_s + t_c} \alpha c_s^2 \rho_d \left( 2 \frac{t_s}{t_c} \frac{t_s + 2t_c}{t_s + t_c} \bar{\nabla}_{(\alpha_d} U_{\beta_d)} + \bar{\nabla}_{\alpha_d} U_{\beta_d^*} + \bar{\nabla}_{\beta_d} U_{\alpha_d^*} \right) \quad (B 19) $$

$$ T_{\alpha_d \beta_g}^{(1)} = -\frac{1}{2} \frac{t_s t_c^2}{(t_s + t_c)^2} \alpha c_s^2 \rho_d \left[ \left(2 + \frac{t_c}{t_s}\right) \bar{\nabla}_{\alpha_d} U_{\beta_g} + 2 \bar{\nabla}_{\beta_g^*} U_{\alpha_d} + \frac{t_c}{t_s} \bar{\nabla}_{\beta_g^*} U_{\alpha_d^*} \right], \quad (B 20) $$

$$ T_{\alpha_g \beta_g}^{(1)} = -\frac{1}{2} \frac{t_c^2}{t_s + t_c} \alpha c_s^2 \rho_d (\bar{\nabla}_{\alpha_g} U_{\beta_g} + \bar{\nabla}_{\beta_g} U_{\alpha_g}). \quad (B 21) $$

This results in a rheological stress tensor of the form,

$$ T_{\alpha\beta} = p_d \left( 1 + \frac{t_s}{t_c} \Theta_{\alpha\beta}^g \right) g_{\alpha\beta} + \frac{1}{2} p_x (g_{\alpha\beta^*} + g_{\alpha^*\beta}) - 2 \mu_{\alpha\beta}^{\mu\nu} \bar{\nabla}_{\mu} U_{\nu} + O(\mathrm{De}^2) \quad (B 22) $$

where we have introduced the dust pressure, $p_d$, and cross pressure, $p_x$, with $p_d = p_x = \alpha c_s^2 \frac{t_c}{t_s+t_c} \rho_d$; and for convenience we have defined $\Theta_{\alpha\beta}^g$, with $\Theta_{\alpha_d\beta_d}^g = \Theta_{\alpha_d\beta_g}^g = \Theta_{\alpha_g\beta_d}^g = 0$ and $\Theta_{\alpha_g\beta_g}^g = 1$. We have also introduced an anisotropic viscosity tensor, $\mu_{\alpha\beta}^{\mu\nu}$, given by,

$$ \mu_{\alpha\beta}^{\mu\nu} = \mu_d \delta_{(\alpha}^{\mu} \delta_{\beta)}^{\nu} + \mu_x \left( \delta_{(\alpha}^{\mu} \delta_{\beta^*)}^{\nu} + \delta_{(\alpha^*}^{\mu} \delta_{\beta)}^{\nu} \right) + \eta_{\alpha\beta}^{\mu\nu}, \quad (B 23) $$

with

$$ \mu_d = \frac{1}{2} \frac{t_s t_c (t_s + 2t_c)}{(t_s + t_c)^2} \alpha c_s^2 \rho_d, \quad (B 24) $$

$$ \mu_x = \frac{1}{2} \frac{t_c^2}{t_s + t_c} \alpha c_s^2 \rho_d, \quad (B 25) $$

and

$$ \eta_{\alpha_d \beta_d}^{\mu \nu} = \eta_{\alpha_g \beta_g}^{\mu \nu} = 0, \quad (B 26) $$

$$ \eta_{\alpha_d \beta_g}^{\mu \nu} = \frac{1}{2} \mu_x \left[ \frac{t_c^2 - t_s^2}{t_c} \delta_{\alpha_d}^{\mu} \delta_{\beta_g}^{\nu} + (t_s - t_c) \delta_{\beta_g^*}^{\mu} \delta_{\alpha_d}^{\nu} - (t_s + t_c) \delta_{\alpha_d}^{\mu} \delta_{\beta_g^*}^{\nu} + \frac{t_c}{t_s} \delta_{\beta_g^*}^{\mu} \delta_{\alpha_g}^{\nu} \right]. \quad (B 27) $$

### B.3. Dust velocity correlations in the short stopping time limit.

In this section we consider the short stopping time behaviour of the rheological stress tensor. Naively one might expect the dust velocity correlations to match the gas velocity correlations due to the tight coupling between the gas and dust. This would mean the dust stress tensor would be given by $T_{\alpha_d \beta_d} = f_d R_{\alpha_d^* \beta_d^*}$ where $f_d$ is the dust to gas ratio. We shall show that this is only the case when the dust experiences the turbulence as a continuum where the dust
interacts with many turbulent eddies over the lengthscale on which the dust fluid varies. If, however an individual eddy transports a dust particle a significant distance in the fluid then the dust velocity correlations can deviate strongly from those of the gas.

We wish to compare the evolutionary equations for the gas Reynolds stress to that of the dust stress tensor. The gas Reynolds stress evolves according to

$$ \tilde{D}R_{ij} + 2R_{k(i}\nabla^k u^g_{j)} + R_{ij}\nabla_k u^k_g = -\frac{2}{t_c}(R_{ij} - \rho_g D_{ij}), \quad (B 28) $$

where $\tilde{D}$ is the Lagrangian time derivative with respect to the mean gas flow, $u^g$. The dust stress tensor evolves according to

$$ \overline{D}T_{\alpha\beta} + 2T_{\gamma(\alpha}\overline{\nabla}^{\gamma}U_{\beta)} + T_{\alpha\beta}\overline{\nabla}_{\gamma}U^{\gamma} = -2(T^{\gamma}_{(\alpha}C_{\beta)\gamma} - \rho_d D_{\alpha\beta}), \quad (B 29) $$

where $\overline{D}$ is the Lagrangian time derivative with respect to the mean dust flow.

Because of the factor of the dust to gas ratio between the dust rheological and the gas Reynolds stress, it is more convenient to work with the respective velocity correlation tensors. The dust velocity correlation tensor $W_{\alpha\beta} = T_{\alpha\beta}/\rho_d$ which evolves according to

$$ \overline{D}W_{\alpha\beta} + 2W_{\gamma(\alpha}\overline{\nabla}^{\gamma}U_{\beta)} = -2(W^{\gamma}_{(\alpha}C_{\beta)\gamma} - D_{\alpha\beta}), \quad (B 30) $$

We can also write the evolutionary equation for the gas velocity correlation tensor, $R_{ij}/\rho_g$, in the form

$$ \tilde{D}_2(R_{ij}/\rho_g) = -\frac{2}{t_c}((R_{ij}/\rho_g) - D_{ij}), \quad (B 31) $$

where we have introduced the differential operator $\tilde{D}_2$ which, when acting on $R_{ij}/\rho_g$, is given by

$$ \tilde{D}_2(R_{ij}/\rho_g) = \tilde{D}(R_{ij}/\rho_g) + (2/\rho_g)R_{k(i}\nabla^k u^g_{j)} \quad (B 32) $$

We now consider small dust grains ($t_s \to 0$) embedded in a gas flow with mean velocity $u^g_i$. The ‘dust’ components of the dust fluid momentum equation, in the limit $t_s \to 0$, simplify to

$$ U^{\alpha d} = U^{\alpha_d}, \quad (B 33) $$

while the ‘dummy gas’ components are

$$ \rho_d \overline{D} U_{\alpha_g} = \rho_d F_{\alpha_g} - \overline{\nabla}^{\beta_d} T_{\alpha_g \beta_d} - \frac{1}{t_c} \rho_d (U_{\alpha_g} - U^g_{\alpha_g}). \quad (B 34) $$

While the mean dust velocity is tightly coupled to the mean gas velocity ($U^{\alpha_d} = U^{\alpha_g}$) the mean gas velocity as experience by the dust is not generally the same as the mean gas velocity experience by the gas, $u^g_i$. This is because the dust is effectively a subsample of the gas velocity field and can experience a mean gas velocity relative to the gas frame due to correlation in the gas turbulence. This distinction is vital for allowing zero stopping time particles to diffuse in gas turbulence.

We can write the 6D dust velocity as follows

$$ U^{\alpha_d} = u^{\alpha_d} + \Delta U^{\alpha_d}, \quad U^{\alpha_g} = u^{\alpha_g} + \Delta U^{\alpha_g}, \quad (B 35) $$
where $\Delta U^i$ is the relative velocity with respect to the mean gas flow experienced by the gas, which need not be small. With this velocity for the dust flow the Lagrangian time derivative, $\tilde{D}$, can be related to $\tilde{D}$ by

$$ \tilde{D} = \tilde{D} + \Delta U^i \bar{\nabla}_y. \qquad (B 36) $$

Substituting this velocity into Equation B 30 and separating the dust and dummy gas components of the dust constitutive relation we get

$$ \tilde{D} W_{\alpha_d \beta_d} + 2 W_{\gamma(\alpha_d} \bar{\nabla}^\gamma u_{\beta_d)}^g + \Delta U^\gamma \bar{\nabla}_\gamma W_{\alpha_d \beta_d} + 2 W_{\gamma(\alpha_d} \bar{\nabla}^\gamma \Delta U_{\beta_d)} = - \frac{1}{t_s} (2 W_{\alpha_d \beta_d} - W_{\alpha_d}^* \beta_d - W_{\alpha_d \beta_d}^*), \qquad (B 37) $$

$$ \tilde{D} W_{\alpha_d \beta_g} + W_{\gamma \alpha_d} \bar{\nabla}^\gamma (u_{\beta_g}^g + \Delta U \beta_g^*) + W_{\gamma \beta_g} \bar{\nabla}^\gamma (u_{\alpha_d}^g + \Delta U \alpha_d) + \Delta U^\gamma \bar{\nabla}_\gamma W_{\alpha_d \beta_g} = - \left(\frac{1}{t_s} + \frac{1}{t_c}\right) W_{\alpha_d \beta_g} + \frac{1}{t_s} W_{\alpha_d}^* \beta_g, \qquad (B 38) $$

$$ \tilde{D} W_{\alpha_g \beta_g} + 2 W_{\gamma(\alpha_g} \bar{\nabla}^\gamma u_{\beta_g)}^g + \Delta U^\gamma \bar{\nabla}_\gamma W_{\alpha_g \beta_g} + 2 W_{\gamma(\alpha_g} \bar{\nabla}^\gamma \Delta U_{\beta_g)}^* = - \frac{2}{t_c} (W_{\alpha_g \beta_g} - \alpha c_s^2 g_{\alpha_g \beta_g}). \qquad (B 39) $$

Taking the short stopping time limit of Equations B 37 and B 38 leads to

$$ W_{\alpha_d \beta_d} = \frac{1}{2} (W_{\alpha_d}^* \beta_d + W_{\alpha_d \beta_d}^*), \qquad (B 40) $$

and

$$ W_{\alpha_d \beta_g} = W_{\alpha_d}^* \beta_g. \qquad (B 41) $$

We can use these relations to simplify the ‘dummy gas’ components (Equation B 39), which can be rearranged to obtain

$$ t_c \tilde{D}_2 W_{\alpha_g \beta_g} + 2(W_{\alpha_g \beta_g} - \alpha c_s^2 g_{\alpha_g \beta_g}) = -t_c (\Delta U^\gamma \bar{\nabla}_\gamma W_{\alpha_g \beta_g} + 2W_{\gamma(\alpha_g} \bar{\nabla}^\gamma \Delta U_{\beta_g)}). \qquad (B 42) $$

If $\Delta U^*$ is the characteristic scale of the relative velocity $\Delta U^\alpha$ and $L$ is a characteristic lengthscale of variations in the fluid flow then we can introduce an Eddy-Knudsen number,

$$ \text{Kn}_e = \frac{\lambda}{L} = \frac{t_c \Delta U^*}{L}, \qquad (B 43) $$

where $\lambda = t_c \Delta U^*$ is the mean free path of a dust grain in the turbulent flow representing the lengthscale a dust grain is transported by a single eddy. When $\text{Kn}_e \ll 1$ the dust experiences the gas turbulence as a continuum, interacting with a large number of turbulent eddies on the lengthscale of the fluid. When $\text{Kn}_e \gtrsim 1$ the dynamics of a dust grain is dominated by the last eddy that it interacted with - in a similar manor to the effects of individual particle collisions in weakly collisional gasses. Rescaling the right-hand side of Equation B 42 and making use of the Eddy-Knudsen number we arrive at
$$t_c \tilde{D}_2 W_{\alpha_g \beta_g} + 2(W_{\alpha_g \beta_g} - \alpha c_s^2 g_{\alpha_g \beta_g}) = -K n_e \left( \frac{\Delta U^\gamma}{\Delta U^*} L \bar{\nabla}_\gamma W_{\alpha_g \beta_g} + 2 W_\gamma (\alpha_g L \bar{\nabla}^\gamma \frac{\Delta U_{\beta_g}}{\Delta U^*}) \right). \quad (B44)$$

In the limit $K n_e \to 0$ this matches Equation B 31 for the turbulent gas velocity correlations. Thus we conclude that in limit $t_s$, $K n_e \to 0$ the dust velocity correlations are set by the gas velocity correlations. However when $K n_e \gtrsim 1$ the dust velocity correlations no longer match those of the gas as the dust velocity correlations are strongly affected by individual eddies.

## Appendix C. Higher moments of the Fokker-Planck equation terms

The individual terms in the Fokker-Planck equation satisfy the follow relations, which are important for deriving the evolutionary equations for the higher velocity moments,

$$\int (V_{\alpha_1} - U_{\alpha_1}) \cdots (V_{\alpha_k} - U_{\alpha_k}) \frac{\partial p}{\partial t} d^6 \mathbf{V} = \frac{\partial \Pi_{\alpha_1 \cdots \alpha_k}}{\partial t} + k \Pi_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial t} U_{\alpha_k}, \quad (C1)$$

$$\int (V_{\alpha_1} - U_{\alpha_1}) \cdots (V_{\alpha_k} - U_{\alpha_k}) \frac{\partial}{\partial X_{\sigma}} [V_{\sigma} p] d^6 \mathbf{V} \\ = \frac{\partial}{\partial X_{\sigma}} (\Pi_{\alpha_1 \cdots \alpha_k} \sigma + U_{\sigma} \Pi_{\alpha_1 \cdots \alpha_k}) + k \Pi_{(\alpha_1 \cdots \alpha_{k-1})}^{\sigma} \frac{\partial}{\partial X_{\sigma}} U_{\alpha_k} + k U^{\sigma} \Pi_{(\alpha_1 \cdots \alpha_{k-1})} \frac{\partial}{\partial X_{\sigma}} U_{\alpha_k}, \quad (C2)$$

$$\int (V_{\alpha_1} - U_{\alpha_1}) \cdots (V_{\alpha_k} - U_{\alpha_k}) \frac{\partial}{\partial V_{\sigma}} [p \nabla_{\sigma} \Phi + p C_{\sigma \gamma} (V^{\gamma} - U_{\sigma}^{\gamma})] d^6 \mathbf{V} \\ = -k \Pi_{(\alpha_1 \cdots \alpha_{k-1})} [\nabla_{\alpha_k} \Phi + C_{\alpha_n}^{\gamma} (U_{\gamma} - U_{\alpha_k}^{\gamma})] - k \Pi_{(\alpha_1 \cdots \alpha_{k-1})} C_{\alpha_n}^{\gamma} \quad (C3)$$

$$\int (V_{\alpha_1} - U_{\alpha_1}) \cdots (V_{\alpha_k} - U_{\alpha_k}) D_{\mu\nu} \frac{\partial^2 p}{\partial V_{\mu} \partial V_{\nu}} d^6 \mathbf{V} = k(k-1) \Pi_{(\alpha_1 \cdots \alpha_{k-2})} D_{\alpha_{k-1} \alpha_{k-2}}. \quad (C4)$$

## REFERENCES

* ARAKI, S. & TREMAINE, S. 1986 The dynamics of dense particle disks. *Icarus* **65** (1), 83–109.

* ARMITAGE, PHILIP J. 2020 *Astrophysics of Planet Formation*, 2nd edn. Cambridge University Press.

* BAI, XUE-NING & STONE, JAMES M. 2010 Particle-gas Dynamics with Athena: Method and Convergence. *ApJS* **190** (2), 297–310, arXiv: 1005.4980.

* BAINES, M. J., WILLIAMS, I. P. & ASEBIOMO, A. S. 1965 Resistance to the motion of a small sphere moving through a gas. *MNRAS* **130**, 63.

* BALBUS, S. A. & HAWLEY, J. F. 1991 A powerful local shear instability in weakly magnetized disks. I - Linear analysis. II - Nonlinear evolution. *ApJ* **376**, 214–233.

* BALSARA, DINSHAW S., TILLEY, DAVID A., RETTIG, TERRENCE & BRITTAIN, SEAN D. 2009 Dust settling in magnetorotationally driven turbulent discs - I. Numerical methods and evidence for a vigorous streaming instability. *MNRAS* **397** (1), 24–43, arXiv: 0810.0246.

* BARKER, A. J. & OGILVIE, G. I. 2014 Hydrodynamic instability in eccentric astrophysical discs. *MNRAS* **445**, 2637–2654, arXiv: 1409.6488.

* BARRIÈRE-FOUCHET, L., GONZALEZ, J. F., MURRAY, J. R., HUMBLE, R. J. & MADDISON, S. T. 2005 Dust distribution in protoplanetary disks. Vertical settling and radial migration. *A&A* **443** (1), 185–194, arXiv: astro-ph/0508452.
BENÍTEZ-LLAMBAY, PABLO & MASSET, FRÉDÉRIC S. 2016 FARGO3D: A New GPU-oriented MHD Code. *ApJS* **223** (1), 11, arXiv: 1602.02359.

BI, JIAQING, LIN, MIN-KAI & DONG, RUOBING 2021 Puffed-up Edges of Planet-opened Gaps in Protoplanetary Disks. I. Hydrodynamic Simulations. *ApJ* **912** (2), 107, arXiv: 2103.09254.

BINKERT, FABIAN 2023 Beyond diffusion: a generalized mean-field theory of turbulent dust transport in protoplanetary discs. *MNRAS* **525** (3), 4299–4320, arXiv: 2306.06103.

BOYLEV, ALEKSANDR VASIL'EVICH 1982 The chapman-enskog and grad methods for solving the boltzmann equation. In *Akademiia Nauk SSSR Doklady*, **v**, vol. 262, pp. 71–75.

BOYLEV, A. V. 2018 Boltzmann equation and hydrodynamics beyond navier–stokes. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **376** (2118), 20170227, arXiv: https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2017.0227.

BOOTH, RICHARD A. & CLARKE, CATHIE J. 2021 Modelling the delivery of dust from discs to ionized winds. *MNRAS* **502** (2), 1569–1578, arXiv: 2101.04121.

BORDERIES, N., GOLDREICH, P. & TREMAINE, S. 1985 A granular flow model for dense planetary rings. *Icarus* **63** (3), 406–420.

CAPECELATRO, JESSE, DESJARDINS, OLIVIER & FOX, RODNEY O 2016a Strongly coupled fluid-particle flows in vertical channels. i. reynolds-averaged two-phase turbulence statistics. *Physics of Fluids* **28** (3).

CAPECELATRO, JESSE, DESJARDINS, OLIVIER & FOX, RODNEY O 2016b Strongly coupled fluid-particle flows in vertical channels. ii. turbulence modeling. *Physics of Fluids* **28** (3).

CARBALLIDO, AUGUSTO, FROMANG, SÉBASTIEN & PAPALOIZOU, JOHN 2006 Mid-plane sedimentation of large solid bodies in turbulent protoplanetary discs. *MNRAS* **373** (4), 1633–1640, arXiv: astro-ph/0610075.

CHAPMAN, SYDNEY & COWLING, THOMAS GEORGE 1990 *The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases*. Cambridge university press.

COMMERÇON, B., LEBREUILLY, U., PRICE, D. J., LOVASCIO, F., LAIBE, G. & HENNEBELLE, P. 2023 Dynamics of dust grains in turbulent molecular clouds. Conditions for decoupling and limits of different numerical implementations. *A&A* **671**, A128, arXiv: 2301.04946.

CSANADY, GT 1963 Turbulent diffusion of heavy particles in the atmosphere. *Journal of Atmospheric Sciences* **20** (3), 201–208.

DIAMOND, P. H., LIANG, Y.-M., CARRERAS, B. A. & TERRY, P. W. 1994 Self-regulating shear flow turbulence: A paradigm for the l to h transition. *Phys. Rev. Lett.* **72**, 2565–2568.

DUBRULLE, B., MORFILL, G. & STERZIK, M. 1995 The dust subdisk in the protoplanetary nebula. *Icarus* **114** (2), 237–246.

ELEUTERIO, F TORO 1999 Riemann solvers and numerical methods for fluid dynamics.

EPSTEIN, PAUL S. 1924 On the Resistance Experienced by Spheres in their Motion through Gases. *Physical Review* **23** (6), 710–733.

FLOCK, MARIO, NELSON, RICHARD P., TURNER, NEAL J., BERTRANG, GESA H. M., CARRASCO-GONZÁLEZ, CARLOS, HENNING, THOMAS, LYRA, WLADIMIR & TEAGUE, RICHARD 2017 Radiation Hydrodynamical Turbulence in Protoplanetary Disks: Numerical Models and Observational Constraints. *ApJ* **850** (2), 131, arXiv: 1710.06007.

FOX, RODNEY O. 2014 On multiphase turbulence models for collisional fluid-particle flows. *Journal of Fluid Mechanics* **742**, 368–424.

FROMANG, S. & PAPALOIZOU, J. 2006 Dust settling in local simulations of turbulent protoplanetary disks. *A&A* **452** (3), 751–762, arXiv: astro-ph/0603153.

GAVRILYUK, S. & GOUIN, H. 2012 Geometric evolution of the reynolds stress tensor. *International Journal of Engineering Science* **59**, 65–73, the Special Issue in honor of VICTOR L. BERDICHEVSKY.

GOLDREICH, P. & TREMAINE, S. 1978 The excitation and evolution of density waves. *ApJ* **222**, 850–858.

GRAD, HAROLD 1948 *Approximation to the Boltzmann equation by moments*. New York University.

GRAD, HAROLD 1949 On the kinetic theory of rarefied gases. *Communications on pure and applied mathematics* **2** (4), 331–407.

HARTEN, AMIRAM, LAX, PETER D & LEER, BRAM VAN 1983 On upstream differencing and godunov-type schemes for hyperbolic conservation laws. *SIAM review* **25** (1), 35–61.

HAWLEY, JOHN F & BALBUS, STEVEN A 1991 A powerful local shear instability in weakly magnetized disks. ii. nonlinear evolution. *The Astrophysical Journal* **376**, 223.

HAWLEY, J. F., GAMMIE, C. F. & BALBUS, S. A. 1995 Local Three-dimensional Magnetohydrodynamic Simulations of Accretion Disks. *ApJ* **440**, 742.
HOBSON, MICHAEL PAUL, EFSTATHIOU, GEORGE P & LASENBY, ANTHONY N 2006 *General relativity: an introduction for physicists*. Cambridge University Press.

INNOCENTI, ALESSIO, FOX, R. & CHIBBARO, SERGIO 2019 A lagrangian probability-density-function model for collisional turbulent fluid-particle flows. *Journal of Fluid Mechanics* **862**, 449.

INNOCENTI, ALESSIO, FOX, RODNEY O. & CHIBBARO, SERGIO 2021 A Lagrangian probability-density-function model for turbulent particle-laden channel flow in the dense regime. *Physics of Fluids* **33** (5), 053308, arXiv: https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/5.0045690/16085748/053308_1_online.pdf.

LAIBE, GUILLAUME, BRÉHIER, CHARLES-EDOUARD & LOMBART, MAXIME 2020 On the settling of small grains in dusty discs: analysis and formulae. *MNRAS* **494** (4), 5134–5147, arXiv: 2004.03689.

LAIBE, GUILLAUME & PRICE, DANIEL J. 2012a Dusty gas with smoothed particle hydrodynamics - I. Algorithm and test suite. *MNRAS* **420** (3), 2345–2364, arXiv: 1111.3090.

LAIBE, GUILLAUME & PRICE, DANIEL J. 2012b Dusty gas with smoothed particle hydrodynamics - II. Implicit timestepping and astrophysical drag regimes. *MNRAS* **420** (3), 2365–2376, arXiv: 1111.3089.

LAIBE, GUILLAUME & PRICE, DANIEL J. 2014 Dusty gas with one fluid. *MNRAS* **440** (3), 2136–2146, arXiv: 1402.5248.

LARUE, RÉMY, LATTER, HENRIK & REIN, HANNO 2023 Thermal hysteresis and front propagation in dense planetary rings. *MNRAS* **520** (1), 1128–1145, arXiv: 2301.03289.

LATTER, HENRIK N. & OGILVIE, GORDON I. 2008 Dense planetary rings and the viscous overstability. *Icarus* **195** (2), 725–751, arXiv: 0803.4123.

LATTER, HENRIK N. & PAPALOIZOU, JOHN 2017 Local models of astrophysical discs. *MNRAS* **472** (2), 1432–1446, arXiv: 1708.09285.

LESUR, G., ERCOLANO, B., FLOCK, M., LIN, M. K., YANG, C. C., BARRANCO, J. A., BENITEZ-LLAMBAY, P., GOODMAN, J., JOHANSEN, A., KLAHR, H., LAIBE, G., LYRA, W., MARCUS, P., NELSON, R. P., SQUIRE, J., SIMON, J. B., TURNER, N., UMURHAN, O. M. & YOUDIN, A. N. 2022 Hydro-, Magnetohydro-, and Dust-Gas Dynamics of Protoplanetary Disks. *arXiv e-prints* p. arXiv:2203.09821, arXiv: 2203.09821.

LEVERMORE, C DAVID 1996 Moment closure hierarchies for kinetic theories. *Journal of statistical Physics* **83**, 1021–1065.

LIN, MIN-KAI 2019 Dust settling against hydrodynamic turbulence in protoplanetary discs. *MNRAS* **485** (4), 5221–5234, arXiv: 1903.03620.

LIN, MIN-KAI & YOUDIN, ANDREW N. 2015 Cooling Requirements for the Vertical Shear Instability in Protoplanetary Disks. *ApJ* **811** (1), 17, arXiv: 1505.02163.

LIN, MIN-KAI & YOUDIN, ANDREW N. 2017 A Thermodynamic View of Dusty Protoplanetary Disks. *ApJ* **849** (2), 129, arXiv: 1708.02945.

LYNCH, ELLIOT M. & OGILVIE, GORDON I. 2021 Importance of magnetic fields in highly eccentric discs with applications to tidal disruption events. *MNRAS* **501** (4), 5500–5516, arXiv: 2101.01221.

MASSET, F. S. 2000 FARGO: A Fast Eulerian Transport Algorithm for Differentially Rotating Disks. In *Disks, Planetesimals, and Planets* (ed. G. Garzón, C. Eiroa, D. de Winter & T. J. Mahoney), *Astronomical Society of the Pacific Conference Series*, vol. 219, p. 75.

MIGNONE, A., FLOCK, M. & VAIDYA, B. 2019 A Particle Module for the PLUTO Code. III. Dust. *ApJS* **244** (2), 38, arXiv: 1908.10793.

MINIER, JEAN-PIERRE 2001 Probabilistic approach to turbulent two-phase flows modelling and simulation: theoretical and numerical issues 7 (3-4), 295–310.

MINIER, JEAN-PIERRE 2015 On lagrangian stochastic methods for turbulent polydisperse two-phase reactive flows. *Progress in Energy and Combustion Science* **50**, 1–62.

MINIER, JEAN-PIERRE 2016 Statistical descriptions of polydisperse turbulent two-phase flows. *Physics Reports* **665** (2016), 1–122.

MINIER, JEAN-PIERRE, CHIBBARO, SERGIO & POPE, STEPHEN B 2014 Guidelines for the formulation of lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent flows. *Physics of Fluids* **26** (11), 113303.

MINIER, JEAN-PIERRE & HENRY, CHRISTOPHE 2023 The dynamics of discrete particles in turbulent flows: open issues and current challenges in statistical modeling. *arXiv e-prints* p. arXiv:2311.01921, arXiv: 2311.01921.

MINIER, JEAN-PIERRE & PEIRANO, ERIC 2001 The pdf approach to turbulent polydispersed two-phase flows. *Physics Reports* **352** (1), 1–214.

MINIER, JEAN-PIERRE, PEIRANO, ERIC & CHIBBARO, SERGIO 2004 Pdf model based on langevin equation
for polydispersed two-phase flows applied to a bluff-body gas-solid flow. *Physics of fluids* **16** (7), 2419–2431.

NELSON, RICHARD P., GRESSEL, OLIVER & UMURHAN, ORKAN M. 2013 Linear and non-linear evolution of the vertical shear instability in accretion discs. *MNRAS* **435** (3), 2610–2632, arXiv: 1209.2753.

OGILVIE, G. I. 2003 On the dynamics of magnetorotational turbulent stresses. MNRAS **340** (3), 969–982, arXiv: astro-ph/0212442.

OGILVIE, G. I. & BARKER, A. J. 2014 Local and global dynamics of eccentric astrophysical discs. MNRAS **445**, 2621–2636, arXiv: 1409.6487.

OGILVIE, G. I. & LATTER, H. N. 2013a Hydrodynamic instability in warped astrophysical discs. MNRAS **433**, 2420–2435, arXiv: 1303.0264.

OGILVIE, G. I. & LATTER, H. N. 2013b Local and global dynamics of warped astrophysical discs. MNRAS **433**, 2403–2419, arXiv: 1303.0263.

OLDROYD, J. G. 1950 On the Formulation of Rheological Equations of State. *Proceedings of the Royal Society of London Series A* **200** (1063), 523–541.

ORMEL, C. W. & CUZZI, J. N. 2007 Closed-form expressions for particle relative velocities induced by turbulence. A&A **466** (2), 413–420, arXiv: astro-ph/0702303.

ORMEL, CHRIS W. & LIU, BEIBEI 2018 Catching drifting pebbles. II. A stochastic equation of motion for pebbles. A&A **615**, A178, arXiv: 1803.06150.

PAPALOIZOU, J. C. B. 2005a Global numerical simulations of differentially rotating disks with free eccentricity. A&A **432** (3), 757–769.

PAPALOIZOU, J. C. B. 2005b The local instability of steady astrophysical flows with non circular streamlines with application to differentially rotating disks with free eccentricity. A&A **432** (3), 743–755, arXiv: astro-ph/0412587.

PEIRANO, E., CHIBBARO, S., POZORSKI, J. & MINIER, J.-P. 2006 Mean-field/pdf numerical approach for polydispersed turbulent two-phase flows. *Progress in Energy and Combustion Science* **32** (3), 315–371.

POPE, S.B. 1985 Pdf methods for turbulent reactive flows. *Progress in Energy and Combustion Science* **11** (2), 119–192.

POPE, SB 1987 Consistency conditions for random-walk models of turbulent dispersion. *The Physics of fluids* **30** (8), 2374–2379.

POPE, STEPHEN B. 2000 *Turbulent Flows*. Cambridge University Press.

POPE, STEPHEN B 2002 A stochastic lagrangian model for acceleration in turbulent flows. *Physics of Fluids* **14** (7), 2360–2375.

ROE, PHILIP L 1981 Approximate riemann solvers, parameter vectors, and difference schemes. *Journal of computational physics* **43** (2), 357–372.

SAWFORD, BL 1991 Reynolds number effects in lagrangian stochastic models of turbulent dispersion. *Physics of Fluids A: Fluid Dynamics* **3** (6), 1577–1586.

SVANBERG, ELEONORA, CUI, CAN & LATTER, HENRIK N. 2022 Wavelike nature of the vertical shear instability in global protoplanetary discs. MNRAS **514** (3), 4581–4587, arXiv: 2206.03840.

TESTI, L., BIRNSTIEL, T., RICCI, L., ANDREWS, S., BLUM, J., CARPENTER, J., DOMINIK, C., ISELLA, A., NATTA, A., WILLIAMS, J. P. & WILNER, D. J. 2014 Dust Evolution in Protoplanetary Disks. In *Protostars and Planets VI* (ed. Henrik Beuther, Ralf S. Klessen, Cornelis P. Dullemond & Thomas Henning), pp. 339–361, arXiv: 1402.1354.

THOMSEN, LEON 1986 Weak elastic anisotropy. *Geophysics* **51** (10), 1954–1966.

THOMSON, DJ 1987 Criteria for the selection of stochastic models of particle trajectories in turbulent flows. *Journal of fluid mechanics* **180**, 529–556.

TORO, ELEUTERIO F, SPRUCE, MICHAEL & SPEARES, WILLIAM 1994 Restoration of the contact surface in the hll-riemann solver. *Shock waves* **4**, 25–34.

VAN LEER, BRAM 2006 Upwind and high-resolution methods for compressible flow: From donor cell to residual-distribution schemes. In *16th aiaa computational fluid dynamics conference*, p. 3559.

WHIPPLE, F. L. 1972 On certain aerodynamic processes for asteroids and comets. In *From Plasma to Planet* (ed. Aina Elvius), p. 211.

WITHERS, C.S. 1985 The moments of the multivariate normal. *Bulletin of the Australian Mathematical Society* **32** (1), 103–107.

YANG, CHAO-CHIN & JOHANSEN, ANDERS 2016 Integration of Particle-gas Systems with Stiff Mutual Drag Interaction. *ApJS* **224** (2), 39, arXiv: 1603.08523.
YOUDIN, A. & JOHANSEN, A. 2007 Protoplanetary Disk Turbulence Driven by the Streaming Instability: Linear Evolution and Numerical Methods. *ApJ* **662** (1), 613–626, arXiv: astro-ph/0702625.

YOUDIN, ANDREW N. & LITHWICK, YORAM 2007 Particle stirring in turbulent gas disks: Including orbital oscillations. *Icarus* **192** (2), 588–604, arXiv: 0707.2975.

ZHU, ZHAOHUAN, STONE, JAMES M., RAFIKOV, ROMAN R. & BAI, XUE-NING 2014 Particle Concentration at Planet-induced Gap Edges and Vortices. I. Inviscid Three-dimensional Hydro Disks. *ApJ* **785** (2), 122, arXiv: 1308.0648.