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

Abstract

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.

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 that 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 gases effectively.

The most physically accurate method of evolving dust grains in fluids is an N-body 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 length scales 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, Flock & Vaidya 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 the order of 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, Lin & Dong 2021). In this paper we derive such a fluid model, starting from a stochastic differential equation (SDE) 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) that is dominated by gravity and the dust–gas interactions. 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 (Grad 1948, 1949; Bobylev 1982, 2018; Chapman & Cowling 1990). 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 mixtures there is an additional dust Reynolds stress from the turbulent motion.

Stochastic differential equations have been used to model turbulent motion in fluids (e.g. Pope 1987; Thomson 1987; Sawford 1991; Minier, Peirano & Chibbaro 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, Morfill & Sterzik 1995; Minier 2001; Carballido, Fromang & Papaloizou 2006; Youdin & Lithwick 2007; Minier, Chibbaro & Pope 2014; Minier 2015; Ormel & Liu 2018; Laibe, Bréhier & Lombart 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 SDEs describing the motion of a single dust grain in a turbulent gas. To do this, we preform a moment expansion of the Fokker–Planck equation associated with the SDEs, similar to that preformed by Youdin & Lithwick (2007) but without the restrictive assumption that the correlation time is the shortest time scale 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 six-dimensional (6-D) 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-D 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 § 3 we consider a SDE for motion of a single dust grain in a turbulent gas disc. In § 4 we derive the dust fluid equations by performing a moment expansion of the Fokker–Planck equation associated with the SDE introduced in § 3 and discuss our closure scheme. Sections 5–7 describe 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. Section 7 looks at the behaviour of the dust rheological (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) stress tensor in rotating shear flows. In § 8 we suggest possible refinements that could be made to the model. We present our conclusions in § 9 and further mathematical derivations are given in the appendices.

2. Overview of astrophysical flows

In this section we 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 that dominate the gravitational field. The gas in such a system has the following properties.

1. (i) 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.

2. (ii) 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.

3. (iii) The discs are Rayleigh stable, however, they can exhibit subsonic hydrodynamical or magnetohydrodynamical turbulence. Magnetohydrodynamical turbulence in discs  – due to the magnetorotational instability (MRI) (Balbus & Hawley 1991; Hawley & Balbus 1991; Hawley, Gammie & Balbus 1995) – is much stronger than hydrodynamical turbulence; see, e.g. vertical shear instability (Nelson, Gressel & Umurhan 2013; Lin & Youdin 2015; Flock et al. 2017; Svanberg, Cui & Latter 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.

4. (iv) The disc is stratified with a pressure scale height $H \thickapprox 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 a vertical extent approximately equal to the scale height.

5. (v) Characteristic time scales are the orbital period of the order of $1$ day–$10^3$ years (depending on the position in the disc). Characteristic length scales are the scale height $H \lesssim 0.1 R$ and cylindrical radius $R \sim 0.1\unicode{x2013}100\,\mathrm {AU}$ ${\sim }10^7\unicode{x2013}10^{10}\,\mathrm {km}$.

6. (vi) Molecular viscosity is typically sufficiently low that it can be neglected – although the Kolmogorov scale is the order of $10\,\mathrm {m}$ (Armitage 2020).

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

1. (i) The dust is polydispersed with sizes between the order of micrometres and $10\,\mathrm {cm}$ and forms a near continuous distribution in size space; however, we 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 that is dependent on the observational wavelength.

2. (ii) Dust to gas mass ratio is typically ${\gtrsim }0.01$, with the vast majority of the mass in the largest grains (Testi et al. 2014).

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

4. (iv) The mean free path for dust–dust collisions is ${\sim }10^{5}\,\mathrm {km}$, with the collision time scale 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 $\boldsymbol {v}^{g}$. The position $\boldsymbol {x}$ and velocity $\boldsymbol {v}$ for a dust particle, subject to force per unit mass, $\boldsymbol {f}$, and gas drag are given by the following set of differential equations:

(3.1)$$\begin{gather} {\rm d}{\kern0.9pt}x_i = v_i \,{\rm d} t, \end{gather}$$

(3.2)$$\begin{gather}{\rm d} v_i = f_i \,{\rm d} t - \frac{1}{t_s} (v_i - v_i^{g})\,{\rm d} t. \end{gather}$$

Here $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 = - \boldsymbol {\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 $\boldsymbol {x}$, $\boldsymbol {v}$, $\boldsymbol {v}^{g}$ and $\boldsymbol {f}$, respectively. These are related to the contravariant components $x^i$, $v^i$, $v^i_{g}$ and $f^i$ via the metric tensor $\gamma _{i j}$, where $x_i = \gamma _{i j} x^{j}$, $v_i = \gamma _{i j} v^{j}$, $v_i^{g} = \gamma _{i j} v_{g}^{j}$ and $f_i = \gamma _{i j} 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, Efstathiou & Lasenby (2006) for details).

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

(3.3)\begin{equation} t_{s} = \frac{\rho_{grain} s}{\rho_{g} c_s} \sqrt{\frac{{\rm \pi} \gamma}{8}}, \end{equation}

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

(3.4)\begin{equation} {\rm d} v_i^{g} ={-}\frac{1}{t_c} v_i^{g} \,{\rm d} t + \sqrt{\frac{2 \alpha}{t_c}} c_s \,{\rm d} W_{i}, \end{equation}

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 time scale.

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

(3.5)\begin{equation} {\rm d} v_i^{g} = f_i^{g} \,{\rm d} t - \frac{1}{t_c} (v_i^{g} - u_i^{g}) \,{\rm d} t + \sqrt{\frac{2 \alpha}{t_c}} c_s \,{\rm d} W_i, \end{equation}

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 (A10)–(A12) 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} = -\boldsymbol {\nabla }_i \phi - \rho _g^{-1} \boldsymbol {\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}$, (3.5) amounts to modelling the pressure fluctuation and dissipation terms as being responsible for the Ornstein–Uhlenbeck terms present above (Pope 2000). Here $f_i^g$, $\alpha$, $t_c$, $c_s$ and $\boldsymbol {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 sound speed 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$ and we give a brief discussion of this possibility in § 8.

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

(3.6)$$\begin{gather} {\rm d}{\kern0.9pt}x_i = v_i\,{\rm d} t, \end{gather}$$

(3.7)$$\begin{gather}{\rm d} v_i = f_i \,{\rm d} t - \frac{1}{t_s} (v_i - v_i^{g}) \,{\rm d} t, \end{gather}$$

(3.8)$$\begin{gather}D_d v_i^{g} = f_i^{g} \,{\rm d} t - \frac{1}{t_c} (v_i^{g} - u_i^{g}) \,{\rm d} t + \sqrt{\frac{2 \alpha}{t_c}} c_s \,{\rm d} W_i. \end{gather}$$

Now one can regard 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

(3.9)\begin{equation} D_d v_i^{g} = {\rm d} v_i^{g} - (u^k - u^{k}_{g}) \boldsymbol{\nabla}_{k} u_{i}^{g} \,{\rm d} t . \end{equation}

This can be thought of as a separate ‘advection’ step that, 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 six 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^{k}_{g}) \boldsymbol {\nabla }_{k} u_{i}^{g}$):

(3.10)\begin{equation} \left.\begin{gathered} {\rm d} X_{\alpha} = V_{\alpha} \,{\rm d} t, \\ {\rm d} V_{\alpha} = F_{\alpha} \,{\rm d} t - C_{\alpha \beta} (V^{\beta} - U_g^{\beta}) \,{\rm d} t + \sigma_{\alpha \beta} \,{\rm d} W^{\beta}. \end{gathered}\right\} \end{equation}

Here we have adopted the convention that Greek indices are over the 6-D space and Latin indices are taken over the three-dimensional (3-D) space. These 6-D indices are raised and lowered with a 6-D metric tensor $g_{\alpha \beta }$, constructed from $\gamma _{i j}$, 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 \times 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 6-D 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 is the 6-D form of the last term in (3.8). In addition to simplifying the subsequent derivations, (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 particles in rotating shear flows experience the same turbulence as particles 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-D space

The three 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 6-D system agrees with the 3-D system, the 6-D system must posses 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) 3-D coordinate system.

