Dynamics of Dusty Vortices I: Extensions and limitations of the terminal velocity approximation

TL;DR

This paper investigates the terminal velocity approximation in dusty vortices within protoplanetary discs, highlighting its limitations near shocks and providing a numerical solver implementation in FARGO3D.

Contribution

It introduces a numerical solver for the terminal velocity equations in FARGO3D, demonstrating its application to viscous dusty protoplanetary discs and analyzing its limitations.

Findings

The model works well for viscous discs using baricenter velocity.

The approximation breaks down near shocks, conflicting with the two-fluid model.

The implemented solver is fast, stable, and second-order accurate for smooth problems.

Abstract

Motivated by the stability of dust laden vortices, in this paper we study the terminal velocity approximation equations for a gas coupled to a pressureless dust fluid and present a numerical solver for the equations embedded in the FARGO3D hydrodynamics code. We show that for protoplanetary discs it is possible to use the baricenter velocity in the viscous stress tensor, making it trivial to simulate viscous dusty protoplanetary discs with this model. We also show that the terminal velocity model breaks down around shocks, becoming incompatible with the two fluid model it is derived from. Finally we produce a set of test cases for numerical schemes and demonstrate the performance of our code on these tests. Our implementation embedded in FARGO3D using an unconditionally stable explicit integrator is fast, and exhibits the desired second order spatial convergence for smooth problems.

Equations

$$\begin{array}{rl}& {\partial }_t{\rho }_g+\nabla \cdot \left({\rho }_g{\mathbf{u}}_g\right)=0\end{array}$$

$$\begin{array}{rl}& {\partial }_t\left({\rho }_g{\mathbf{u}}_g\right)+\nabla \cdot \left({\rho }_g{\mathbf{u}}_g{\mathbf{u}}_g\right)+c_s^2\nabla {\rho }_g=\nabla \cdot {\mathsf{T}}_g-{\rho }_g\nabla \Phi -\\ & \frac{{\rho }_d{\rho }_g}{t_s\rho }\left({\mathbf{u}}_g-{\mathbf{u}}_d\right)\end{array}$$

$$\begin{array}{rl}& {\partial }_t{\rho }_d+\nabla \cdot \left({\rho }_d{\mathbf{u}}_d\right)=0\end{array}$$

$$\begin{array}{rl}& {\partial }_t\left({\rho }_d{\mathbf{u}}_d\right)+\nabla \cdot \left({\rho }_d{\mathbf{u}}_d{\mathbf{u}}_d\right)=-{\rho }_g\nabla \Phi -\frac{{\rho }_d{\rho }_g}{t_s\rho }\left({\mathbf{u}}_d-{\mathbf{u}}_g\right)\end{array}$$

$$\begin{array}{rl}& {\partial }_t\rho +\nabla \cdot \left(\rho \mathbf{u}\right)=0\end{array}$$

$$\begin{array}{rl}& {\partial }_t\left(\rho \mathbf{u}\right)+\nabla \cdot \left(\rho \mathbf{uu}+\frac{P}{c_s^2}\left(1-\frac{P}{c_s^2\rho }\right)\Delta \mathbf{u}\Delta \mathbf{u}\right)+\nabla P=\\ & \nabla \cdot {\mathsf{T}}_{\mathsf{g}}-\rho \nabla \Phi \end{array}$$

$$\begin{array}{rl}& {\partial }_{t}P+\nabla \cdot \left(P\mathbf{u}-P\left(1-\frac{P}{c_s^2\rho }\right)\Delta \mathbf{u}\right)=0\end{array}$$

$$\begin{array}{rl}& {\partial }_t\Delta \mathbf{u}+\nabla \cdot \left(\mathbf{u}\Delta \mathbf{u}-\left(\frac{1}{2}-\frac{P}{c_s^2\rho }\right)\Delta \mathbf{u}\Delta \mathbf{u}\right)=\\ & -\frac{\Delta \mathbf{u}}{t_s}+c_s^2\frac{\nabla P}{P}-c_s^2\frac{\nabla \cdot {\mathsf{T}}_{\mathsf{g}}}{P}.\end{array}$$

$$\begin{array}{rl}& {\partial }_t\tilde{\rho }+\nabla \cdot \left(\tilde{\rho }\tilde{\mathbf{u}}\right)=0\end{array}$$

$$\begin{array}{rl}& {\partial }_t\left(\tilde{\rho }\tilde{\mathbf{u}}\right)+\nabla \cdot \left(\tilde{\rho }\tilde{\mathbf{u}}\tilde{\mathbf{u}}+\tilde{P}\overline{f_0}\left(1-\frac{\tilde{P}\overline{f_0}}{\tilde{\rho }}\right)\Delta \tilde{\mathbf{u}}\Delta \tilde{\mathbf{u}}\right)+\frac{\overline{f_0}}{\mathrm{M}{\mathrm{a}}^2}\nabla \tilde{P}=\\ & \frac{\overline{f_0}\nabla \cdot {\tilde{\mathsf{T}}}_g}{\text{Re}}-\frac{\tilde{\rho }\nabla \Phi }{\mathrm{F}{\mathrm{r}}^2}\end{array}$$

$$\begin{array}{rl}& {\partial }_t\tilde{P}+\nabla \cdot \left(\tilde{P}\tilde{\mathbf{u}}-\tilde{P}\left(1-\frac{\tilde{P}\overline{f_0}}{\tilde{\rho }}\right)\Delta \tilde{\mathbf{u}}\right)=0\end{array}$$

$$\begin{array}{rl}& {\partial }_t\Delta \tilde{\mathbf{u}}+\nabla \cdot \left(\tilde{\mathbf{u}}\Delta \tilde{\mathbf{u}}-\left(\frac{1}{2}-\frac{\tilde{P}\overline{f_0}}{\tilde{\rho }}\right)\Delta \tilde{\mathbf{u}}\Delta \tilde{\mathbf{u}}\right)=\\ & -\frac{\Delta \tilde{\mathbf{u}}}{\widetilde{t_s}\text{St}}+\frac{\nabla \tilde{P}}{\mathrm{M}{\mathrm{a}}^2\tilde{P}}-\frac{\nabla \cdot {\tilde{\mathsf{T}}}_g}{\tilde{P}\text{Re}}.\end{array}$$

