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.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 Reference Pope1985, Reference Pope2000; Minier & Peirano Reference Minier and Peirano2001). 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 Reference Pope1985, Reference Pope2000; Minier & Peirano Reference Minier and Peirano2001):

(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 (Reference Youdin and Lithwick2007). 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 Reference Fox2014; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2016b).) 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 (Reference Ogilvie2003) and Lynch & Ogilvie (Reference Lynch and Ogilvie2021), $\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 Reference Withers1985) 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 (Reference Minier and Henry2023), 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.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 Reference Minier2001; Minier et al. Reference Minier, Peirano and Chibbaro2004, Reference Minier, Chibbaro and Pope2014) 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 Reference Gavrilyuk and Gouin2012)) – 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 Reference Minier2016). 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.

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 Reference Harten, Lax and van Leer1983) and Roe (Roe Reference Roe1981) approximate-Reimann solvers, and constant reconstruction. We use an operator-split Van-Leer integrator (Eleuterio Reference Eleuterio1999; Van Leer Reference Van Leer2006) 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 Reference Toro, Spruce and Speares1994) 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 Reference Masset2000) 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 Reference Benítez-Llambay and Masset2016). 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.