To construct this coordinate system, we first consider the coordinates of the underlying 3-D system with metric tensor $\gamma _{i j}$ and associated Christoffel symbols $\mathcal {T}_{i j}^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 ${\cdot }^{*} : \{1\cdots 6\} \rightarrow \{1\cdots 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 make use of symmetrising/antisymmetrising operations on the tensor indices with $E_{(\alpha _1 \cdots \alpha _n)}$ and $E_{[\alpha _1 \cdots \alpha _n]}$, for some tensor $\boldsymbol{\mathsf{E}}$, being symmetrisation and antisymmetrisation 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/antisymmetrisation operations, but instead follows the obvious order of operations such that

(3.11)$$\begin{gather} E_{(\alpha, \beta^{*})} = \tfrac{1}{2} ( E_{\alpha \beta^{*}} + E_{\beta^{*} \alpha}), \end{gather}$$

(3.12)$$\begin{gather}E_{(\alpha, \beta)^{*}} = \tfrac{1}{2} ( E_{\alpha \beta^{*}} + E_{\beta \alpha^{*}}), \end{gather}$$

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,

(3.13)\begin{equation} \bar{\cdot} := \int {\cdot} J_{g} \,{\rm d}^3 x_g, \end{equation}

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

(3.14)\begin{equation} \overline{\boldsymbol{\nabla}_{\alpha_g} \boldsymbol{\mathsf{E}}} = 0. \end{equation}

For Cartesian gas displacement coordinates, this integrating out of the non-physical space is straightforward. 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 3-D coordinate system (e.g. when calculating the gas drag). It is instead useful to set-up the geometry of our 6-D space so that the rotation happens automatically. To do this, we first introduce the metric tensor of the 6-D coordinate system:

(3.15)\begin{equation} g_{\alpha \beta} = \begin{cases} \gamma_{\alpha \beta} , & \alpha, \beta \in \{1,2,3\} , \\ \gamma_{\alpha^{*} \beta^{*}} , & \alpha, \beta \in \{4,5,6\} , \\ 0 , & \mathrm{otherwise}. \end{cases}\end{equation}

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

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

2. (ii) Translational invariance with respect to the gas displacement such that $\bar {\boldsymbol {\nabla }}_{\alpha _g} \boldsymbol{\mathsf{E}} (\boldsymbol {x}_d) = 0$ for tensoral quantity $\boldsymbol{\mathsf{E}}$.

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

Property (ii) ensures that $\bar {\boldsymbol {\nabla }}_{\alpha _g} \boldsymbol{\mathsf{E}} (\boldsymbol {x}_d) = \overline {\boldsymbol {\nabla }_{\alpha _g} \boldsymbol{\mathsf{E}}}$, where $\boldsymbol {\nabla }_i$ is the covariant derivative, and allows us to carry out the integral over the dummy gas directions by replacing covariant derivatives with $\bar {\boldsymbol {\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 6-D vectors. This can be seen considering $\boldsymbol {A} = \boldsymbol {U}$ and considering the action of the Lagrangian time derivative, $D = \partial _{t} + \boldsymbol {U} \boldsymbol {\cdot } \bar {\boldsymbol {\nabla }}$, on the vectors $\boldsymbol {B}$ and $\tilde {\boldsymbol {B}}$. As $B^{\alpha } = \tilde {B}^{\alpha ^{*}}$, one requires that $(D \boldsymbol {B})^{\alpha } = (D \tilde {\boldsymbol {B}})^{\alpha ^{*}}$, which requires condition (iii) as $\boldsymbol {U}$ is arbitrary and $(\partial _t \boldsymbol {B})^{\alpha } = (\partial _t \tilde {\boldsymbol {B}})^{\alpha ^{*}}$.

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

(3.16a,b)\begin{equation} \bar{\boldsymbol{\nabla}}_{\alpha_d} \hat{e}_{\beta_d} = \mathcal{T}_{\alpha_d \beta_d}^{\gamma_d} \hat{e}_{\gamma_d},\quad \bar{\boldsymbol{\nabla}}_{\alpha_d} \hat{e}_{\beta_g} = \mathcal{T}_{\alpha_d \beta_g^{*}}^{\gamma_g^{*}} \hat{e}_{\gamma_g}, \end{equation}

with $\bar {\boldsymbol {\nabla }}_{\alpha _g} \hat {e}_{\beta } = 0$. As $\mathcal {T}_{i j}^{k}$ are the Christoffel symbol components for the 3-D coordinate system associated with the metric $\gamma _{i j}$, it is straightforward to show that this connection satisfies property (i). Property (ii) follows from $\bar {\boldsymbol {\nabla }}_{\alpha _g} \hat {e}_{\beta } = 0$. Finally, for property (iii),

(3.17)\begin{gather} (\boldsymbol{A} \boldsymbol{\cdot} \bar{\boldsymbol{\nabla}} \boldsymbol{B})^{\beta} = A^{\alpha} \bar{\boldsymbol{\nabla}}_{\alpha} B^{\beta} = A^{\alpha} \partial_{\alpha} B^{\beta} + A^{\alpha} \mathcal{T}_{\alpha \gamma}^{\beta} B^{\gamma}, \end{gather}

(3.18)\begin{gather}\begin{aligned}[b] (\boldsymbol{A} \boldsymbol{\cdot} \bar{\boldsymbol{\nabla}} \tilde{\boldsymbol{B}})^{\beta} = A^{\alpha} \bar{\boldsymbol{\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^{*}} = (\boldsymbol{A} \boldsymbol{\cdot} \bar{\boldsymbol{\nabla}} \boldsymbol{B})^{\beta*} .\end{aligned} \end{gather}

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 is not the Levi-Civita connection). This torsion arises when $\bar {\boldsymbol {\nabla }}_{\alpha _d} \hat {e}_{\beta _g} \ne 0$ as $\bar {\boldsymbol {\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

(3.19)\begin{equation} S_{\alpha \beta}^{\gamma} \hat{e}_{\gamma} = \bar{\boldsymbol{\nabla}}_{\alpha} \hat{e}_{\beta} - \bar{\boldsymbol{\nabla}}_{\beta} \hat{e}_{\alpha}, \end{equation}

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

(3.20)\begin{equation} S_{\alpha_d \beta_g}^{\gamma_g} ={-} S_{\beta_g \alpha_d}^{\gamma_g} = \mathcal{T}_{\alpha_d \beta_g^{*}}^{\gamma_g^{*}}, \end{equation}

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}$,

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

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

(3.22)\begin{equation} F_{\alpha} = \begin{cases} f_{\alpha},\quad \beta \in \{1, 2, 3 \}, \\[3pt] f^g_{\alpha^{*}} + (U^{\beta} - U^{\beta}_{g}) \boldsymbol{\nabla}_{\beta} U_{\alpha}^{g},\quad \alpha \in \{4, 5, 6 \}, \end{cases} \end{equation}

while the drag tensor is

(3.23)\begin{equation} C_{\alpha \beta} = \begin{cases} \dfrac{1}{t_s} g_{\alpha \beta} , & \alpha, \beta \in \{1,2,3\} , \\[10pt] -\dfrac{1}{t_s} g_{\alpha \beta^{*}} , & \alpha \in \{1,2,3\},\quad \beta \in \{4,5,6 \} , \\[6pt] 0 , & \alpha \in \{4,5,6\} , \quad \beta \in \{1,2,3 \} , \\[5pt] \dfrac{1}{t_c} g_{\alpha \beta} , & \alpha, \beta \in \{4,5,6\}, \end{cases} \end{equation}

while the diffusion tensor is

(3.24)\begin{equation} D_{\alpha \beta} = \begin{cases} \dfrac{\alpha c_s^2}{t_c} g_{\alpha \beta} , & \alpha, \beta \in \{4,5,6\} , \\[8pt] 0 , & \mathrm{otherwise} . \end{cases} \end{equation}

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

(3.25)\begin{equation} D_{\alpha \beta} = \begin{cases} \dfrac{c_s^2}{t_c} a_{\alpha \beta} , & \alpha, \beta \in \{4,5,6\} , \\[8pt] 0 , & \mathrm{otherwise}. \end{cases} \end{equation}

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 associated with (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(\boldsymbol {X},\boldsymbol {V})$. By use of Ito's lemma this evolves according to

(4.1)\begin{equation} {\rm d} A = \frac{\partial A}{\partial X_\alpha} {\rm d} X_\alpha + \frac{\partial A}{\partial V_{\alpha}} {\rm d} V_{\alpha} + \frac{1}{2} \frac{\partial^2 A}{\partial V_\alpha \partial V_{\beta}} \langle {\rm d} V_{\alpha}, {\rm d} V_{\beta} \rangle, \end{equation}

where the angle bracket $\langle {\cdot }, {\cdot } \rangle$ denotes the covariance. The covariance of a Wiener process $d W_{\alpha }$ is given by

(4.2)\begin{equation} \langle {\rm d} W^{\alpha}, {\rm d} W^{\beta} \rangle = g^{\alpha \beta} \,{\rm d} t . \end{equation}

This leads to the following covariance of velocity:

(4.3)\begin{equation} \langle {\rm d} V_{\alpha}, {\rm d} V_{\beta} \rangle = \sigma_{\alpha \mu} \sigma_{\beta \nu} \langle {\rm d} W^{\mu} {\rm d} W^{\nu} \rangle = 2 D_{\alpha \beta} \,{\rm d} t. \end{equation}

Here we have introduced the diffusion tensor, $D_{\alpha \beta } = \frac {1}{2} g^{\mu \nu } \sigma _{\alpha \mu } \sigma _{\beta \nu }$. Substituting (3.10) into (4.1), we arrive at

(4.4)\begin{align} {\rm d} A &= \frac{\partial A}{\partial X_\alpha} V_{\alpha} \,{\rm d} t + \frac{\partial A}{\partial V_{\alpha}} [F_{\alpha} - C_{\alpha \beta} (U^{\beta} - U_g^{\beta})] \,{\rm d} t \nonumber\\ &\quad +\frac{\partial A}{\partial V_{\alpha}} \sigma_{\alpha \beta} \,{\rm d} W^{\beta} + D_{\alpha \beta} \frac{\partial^2 A}{\partial V_\alpha \partial V_{\beta}}{\rm d} t . \end{align}

The expectation of $A$ is given by

(4.5)\begin{equation} \mathbb{E} [A] = \int p^{L} (\boldsymbol{X},\boldsymbol{V},t,\boldsymbol{X}_0,\boldsymbol{V}_0,t_0) A(\boldsymbol{X},\boldsymbol{V}) \,{\rm d}^6 \boldsymbol{X} \,{\rm d}^6 \boldsymbol{V}, \end{equation}

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

(4.6)\begin{equation} \mathbb{E} [{\rm d} A] = \int {\rm d} p^{L} (\boldsymbol{X},\boldsymbol{V},t,\boldsymbol{X}_0,\boldsymbol{V}_0,t_0) A(\boldsymbol{X},\boldsymbol{V}) \,{\rm d}^6 \boldsymbol{X} \,{\rm d}^6 \boldsymbol{V} . \end{equation}

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

(4.7)\begin{align} & \int A \left\{{\rm d} p^{L} + \frac{\partial}{\partial X_{\alpha}} (p^{L} V_{\alpha}) \,{\rm d} t + \frac{\partial}{\partial V_{\alpha}} [(F_{\alpha} - C_{\alpha \beta} (U^{\beta} - U_g^{\beta})) p^{L}]\,{\rm d} t \right. \nonumber\\ &\quad \left.- \,D_{\alpha \beta}\frac{\partial^2 p^{L}}{\partial V_{\alpha} V_{\beta} } {\rm d} t \right\} {\rm d}^6 \boldsymbol{X} \,{\rm d}^6 \boldsymbol{V}, \end{align}

provided that $\int _{\partial } p^{L} A \boldsymbol {V} \boldsymbol {\cdot } \textrm {d} \boldsymbol {S} \,\textrm {d}^6 \boldsymbol {V} = 0$, where $\int _{\partial }\textrm {d} \boldsymbol {S}$ 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$,

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

This equation gives an evolutionary equation for the Lagrangian transition probability density function (PDF), describing the probability of finding a particle at $\boldsymbol {X}$, $\boldsymbol {V}$ at time $t$ conditional on it being located at $\boldsymbol {X}_0, \boldsymbol {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 (\boldsymbol {X}, \boldsymbol {V}, t)$, which is the expected mass density of particles at $\boldsymbol {X}$, $\boldsymbol {V}$ at time $t$ (Pope 1985, 2000; Minier & Peirano 2001). This will contain contributions from particles with different initial conditions $(\boldsymbol {X}_0, \boldsymbol {V}_0)$, arriving from differing trajectories. This can be obtained from the Eulerian MDF at $t_0$, $p (\boldsymbol {X}_0, \boldsymbol {V}_0, t_0)$, by using the transition PDF and integrating over the initial positions and velocities (Pope 1985, 2000; Minier & Peirano 2001):

(4.9)\begin{equation} p (\boldsymbol{X}, \boldsymbol{V}, t) = \int p^{L} (\boldsymbol{X},\boldsymbol{V},t,\boldsymbol{X}_0,\boldsymbol{V}_0,t_0) p (\boldsymbol{X}_0, \boldsymbol{V}_0, t_0) \,{\rm d}^6 \boldsymbol{X}_0 \,{\rm d}^6 \boldsymbol{V}. \end{equation}

We can obtain the Fokker–Planck equation for $p$ by multiplying (4.8) by $p (\boldsymbol {X}_0, \boldsymbol {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

(4.10)\begin{equation} \frac{\partial p}{\partial t} + \frac{\partial}{\partial X^{\alpha}} (p V^{\alpha}) + \frac{\partial}{\partial V_{\alpha}} [(F_{\alpha} - C_{\alpha \beta} (V^{\beta} - U_g^{\beta})) p] = D_{\alpha \beta} \frac{\partial^2 p}{\partial V_{\alpha} V_{\beta} }. \end{equation}

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 manor to that done in kinetic theory. In performing this moment expansion we wish to arrive at a set of partial differential equations (PDEs) in space and time from the initial PDE in $(t,\boldsymbol {X},\boldsymbol {V})$. This means we need to compute a moment expansion in $\boldsymbol {V}$. A similar procedure was carried out by Youdin & Lithwick (2007). Defining the velocity moments of $p$ as follows:

(4.11)$$\begin{gather} \rho_6 := \int p \,{\rm d}^6 \boldsymbol{V}, \end{gather}$$

(4.12)$$\begin{gather}\rho_6 U_{\alpha} := \int p V_{\alpha} \,{\rm d}^6 \boldsymbol{V}, \end{gather}$$

(4.13)$$\begin{gather}\varPi_{\alpha \beta} := \int p (V_{\alpha} - U_{\alpha}) (V_{\beta} - U_{\beta}) \,{\rm d}^6 \boldsymbol{V}, \end{gather}$$

(4.14)$$\begin{gather}\varPi_{\alpha_1 \cdots \alpha_k} := \int p (V_{\alpha} - U_{\alpha}) \cdots (V_{\alpha_k} - U_{\alpha_k}) \,{\rm d}^6 \boldsymbol{V}. \end{gather}$$

Note that this moment expansion is in the 6-D space, so that $\rho _6$ is the 6-D mass density and $\varPi _{\alpha \beta }$ is the 6-D rheological stress tensor. (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, Desjardins & Fox 2016b).) Furthermore, we have chosen a normalisation such that

(4.15)\begin{equation} \int p \,{\rm d}^6 \boldsymbol{V} \,{\rm d}^3 \boldsymbol{x}_{g} = \rho_d, \end{equation}

where $\rho _d$ is the dust density (i.e. the density of the dust phase, this is equal to the grain density, $\rho _{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 \varPi _{\alpha \beta } \,\textrm {d}^3 \boldsymbol {x}_{g}$ has the same units as the gas pressure.

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

(4.16)\begin{equation} \dot{\rho}_6 + \boldsymbol{\nabla}_{\alpha} [\rho_6 U^{\alpha}] = 0. \end{equation}

The first $\boldsymbol {V}$ moment of (4.10) leads to the (6-D) dust momentum equation,

(4.17)\begin{equation} \frac{\partial}{\partial t} [\rho_6 U_{\alpha}] + \boldsymbol{\nabla}^{\beta} [\varPi_{\alpha \beta} + \rho_6 U_{\alpha} U_{\beta}] - \rho_6 F_{\alpha} + \rho_6 C_{\alpha \beta} (U^{\beta} - U_g^{\beta}) = 0 . \end{equation}

Taking the second $\boldsymbol {V}$ moment yields a constitutive relation for the (6-D) dust stress tensor,

(4.18)\begin{align} & \frac{\partial}{\partial t} [\varPi_{\alpha \beta} + \rho_6 U_{\alpha} U_{\beta}] + \boldsymbol{\nabla}^{\gamma} [\varPi_{\alpha \beta \gamma} + 3 U_{(\alpha} \varPi_{\beta \gamma)} + \rho_6 U_{\alpha} U_{\beta} U_{\gamma}] \nonumber\\ &\quad -2 \rho_6 U_{(\alpha} [F_{\beta)} - C_{\beta) \gamma} (U^{\gamma} - U_g^{\gamma})] + 2 \varPi^{\gamma}_{(\alpha} C_{\beta) \gamma}= 2 \rho_6 D_{(\alpha \beta)}. \end{align}

Here we have made use of the notation for the symmetrisation of the tensor indices. As we make extensive use of this notation, we give explicit expressions for the symmetrised terms in the above equation as a illustrative example, $U_{(\alpha } \varPi _{\beta \gamma )} = \frac {1}{3} (U_{\alpha } \varPi _{\beta \gamma } + U_{\beta } \varPi _{\alpha \gamma } + U_{\gamma } \varPi _{\alpha \beta } )$ and $2 U_{(\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 manor. 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

(4.19) \begin{align} & \frac{\partial \varPi_{\alpha_1 \cdots \alpha_k}}{\partial t} + k \varPi_{(\alpha_1 \cdots \alpha_{k-1}} [D U_{\alpha_k} + \boldsymbol{\nabla}_{\alpha_k )} \phi + C_{\alpha_k)}^{\gamma} (U_{\gamma} - U_{\gamma}^{g})] \nonumber\\ &\quad +\boldsymbol{\nabla}^{\sigma} [\varPi_{\alpha_1 \cdots \alpha_{k} \sigma} + U_{\sigma} \varPi_{\alpha_1 \cdots \alpha_k}] \nonumber\\ &\quad + k \varPi_{\sigma (\alpha_1 \cdots \alpha_{k-1}} [\boldsymbol{\nabla}^{\sigma} U_{\alpha_k)} + C_{\alpha_k)}^{\sigma}] = k (k - 1) \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}. \end{align}

Making use of the dust momentum equation this simplifies to

(4.20)\begin{align} & (D + \boldsymbol{\nabla}_{\sigma} U^{\sigma}) \varPi_{\alpha_1 \cdots \alpha_{k}} + k \varPi_{\sigma ( \alpha_1 \cdots \alpha_{k-1}} \boldsymbol{\nabla}^{\sigma}U_{\alpha_k)} + \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_1 \cdots \alpha_k \sigma} \nonumber\\ &\quad ={-} k \left[\varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} C_{\alpha_k) \sigma} - \frac{1}{\rho} \varPi_{(\alpha_1 \cdots \alpha_{k-1}} \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_{k } ) \sigma} - (k - 1) \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}\right]. \end{align}

Taking $k=2$ in the above equation we recover the constitutive relation for $\varPi _{\alpha \beta }$ (to obtain this, we note that $\varPi _{\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

(4.21)$$\begin{gather} \mathcal{D} \varPi_{\alpha_1 \cdots \alpha_k} = D \varPi_{\alpha_1 \cdots \alpha_k} +k \varPi_{\gamma (\alpha_1 \cdots \alpha_{k - 1}} \boldsymbol{\nabla}_{\alpha_k)} U^{\gamma} + \varPi_{\alpha_1 \cdots \alpha_k} \boldsymbol{\nabla}_{\gamma} U^{\gamma}, \end{gather}$$

(4.22)$$\begin{gather}\mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_k} = D \varPi_{\alpha_1 \cdots \alpha_k} + \varPi_{\alpha_1 \cdots \alpha_k} \boldsymbol{\nabla}_{\gamma} U^{\gamma}, \end{gather}$$

(4.23)$$\begin{gather}\mathcal{D}_2 \varPi_{\alpha_1 \cdots \alpha_k} = D \varPi_{\alpha_1 \cdots \alpha_k} + k \varPi_{\gamma (\alpha_1 \cdots \alpha_{k - 1}} \boldsymbol{\nabla}^{\gamma} U_{\alpha_k)} + \varPi_{\alpha_1 \cdots \alpha_k} \boldsymbol{\nabla}_{\gamma} U^{\gamma}. \end{gather}$$

The first of these is closely related to the convective Maxwell derivative, with $\mathcal {D} \varPi _{\alpha _1 \cdots \alpha _k = 0}$ implying that the tensoral quantity $\rho _6^{-1} \varPi _{\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 6-D flow.

Rearranging the continuity, momentum and constitutive equations, and making use of the operator $\mathcal {D}_2$, we obtain

(4.24)$$\begin{gather} D \rho_6 ={-} \rho_6 \boldsymbol{\nabla}_{\alpha} u^{\alpha}, \end{gather}$$

(4.25)$$\begin{gather}\rho_6 D U_{\alpha} = \rho_6 F_{\alpha} - \boldsymbol{\nabla}^{\beta} \varPi_{\alpha \beta} - \rho_6 C_{\alpha \beta} (U^{\beta} - U_g^{\beta}), \end{gather}$$

(4.26)$$\begin{gather}\mathcal{D}_2 \varPi_{\alpha \beta} ={-}\boldsymbol{\nabla}^{\gamma} \varPi_{\alpha \beta \gamma} -2 (\varPi^{\gamma}_{(\alpha} C_{\beta) \gamma} - \rho_6 D_{(\alpha \beta)}). \end{gather}$$

As the right-hand side of (4.26) is symmetrised, this ensures that $\varPi _{\alpha \beta }$ remains symmetric for symmetric initial conditions. Using a similar argument to that advanced in Ogilvie (2003) and Lynch & Ogilvie (2021), $\varPi _{\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

(4.27)\begin{align} \mathcal{D}_2 \varPi_{\alpha_1 \cdots \alpha_k} &={-}\boldsymbol{\nabla}^{\gamma}\varPi_{\gamma \alpha_1 \cdots \alpha_k} - k \left[\varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} C_{\alpha_k) \sigma } - \frac{1}{\rho} \varPi_{(\alpha_1 \cdots \alpha_{k-1}} \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_{k }) \sigma} \right. \nonumber\\ &\quad \left.\vphantom{\frac{1}{\rho}} - (k - 1) \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}\right]. \end{align}

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

(4.28)\begin{equation} \mathcal{D} \varPi_{\alpha \beta} ={-}2 (\varPi^{\gamma}_{(\alpha} A_{\beta) \gamma} - \rho_6 D_{(\alpha \beta)}). \end{equation}

Here we have defined

(4.29)\begin{equation} A_{\alpha \beta} = C_{\alpha \beta} - \omega^{\gamma} \varepsilon_{\alpha \beta \gamma}, \end{equation}

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

(4.30)\begin{equation} \omega^{\gamma} \varepsilon_{\gamma \alpha \beta} = \boldsymbol{\nabla}_{\alpha} U_{\beta} - \boldsymbol{\nabla}_{\beta} U_{\alpha} . \end{equation}

The evolutionary equation for the $k$th velocity moment can be similarly rewritten. In the full 6-D model, with the Levi-Civita connection, (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 (4.28) act like sources/sinks for the stress tensor. Unfortunately, in the presence of torsion the constitutive equation based on (4.28) ends up more complicated to manipulate than that based on (4.26) owing to the addition of terms involving the torsion tensor. As such, we stick to (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$,

(4.31)\begin{align} \mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_{k}} &={-}\boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_1 \cdots \alpha_k \sigma} -k \left[ \varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma } \vphantom{\frac{1}{\rho}}\right.\nonumber\\ &\quad \left.-\frac{1}{\rho} \varPi_{(\alpha_1 \cdots \alpha_{k-1}} \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_{k } ) \sigma} - (k - 1) \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}\right], \end{align}