$${\tilde{\mathsf{T}}}_g=\tilde{\nu }\tilde{P}\left[\nabla {\tilde{\mathbf{u}}}_g+{\left(\nabla {\tilde{\mathbf{u}}}_g\right)}^T-\frac{2}{3}\left(\nabla \cdot {\tilde{\mathbf{u}}}_g\right)\mathrm{I}\right],$$

$$\xi (x,t,\mathrm{St})={\xi }^{(0)}(x,t)+\mathrm{St}{\xi }^{(1)}(x,t)+\mathcal{O}(\mathrm{S}{\mathrm{t}}^2),$$

$$\underset{\mathrm{St}\to 0}{\lim }(\Delta \tilde{\mathbf{u}})=0,$$

$$\Delta {\tilde{\mathbf{u}}}^{(1)}=\frac{\widetilde{t_s}}{\tilde{P}}\left(\frac{\nabla \tilde{P}}{\mathrm{M}{\mathrm{a}}^2}-\frac{\nabla \cdot {\tilde{\mathsf{T}}}_g}{\mathrm{Re}}\right)+\mathcal{O}(\mathrm{St})$$

$$\Delta \tilde{\mathbf{u}}=\mathrm{St}\frac{\widetilde{t_s}}{\tilde{P}}\left(\frac{\nabla \tilde{P}}{\mathrm{M}{\mathrm{a}}^2}-\frac{\nabla \cdot {\tilde{\mathsf{T}}}_g}{\mathrm{Re}}\right)+\mathcal{O}({\mathrm{St}}^2).$$

$$\mathrm{R}{\mathrm{e}}^{-1}=\mathcal{O}(\mathrm{St}),$$

$${\tilde{\mathbf{u}}}_g=\tilde{\mathbf{u}}+\mathcal{O}(\mathrm{St});$$

$$\frac{{\tilde{\mathbf{u}}}_g}{\mathrm{Re}}=\frac{\tilde{\mathbf{u}}}{\mathrm{Re}}+\mathcal{O}\left(\frac{\mathrm{St}}{\mathrm{Re}}\right)=\frac{\tilde{\mathbf{u}}}{\mathrm{Re}}+\mathcal{O}\left(\mathrm{S}{\mathrm{t}}^2\right),$$

$$\frac{{\tilde{\mathsf{T}}}_g}{\mathrm{Re}}=\frac{1}{\mathrm{Re}}\tilde{\nu }P\left[\nabla \tilde{\mathbf{u}}+{\left(\nabla \tilde{\mathbf{u}}\right)}^T-\frac{2}{3}\left(\nabla \cdot \tilde{\mathbf{u}}\right)\mathrm{I}\right]+\mathcal{O}(\mathrm{S}{\mathrm{t}}^2)$$

$$=\frac{\tilde{\mathsf{T}}}{\mathrm{Re}}+\mathcal{O}(\mathrm{S}{\mathrm{t}}^2).$$

$${\partial }_t\tilde{\rho }+\nabla \cdot \left(\tilde{\rho }\tilde{\mathbf{u}}\right)=0$$

$${\partial }_t\left(\tilde{\rho }\tilde{\mathbf{u}}\right)+\nabla \cdot \left(\tilde{\rho }\tilde{\mathbf{u}}\tilde{\mathbf{u}}\right)+\frac{\overline{f_0}}{\mathrm{M}{\mathrm{a}}^2}\nabla \tilde{P}=\frac{\overline{f_0}\nabla \cdot \tilde{\mathsf{T}}}{\text{Re}}-\frac{\tilde{\rho }\nabla \Phi }{\mathrm{F}{\mathrm{r}}^2}$$

$${\partial }_t\tilde{P}+\nabla \cdot \left(\tilde{P}\tilde{\mathbf{u}}\right)=\frac{\mathrm{St}}{\mathrm{M}{\mathrm{a}}^2}\nabla \cdot \left[\widetilde{t_s}\left(1-\frac{\tilde{P}}{c_s^2\tilde{\rho }}\right)\nabla \tilde{P}\right].$$

$${\partial }_t\rho +\nabla \cdot \left(\rho \mathbf{u}\right)=0$$

$${\partial }_t\left(\rho \mathbf{u}\right)+\nabla \cdot \left(\rho \mathbf{uu}\right)+\nabla P=\nabla \cdot \mathsf{T}-\rho \nabla \Phi$$

$${\partial }_{t}P+\nabla \cdot \left(P\mathbf{u}\right)=\nabla \cdot \left[t_s\left(1-\frac{P}{c_s^2\rho }\right)\nabla P\right].$$

$$i\omega \delta \left(\rho \right)+ik{\rho }_0\delta \left(\mathbf{u}\right)=0$$

$$i\omega {\rho }_0\delta \left(\mathbf{u}\right)=ik\delta \left(P\right)$$

$$i\omega \delta \left(P\right)+ik{P}_0\delta \left(\mathbf{u}\right)=-c_s^2{t}_s\left(1-\frac{P_0}{c_s^2{\rho }_0}\right)k^2\delta \left(P\right)$$

$$\omega =\frac{i{k}^2{D}_c}{2}\pm |k|\sqrt{\frac{P_0}{{\rho }_0}-\frac{D_c^2}{4}},$$

$$\omega =\frac{i{c}_s^2{t}_s{f}_d{k}^2\pm \sqrt{-c_s^4{t}_s^2{f}_d^2{k}^4+4k^2{c}_s^2(1-f_d)}}{2}.$$

$${\partial }_t(\rho )+{\partial }_x(\rho u)=0$$

$${\partial }_t(\rho u)+{\partial }_x(\rho u^2+P)=0$$

$${\partial }_t(P)+{\partial }_x(Pu)=C(P,\rho )$$

$$C(P,\rho )=c_s^2{\partial }_x\left[t_s\left(1-\frac{P}{c_s^2\rho }\right){\partial }_{x}P\right].$$

$${\partial }_x(\rho u)=0$$

$${\partial }_x(\rho u^2+P)=0$$

$${\partial }_x(Pu)=C(P,\rho )$$

$$B=$$

Full text

††pubyear: 2019††pagerange: Dynamics of Dusty Vortices I: Extensions and limitations of the terminal velocity approximation–A

1 Introduction

Dust is ubiquitous in space. In most astrophysical scenarios dust is coupled to a gaseous phase, such that its evolution is determined by an interplay of gas dynamics, external forces, and self gravity. Gas-dust dynamics plays an important role in AGB star winds (e.g. Dominik et al., 1989; Höfner & Olofsson, 2018), and planetary and substellar atmospheres (e.g. Juncher et al., 2017; Johnas et al., 2008; Taylor et al., 2007) as well as in protoplanetary discs. The behaviour of dust in such a mixture can be understood as an ensemble of dust particles feeling the external fields, as well as a drag force from the surrounding gas. Many astrophysical flows, especially protoplanetary discs tend to have very well but not perfectly coupled dust and gas (e.g. Testi et al., 2014). Effectively modelling the behaviour of these dusty gases is key to understanding many astrophysical phenomena, from planet formation to the formation of complex molecules in the inter stellar medium (e.g. Wakelam et al., 2017). This is especially true when tying the models to infrared (IR) observations from telescopes like the Atacama Large Millimeter Array (ALMA) as observed IR emission from gas clouds and discs is largely thermal emission from millimeter and sub millimeter dust particles in these objects. The brightness of dust continuum emission makes it in general easier to observe than gas. This often makes dust the target for observations of protoplanetary discs, (see e.g Zurlo et al., 2018).

Dust comprises only a small fraction of interstellar matter (~1% by mass), but its effects on the dynamics of flows can be important such as with the damping of vortices (Inaba & Barge, 2006). In other cases the dynamics of the dust phase in a gas flow can be quite different than the underlying gas flow. For example in the case of protoplanetary discs it is much easier to clear gaps in the dust than it is in the gas, leading to much deeper gaps in the dust density than in the gas density (e.g. Paardekooper & Mellema, 2004). For protoplanetary discs, dust is both important in the evolution of the disc, among other things playing a fundamental role in planet formation, and key to observations. The highest resolution data for protoplanetary discs has been produced by ALMA (e.g. Andrews et al., 2016) taking observations in sub-millimeter continuum emission. This fine dust is in the very well, but not perfectly coupled regime. This makes understanding the behaviour of well coupled dust in a gas crucial also in the interpretation of observations for protoplanetary discs.

One area where dust-gas interactions can play an important role is in the evolution of dust-laden vortices. These dusty vortices are understood to be important in the dynamics of protoplanetary discs. In gas only discs, vortices arising from instabilities may be long lived (e.g. Fu et al., 2014a), or even self sustaining (e.g. Lesur & Papaloizou, 2010). Protoplanetary discs however, contain dust, which has been shown to collect in vortices (Barge & Sommeria, 1995) and may therefore aid the process of planet formation (Lyra, W. et al., 2009). These dusty vortices do not behave like their gas only counterparts, with dust damping the vortex and driving it unstable as shown in Fu et al. (2014b), Chang & Oishi (2010) and Railton & Papaloizou (2014).

To fully capture the physics of the gas-dust mixture, a model that solves the motion of the individual dust particles in the gas flow would be required, studying this kind of model though is currently not computationally viable for large systems, and is virtually intractable analytically. Further study of the driving mechanisms behind these instabilities requires high resolution simulations as well as analytical study of the governing equations; both are good reasons to apply simplified models to the problem.

Several approaches have been taken to model dusty astrophysical flows. The cheapest methods use a small number of test particles feeling a drag force due to the fluid around them and an additional force due to any external potential field they may be in, like was done by McNally et al. (2018) or Clarke et al. (2015). It is more challenging to simulate the back reaction of the dust on the gas, many methods attempting to do this also model the dust as a fluid. This kind of approximation is possible when the dust particles’ mean free path through the fluid is short compared to the typical length scale of the flow, allowing for the description of the dust through bulk properties (Garaud et al., 2004). Examples of this kind of model are the Lagrangian dust model described in Osiptsov (2000) or two fluid Eulerian models like Johansen et al. (2004). This kind of problem can also be quite intuitively approached using an SPH (smoothed particle hydrodynamics) approach, where some of the fluid particles are dust, and others are gas, like is done in Price & Laibe (2015).

As a possible further simplification to modelling back reacting dust in a gas, one may consider one fluid terminal velocity approximations like the ones from Laibe & Price (2014) and Lin & Youdin (2017). These approximations require for less equations to be solved, while still capturing the effects of the dust on the gas. With this kind of approximation it becomes possible to more easily run simulations at resolutions where with two fluids it would have required large amounts of computing power, or allow for parameter studies to be done where with two fluids only a single run could have been done. This is especially true for very well coupled dust gas mixtures, as the two fluid equations become stiffer with stronger coupling between the dust and gas, while the opposite is true for the terminal velocity approximation. These one fluid approximations also have the added benefit of being more tractable with analytic methods. As one fluid terminal velocity approximations have become more popular, several numerical codes have been written to solve this kind of problem like in Lin & Youdin (2017), Hutchison et al. (2018) for multiple grain sizes, or our own variation on FARGO3D described in section 5.1.

Motivated by our interest in the stability of dusty vortices, in this paper we try to extend understanding of the terminal velocity approximation, formulating a viscous version of the equations, and study the limitations and breakdown of the approximation around sharp pressure gradients and shocks. We also propose some analytic solutions to possibly be used as test cases for one fluid terminal velocity approximation solvers, similarly to what was recently done by Lehmann & Wardle (2018) for two fluid dusty gas hydrodynamics solvers.

The rest of this paper is structured as follows. In section 2 we re-derive the terminal velocity approximation, showing how it can consistently be extended to viscous discs without any additional source terms. We then in section 3 detail several analytic solutions of the equations, a small amplitude wave solution, a solution for shock profiles, and a purely diffusive solution. In section 4 we compare the one fluid and two fluid shock solutions and show that a first order error arises in the terminal velocity approximation around shocks. In section 5.1 we detail an implementation embedded in FARGO3D and its performance. In section 6 we show the performance of our implementation of FARGO3D on the analytical test problems derived in section 3. In sections 7 and 8, we discuss and summarise the results of the paper.

2 Dusty Gases

A dusty fluid flow can be modelled as two interacting fluids, with a gas phase and a pressureless cold dust phase. The description of dust as a fluid is possible provided that the dust is coupled to the gas, allowing the dust particles to exhibit collective behaviour (Garaud et al., 2004). A gas and cold dust phase can be understood as only coupled by a drag force. The drag force is modelled as an Epstein drag. This is applicable high Knudsen number (the dimensionless relation between characteristic length and mean free path in a system) regime. This means that the dust particles are smaller than the gas molecules’ mean free path. The two fluid model however is complex to deal with, having at least two additional fields and equations than a pure gas in the form of the continuity and momentum equation for the dust, which require the addition of a dust density and dust velocity field. This makes the problem both computationally expensive to solve numerically, and difficult to study analytically. Several approaches have been taken to reduce the complexity of the problem, the one we are mainly interested in is the terminal velocity approximation. The terminal velocity approximation is a first order in Stokes number (a dimensionless number describing the relation between drag and momentum of a body in a flow) approximation on the one fluid formulation of the dusty gas equations and hence works best for small particles. Here we run through the derivation of the one fluid approximation for a pressureless dust fluid coupled to an isothermal gas. We use some arguments from dimensional analysis to demonstrate how the equations can be extended to include viscous terms in the case of protoplanetary discs, at no additional cost.