where we have introduced $B_{\alpha \beta } = C_{\alpha \beta } + \boldsymbol {\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 show that when the fluid is thermally stable, and the turbulent velocity small relative to the fluid velocity, the third velocity moment typically decays until it is asymptotically small relative to the stress tensor, we can therefore drop the $\boldsymbol {\nabla }^{\gamma } \varPi _{\alpha \beta \gamma }$ in the constitutive relation and close the moment expansion at the second velocity moment.

4.3.1. Well-coupled ordering scheme

Previous authors have noted that when the dust is well coupled to the gas $({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 ${St} = O(\epsilon )$. We consider units such that $U^{\alpha } = O(1)$, $\mathcal {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 $\boldsymbol {\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 the rescaled velocity moment $\tilde {\varPi }_{\alpha _1 \cdots \alpha _k}$, such that $\varPi _{\alpha _1 \cdots \alpha _k} = \epsilon ^{\delta _k} \tilde {\varPi }_{\alpha _1 \cdots \alpha _k}$, and the stretched/rescaled variable $\tilde {X} = X/\epsilon$, such that $\boldsymbol {\nabla } = \epsilon ^{-1} \tilde {\boldsymbol {\nabla }}$, then we arrive at a rescaled equation for the $k$th velocity moment:

(4.32)\begin{align} \epsilon^{\delta_k} \mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_{k}} &={-} \epsilon^{\delta_{k+1} - 1} \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_1 \cdots \alpha_k \sigma} - k \left[\epsilon^{\delta_k - 1} \varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma } \vphantom{\frac{\epsilon^{\delta_{k-1} + \delta_2 -1}}{\rho}}\right. \nonumber\\ &\quad \left.-\frac{\epsilon^{\delta_{k-1} + \delta_2 -1}}{\rho} \varPi_{(\alpha_1 \cdots \alpha_{k-1}} \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_{k } ) \sigma} - (k - 1) \epsilon^{\delta_{k-2} + 2} \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}\right] . \end{align}

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