2.1 From two fluids to one fluid approximation

The two fluid dust equations describe a dusty fluid in a potential field $\Phi$ as two fluids coupled by a drag force parametrised as a stopping time $t_{s}$. One of the fluids is the pressure-less ”dust fluid” with velocity $\mathbf{u}_{d}$ and density $\rho_{d}$, the other is the isothermal gas described by its velocity $\mathbf{u}_{g}$, density $\rho_{g}$, and a sound speed $c_{s}$, the pressure in this fluid is given by $P=c_{s}^{2}\rho_{g}$. For dust coupled to a viscous gas we also add the stress tensor $\mathsf{T}_{g}$ in the gas equations.

[TABLE]

The velocity can then be re-expressed as a bulk velocity $\mathbf{u}$ such that: $\mathbf{u}=(\rho_{g}\mathbf{u}_{g}+\rho_{d}\mathbf{u}_{d})/\rho$ where $\rho=\rho_{d}+\rho_{g}$, the total density, and the velocity difference between gas and dust $\Delta\mathbf{u}=\mathbf{u}_{d}-\mathbf{u}_{g}$. To describe the gas to dust ratio we choose a gas to total mass ratio $f_{0}=\rho_{g}/\rho$, for an isothermal gas $f_{0}=P/c_{s}^{2}\rho$, as the pressure $P=c_{s}^{2}\rho_{g}$. This allows for the reconstruction of the two fluid equations as a set of one fluid equations without need for any approximations.

[TABLE]

The non dimensional formulation of the 1 fluid gas dust equations is constructed by choosing a typical length scale $L$, velocity $\bar{u}$, density $\bar{\rho}$ and gas fraction $\bar{f_{0}}=\bar{P}/(c_{s}^{2}\bar{\rho})$, and decomposing the flow fields, $\xi=\bar{\xi}\tilde{\xi}$ where $\tilde{\xi}$ is the dimensionless field,

[TABLE]

The dimensionless numbers produced are the Froude number $\mathrm{Fr}=\bar{u}/\sqrt{gL}$, the Stokes number $\mathrm{St}=t_{s}\bar{u}/L$, the Mach number $\mathrm{Ma}=\bar{u}/c_{s}$, and the Reynolds number $\mathrm{Re}=\bar{u}L/\bar{\nu}$. We keep a non-dimensional stopping time in the equations for clarity, which is understood to be $\mathcal{O}(1)$. The nondimensional viscous stress tensor $\tilde{\mathsf{T}}_{g}$ is given by,

[TABLE]

where we keep a viscosity $\tilde{\nu}$ which is understood to be $\mathcal{O}(1)$.

2.2 Asymptotic expansions and the $\mathrm{St}\ll 1$ approximation

Using a power series expansion in $\mathrm{St}$ for the flow fields such that a flow variable $\xi(x,t,\mathrm{St})$ becomes

[TABLE]

we thus exploit the assumption that we are in the $\mathrm{St}\ll 1$ regime. The only flow field that depends explicitly on $\mathrm{St}$ is $\Delta\tilde{\mathbf{u}}$, since equation (LABEL:eq:twofluidnondim.4) is the only one with a term of $\mathrm{St}$. As the coupling approaches perfect coupling we expect the velocity difference between the dust and gas phases to go to zero so,

[TABLE]

which implies that $\Delta\tilde{\mathbf{u}}^{(0)}=0$. Substituting the expansion for $\Delta\tilde{\mathbf{u}}$ into equation (LABEL:eq:twofluidnondim.4), we obtain,

[TABLE]

this expression can be substituted back into the expansion of $\Delta\tilde{\mathbf{u}}$ to give the expression

[TABLE]

The form of the viscous terminal velocity approximation depends on the relative magnitudes of Re, St and Ma. For it to be possible to extract a simple terminal velocity approximation from these equations we have to be in a regime where $\mathrm{Re}\gg\mathrm{Ma}^{2}$ which can be formalised by setting $\mathrm{Re}=\mathcal{O}\left(\mathrm{St}^{-1}\right)$. The viscous term in (17) then gets pushed to higher order and taking this first order approximation in $\mathrm{St}$ for an inviscid gas one then retrieves a terminal velocity approximation as derived by Lin & Youdin (2017).

Having taken the terminal velocity approximation, (LABEL:eq:twofluidnondim.4) vanishes, so the viscous stress tensor $\mathrm{T}_{g}$ (13) only appears in equation (LABEL:eq:twofluidnondim.2), with factor of $\mathrm{Re^{-1}}$. As shown in (17), for the terminal velocity approximation to hold $\mathrm{Re}$ must be larger than $\mathcal{O}\left(\mathrm{St}^{-1}\right)$, so

[TABLE]

such that

[TABLE]

the stress tensor can be rewritten as

[TABLE]

It can be shown that in fact, for protoplanetary discs $\mathrm{Re}\gg\mathrm{Ma^{2}}$ for $\alpha$ viscosity. Considering a protoplanetary disc with Reynolds number $\mathrm{Re}={\bar{u}\bar{L}}/{\bar{\nu}}$ one can take the characteristic velocity $\bar{u}$ to be the Keplerian velocity $u_{K}=c_{s}\mathrm{Ma}$, such that $\mathrm{Re}=c_{s}\mathrm{Ma}\bar{L}/\bar{\nu}$. In the alpha viscosity formalism $\nu=\alpha c_{s}H$ (Shakura & Sunyaev, 1973), taking $\bar{L}$ as the orbital radius makes $\bar{L}/H=\mathrm{Ma}$ so the Reynolds number becomes $\mathrm{Re}=\mathrm{Ma^{2}}/\alpha$, given that $\alpha\ll 1$ (Rafikov, 2017); $\mathrm{Re}\gg\mathrm{Ma^{2}}$. It is therefore correct to first order in $\mathrm{St}$ to use the gas viscous stress tensor for protoplanetary discs. The nondimensional equations can therefore be rewritten in one fluid form to include a viscous term as well,