(4.33)\begin{align} & \epsilon \mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_{k}} + k [\varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma } - (k - 1) \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}] \nonumber\\ &\quad =\epsilon^{3} \left[ -\boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_1 \cdots \alpha_k \sigma} + \frac{k}{\rho} \varPi_{(\alpha_1 \cdots \alpha_{k-1}}\boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_{k } ) \sigma} \right]. \end{align}

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

(4.34)\begin{align} \epsilon \mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_{k}} &={-} \boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_1 \cdots \alpha_k \sigma} - k \left[\varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma } \vphantom{\frac{1}{\rho}}\right. \nonumber\\ &\quad \left.-\frac{1}{\rho} \varPi_{(\alpha_1 \cdots \alpha_{k-1}}\boldsymbol{\nabla}^{\sigma} \varPi_{\alpha_{k } ) \sigma} - (k - 1) \varPi_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})} \right]. \end{align}

Thus, we find that the correction to the evolutionary equation for the second velocity moment $\varPi _{\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 \varPi _{\alpha \beta }$) that will no longer be negligible.

According to (4.34) the evolution of the third velocity moment will depend on gradients of the fourth 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,

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

with second velocity moment

(4.36)\begin{equation} W_{\alpha \beta} = \int (V_{\alpha} - U_{\alpha}) (V_{\beta} - U_{\beta}) f \,{\rm d}^n V . \end{equation}

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

(4.37)\begin{equation} W_{\alpha_1 \cdots \alpha_k} = \int (V_{\alpha_1} - U_{\alpha_1}) \cdots (V_{\alpha_k} - U_{\alpha_k}) f \,{\rm d}^n V . \end{equation}

For odd $k$, $W_{\alpha _1 \cdots \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:

(4.38)\begin{equation} W_{\alpha_1 \cdots \alpha_k} = (k-1) W_{(\alpha_1 \cdots \alpha_{k-2}} W_{\alpha_{k-1}) \alpha_k} . \end{equation}

By symmetry of the velocity moments we also have $W_{\alpha _1 \cdots \alpha _k} = W_{(\alpha _1 \cdots \alpha _k)} = (k-1) W_{(\alpha _1 \cdots \alpha _{k-2}} W_{\alpha _{k-1} \alpha _k)}$.

Starting from the assumption that the second velocity moment evolves according to

(4.39)\begin{equation} D W_{\alpha_1 \alpha_2} ={-} 2 [W^{\sigma}_{(\alpha_1} B_{\alpha_2) \sigma} - D_{\alpha_{1} \alpha_{2}}], \end{equation}

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

(4.40)\begin{equation} D W_{\alpha_1 \cdots \alpha_k} ={-} k [W^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma} - (k-1) W_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k)}}]. \end{equation}

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

(4.41)\begin{align} D W_{\alpha_1 \cdots \alpha_k} &= (k-1) W_{(\alpha_{1} \alpha_2} D W_{\alpha_3 \cdots \alpha_{k})} + (k-1) W_{(\alpha_1 \cdots \alpha_{k-2}} D W_{\alpha_{k-1} \alpha_k)} \nonumber\\ &={-}(k-1) (k-2) W_{(\alpha_1 \alpha_2} [W^{\sigma}_{\alpha_3 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma} - (k-3) W_{\alpha_3 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}] \nonumber\\ &\quad -2 (k-1) W_{(\alpha_1 \cdots \alpha_{k-2}} [W^{\sigma}_{\alpha_{k-1}} B_{\alpha_k) \sigma} - D_{\alpha_{k-1} \alpha_{k})}] \nonumber\\ &={-}(k-1) [(k-2) W_{\sigma (\alpha_1 \cdots \alpha_{k-3}} B_{\alpha_{k-2}}^{\sigma} W_{\alpha_{k-1} \alpha_k)} + 2 W_{(\alpha_1 \cdots \alpha_{k-2} } W^{\sigma}_{\alpha_{k-1}} B_{\alpha_k) \sigma}] \nonumber\\ &\quad +(k-1) (k-3) [(k-2) W_{(\alpha_1 \cdots \alpha_{k-4}} D_{\alpha_{k-3} \alpha_{k-2}} W_{\alpha_{k-1} \alpha_k)} \nonumber\\ &\quad + 2 W_{(\alpha_1 \cdots \alpha_{k-4}} W_{\alpha_{k-3} \alpha_{k-2}} D_{\alpha_{k-1} \alpha_k)}] \nonumber\\ &={-} k [W^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma} - (k-1) W_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k)}}]. \end{align}

Thus, we see that if the $(k-2)$th velocity moment evolves according to (4.40) and the second velocity moment evolves according to (4.39), then the $k$th velocity moment also evolves according to (4.40). Starting with the fourth velocity moment we see that, given (4.39) and (4.41), it evolves according to (4.40). We can thus proceed by induction to arbitrary $k$, and conclude that $W_{\alpha _1 \cdots \alpha _k}$ evolve according to (4.40).

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

(4.42)$$\begin{gather} \varPi_{\alpha_1 \cdots \alpha_k} = \epsilon^{k} \rho (\boldsymbol{\xi}, \boldsymbol{x}) W_{\alpha_1 \cdots \alpha_k} (\boldsymbol{x}) + \epsilon^{k + 1} \varSigma_{\alpha_1 \cdots \alpha_k} (\boldsymbol{\xi}, \boldsymbol{x}), \end{gather}$$

(4.43)$$\begin{gather}U = U_0 (\boldsymbol{x}) + \epsilon^{\kappa} u_0 (\boldsymbol{\xi}, \boldsymbol{x}), \end{gather}$$

(4.44)$$\begin{gather}\boldsymbol{\nabla} = \epsilon^{{-}1} \frac{\partial}{\partial \boldsymbol{\xi}}+ \frac{\partial}{\partial \boldsymbol{x}}. \end{gather}$$

Here $\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 $\varSigma _{\alpha _1 \cdots \alpha _k}$. While we can often treat $\kappa = 2$ (i.e. the part of the mean velocity that varies on the dust length scale is $O(\epsilon ^2)$), we assume that $\kappa = 1$ throughout as this will allow for a wider range of dust flows.

Substituting (4.42)–(4.44) into (4.31) and making use of $\mathcal {D}_1 \rho = 0$, the evolutionary equation for the $k$th velocity moment becomes

(4.45)\begin{align} \mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_k} &= \epsilon^{k} \rho D W_{\alpha_1 \cdots \alpha_k} + \epsilon^{k+1} u_0^{\alpha} \frac{\partial}{\partial x^{\alpha}} W_{\alpha_1 \cdots \alpha_k} + \epsilon^{k+1} \mathcal{D}_1 \varSigma_{\alpha_1 \cdots \alpha_k} \nonumber\\ &={-}\epsilon^{k} \left( \frac{\partial}{\partial \xi^{\sigma}} + \epsilon \frac{\partial}{\partial x^{\sigma}} \right) (\rho W_{\alpha_1 \cdots \alpha_k \sigma}) - \epsilon^{k+1} \left( \frac{\partial}{\partial \xi^{\sigma}} + \epsilon \frac{\partial}{\partial x^{\sigma}} \right) \varSigma_{\alpha_1 \cdots \alpha_k \sigma} \nonumber\\ &\qquad - k \epsilon^{k} \left[(\rho W^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} + \epsilon \varSigma^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}}) B_{\alpha_k) \sigma} \vphantom{\left(\frac{\partial}{\partial \xi^{\sigma}} + \epsilon \frac{\partial}{\partial x^{\sigma}} \right)}\right.\nonumber\\ &\qquad - (\rho W_{(\alpha_1 \cdots \alpha_{k-1}} + \epsilon \varSigma_{(\alpha_1 \cdots \alpha_{k-1}}) \left(\frac{\partial}{\partial \xi^{\sigma}} + \epsilon \frac{\partial}{\partial x^{\sigma}} \right) (\rho W_{\alpha_k) \sigma} + \epsilon \varSigma_{\alpha_k) \sigma}) \nonumber\\ &\qquad \left.\vphantom{\left(\frac{\partial}{\partial \xi^{\sigma}} + \epsilon \frac{\partial}{\partial x^{\sigma}} \right)} -(k-1) (\rho W_{(\alpha_1 \cdots \alpha_{k-2}} + \epsilon \varSigma_{(\alpha_1 \cdots \alpha_{k-2}}) D_{\alpha_{k-1} \alpha_{k})} \right], \end{align}

where, here, $D = \partial _t + U_0^{\alpha } \boldsymbol {\nabla }_{\alpha }$ is the Lagrangian time derivative with respect to the leading-order flow described by $U_0$.

Making use of (4.40) for the evolution of $W_{\alpha _1 \cdots \alpha _k}$, along with the recurrence relation for $W_{\alpha _1 \cdots \alpha _k}$ (4.38) and rearranging we obtain an equation for the evolution of the non-Maxwellian part of the velocity moment:

(4.46)\begin{align} & \mathcal{D}_1 \varSigma_{\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}} \varSigma_{\alpha_1 \cdots \alpha_k \sigma} \nonumber\\ &\qquad + k \left[\varSigma^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma} - \rho W_{(\alpha_1 \cdots \alpha_{k-1}} \frac{\partial}{\partial \xi^{\sigma}} \varSigma_{\alpha_k) \sigma} \right.\nonumber\\ &\qquad \left.-\varSigma_{(\alpha_1 \cdots \alpha_{k-1}} \frac{\partial}{\partial \xi^{\sigma}} \rho W_{\alpha_k) \sigma} - (k-1) \varSigma_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1}\alpha_{k})}\right] \nonumber\\ &\quad = \epsilon \left[k \rho W_{(\alpha_1 \cdots \alpha_{k-1}} \frac{\partial}{\partial x^{\sigma}} \varSigma_{\alpha_k) \sigma} + k \varSigma_{(\alpha_1 \cdots \alpha_{k-1}} \frac{\partial}{\partial x^{\sigma}} (\rho W_{\alpha_k) \sigma} + \epsilon \varSigma_{\alpha_k) \sigma}) \right.\nonumber\\ &\qquad \left.+\,k \varSigma_{(\alpha_1 \cdots \alpha_{k-1}} \frac{\partial}{\partial \xi^{\sigma}} \varSigma_{\alpha_k) \sigma} - \frac{\partial}{\partial x^{\sigma}} \varSigma_{\alpha_1 \cdots \alpha_k \sigma}\right]. \end{align}

Here the terms on the right-hand side are all subleading. Dropping these subleading terms we obtain

(4.47)\begin{align} 0 &= \mathcal{D}_1 \varSigma_{\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}} \varSigma_{\alpha_1 \cdots \alpha_k \sigma} \nonumber\\ &\quad + k \left[\varSigma^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma} - \rho W_{(\alpha_1 \cdots \alpha_{k-1}} \frac{\partial}{\partial \xi^{\sigma}} \varSigma_{\alpha_k) \sigma} \right.\nonumber\\ &\quad \left.-\varSigma_{(\alpha_1 \cdots \alpha_{k-1}} W_{\alpha_k) \sigma} \frac{\partial}{\partial \xi^{\sigma}} \rho - (k-1) \varSigma_{(\alpha_1 \cdots \alpha_{k-2}} D_{\alpha_{k-1} \alpha_{k})}\right]. \end{align}

This confirms that the asymptotic ordering scheme (4.42)–(4.44) is self-consistent and the non-Maxwellian terms are suppressed by a factor of $\epsilon \sim L_{dust}/L_{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

(4.48)\begin{align} \boldsymbol{\nabla}_{\beta} \varPi^{\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) \varSigma_{\beta \alpha} \nonumber\\ &= \epsilon W^{\beta \alpha} \frac{\partial}{\partial \xi^{\beta}} \rho + O(\epsilon^{2}) . \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 ${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_{dust} \ll L_{gas}$ (with the length scale 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 that are not strongly stirred, in the presence of vertical gravity, as these would be expected to settle into a hydrostatically supported dust layer that 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, allows for a fluid treatment of the dust, is beyond the scope of this work; in this section we present an argument showing that velocity moments that 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 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 the $k$th velocity moment are those involving $\varPi _{\alpha _1 \cdots \alpha _k}$. Thus, evolution of the $k$th velocity moment is approximately described by

(4.49)\begin{equation} \mathcal{D}_1 \varPi_{\alpha_1 \cdots \alpha_k} ={-}k \varPi^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma}. \end{equation}

Defining $W_{\alpha _1 \cdots \alpha _k} = \rho ^{-1} \varPi _{\alpha _1 \cdots \alpha _k}$, (4.49) simplifies to

(4.50)\begin{equation} D W_{\alpha_1 \cdots \alpha_k} ={-}k W^{\sigma}_{(\alpha_1 \cdots \alpha_{k-1}} B_{\alpha_k) \sigma} . \end{equation}

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,

(4.51)\begin{equation} D Y^{\alpha} = B_{\beta}^{\alpha} Y^{\beta}, \end{equation}

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 (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

(4.52)\begin{align} D \chi &= Y^{\alpha_{n+1}} \cdots Y^{\alpha_{k}} D Q_{\alpha_{n+1} \cdots \alpha_{k}} + (k - n) Q_{\alpha_{n+1} \cdots \alpha_{k}^{d}}Y^{\alpha_{n+1}} \cdots Y^{\alpha_{k-1}} D Y^{\alpha_{k}} \nonumber\\ &={-} (k - n) Y^{\alpha_{n+1}} \cdots Y^{\alpha_{k}} Q_{\sigma \alpha_{n+1} \cdots \alpha_{k-1} } M_{\alpha_{k}}^{\sigma} \nonumber\\ &\quad + (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} \nonumber\\ &= 0. \end{align}

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 $|\boldsymbol {y}_0| = 1$, such that

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

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

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

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

(4.55)\begin{equation} \mathcal{Y}^2 = \frac{Q_{\alpha \beta} y_0^{\alpha} y_0^{\beta}} { Q_{\alpha \beta} y^{\alpha} y^{\beta}}. \end{equation}

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 $\varPi _{\alpha \beta }$ where heating by the disc turbulence is not balanced by cooling from the $\varPi ^{\sigma }_{(\alpha } B_{\beta ) \sigma }$ term, and would thus experience thermal runaway. Thus, for $\delta > 0$, there exists a $t = t_{cool} > t_0$ such that the components of $Q_{\alpha \beta }$ at $t=t_\textrm {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

(4.56)\begin{equation} 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} \end{equation}

for $t > t_{cool}$. From this we can conclude that $\mathcal {Y} > 1$ for $t > t_{cool}$. By choosing $t$ large enough we can make $\mathcal {Y}$ arbitrarily large. Typically, one expects $t_{cool}$ to be of the order of 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

(4.57)\begin{align} Q_{\alpha_{n+1} \cdots \alpha_{k}}|_{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} \nonumber\\ &= ( Q_{\alpha_{n+1} \cdots \alpha_{k}} Y^{\alpha_{n+1}} \cdots Y^{\alpha_{k}})_{t=t_1} \nonumber\\ &= \mathcal{Y}^{k} Q_{\alpha_{n+1} \cdots \alpha_{k}}|_{t=t_1} y^{\alpha_{n+1}} \cdots y^{\alpha_{k}}. \end{align}

Rearranging this we obtain

(4.58)\begin{align} & Q_{\alpha_{n+1} \cdots \alpha_{k}}|_{t=t_1} y^{\alpha_{n+1}} \cdots y^{\alpha_{k}} \nonumber\\ &\quad =\mathcal{Y}^{{-}k} Q_{\alpha_{n+1} \cdots \alpha_{k}}|_{t=t_0} y_0^{\alpha_{n+1}} \cdots y_0^{\alpha_{k}} \le \mathcal{Y}^{{-}k} \sum_{\alpha_{n+1} \cdots \alpha_{k}} |Q_{\alpha_{n+1} \cdots \alpha_{k}}|_{t=t_0}|. \end{align}

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}}|_{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, $\varPi _{\alpha _1 \cdots \alpha _k}$, decays when its evolution can be well approximated by (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 moments simultaneously that can be sustained and will not damp towards the fluid ordering scheme. For now, we 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 $\varPi _{\alpha \beta \gamma } = 0$ and only consider the first three moments of the Fokker–Planck equation. This yields a continuity, momentum and constitutive relation for a 6-D dust fluid. This dust fluid has a high degree of symmetry as physical properties must be independent of the gas displacement $\{\boldsymbol {x}_g\}$. (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.) One can, therefore, integrate out these redundant degrees of freedom.

In integrating out the dummy gas degrees of freedom, we replace the connections $\boldsymbol {\nabla }_{\alpha }$ with the connections $\bar {\boldsymbol {\nabla }}_{\alpha }$ that 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 {\varPi }_{\alpha \beta }$ and we can always choose the size of the dummy gas dimensions such that $U^{\alpha } = \overline {U^{\alpha }}$. With these choices the dust fluid equations are

(4.59)$$\begin{gather} \bar{D} \rho_d ={-} \rho_d \bar{\boldsymbol{\nabla}}_{\alpha} U^{\alpha}, \end{gather}$$

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

(4.61)$$\begin{gather}\bar{\mathcal{D}}_2 T_{\alpha \beta} ={-}2 (T^{\gamma}_{(\alpha} C_{\beta) \gamma} - \rho_d D_{(\alpha \beta)}), \end{gather}$$

where

(4.62)\begin{equation} \bar{D} = \partial_t + U^{\alpha} \bar{\boldsymbol{\nabla}}_{\alpha} = \begin{pmatrix} \partial_t + u^i_d \boldsymbol{\nabla}_i & \boldsymbol{0} \\ \boldsymbol{0} & \partial_t + u^i_d \boldsymbol{\nabla}_i \end{pmatrix}, \end{equation}

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

(4.63)\begin{equation} \bar{\mathcal{D}}_2 T_{\alpha \beta} = \bar{D} T_{\alpha \beta} + 2 T_{\gamma_d (\alpha} \bar{\boldsymbol{\nabla}}^{\gamma_d} U_{\beta)} + T_{\alpha \beta} \bar{\boldsymbol{\nabla}}_{\gamma_d} U^{\gamma_d}. \end{equation}

Finally, to highlight the effects of the torsion, we consider it's contribution to the dust fluid vorticity,

(4.64)\begin{equation} \omega^{\gamma} \epsilon_{\gamma \alpha \beta} = 2 \bar{\partial}_{[\alpha} U_{\beta]} + S^{\gamma_g}_{\alpha \beta} U_{\gamma_g}, \end{equation}

meaning that

(4.65)\begin{equation} \omega^{\gamma} \epsilon_{\gamma \alpha_g \beta_d} ={-} \bar{\partial}_{\beta_d} U_{\alpha_g} - \mathcal{T}_{\beta_d \alpha_g^{*}}^{\gamma_g^*} U_{\gamma_g} . \end{equation}

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 }$ given by (3.23) and (3.24), then (4.59)–(4.61) are explicitly

(4.66)$$\begin{gather} (\partial_t + U^{\alpha_d} \partial_{\alpha_d} ) \rho_d ={-} \rho_d \partial_{\alpha_d} U^{\alpha_d}, \end{gather}$$

(4.67)$$\begin{gather}\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}), \end{gather}$$