[TABLE]

In dimensional form these become,

[TABLE]

This is the same set of equations derived in Lin & Youdin (2017), the locally isothermal terminal velocity approximation (LITVA), with the addition of a viscous term. The equations (25-27) can also be obtained through assuming an isothermal gas in the Laibe & Price (2014) terminal velocity approximation. The dust then behaves as a cooling term acting on the gas, removing the need to solve a separate dust advection equation.

3 Test Cases

To test our implementation of a one fluid gas dust solver we consider some analytic test cases. We produced analytic solutions for three initial conditions, a small amplitude wave, a steady state shock, and a pure diffusion problem. The latter two could be especially useful as they complement the existing studies of analytic solutions already found in literature like two fluid shock solutions (Lehmann & Wardle, 2018) and a dusty wave dispersion analogous to our first test case (Laibe & Price, 2011). Analytic solutions like these can not only be used as test cases but also give some insight to the behaviour of the equations and better understand simulation results.

3.1 Dusty Wave Solutions

Using linear perturbations it is possible to obtain a dispersion relation for small amplitude sound waves in a dusty gas similarly to what Laibe & Price (2016) did for their one fluid formulation. Starting from equations (25-27) in the inviscid regime and linearising by applying a small wavelike perturbation $\delta(F)=F^{\prime}\psi$ to a background $F_{0}$, where $\psi=e^{i(kx+\omega t)}$ and $F^{\prime}\ll F_{0}$. Given that $F^{\prime}\ll F_{0}$ terms which are of order $\mathcal{O}(F^{\prime 2})$ can be neglected. This leaves us with

[TABLE]

from which a dispersion relation can be calculated. This dispersion relation is complex, as there are both wavelike and damping like parts to the solution.

[TABLE]

where $D_{c}=c_{s}^{1}t_{s}\left(1-P_{0}/c_{s}^{2}\rho_{0}\right)$ or,

[TABLE]

Having a complex dispersion relation requires the perturbations $F^{\prime}$ to be complex too, with the consequence that for dusty sound waves the wave in pressure lags the wave in density by a phase proportional to $f_{d}$, the gas to dust ratio, defined as $f_{d}=1-P/c_{s}^{2}\rho$. This relation provides an analytic solution of the gas dust equations in the linear regime, which we used to test the convergence properties of our scheme. It should be noted that this solution is a wave damping in time where each branch corresponds to a left travelling or right travelling wave, the superposition of the two producing a damped standing wave. It is important to note that the damping in this solution is entirely due to the dust cooling term and not some viscous term.

3.2 Dusty Shock Solutions

A solution involving steep gradients and large amplitudes to the one fluid dust gas approximation is the solution around a shock front. While there exists an analytic solution for this case, it should in fact be noted that the terminal velocity approximation breaks down around shocks. This solution is unphysical, as will be discussed more rigorously in section 4. Here we construct a steady state shock solution; this solution holds true only when the shock has reached steady state as $t\rightarrow\infty$, and the diffusion term is balanced by the advective terms. We could not construct a full solution like the Sod (1978) solutions, as the diffusive term breaks scale invariance for the system of equations, making a solution especially for the rarefaction wave very hard to find. Starting with the system of equation (25-27) in the inviscid regime, constructed in 1D

[TABLE]

Considering the equations(31-33) in the frame co-moving with the shock and allowing for a steady state to be reached,

[TABLE]

the invariants for these equations (35) are,

[TABLE]

where $\rho_{L},P_{L},u_{L}$ are the "left" preshock values, and $\rho_{R},P_{R},u_{R}$ are the "right" postshock values, in a shock travelling from right to left. Additionally $f_{dR}=f_{dL}$, however $f_{d}$ is not constant over the shock, as that would require $f_{d}=0$ or $t_{s}=0$, i.e. either the no dust or perfect coupling solutions respectively. The jump conditions determining $u_{L}$ and $u_{R}$ can be calculated the same way as they would for the gas only case. The jump depends only on the Mach number of the shock and the soundspeed,

[TABLE]

In the dusty case however, the soundspeed in the dusty fluid is the modified soundspeed $c_{s}^{\prime}=f_{0}c_{s}$ The system of equations can be thus reduced to a single ODE (41), which can then be solved numerically to produce the profile of a steady state dusty shock.

[TABLE]

The ODE (41) describing the pressure profile of a shock is first order and therefore has some form of analytic solution (even if not necessarily integrable, as is the case here), which can be found through substitution to be,

[TABLE]

where $\xi=0.5\left(\frac{u_{L}^{2}-a^{2}}{c^{2}-a^{2}}\right)$ and $z=\frac{B-2P}{B-2P_{L}}$, for shocks $u_{L}<a<c)$. $P_{c}$ is the central value of $P$ at the shock; to obtain this value we have to assume that the shock is symmetric, which from the $\tanh(x)$ like solution is a fair assumption (Fig. 1). Equation (42) can be simplified to

[TABLE]

This expression also has no explicit roots for $z$ and therefore $P_{x}$. Thus, a numerical root finder like Newton-Rhapson must be used, producing the same results as those obtained through the numerical integration of equation (35). For all $\xi<1$ this solution is smooth, however on a large scale, this still can be approximated as a jump described by the perfectly coupled jump conditions. As mentioned before, this solution is only an asymptotic solution, that the fluid relaxes into. The rate at which this solution is approached does not behave intuitively; the diffusion timescale actually grows with $t_{s}$. It can be shown that this is the case by changing the time coordinate to a ”diffusion time” $\tau=t\cdot t_{s}$, which removes $t_{s}$ from the cooling term, given that $\tau\propto t_{s}$ the rate at which a time dependant solution approaches the stationary solution; larger stopping times lead to solutions which which evolve slower towards the steady state solution than smaller stopping times. The most obvious case of this is the $t_{s}=0$ perfect coupling case described by Laibe & Price (2014), where the steady state solution is reached immediately. As with the wave solution from section 3.1 the shock smoothing is entirely due to dust gas interactions.

3.3 A purely diffusive solution, and its relation to the porous medium equation

It is also possible to test the behaviour of a scheme for just the diffusion problem. This is actually a useful test, as it is not obvious that a diffusion solver would be well behaved for nonlinear diffusion problems such as this one. Assuming that the evolution is entirely dominated by the dust source term from equation (27), the only equation that remains is the anomalous diffusion equation described by the time evolution of the cooling term,