(4.68)\begin{gather}\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^{\beta_d}_{g}) \partial_{\beta_d} U_{\alpha_d}^{g} \nonumber\\ &\quad - \partial_{\beta_d} T_{\alpha_g}^{\beta_d} - \frac{\rho_d}{t_c} (U_{\alpha_g} - U^g_{\alpha_g}) ,\end{aligned} \end{gather}

(4.69)\begin{gather}\begin{aligned} & (\partial_t + U^{\alpha_d} \partial_{\alpha_d} ) T_{\alpha_d \beta_d} + T^{\gamma_d}_{\alpha_d} \partial_{\gamma_d} U_{\beta_d} + T^{\gamma_d}_{\beta_d} \partial_{\gamma_d} U_{\alpha_d} + T_{\alpha_d \beta_d} \partial_{\gamma_d} U^{\gamma_d} \nonumber\\ &\quad ={-}\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} \end{gather}

(4.70)\begin{gather}\begin{aligned} & (\partial_t + U^{\alpha_d} \partial_{\alpha_d} ) T_{\alpha_d \beta_g} + T^{\gamma_d}_{\alpha_d} \partial_{\gamma_d} U_{\beta_g} + T^{\gamma_d}_{\beta_g} \partial_{\gamma_d} U_{\alpha_d} + T_{\alpha_d \beta_g} \partial_{\gamma_d} U^{\gamma_d} \nonumber\\ &\quad ={-} \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} \end{gather}

(4.71)\begin{gather}\begin{aligned} & (\partial_t + U^{\alpha_d} \partial_{\alpha_d} ) T_{\alpha_g \beta_g} + T^{\gamma_d}_{\alpha_g} \partial_{\gamma_d} U_{\beta_g} + T^{\gamma_d}_{\beta_g} \partial_{\gamma_d} U_{\alpha_g} + T_{\alpha_g \beta_g} \partial_{\gamma_d} U^{\gamma_d} \nonumber\\ &\quad ={-} \frac{2}{t_c} (T_{\alpha_g \beta_g} - \alpha c_s^2 \rho_d \delta_{\alpha_g \beta_g}).\end{aligned} \end{gather}

In the 3-D picture we have $u_d^{i} = U^{i}$, $u_s^{i} = U^{i^{*}}$ are the (3-D) dust velocity and fluid seen, respectively, and $p_{i j} = T_{i j}$, $\tau _{i j} = T_{i j^{*}}$, $\sigma _{i j} = 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

(4.72)$$\begin{gather} (\partial_t + u_d^{i} \partial_{i} ) \rho_d ={-} \rho_d \partial_{i} u_d^{i} , \end{gather}$$

(4.73)$$\begin{gather}\rho_d (\partial_t + u_d^{j} \partial_{j} ) u^d_i ={-}\rho_d \partial_{i} \phi - \partial_{j} p^{j}_{i} - \frac{\rho_d}{t_s} (u^d_i - u^s_i), \end{gather}$$

(4.74)$$\begin{gather}\rho_d (\partial_t + u_d^{j} \partial_{j} ) u^s_i ={-}\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^{j}_{i} - \frac{\rho_d}{t_c} (u^s_i - u^g_{i}), \end{gather}$$

(4.75)$$\begin{gather}(\partial_t +u_d^{k} \partial_{k} ) p_{i j} + p^{k}_{i} \partial_{k} u^d_{j} + p^{k}_{j} \partial_{k} u^d_{i} + p_{i j} \partial_{k} u_d^{k} ={-}\frac{2}{t_s} p_{i j} + \frac{1}{t_s} \tau_{i j} + \frac{1}{t_s} \tau_{j i}, \end{gather}$$

(4.76)$$\begin{gather}(\partial_t + u_d^{k} \partial_{k} ) \tau_{i j} + p^{k}_{i} \partial_{k} u^s_{j} + \tau^{k}_{j} \partial_{k} u^d_{i} + \tau_{i j} \partial_{k} u_d^{k} ={-} \left(\frac{1}{t_c} + \frac{1}{t_s}\right) \tau_{i j} + \frac{1}{t_s} \sigma_{i j}, \end{gather}$$

(4.77)$$\begin{gather}(\partial_t + u_d^{k} \partial_{k} ) \sigma_{i j} + \tau^{k}_{i} \partial_{k} u^s_{j} + \tau^{k}_{j} \partial_{k} u^s_{i} + \sigma_{i j} \partial_{k} u_d^{k} ={-} \frac{2}{t_c} (\sigma_{i j} - \alpha c_s^2 \rho_d \delta_{i j}). \end{gather}$$

5. Properties of the dust fluid model

We 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 experiences 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 3-D picture), in such a set-up, evolves according to

(5.1)\begin{align} \partial_t (\rho_d/\rho_g) &={-} \boldsymbol{\nabla}_i (\rho_d u_d^i/\rho_g) \end{align}

(5.2)\begin{align} &={-} \boldsymbol{\nabla}_i (\rho_d u_s^i/\rho_g). \end{align}

Thus, in order that the gas turbulence disperses the dust, we require the velocity of the fluid seen $u_s^i \ne 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 that 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 time scales 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 § 7 we explore this effect further by considering the dust stress tensor in a rotating shear flow.

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

5.3. Viscoelasticity

The dust fluid exhibits viscoelastic behaviour (see Appendix B.2) that is easiest to see when $t_s \sim t_c = O({De})$, where ${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 time scale $t_f$. When ${De} \gg 1$, the dust stress tensor evolves according to

(5.3)\begin{equation} \bar{\mathcal{D}}_2 T_{\alpha \beta} = \bar{\mathcal{D}} T_{\alpha \beta} - 2 T^{\gamma}_{(\alpha} \varepsilon_{\beta) \gamma \sigma} \omega^{\sigma} = 0 . \end{equation}

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 ${De} \ll 1$, the stress tensor is approximately

(5.4)\begin{equation} T_{\alpha \beta} = p_d \left(1 + \frac{t_s}{t_c}\varTheta^g_{\alpha \beta}\right) g_{\alpha \beta} + \frac{1}{2} p_{x} (g_{\alpha \beta^*} + g_{\alpha^* \beta}) - 2 \mu_{\alpha \beta}^{\mu \nu} \bar{\boldsymbol{\nabla}}_{\mu} U_{\nu} + O({De}^2). \end{equation}

At leading order this consists of an isotropic, isothermal, effective, dust pressure with sound speed $\sqrt {{\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 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 viscous terms, for small dust grains, are typically negligible and the dust primarily behaves like an inviscid isothermal gas with a lower temperature than the gas phase. The difference between $U^g_{\alpha _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 (5.4) arises due to the fact that in the $\mathrm {De} \ll 1$ limit $\boldsymbol {v} - \boldsymbol {v}^g$ and $\boldsymbol {v}^g - \boldsymbol {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 $\mathrm {St} \sim 1$, only $\boldsymbol {v}^g - \boldsymbol {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 the order of the stopping time. In the large-Deborah-number limit the fast terms in (3.8) are negligible, resulting in a fully non-local behaviour for $T_{\alpha \beta }$ (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 sees the turbulence as a continuum or is sensitive to individual eddies is determined by a form of ‘eddy-Knudsen number’:

(5.5)\begin{equation} {Kn}_{e} = \frac{\lambda}{L} = \frac{t_c \Delta U^{*}}{L}. \end{equation}

Here $\lambda = t_c \Delta U^{*}$ is the mean free path of a dust grain in the turbulent flow representing the length scale 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 length scale of variations in the fluid flow.

When ${Kn}_{e} \ll 1$, a dust grain interacts with many turbulent eddies over the length scale on which the dust fluid varies, meaning the dust experiences the turbulence as a continuum of stochastic perturbations. When ${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 simplifies to (see Appendix B.3)

(5.6)\begin{align} & 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} \boldsymbol{\nabla}^{k} u^{g}_{\beta_g)} + 2 \left(\frac{T_{\alpha_g \beta_g}}{\rho_d} - \alpha c_s^2 g_{\alpha_g \beta_g}\right) \nonumber\\ &\quad ={-} {Kn}_{e} \left[\frac{\Delta U^{\gamma}}{\Delta U^{*}} L \bar{\boldsymbol{\nabla}}_{\gamma} \left(\frac{T_{\alpha_g \beta_g}}{\rho_d}\right) + 2 \frac{T_{\gamma (\alpha_g}}{\rho_d} L \bar{\boldsymbol{\nabla}}^{\gamma} \frac{\Delta U_{\beta_g)}}{\Delta U^{*}}\right]. \end{align}

When ${Kn}_{e} \rightarrow 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 ${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 rearrange the equations into hyperbolic form, which is useful for some types of numerical solver and for understanding the wave modes in the system. We wish to find a state vector $\boldsymbol {W}$, matrices $\boldsymbol{\mathsf{A}}_i$ and source vector such that the dust fluid equations take the form

(6.1)\begin{equation} \frac{\partial \boldsymbol{W}}{\partial t} + \boldsymbol{\mathsf{A}}^{\alpha} \bar{\boldsymbol{\nabla}}_{\alpha} \boldsymbol{W} = \boldsymbol{S}. \end{equation}

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

(6.2)$$\begin{gather} \dot{\rho}_d + U^{\alpha} \bar{\boldsymbol{\nabla}}_{\alpha} \rho_d + \rho_d \bar{\boldsymbol{\nabla}}_{\alpha} U^{\alpha} = 0, \end{gather}$$

(6.3)$$\begin{gather}\dot{T}_{\alpha \beta} + U^{\gamma} \bar{\boldsymbol{\nabla}}_{\gamma} T_{\alpha \beta} + 2 T_{\gamma (\alpha} \bar{\boldsymbol{\nabla}}^{\gamma} U_{\beta)} + T_{\alpha \beta} \bar{\boldsymbol{\nabla}}_{\gamma} U^{\gamma} ={-}2 (T^{\gamma}_{(\alpha} C_{\beta) \gamma} - \rho_d D_{\alpha \beta}), \end{gather}$$

(6.4)$$\begin{gather}\dot{U}_{\alpha} + U^{\beta} \bar{\boldsymbol{\nabla}}_{\beta} U_{\alpha} + \frac{1}{\rho_d} \bar{\boldsymbol{\nabla}}^{\beta} T_{\beta \alpha} = F_{\alpha} - C_{\alpha \beta} (U^{\beta} - U_g^{\beta}), \end{gather}$$

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

(6.5)\begin{equation} \boldsymbol{W} = \begin{pmatrix} \rho_d \\ \boldsymbol{T}\\ \boldsymbol{U} \end{pmatrix}. \end{equation}

The source vector is

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

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

(6.7) \begin{equation} \boldsymbol{\mathsf{A}}^{\alpha} = \begin{pmatrix} U^{\alpha} & \boldsymbol{0} & \rho_d \hat{\boldsymbol{e}}^{\alpha} \\ 0 & U^{\alpha} \boldsymbol{\mathsf{I}} & \boldsymbol{M}^{\alpha} \\ 0 & \dfrac{\boldsymbol{\mathsf{I}}}{\rho_d} \hat{\boldsymbol{e}}^{\alpha} & U^{\alpha} \boldsymbol{\mathsf{I}} \end{pmatrix}, \end{equation}

where $\boldsymbol{\mathsf{I}}$ denote the identity matrix and

(6.8)\begin{equation} (\boldsymbol{M}^{\alpha} )_{\sigma \beta \gamma} = T_{\sigma \beta} \delta_{\gamma}^{\alpha} + 2 T^{\alpha}_{(\sigma} \delta_{\beta)}^{\gamma}, \end{equation}

such that

(6.9)\begin{equation} (\boldsymbol{M}^{\alpha} \bar{\boldsymbol{\nabla}}_{\alpha} \boldsymbol{U})_{\sigma \beta} = (\boldsymbol{M}^{\alpha} )_{\sigma \beta \gamma} \bar{\boldsymbol{\nabla}}_{\alpha} U^{\gamma}. \end{equation}

For the system to be hyperbolic, we must show that all the eigenvalues of $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 {\boldsymbol {n}} = \hat {\boldsymbol {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 $\boldsymbol {n}$ has non-zero components in the dummy gas directions as we require physical quantities to be independent of $\boldsymbol {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,\ldots,6$). We similarly separate out the stress tensor into compression along the $x$ direction, $T_{1 1}$, shear components in the $x$ direction, $T_{1 \alpha }$, and the components of the stress in other directions, $T_{\alpha \beta }$. We split the momentum and constitutive relation in a similar manor, which results in the following state vector:

(6.10)\begin{equation} \boldsymbol{W} = \begin{pmatrix} \rho \\ T_{1 1} \\ T_{1 \alpha} \\ T_{\alpha \beta} \\ U_{1} \\ U_{\alpha} \end{pmatrix}. \end{equation}

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

(6.11) \begin{equation} \boldsymbol{A}^{1} - v \boldsymbol{\mathsf{I}} = \begin{pmatrix} u^{x} - v & 0 & \boldsymbol{0} & \boldsymbol{0} & \rho_d & \boldsymbol{0} \\ 0 & u^{x} - v & \boldsymbol{0} & \boldsymbol{0} & 3 T_{1 1} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & (u^{x} - v) \boldsymbol{\mathsf{I}} & \boldsymbol{0} & 2 T_{1 \alpha} & T_{1 1} \boldsymbol{\mathsf{I}} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & (u^{x} - v) \boldsymbol{\mathsf{I}} & T_{\alpha \beta} & 2 T_{1 (\alpha} \hat{\boldsymbol{e}}^{\boldsymbol{T}}_{\beta)} \\ 0 & \dfrac{1}{\rho_d} & \boldsymbol{0} & \boldsymbol{0} & (u^{x} - v) & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \dfrac{1}{\rho_d} \boldsymbol{\mathsf{I}} & \boldsymbol{0} & \boldsymbol{0} & (u^{x} - v) \boldsymbol{\mathsf{I}} \end{pmatrix}, \end{equation}

which has a characteristic equation

(6.12)\begin{equation} (v - u^{x})^{16} \left((v - u^{x})^2 - 3 \frac{T_{1 1}}{\rho_d} \right) \left( (v - u^{x})^2 - \frac{T_{1 1}}{\rho_d} \right)^5 = 0 . \end{equation}

This results in 16 non-propagating (in the fluid frame) wavemodes, with wavespeed $v = u^{x}$. These consist of the entropy wave with eigenvector $\left (\begin{smallmatrix} 1 \\ 0_{27} \end{smallmatrix}\right )$ and 15 ‘stress’ waves with eigenvectors $\left (\begin{smallmatrix} 0_7 \\ \hat {\boldsymbol {e}}_{\alpha \beta } \\ 0_{6} \end{smallmatrix}\right )$, where we have introduced the notation $0_n = \begin{smallmatrix} 0 \\ \vdots \\ 0\end{smallmatrix}$, 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 {{3 T_{1 1}}/{\rho _d}}$. The P waves are analogous to sound waves, but with an anisotropic sound speed, with seismic wavespeed anistropy being a well-known phenomena in geophysics (Thomsen 1986). These wavemodes have eigenvectors

(6.13)\begin{equation} \begin{pmatrix} \pm \rho_d \\ \pm 3 T_{1 1} \\ \pm 3 T_{1 \alpha} \\ \pm \left(T_{\alpha \beta} + \dfrac{2}{T_{1 1}} T_{1 (\alpha} T_{\beta) 1}\right) \\ - \sqrt{\dfrac{3 T_{1 1}}{\rho_d}} \\ -\sqrt{\dfrac{3}{\rho_d T_{1 1}}} T_{1 \alpha} \end{pmatrix}. \end{equation}

Finally, there are 10 propagating waves that can be identified as S waves, with wavespeed $v = u^{x} \pm \sqrt {{T_{1 1}}/{\rho _d}}$. As is typical for elastic media, the S waves have slower wavespeeds than the P waves. These wavemodes have eigenvectors

(6.14)\begin{equation} \begin{pmatrix} 0 \\ 0 \\ \pm T_{1 1} \hat{\boldsymbol{e}}_{\gamma} \\ \pm 2 T_{1 (\alpha} \delta^{\gamma}_{\beta)} \\ 0 \\ -\sqrt{\dfrac{T_{1 1}}{\rho_d}} \hat{\boldsymbol{e}}_{\gamma} \end{pmatrix} . \end{equation}

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

(6.15)\begin{equation} \boldsymbol{S} = \begin{pmatrix} 0 \\ - 2 T^{1}_{1} C_{1 1} - 2 T^{\gamma}_{1} C_{1 \gamma} \\ - 2 T^{1}_{(1} C_{\alpha) 1} - 2 T^{\gamma}_{(1} C_{\alpha) \gamma} \\ - 2 T^{1}_{(\alpha} C_{\beta) 1} - 2 T^{\gamma}_{(\alpha} C_{\beta) \gamma} + 2 \rho_d D_{\alpha \beta} \\ F_{1} - C_{1 1} (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}. \end{equation}

Thus, we see that there are no source terms for the entropy wave. Turbulent diffusion ($\boldsymbol {D}$) and the drag-dependent coupling between pressure tensor components are sources/sinks of the stress waves. Finally, the force per unit mass $\boldsymbol {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 structures. In such a situation the disc possesses 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 set-up in experimental fluid dynamics (e.g. Taylor–Couette flows). To study this problem, we adopt (6-D) cylindrical polar coordinates $(R,\phi,z,R_g,\phi _g,z_g)$ with metric tensor components

(7.1a,b)\begin{equation} g_{R R} = g_{R_g R_g} = g_{z z} = g_{z_g z_g} = 1,\quad g_{\phi \phi} = g_{\phi_g \phi_g} = R^2 . \end{equation}

and connection coefficients

(7.2)$$\begin{gather} \varGamma_{\phi \phi}^{R} = \varGamma_{\phi \phi_g}^{R_g} ={-}R, \end{gather}$$

(7.3)$$\begin{gather}\varGamma_{\phi R}^{\phi} = \varGamma_{R \phi}^{\phi} = \varGamma_{\phi R_g}^{\phi_g} = \varGamma_{R \phi_g}^{\phi_g} = 1/R, \end{gather}$$

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 $\varOmega = \varOmega (R)$. This leads to the 6-D mean velocity of the dust fluid of

(7.4)\begin{equation} U^{\gamma} = \varOmega (R) (\hat{e}^{\alpha}_{\phi} + \hat{e}^{\alpha}_{\phi_g}). \end{equation}

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 $\bar {\mathcal {D}}_2$ (when acting on $T_{\alpha \beta }$) is

(7.5)\begin{align} \bar{\mathcal{D}}_2 T_{\alpha \beta} = \bar{D} T_{\alpha \beta} + 2 T^{R}_{(\alpha} (\hat{e}_{\beta)}^{\phi} + \hat{e}_{\beta)}^{\phi_g}) \partial_{R} (R^2 \varOmega) -2 \varOmega \varGamma_{\phi (\alpha}^{\gamma} T_{\beta) \gamma} - 2 R^2 \varOmega (\varGamma_{\gamma (\alpha}^{\phi} + \varGamma_{\gamma (\alpha}^{\phi_g}) T_{\beta)}^{\gamma}. \end{align}

We are interested in the steady-state solution to the stress tensor with $\bar {D} T_{\alpha \beta } = 0$. We thus have the following for the constitutive relation in the steady rotating shear flow:

(7.6)\begin{align} & 4 R (\varOmega - A) T^{R}_{(\alpha} (\hat{e}_{\beta)}^{\phi} + \hat{e}_{\beta)}^{\phi_g}) -2 \varOmega \varGamma_{\phi (\alpha}^{\gamma} T_{\beta) \gamma} - 2 R^2 \varOmega (\varGamma_{\gamma (\alpha}^{\phi} + \varGamma_{\gamma (\alpha}^{\phi_g}) T_{\beta)}^{\gamma} \nonumber\\ &\quad ={-}2 (T^{\gamma}_{(\alpha} C_{\beta) \gamma} - \rho_d D_{\alpha \beta}). \end{align}

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

Explicitly, the ‘dust–dust’ components of (7.6), which can be thought of as the equations governing the behaviour of the 3-D dust stress, are

(7.7)$$\begin{gather} -\frac{4 \varOmega}{R} T_{R \phi} ={-}\frac{2}{t_s} (T_{R R} - T_{R R_g}) , \end{gather}$$

(7.8)$$\begin{gather}- \frac{2 \varOmega}{R} T_{\phi \phi} + 2 R (\varOmega - A) T_{R R} ={-}\frac{2}{t_s} T_{R \phi} + \frac{1}{t_s} (T_{R \phi_g} + T_{R_g \phi}), \end{gather}$$

(7.9)$$\begin{gather}4 R (\varOmega - A) T_{R \phi} ={-} \frac{2}{t_s} (T_{\phi \phi} - T_{\phi \phi_g}) . \end{gather}$$

The ‘dust–gas’ components equation (7.6), which governs the behaviour of the ‘cross-stress’ for the cross-correlation between the dust and gas velocities are

(7.10)$$\begin{gather} -\frac{2 \varOmega}{R} T_{R_g \phi} - \frac{\varOmega}{R} (T_{R \phi} + T_{R \phi_g}) ={-} \left(\frac{1}{t_s} + \frac{1}{t_c} \right) T_{R R_g} + \frac{1}{t_s} T_{R_g R_g}, \end{gather}$$

(7.11)$$\begin{gather}-\frac{2 \varOmega}{R} T_{\phi \phi_g} - \varOmega R (T_{R R} - T_{R R_g}) + 2 R (\varOmega - A) T_{R R} ={-} \left(\frac{1}{t_s} + \frac{1}{t_c} \right) T_{R \phi_g} + \frac{1}{t_s} T_{R_g \phi_g}, \end{gather}$$

(7.12)$$\begin{gather}-\frac{\varOmega}{R} (T_{\phi \phi} + T_{\phi \phi_g}) + 2 R (\varOmega - A) T_{R R_g} ={-} \left(\frac{1}{t_s} + \frac{1}{t_c} \right) T_{R_g \phi} + \frac{1}{t_s} T_{R_g \phi_g}, \end{gather}$$

(7.13)$$\begin{gather}2 R (\varOmega - A) (T_{R \phi} + T_{R \phi_g}) + R \varOmega (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} . \end{gather}$$

Finally, the ‘gas–gas’ components of (7.6), which describe the behaviour the (3-D) gas Reynolds stress along the dust trajectory, are

(7.14)$$\begin{gather} -\frac{2 \varOmega}{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), \end{gather}$$

(7.15)$$\begin{gather}-\frac{\varOmega}{R} (T_{\phi \phi_g} + T_{\phi_g \phi_g}) - \varOmega R (T_{R R_g} - T_{R_g R_g}) + 2 R (\varOmega - A) T_{R R_g} ={-}\frac{2}{t_c} T_{R_g \phi_g} , \end{gather}$$

(7.16)$$\begin{gather}2 R \varOmega T_{R_g \phi_g} + 2 R (\varOmega - 2 A) T_{R \phi_g} ={-}\frac{2}{t_c} (T_{\phi_g \phi_g} - \alpha c_s^2 \rho_d R^2) . \end{gather}$$

It is straightforward, if rather laborious, to invert the above equations. However, the resulting expressions are somewhat cumbersome and not particularly informative. We 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 available at https://doi.org/10.1017/jfm.2024.1088). We can then numerically explore the behaviour of the dust pressure tensor.

Figure 1 shows how the horizontal stress tensor components change with $A/\varOmega$ 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 possesses a steady solution. Together these show the general behaviour of the dust stress tensor. This shows that the dust stress tends to get weaker at larger Stokes numbers and tends to isotropy in the absence of shear $A/\varOmega \rightarrow 0$. There are singularities/zeros of the pressure tensor at large $A/\varOmega$ and ${St}$ associated with the breakdown of the fluid dust description. For $\tau _c \lesssim 1$, these asymptote to the Rayleigh stability criterion for large ${St}$, for small Stokes numbers, dust drag/cooling helps to regularise the behaviour of the pressure tensor allowing for steady solutions at higher $A/\varOmega$. 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/\varOmega \sim 2$ for small Stokes numbers. This occurs because of a breakdown of the turbulence model for large $\tau _c$ and $A/\varOmega$ (see Appendices A.2 and A.3).

Figure 1. Horizontal stress tensor components in a rotating shear flow, as a function of ${St}$ and $A/\varOmega$, with different dimensionless correlation times.

Figure 2. Locations in the parameter space of zeros of the stress tensor components, indicating a lack of steady solutions. Black solid line: $\tau _c=10^{-3}$, black dashed line: $\tau _c=10^{-2}$, black dotted line: $\tau _c=10^{-2}$, blue solid line: $\tau _c=1$, blue dashed line: $\tau _c=10^{3}$. The region of the parameter space above these lines contains no (physical) steady solutions.

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} \varOmega$. (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.

Figure 3. Horizontal stress tensor components for a fluid with a dimensionless correlation time, $\tau _c = 1$. (a) A Rayleigh stable Keplerian shear flow where $\varOmega \propto R^{-3/2}$ ($A = \frac {3}{4} \varOmega$). (b) Rayleigh unstable shear flow with $A = 1.1 \varOmega$. 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.

The right-hand plot shows a Rayleigh unstable shear flow with $A = 1.1 \varOmega$. (The Reynolds stress of the gas is shown in the right-hand plot of figure 7.) Again at low Stokes numbers the dust stress is set by the gas velocity correlations. The stress tensor components diverge as the Stokes number increases, and one approaches breakdown 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. The S-wave 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.

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.

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 one-dimensional 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 (A10)–(A12) perturbatively, with the leading-order terms being that due to gravity and circular Keplerian motion. We then consider the first-order correction to the gas velocity due to the gas pressure and turbulence. This has the 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 time scale 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 of the HLL (Harten, Lax & van Leer 1983) and Roe (Roe 1981) approximate-Reimann solvers, and constant reconstruction. We use an operator-split Van-Leer integrator (Eleuterio 1999; Van Leer 2006) 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 ensuring 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,

(7.17)\begin{equation} \rho_d = 0.1 + 0.9 [\tanh(2 R + 1) + 1 ] . \end{equation}

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 § 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 starts 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 § 7.1. The 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.

Figure 5. Density plots for the rotating shear flow, (a) constant density flow, (b) 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 (a) and every 200 inner orbits (b).

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

The simulations exhibit several of the expected features 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 absence 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-type (Toro, Spruce & Speares 1994) 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\textrm{--} 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 }$. One obvious improvement would be to implement a FARGO-like(Masset 2000) advection scheme and subtract off the orbital motion, to reduce the error on $T_{\phi \phi }$. It is 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}$,

(8.1)\begin{equation} f_i^{g} ={-} \boldsymbol{\nabla}_i \phi - \frac{1}{\rho_g} \boldsymbol{\nabla}_i p_g - \frac{1}{t_s} f_d (u_i^{g} - u_i^{d}), \end{equation}

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 (8.1) by replacing the mean flow back reaction term with $- ({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, $\boldsymbol {v}^{g}$ is the velocity of a fluid element (seen), while $\boldsymbol {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 to make $t_c$ and $\alpha$ dependent 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 that the total energy in the gas plus dust system is conserved (particularly to ensure that the turbulence does not 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 that 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, Goldreich & Tremaine 1985; Araki & Tremaine 1986; Latter & Ogilvie 2008; Larue, Latter & Rein 2023). The inclusion of dust collisions in a SDE model for particle laden flows, and associate moment closure, has been studied in Innocenti, Fox & Chibbaro (2019, 2021), Capecelatro, Desjardins & Fox (2016a) and Capecelatro et al. (2016b). 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, which are an additional source of stochasticity on top of the gas turbulence that 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 (e.g. Sawford 1991; Pope 2002) include two time scales (correlation time and Kolmogorov time) and may be 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) gases are very large, this effect is likely only important in a limited region of parameter space – but may be 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 SDEs 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-D anisotropic Maxwell fluid with a dynamically important rheological stress tensor. The 6-D 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-D dust stress tensor by the fluid flow.

In summary our conclusions are as follows.

1. (i) 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.

2. (ii) Collisionless dust has a non-zero anisotropic rheological stress that 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.

3. (iii) The dust can support seismic (P and S) waves

4. (iv) 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 extensions to our model, some of which we intend to pursue in future work.