[TABLE]

There were no real initial conditions in the advection diffusion problem for which this would be the case, so this is probably best thought of as mainly a test problem. Equation 44 is a nonlinear diffusion equation as the diffusion operator can be written in the form $\nabla\cdot\left(D_{P}\nabla P\right)$, which is in fact a sub diffusion equation as $D_{P}=D^{\prime}P^{n}$ where $n>0$ and $D^{\prime}$ is constant in $P$ (Vázquez, 2017). Through the transformation $f_{d}=1-\frac{P}{c_{s}^{2}\rho}$, assuming $\rho$ is constant, we can obtain the equation

[TABLE]

a porous medium equation of the form

[TABLE]

The $n=1$ case for the porous medium equation is also called the Boussinesq equation, which is used to describe the filtration of water underground (Boussinesq, 1903). From a physical standpoint it is interesting that one would recover something of the form of the Boussinesq equation from the set of equations (25-27). The Boussinesq equation describes a filtration of a fluid (water) through a fixed matrix of particles (soil), which for a case when the stopping time is quite long compared to the characteristic time of the system, is actually quite similar to a dusty gas.

There exists a well studied self similar solution to equation (46), the Barenblatt profile Barenblatt (1996), which for the $\mathbb{R}^{1},n=1$ case is described by the expression

[TABLE]

This profile, shown in figure 2, only constitutes a weak solution to the equation (46) as $\partial_{x}B$ is not continuous over all $x$ (Huang, 2013). Due to this, a numerical scheme may misbehave around the boundaries, therefore for a convergence study it might be preferable to use a part of the profile far from the non-smooth points at $u=0$ for any numerical tests. This smooth region can be defined as

[TABLE]

The boundaries in this solution are time dependant and the solution is only stable for $u>0$.

There also exists a separable solution to the Boussinesq equation. This solution takes the form

[TABLE]

The separable solution has the advantage of being smooth and continuous everywhere except at $t=0$, but does have time dependent boundary conditions when solved numerically. Additionally, this solution is not stable for any $t>0$, as the diffusion coefficient in equation (45) becomes negative, making it behave as an anti-diffusion equation. There is a further issue with all these solutions, that being that it is not useful for testing spatial convergence in any differencing scheme higher than the first order. Taking the case of a central difference scheme the spatial derivative can be written as,

[TABLE]

with error terms depending on powers $h^{2}$ and higher. These error terms however also depend on the third spatial derivative of $f_{i}$ which is [math] for this function, meaning that, at least at the first step, there is no spatial discretisation error. It can further be shown that provided the analytic solution can be discretised exactly in space, then the error at all times will not depend on the grid spacing (as shown in appendix A). This kind of problem is quite useful as a test of the time marching scheme, this is because no error arises from the spatial discretisation.

4 Problematic shocks

When studying the solutions to terminal velocity approximation and two fluid shocks, we found some discrepancies in the solutions. Already from a first analysis one can see that in the terminal velocity approximation solution is symmetric and smooth. This does not resemble the two fluid model shocks, where the shock front precedes a region where the velocity difference between the gas and dust decays.

4.1 Two fluid shocks

Solutions describing a two fluid shock can be found, this is not limited to perfectly coupled cases like the ones shown by Laibe & Price (2014). For initial conditions with a jump in a gas there exists a well studied solution in the form of the Sod shock tube (Sod, 1978), such a solution has not been found for two fluid approximation dusty gases. It is however possible to find steady state solutions for shock fronts in this kind of setup (Lehmann & Wardle, 2018). Starting from the two fluid equations for an isothermal gas coupled with a pressureless gas equation (LABEL:eq:twofluid.1) expressed in one dimension,

[TABLE]

Considering a small region around the shock such that the source terms become negligible, we have conserved quantities

[TABLE]

which remain constant over the shock. Assuming that there is no jumps in $\rho_{d}$ and $u_{d}$ the mass fluxes $B_{g}$ and $B_{d}$ must be constant. For a steady state shock

[TABLE]

Considering these invariants it is possible to construct a first order differential equation to describe the profile of the shocks. The momentum equation can be rewritten as an ODE

[TABLE]

and by adding equations (57) and (58) it is possible to eliminate $\rho_{d}$ as

[TABLE]

where $D=\rho_{g}(u_{g}^{2}+c_{s}^{2})+\rho_{d}u_{d}^{2}$ must be constant. This operation yields the ODE

[TABLE]

which can be solved numerically or analytically to produce a density and thus a velocity profile for both the dust and gas. Studying the solution and experiments, it can be seen that when the preshock velocity $u_{R}>c_{s}$ then for a length scale $L_{\mathrm{equil}}\sim c_{s}t_{s}$ the velocity difference between the dust and gas $\Delta u\sim c_{s}$.

4.2 Shocks in the terminal velocity approximation

While the terminal velocity approximation yields well-defined, smooth "shock" profiles (see section 3.2), around shocks the approximation breaks down, rendering these solutions unphysical. The terminal velocity approximation relies on the stopping distance of the dust being very small compared to any other length scale in the fluid. Around a shock this ceases to be true, as the stopping length becomes very long compared to the length scale of the shock. By expanding the equations it is possible to show that in fact the approximation cannot hold around shocks. When describing stationary dusty shocks with a two fluid approximation, there is a finite length over which the gas velocity $u_{g}$ jumps for supersonic to subsonic, while the dust velocity is unchanged. Therefore just behind the shock there is a region where $\Delta u\sim c_{s}$. In the one fluid dusty shock the gas and dust are assumed to be coupled, therefore this region does not appear in this approximation. This results in a first order error around the shock with respect to the two fluid shock. In dimensionless form this error can be written as

[TABLE]

This first order error means that around shocks the terminal velocity approximation breaks down. The only way to approach this problem within the terminal velocity approximation would be to resolve the shock, through higher shock viscosity. Alternatively one could under resolve the shock such that the region where the error appears is much smaller than a single grid cell. Neither of these solutions is fully satisfactory, as increasing the shock viscosity will produce unphysical results, while limiting the resolution might be unsuitable for some problems.

5 Computational methods

5.1 A DUSTY implementation of FARGO3D

In a diffusion solver for a hydrodynamics code there is an inherent tradeoff between stability and speed. Fast and accurate explicit solvers are unstable for long timesteps, as their stability condition is inversely proportional to the spatial resolution squared while the advection step’s stability condition is dependant on just one over the spatial resolution. Implicit integrators allow for longer steps, as they are either unconditionally stable or have a much wider stability region than explicit integrators. Implicit integrators however are less accurate, more computationally expensive, and much more difficult to parallelise. A lot of progress has been made in the field of explicit stiff ODE solvers and these stable solvers can be used to do the timestepping in semi-discretised PDE integrators. We implement a stable integrator to integrate the diffusion term $C(P,\rho)=c_{s}^{2}\nabla\cdot\left[f_{d}\nabla P\right]$ in equation (27). This is the source term from the locally isothermal terminal velocity approximation derived in (Lin & Youdin (2017)), see equation (25-27). Due to FARGO3D being operator split, it is possible to integrate the source terms separately from the advection terms, and produce a solution by adding the two. In our problem, we implemented a stable solver to integrate the diffusion term $C\left(P,\rho\right)$ in a finite difference discretised scheme. For this implementation, we built upon and modified the existing and well tested FARGO3D (Benítez-Llambay & Masset, 2016) operator split code, and added a separate solver to integrate the diffusion term $C(P,\rho)$; the Runge Kutta Legendre second order (RKL2) scheme derived by Meyer et al. (2014). This algorithm provides several benefits with respects to both explicit and implicit schemes. By sacrificing higher order, the scheme gains the stability of an implicit scheme, but at the same time maintains the clarity and ease of paralellisation of an explicit scheme. RKL is not the only unconditionally stable scheme, examples of other explicit unconditionally stable schemes being RKC and ROCK (constructed around Chebyshev polynomials). The stability of these schemes for the integration of parabolic PDEs was analysed in Alexiades et al. (1996), showing that these schemes, when correctly formulated are stable for any step. These RKC methods however have smaller stability regions than RKL methods, meaning that they can produce larger features not present in the true solution. The other implicit approaches are usually inefficient as they require large matrix inversions, as shown in Alexiades et al. (1996) and Caplan et al. (2017). This specific scheme (RKL2) was picked because it was shown to be well behaved with non linear diffusion problems by Meyer et al. (2014), as well as being a very efficient scheme (Caplan et al. (2017)). The RKL2 scheme maintains unconditional stability by being constructed on a Legendre polynomial, which when chosen correctly has a magnitude of $1$, therefore there can be no growing wavemodes, making the scheme stable. For our FARGO3D implementation we used a central finite difference discretisation in space with the previously described RKL2 integrator for the timestepping. This implementation plugs into the FARGO3D machinery without any major changes to the core solver. Some change had to be made to the calculation of the Courant number, or CFL, as the soundspeed in the dusty gas is calculated differently from either an isothermal single gas or an adiabatic gas. We took the gas only sound speed to be the sound speed, as the dust gas sound speed is guaranteed to be smaller than the pure gas sound speed, therefore the CFL step cannot be overestimated, and the code remains stable, at the expense of a slightly shorter timestep than could be possible.

6 Results

We used the solutions from section 3 to test the performance of our dusty version of FARGO3D described in section 5.1. The dusty wave solution was used to run a convergence study for the code. For these runs we used a box size of 10.0 with a wavenumber $k=2\pi/10$. The background values used were $\rho_{0}=50.0$, $P_{0}=25.0$, $v_{0}=0$, $c_{s}=1.0$ and $t_{s}=1.0$. The perturbation was given by taking a density perturbation amplitude of $\rho_{1}=0.01$, all the other perturbation amplitudes being constrained by setting one. The code was run until $t=10.0$, then an $L_{1}$ norm was taken to evaluate the numerical error of the scheme. The RKL2 implementation of the code performed well showing comparable spatial convergence to an RK2 scheme (Fig. 4). The absolute error in the RKL2 scheme is larger than that exhibited by the RK2 scheme, however the RKL2 scheme is much faster than its fully explicit counterpart for any high resolution problem, as the number of calculations scales with $\frac{1}{h^{2}}$ not $\frac{1}{h^{3}}$ as for classic explicit schemes (like RK2), this difference can be clearly seen in fig. 3.

For the diffusion test we ran the diffusion problem from section 3, using the central portion of a Barenblatt solution with $A=36$ in the region $-1<x<1$, far from the non-smooth points. We ran the problem for a time of 0.1 from $t=1.0$ to $t=1.1$. To run this problem we turned off the transport portion of FARGO3D and allowed the code to take the maximum stable time-step (RKL2 being unconditionally stable only ever takes 1 step with a varying number $s$ of supersteps). This experiment shows that increasing the number of RKL stages does not increase the precision of the scheme, only its stability as can be seen in figure 5. In fact the accuracy of an $s$ stage RKL scheme slightly decreases with increasing $s$. This becomes especially obvious in the diffusion test, where we advance only the diffusion term in time. By increasing the spatial resolution in this test we do not reduce the spatial error, but we force smaller timesteps in the RK2 scheme to maintain stability, increasing the precision. In the RKL2 scheme however, the timestep is not reduced, as the number of supersteps is increased to maintain stability instead, leading to the error linearly increasing with the number of supersteps. This in turn means that the RKL2 method is faster, especially at high resolution, being about four times as fast for the $N=3000$ test in terms of total run time (fig. 3).

Regarding the stability of the scheme as a whole, the code is able to remain stable running a discontinuous initial condition, with the shock produced converging on the steady state shock profile calculated in section 3.2, as expected. For this test we used a resolution of 100 cells per unit length, an initial density and pressure jump of a factor of 2 from $\rho_{R}=5.0$ to $\rho_{L}=10.0$ the other parameters being: $t_{s}=0.1$, $c_{s}=1.0$, and $f_{d}=0.5$. The RKL2 scheme under diffused the solution compared to the RK2 scheme as can be seen in figure 6. No real statement can be made here about the relative accuracy of the schemes, as we do not have a time dependent solution for the shock wave. The analytic solution holds as $t\rightarrow\infty$ and we expect the time dependent solution to tend to our analytic solution in the large t limit, but we do not know how a discontinuity evolves towards said solution.

7 Discussion

7.1 Dusty shocks

For an astrophysical flow model it can be argued that the inability to describe shocks is an important flaw. Most astrophysical flows are largely very low density and dominated by external fields, this means that supersonic flows and shocks are quite common. The fact that shocks are common does not mean that there are no scenarios where shocks are not important. In protoplanetary discs, for example the dynamics of many phenomena like the evolution of vortices, and the linear growth of certain instabilities like the streaming instability, or the evolution of the dust distribution in the disc can be studied without needing to model shocks. In these cases the pressure profile is smooth and shocks are not expected to arise.

As we have shown the terminal velocity approximation breaks down around shocks. This means that the local behaviour around shocks is not to be trusted in terminal velocity approximations. If the shock is important to the evolution of the system, one should not expect to obtain correct results from a terminal velocity approximation. In essence the terminal velocity approximation is not well suited for studying structure arising because of or at shocks.

It is not clear however to what extent the error around shocks is an issue for the global behaviour of the system. If the shocks are short lived, and only affect a small area of the domain, the final result may be largely unaffected by the error appearing due to the shocks. Additionally, the length scale of the error is quite small, so it might not be a very serious issue if a shock does arise. As the length scale the error occurs on may, but is not guaranteed to, be much smaller than the grid scale; making the terminal velocity approximation nearly indistinguishable from the two fluid solution. This does not change the fact that caution should be taken when shocks arise in a terminal velocity approximation simulation.

With regards to protoplanetary discs, shocks are probably important in the evolution of the dust distribution, as well as the global evolution, contributing to among other things, the vertical mixing in the disc (Boley & Durisen, 2006) and gaps in the dust density. This limits the useful scope for the terminal velocity approximations, as the main advantages of one fluid approximations over two fluid models are the reduced computational cost and better scaling with system size, making terminal velocity approximations very well suited for high resolution global simulations of discs. The inability to model shocks means one has to be careful when using the terminal velocity approximation in cases where shocks may be important.

7.2 Slightly viscous dusty gases

Evidence suggests that while protoplanetary discs probably have a very low turbulent viscosity, the value of this viscosity is probably non zero (Flaherty et al., 2017). Having a model to simulate slightly viscous dusty gases is therefore of great value. The result obtained in section 2 is especially relevant in that the dusty viscous stress tensor is identical to a pure gas stress tensor and would thus require no modification to a code to be implemented. The addition of a viscosity model is also important in that it adds a way to smooth shocks. As discussed in section 4 if the velocity jump in the shock is very small, then the first order error does not appear. In a more viscous discs very strong shocks are less common, this increases the range of problems to which the terminal velocity approximation could be safely applied to.

7.3 Dusty vortices

We propose that this kind of model is very well suited for the simulation of dusty vortices in protoplanetary discs. To fully model the stability of vortices, a high resolution 3D simulation is required, which is made feasible by the use of a terminal velocity approximation. If the problem is then chosen correctly, i.e. the vortex radius is smaller than the local disc scale height (Paardekooper et al., 2010; Goodman & Rafikov, 2001), then no shocks are expected to arise. This kind of setup can be used to investigate the potential of vortices to act as dust traps to facilitate grain growth, or the long term stability of vortices in dusty discs. As mentioned earlier shocks are not present in this setup, so a very non dissipative scheme could also be used, for example a spectral scheme, which allows for the viscosity to be controlled and minimised. This is especially useful in the case of an inviscid disc, where the vortex lifetime would be controlled by the dust density alone.

8 Conclusions

Single fluid terminal velocity approximations for dusty gases are very promising techniques for the study the evolution of a dusty gas in a protoplanetary disc. The computational complexity of this kind of approximation is lower than in a two fluid approach, especially as the correct dust velocity is recovered even in the long timestep limit on one fluid approximations, which is not the case for two fluid solvers. In this paper we show that the approximation is not limited to inviscid fluids, as for an $\alpha$ viscosity model the viscous stress tensor is unchanged at first order in Stokes number. This result implies that there is no need to modify the one fluid description for it to be compatible with $\alpha$ viscosity model, by consequence it is simple to implement the terminal velocity approximation in an $\alpha$ viscosity model, as no additional changes are required in comparison to an inviscid solver. We find however an important limitation in the terminal velocity approximation, being its inability to model shocks correctly. Around the shock an error of order the soundspeed is made in calculating the barycenter velocity, and the correct structure of the shock is not retrieved. The scale of the error is rather small, given by $\mathrm{Ma}_{\mathrm{shock}}c_{s}t_{s}$, the product of stopping time sound speed and shock Mach number, this can however be larger than the gridscale in a high resolution simulation or in cases where $\mathrm{Ma}$ is very large, so care must be taken when simulating initial conditions that may see shocks arising. While the terminal velocity approximation is not suitable for studying shocks in protoplanetary discs, however there are many cases in planet formation and evolution where this is not an issue. A good example of this is small (less than 1 scale height in radius) vortices. These, in a protoplanetary disc are very important to the evolution of dust distribution, while not giving rise to shocks. Because of this, the terminal velocity approximation is especially well suited to study the behaviour of dust laden vortices. This will be the subject of the next paper in this series.

Acknowledgements

We thank the referee, Jeff Oishi, for his thorough response, which much improved the quality and clarity of our paper. FL is funded by an STFC studentship. SJP is funded by a Royal Society University Research Fellowship.

Appendix A Proof of non mixing errors for quadratic in space problems solved with centered differences schemes

Taking a numerical approximation $\tilde{u}_{i}^{n}$ of an analytic function $u_{i}^{n}$ constructed at timestep $n$ at grid node $i$ using a centered in space and $A$th order explicit in time scheme. One can write this as,

[TABLE]

where $\epsilon_{i}^{n}$ is the numerical error at timestep $n$ in cell $i$. The error term $\epsilon_{i}^{n}$ can be expanded such that,

[TABLE]

Postulating that

[TABLE]

meaning that the numerical error at any time depends only the time step and not on the spatial step. One can show that given that 66 is true for $n$ then it is also true for $n+1$. Starting with the general expression for the numerical solution at $n+1$,

[TABLE]

the error terms then become,

[TABLE]

From 66 it can be seen that,

[TABLE]

therefore,

[TABLE]

From our initial assumptions about the analytic function however we know that it can only be quadratic in space all times for which it is well posed, making the third spatial derivative [math]. The error at step $n+1$ can then be rewritten as

[TABLE]

Given that at $n=1$ 66 is trivially true, as $\tilde{u}_{i}^{0}\equiv u_{i}^{0}$ due to $\tilde{u}_{i}^{0}$ being the initial conditions, then by induction, a numerical approximation $\tilde{u}_{i}^{n}$ to a solution $u_{i}^{n}$, analytically discretisable at any time, must obey the identity

[TABLE]