# A Fast second-order solver for stiff multifluid dust and gas hydrodynamics

Leonardo Krapp Department of Astronomy and Steward Observatory, University of
Arizona, Tucson, AZ 85721, USA  Juan Garrido-Deutelmoser Instituto de
Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna
4860, 782-0436 Macul, Santiago, Chile  Pablo Benítez-Llambay Facultad de
Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal las Torres
2640, Peñalolén, Chile  Data Observatory Foundation, ANID Technology Center
No. DO210001  Kaitlin M. Kratter Department of Astronomy and Steward
Observatory, University of Arizona, Tucson, AZ 85721, USA

###### Abstract

We present MDIRK: a Multifluid second-order Diagonally-Implicit Runge-Kutta
method to study momentum transfer between gas and an arbitrary number ($N$) of
dust species. The method integrates the equations of hydrodynamics with an
Implicit Explicit (IMEX) scheme and solves the stiff source term in the
momentum equation with a diagonally-implicit asymptotically stable Runge-Kutta
method (DIRK). In particular, DIRK admits a simple analytical solution that
can be evaluated with $\mathcal{O}(N)$ operations, instead of standard matrix
inversion, which is $\mathcal{O}(N)^{3}$ . Therefore the analytical solution
significantly reduces the computational cost of the multifluid method, making
it suitable for studying the dynamics of systems with particle-size
distributions. We demonstrate that the method conserves momentum to machine
precision and converges to the correct equilibrium solution with constant
external acceleration. To validate our numerical method we present a series of
simple hydrodynamic tests, including damping of sound waves, dusty shocks, a
multi-fluid dusty Jeans instability, and a steady-state gas-dust drift
calculation. The simplicity of MDIRK lays the groundwork to build fast high-
order asymptotically stable multifluid methods.

algorithms - hydrodynamics - accretion disks - protoplanetary disks

##  1 Introduction

Theories of star and planet formation are grounded in the study of
astrophysical fluids. Consequently, many decades of research have been devoted
to the development of algorithms to capture their non-linear evolution
(LeVeque, 1992; Toro, 2009). For a subset of systems, however, it is
insufficient to consider only the evolution of a single species. For example,
in cold systems, dust is often the dominant opacity source (Pollack et al.,
1994; Draine, 2003; Birnstiel et al., 2018), and thus dictates the thermal
evolution of the gas. However, for a range of gas densities and dust particle
sizes, the coupling between molecular gas and dust is imperfect. The density
ratio between gas and different dust species may not be spatially or
temporally constant, leading to non-negligible relative velocities, and thus
momentum exchange (Epstein, 1924; Whipple, 1972; Adachi et al., 1976) These
multi-species systems are susceptible to a wide range of multi-fluid
instabilites that can substantially alter the evolutionary pathway of a system
(Youdin & Goodman, 2005; Squire & Hopkins, 2018; Krapp & Benítez-Llambay,
2020; Hopkins et al., 2020). A corollary of the fact that dust controls the
thermodynamic behavior of these systems is that it also dramatically impacts
the observational signatures of these objects, especially in continuum
emission at millimeter - infrared wavelengths, and scattered light at optical
wavelengths. Thus accurately predicting the non-uniform distribution of dust
grains of various sizes is also vital for predicting and interpreting data
from ALMA (e.g., Andrews et al., 2018), SPHERE (e.g., Avenhaus et al., 2018)
and forecoming JWST observations.

To fully capture the impact of dust on astrophysical systems, we require
numerical methods that can tackle the problem of momentum exchange between $N$
species in an efficient manner. Several approaches have been developed to
study the interaction between dust and gas numerically. So-called “Hybrid”
integration schemes utilize Eulerian grids for the gas, but separate
Lagrangian super particles for dust (e.g., Youdin & Johansen, 2007; Balsara et
al., 2009; Miniati, 2010; Yang & Johansen, 2016; Mignone et al., 2019). In
smoothed particle hydrodynamics (SPH), both gas and dust are treated as
Lagrangian particles, with careful choice of the kernel for the drag terms
(e.g., Laibe & Price, 2012, 2014; Booth et al., 2015; Price et al., 2018;
Stoyanovskaya et al., 2018). Staggered Semi-analytic Methods (e.g., Fung &
Muley, 2019), are best applied to tightly coupled distributions, where a
staggered time step allows for accurate calculation of dust grain terminal
velocity. Finally, there exists a range of Euler-Euler methods where all
species utilize the same grids, (e.g., Johansen et al., 2004; Paardekooper &
Mellema, 2006; Hanasz et al., 2010) including first-order multifluid (Benítez-
Llambay et al., 2019), second-order multifluid (Huang & Bai, 2022), and
single-fluid with modified energy equation (Chen & Lin, 2018).

Solving for the momentum exchange between gas and particles requires special
focus since it introduces source terms that modify the stiffness ratio to the
system of equations of motion. These source terms are usually modeled as drag
forces that depend on the relative velocity between gas and dust (in both the
linear and non-linear regime), the dust-to-gas density ratio, and the
characteristic stopping time (Whipple, 1972; Adachi et al., 1976). In the
limit of strong coupling (small stopping time relative to a characteristic
time scale) and / or strong feedback (large dust-to-gas density ratio) the
time-explicit numerical integration of drag forces requires an excessively
small time step in order to preserve stability. Such small timesteps
dramatically increase the computational cost of a given calculation.

Alternative methods that circumvent the reduction of the time step are time-
implicit methods. However it is not straightforward to couple a time-implicit
drag-force scheme with a time-explicit numerical method for hydrodynamics,
(see e.g. Huang & Bai 2022). In general, this coupling requires a splitting of
the hydrodynamics equations which, if not done properly, downgrades the
precision of the numerical solution. Moreover, combined time-implicit and
time-explicit methods must still convergence to the correct equilibrium
solution. One example of a time-implicit method was presented in Benítez-
Llambay et al. (2019). The method coupled a first-order time-implicit scheme
for the drag force with an explicit first-order hydrodynamics solver. This
multi-species numerical solver is unconditionally and asymptotically stable
and conserves momentum to machine precision. Moreover, the complexity of the
scheme is order $\mathcal{O}(N)$, with $N$ the number of fluids (Krapp &
Benítez-Llambay, 2020).

In this work, we present MDIRK: a novel multifluid second-order method that is
asymptotically stable, conserves momentum to machine precision, and has
complexity $\mathcal{O}(N)$. We solve the hydrodynamics equations a two-stage
Runge-Kutta Implicit-Explicit (IMEX) method. In addition, we highlight the
benefits of a stiff-accurate multifluid solver, and lay the groundwork to
develop even high-order implicit schemes that retain complexity
$\mathcal{O}(N)$.

This work is organized as follows. In Section 2 we introduce the fluid
equations for a multifluid mixture of gas and dust and showcase the numerical
method. In Section 3 we present the results from a test-suite and validate the
properties of our second-order method. In Section 4 we summarize the main
results of our work and suggest possible extensions to achieve higher order
accuracy.

##  2 Numerical Method

We now describe the numerical method for a 1D system. All of its important
properties such as momentum conservation, convergence to the correct
equilibrium, and order of accuracy are independent of the number of
dimensions. Extensions to 2D and 3D are straightforward, because the drag
force does not couple momentum in orthogonal directions.

Consider a system of $N-1$ dust species that exchange momentum with an
inviscid gas through drag forces. Gas density and velocity are defined as
$\rho_{\rm g}$ and $v_{\rm g}$, respectively. Equivalently, for the $i$-th
dust species we define the density and velocity as $\rho_{{\rm d},i}$ and
${v}_{{\rm d},i}$. We neglect the momentum exchange between dust species. We
assume an isothermal equation of state and the gas pressure is defined as
$P=c^{2}_{\rm s}\rho_{\rm g}$, with $c_{\rm s}$ as the sound speed. Under
these assumptions, the continuity and momentum equations for gas and dust
species are

$\displaystyle\partial_{t}\rho_{\rm g}+\partial_{x}\left(\rho_{\rm g}{v}_{\rm
g}\right)$ $\displaystyle=0\,,$ (1) $\displaystyle\partial_{t}\rho_{{\rm
d},i}+\partial_{x}\left(\rho_{{\rm d},i}{v}_{{\rm d},i}\right)$
$\displaystyle=0\,,$ (2) $\displaystyle\partial_{t}(\rho_{\rm g}{v}_{\rm
g})+\partial_{x}\left(\rho_{\rm g}{v}^{2}_{\rm g}+P\right)$
$\displaystyle=-\sum^{N-1}_{i=1}\rho_{{\rm g}}\alpha_{{\rm g},i}\left({v}_{\rm
g}-{v}_{{\rm d},i}\right)+G_{{\rm g}}\,,$ (3)
$\displaystyle\partial_{t}(\rho_{{\rm d},i}{v}_{{\rm
d},i})+\partial_{x}\left(\rho_{{\rm d},i}{v}^{2}_{{\rm d},i}\right)$
$\displaystyle=-\rho_{{\rm d},i}\alpha_{{\rm d},i}\left({v}_{{\rm
d},i}-{v}_{\rm g}\right)+G_{{\rm d},i}\,,$ (4)

with $i=1,\dotsc,N-1$. The right-hand sides of Equations (3) and (4) containt
the drag force terms and any external forces $G_{\rm g}$ and $G_{{\rm d},i}$,
respectively. The rates of momentum exchange are $\alpha_{{\rm g},i}$ and
$\alpha_{{\rm d},i}$. Momentum conservation implies

$\rho_{\rm g}\alpha_{{\rm g},i}=\rho_{{\rm d},i}\alpha_{{\rm d},i}\,,$ (5)

therefore $\alpha_{{\rm g},i}=\epsilon_{i}\alpha_{{\rm d},i}$, with
$\epsilon_{i}=\rho_{{\rm d},i}/\rho_{\rm g}$. The collision rates may depend
on the density, sound speed, and relative velocity between species. To rewrite
the momentum equations in matrix form, we define

$\displaystyle{\bf u}$ $\displaystyle=(\rho_{\rm g}v_{\rm g},\rho_{{\rm
d},1}v_{{\rm d},1},\dotsc,\rho_{{\rm d},N-1}v_{{\rm d},N-1})\,,$ (6)
$\mathcal{F}$ $\displaystyle=(\rho_{\rm g}v^{2}_{\rm g}+P,\rho_{{\rm
d},1}v^{2}_{{\rm d},1},\dotsc,\rho_{{\rm d},N-1}v^{2}_{{\rm d},N-1})\,,$ (7)
$\displaystyle{\bf G}$ $\displaystyle=(G_{\rm g},G_{{\rm d},1},\dotsc,G_{{\rm
d},N-1})\,,$ (8)

now, Eqs. (3)-(4) can be cast into the form

$\displaystyle\partial_{t}{\bf
u}+\partial_{x}\text{\boldmath$\mathcal{F}$}={\bf M}{\bf u}+{\bf G}\,,$ (9)

where ${\bf M}$ is the $N\times N$ matrix.

$\displaystyle{\bf M}$
$\displaystyle=\begin{pmatrix}-\sum^{N-1}_{k=1}\epsilon_{k}\alpha_{{\rm
d},k}&\alpha_{{\rm d},1}&\alpha_{{\rm d},2}&\dotsc&\alpha_{{\rm d},N-1}\\\
\epsilon_{1}\alpha_{{\rm d},1}&-\alpha_{{\rm d},1}&0&\dotsc&0\\\
\epsilon_{2}\alpha_{{\rm d},2}&0&-\alpha_{{\rm d},2}&\dotsc&0\\\
\dotsc&\dotsc&\dotsc&\dotsc&\dotsc\\\ \epsilon_{N-1}\alpha_{{\rm
d},N-1}&0&0&\dotsc&-\alpha_{{\rm d},N-1}\\\ \end{pmatrix}\,.$ (10)

Note that in our case the matrix ${\bf M}$ is defined from the RHS of the
momentum equation in terms of conserved variables, not primitives. This is a
subtle distinction with the similar coupling matrix ${\bf M_{\rm BKP19}}$,
defined in equation (11) of Benítez-Llambay et al. (2019). In this work the
matrix ${\bf M}=-{\bf R}{\bf M}_{\rm BKP19}{\bf R}^{-1}$, where ${\bf R}$ a
diagonal matrix such that $R_{ii}=\rho_{i}$. The relation between these
matrices, their eigenvectors, and eigenvalues will be utilized to explain the
properties of the numerical method. Both matrices are similar to a real
symmetric matrix and therefore diagonalizable. Moreover, if $\lambda_{k}$ is
an eigenvalue of ${\bf M}$, then $-\lambda_{k}$ is an eigenvalue of ${\bf
M}_{\rm BKP19}$. Note that the eigenvalues of ${\bf M}_{\rm BKP19}$ are all
real and non-negative (there is only one zero eigenvalue) (see Benítez-Llambay
et al., 2019). It is also straightforward to show that if ${\bf x}$ is an
eigenvector of ${\bf M}_{\rm BKP19}$, then $\tilde{\bf x}=-{\bf R}{\bf x}$ is
an eigenvector of the matrix ${\bf M}$. Since ${\bf 1}^{\rm T}=(1,1,\dotsc,1)$
is eigenvector of ${\bf M}_{\rm BKP19}$, it follows that $\tilde{\bf 1}^{\rm
T}=-(\rho_{\rm g},\rho_{{\rm d},1},\dotsc,\rho_{{\rm d},N-1})$ is eigenvector
of ${\bf M}$.

###  2.1 Second-order Implicit-Explicit Method (IMEX)

In this section, we describe the Runge-Kutta IMEX numerical method that will
be utilized to solve the hydrodynamics equations. We begin by focusing on the
implementation of the more complex momentum equation; the continuity equation
is shown in Section 2.2. The novelty of this second-order Runge-Kutta scheme,
from (Ascher et al., 1997), is that it incorporates an explicit term and
implicit term together in a single time update. As we will demonstrate, the
method conserves momentum to machine precision, is asymptotically stable
(L-stable), and converges to the correct equilibrium. We can write Eq. (9) as,

$\displaystyle\partial_{t}{\bf u}$ $\displaystyle=\left[\mathcal{H}({\bf
u})\right]_{\rm exp}+\left[{\bf M}{\bf u}\right]_{\rm imp}\,,$ (11)

with $\mathcal{H}({\bf u})=-\partial_{x}\text{\boldmath$\mathcal{F}$}+{\bf
G}$. The brackets $\left[\right]_{\rm exp}$ and $\left[\right]_{\rm imp}$
indicate which terms are time-explicit and time-implicit, respectively. For
our particular system of equations, the term solved with a diagonally implicit
method admits an analytical solution (Krapp & Benítez-Llambay, 2020), which
transforms the complexity of the method from $\mathcal{O}(N^{3})$ to
$\mathcal{O}(N)$, making it suitable for a system with a large number of
species. Consider the momenta vector ${\bf u}(t_{n})={\bf u}^{n}$. The full
update to a new time $t=t_{n}+\Delta t$ is obtained as follows:

$\displaystyle{\bf u}^{n+1}$ $\displaystyle={\bf u}^{n}+\delta\Delta
t\mathcal{H}({\bf u}^{n})+(1-\delta)\Delta t\mathcal{H}({\bf w}_{1})$
$\displaystyle+b_{1}\Delta t{\bf M}{\bf w}_{1}+b_{2}\Delta t{\bf M}{\bf
w}_{2}$ (12)

where ${\bf w}_{1}$ and ${\bf w}_{2}$ are the solution of Eqs. (13)-(14),
written as

$\displaystyle{\bf w}_{1}$ $\displaystyle={\bf u}^{n}+\tilde{\gamma}\Delta
t\mathcal{H}({\bf u}^{n})+\gamma\Delta t{\bf M}{\bf w}_{1}$ (13)
$\displaystyle{\bf w}_{2}$ $\displaystyle={\bf u}^{n}+\delta\Delta
t\mathcal{H}({\bf u}^{n})+(1-\delta)\Delta t\mathcal{H}({\bf w}_{1})$
$\displaystyle+\beta\Delta t{\bf M}{\bf w}_{1}+\gamma\Delta t{\bf M}{\bf
w}_{2}\,.$ (14)

Alternatively, the solution can be also obtained in terms of ${\bf k}_{1}={\bf
M}{\bf w}_{1}$ and ${\bf k}_{2}={\bf M}{\bf w}_{2}$, by applying the ${\bf M}$
matrix to Eqs .(13)-(14).

$\displaystyle\left({\bf I}-\gamma\Delta t{\bf M}\right){\bf k}_{1}$
$\displaystyle={\bf M}\left({\bf u}^{n}+\tilde{\gamma}\Delta t\mathcal{H}({\bf
u}^{n})\right)$ (15) $\displaystyle\left({\bf I}-\gamma\Delta t{\bf
M}\right){\bf k}_{2}$ $\displaystyle={\bf M}\left({\bf u}^{n}+\delta\Delta
t\mathcal{H}({\bf u}^{n})+\beta\Delta t{\bf k}_{1}\right)$ (16)
$\displaystyle+\,\,{\bf M}\left((1-\delta)\Delta t\mathcal{H}({\bf
w}_{1})\right).$

The major advantage of the diagonally implicit method is that one first solves
for ${\bf k}_{1}$ and then for ${\bf k}_{2}$, which simplifies the
calculations. Note that this is not the case in fully implicit methods since
Eq. (15) will also depend on ${\bf k}_{2}$. Moreover, we can find the
analytical solution of the linear system, which has a complexity
$\mathcal{O}(N)$ (c.f. Krapp & Benítez-Llambay, 2020). For instance, consider
the arbitrary vector ${\bf q}^{n}=(q^{n}_{\rm g},q^{n}_{{\rm
d},1},\dotsc,q^{n}_{{\rm d},N-1})$ and the equation

$\left({\bf I}-\gamma\Delta t{\bf M}\right){\bf k}^{n}={\bf M}{\bf q^{n}}\,,$
(17)

the solution ${\bf k}^{n}=(k^{n}_{\rm g},k^{n}_{{\rm d},1},\dotsc,k^{n}_{{\rm
d},N-1})$ has the components

$\displaystyle{k}^{n}_{\rm g}$ $\displaystyle=\frac{\mathcal{A}-q^{n}_{\rm
g}\mathcal{B}}{1+\gamma\Delta t\mathcal{B}}$ (18) $\displaystyle{k}^{n}_{{\rm
d},i}$ $\displaystyle=\frac{\alpha^{n}_{i}}{1+\gamma\Delta
t\alpha^{n}_{i}}\left(\epsilon^{n}_{i}q^{n}_{\rm g}-q^{n}_{{\rm
d},i}+\gamma\Delta t\epsilon^{n}_{i}k^{n}_{\rm g}\right),$ (19)

for $i=1,\dotsc,N-1$, with

$\displaystyle\mathcal{A}$
$\displaystyle=\sum^{N-1}_{k=1}\frac{\alpha^{n}_{{\rm d},k}q^{n}_{{\rm
d},k}}{(1+\gamma\Delta t\alpha^{n}_{{\rm
d},k})}\,,\quad\mathcal{B}=\sum^{N-1}_{k=1}\frac{\epsilon^{n}_{k}\alpha^{n}_{{\rm
d},k}}{(1+\gamma\Delta t\alpha^{n}_{{\rm d},k})}\,.$ (20)

It can be shown that $\gamma=1\pm 1/\sqrt{2}$ are the only two possible values
for $\gamma$ that provide a second-order accurate and asymptotically stable
two-stages diagonally implicit Runge-Kutta (Rosenbrock, 1963; Alexander,
1977). The rest of the coefficients are chosen to satisfy the order condition
and asymptotic behavior of the solution. In this work, we adopt the values
from (Ascher et al., 1997) (see also Appendix B)

$\displaystyle\tilde{\gamma}$
$\displaystyle=\gamma\,,\quad\delta=1-1/(2\gamma)\,,\quad\beta=1-\gamma\,,$
$\displaystyle b_{1}$ $\displaystyle=1-\gamma\,,\quad b_{2}=\gamma\,.$ (21)

An alternative version of the second-order Implicit-Explicit Runge-Kutta can
be found in Pareschi & Russo (2005) with coefficients $\tilde{\gamma}=1$,
$\delta=1/2$, $\beta=1-2\gamma$, $b_{1}=1/2$ and $b_{2}=1/2$. This alternative
version has an explicit part that matches the optimal second-order Runge-Kutta
derived in Shu & Osher (1988), with important Total-Variation-Diminishing
(TVD) properties. This TVD property also holds for the Runge-Kutta utilized in
this work when $\gamma=1+1/\sqrt{2}$, provided a Courant number $C\lesssim
0.58$.

###  2.2 MDIRK

To solve the fluid equations we have reduced the Runge-Kutta described in
Section 2.1 into a more convenient form, taking into account the continuity
equation as well. To describe the algorithm we define the discrete vectors at
time $t_{n}$ as ${\bf U}_{\rm g}^{n}=(\rho^{n}_{\rm g},\rho^{n}_{\rm
g}v^{n}_{\rm g})=(U^{n}_{\rm g0},U^{n}_{\rm g1})$ and ${\bf K}^{n}_{\rm
g}=(0,k^{n}_{\rm g})$ for the gas, and ${\bf U}_{{\rm
d},i}^{n}=(\rho^{n}_{{\rm d},i},\rho^{n}_{{\rm d},i}v^{n}_{{\rm
d},i})=(U^{n}_{\rm d0},U^{n}_{\rm d1})$, ${\bf K}^{n}_{{\rm
d},i}=(0,k^{n}_{{\rm d},i})$ for the $i$-th dust-species, with $k^{n}_{\rm g}$
and $k^{n}_{{\rm d},i}$ defined in Eq. (18) and Eq. (19), respectively.
Therefore, ${\bf K}_{\rm g}$ and ${\bf K}_{{\rm g},i}$ depend on the density
and momentum of all the fluids. The zero components of these vectors ensure
that the drag force is only considered in the momentum equation.

The flux of mass and momentum are included in the $\mathcal{L}$ operator, as
well as external forces, except of course the drag force. To calculate the
fluxes we utilize a standard Riemann solver (see Appendix E).

For each fluid111Since Eqs. (2.2) and (23) are agnostic of the type of fluid
we omitted the subscript that specifies gas or dust species., a full update of
the conserved variables from time step $t_{n}$ to time step
$t_{n+1}=t_{n}+\Delta t$ is obtained as follows

$\displaystyle{\bf U}^{*}$ $\displaystyle={\bf U}^{n}+\tilde{\gamma}\Delta
t\,\mathcal{L}({\bf U}^{n})\,$ $\displaystyle{\bf U}^{**}$ $\displaystyle={\bf
U}^{n}+(1-\delta)\Delta t\,\mathcal{L}({\bf U}^{*}+\gamma\Delta t{\bf
K}^{*})\,$ $\displaystyle{\bf U}^{n+1}$ $\displaystyle={\bf
U}^{p}+(b_{1}-\beta)\Delta t{\bf K}^{*}+b_{2}\Delta t{\bf K}^{p}\,,$ (22)

with

${\bf U}^{p}={\bf U}^{**}+\frac{\delta}{\tilde{\gamma}}({\bf U}^{*}-{\bf
U}^{n})+\beta\Delta t{\bf K}^{*}\,.$ (23)

The coefficients $\tilde{\gamma},\delta,\beta,b_{1}$ and $b_{2}$ are given in
Eq. (2.1). The vectors ${\bf K}^{*}$ and ${\bf K}^{p}$ are obtained from the
analytical solution given in Eqs. (18)-(19). When utilizing these analytical
solutions, the densities and momenta must be taken from the partial updates
${\bf U}^{*}$ and ${\bf U}^{p}$, respectively. For example,

$\displaystyle{\bf K}^{*}_{\rm g}$ $\displaystyle=\begin{bmatrix}0\\\
\frac{\mathcal{A}-U^{*}_{\rm g1}\mathcal{B}}{1+\gamma\Delta
t\mathcal{B}}\end{bmatrix}\,,$ (24) $\displaystyle{\bf K}^{*}_{{\rm d},i}$
$\displaystyle=\frac{\alpha^{*}_{i}}{1+\gamma\Delta
t\alpha^{*}_{i}}\left(\begin{bmatrix}0\\\ \epsilon^{*}_{i}U^{*}_{\rm
g1}-U^{*}_{{\rm d1},i}\end{bmatrix}+\gamma\Delta t\epsilon^{*}_{i}{\bf
K}^{*}_{\rm g}\right)$ (25)

with

$\displaystyle\mathcal{A}=\sum^{N-1}_{k=1}\frac{\alpha^{*}_{{\rm
d},k}U^{*}_{{\rm d1},k}}{(1+\gamma\Delta t\alpha^{*}_{{\rm
d},k})}\,,\quad\mathcal{B}=\sum^{N-1}_{k=1}\frac{\epsilon^{*}_{k}\alpha^{*}_{{\rm
d},k}}{(1+\gamma\Delta t\alpha^{*}_{{\rm d},k})}\,.$ (26)

Note that when the vector ${\bf K}$ is set to zero (no drag forces), second-
order accuracy can be obtained by setting $\tilde{\gamma}=1$ and $\delta=1/2$.

Figure 1: From left to right we plot the evolution of the velocity for the gas
and two dust species. The dashed line correspond to the analytical solution,
given by Equation (30). The unfilled circles were obtain with different
integrators, where the green lines represent the first-order method, red and
purples lines the second-order method with $\gamma=\gamma^{+}=1+1/\sqrt{2}$
and $\gamma=\gamma^{-}=1-1/\sqrt{2}$, respectively. In the rightmost panel, we
plot, the error by Equation (31) as function of the time step $\Delta t$, for
the three cases. The dotted black lines represent the slope for a first and
second order. Upper and lower panels show the results for the tests 1 and 2,
respectively. The time step used are $\Delta t_{1}=0.005$, $\Delta
t_{2}=0.05$, respectively. See Section 3.1.

###  2.3 Momentum conservation

The drag-force between species is an internal force that cannot change the
total momentum of the system. Our second-order solver has this property as
well. Momentum conservation follows from

$\displaystyle U_{{\rm g}1}^{n+1}+\sum_{k}U_{{\rm d}1,k}^{n+1}$
$\displaystyle=U_{{\rm g}1}^{n}+\sum_{k}U_{{\rm d}1,k}^{n}+b_{1}\Delta
t\left(k^{*}_{{\rm g}}+\sum_{k}k^{*}_{{\rm d},k}\right)$
$\displaystyle+b_{2}\Delta t\left(k^{p}_{{\rm g}}+\sum_{k}k^{p}_{{\rm
d},k}\right)$ $\displaystyle=U_{{\rm g}1}^{n}+\sum_{k}U_{{\rm d}1,k}^{n}\,,$
(27)

where for any arbitrary ${\bf q}$ (which for example can be taken as ${\bf
u}^{*}$ or ${\bf u}^{p}$), the vector ${\bf k}$ satisfies

$\displaystyle k_{{\rm g}}+\sum_{k}k_{{\rm d},k}$ $\displaystyle=k_{{\rm
g}}+\sum_{k}\frac{\alpha_{k}}{1+\gamma\Delta
t\alpha_{k}}\left(\epsilon_{k}q_{\rm g}-q_{{\rm d},k}+\gamma\Delta
t\epsilon_{k}k_{{\rm g}}\right)$ $\displaystyle=k_{{\rm
g}}\left(1+\gamma\Delta t\mathcal{B}\right)-\left(\mathcal{A}-q_{\rm
g}\mathcal{B}\right)$ $\displaystyle=\frac{\mathcal{A}-q_{\rm
g}\mathcal{B}}{(1+\gamma\Delta t\mathcal{B})}\,\left(1+\gamma\Delta
t\mathcal{B}\right)-\left(\mathcal{A}-q_{\rm g}\mathcal{B}\right)$
$\displaystyle=0\,,$ (28)

recalling $\mathcal{A}$ and $\mathcal{B}$ as defined in Eq. (20). We have also
utilized the analytical solutions provided in Eqs. (18)-(19) (see Section
2.1).

##  3 Method Test

In this section, we present a test-suite that validates our second-order
method, and illustrates that the computational costs scales linearly with
number of species. We first compare our implicit second-order drag force
solver with the first-order implicit solver of Benítez-Llambay et al. (2019).
Next, we present four test cases using MDIRK method: damped sound waves, a
multi-fluid shock, an idealized, multi-fluid Jeans instability, and steady-
state drift in a 2D shearing box.

###  3.1 Damping

As we show in Appendix A, our drag force solver provides solutions that
converge asymptotically to the velocity of the center of mass. We show this
property and study the accuracy of the method by solving the momentum exchange
between gas and two dust species. The problem reduces to solving the following
equation for the momenta of all fluids

$\partial_{t}{\bf u}={\bf M}{\bf u}\,.$ (29)

for a given initial condition ${\bf u}_{0}$. After solving the eigenvalue
problem ${\bf M}\bar{\bf x}=\lambda\bar{\bf x}$, the solution of Eq. (29) can
be expressed as ${\bf u}(t)=\sum_{i=0}^{N-1}\bar{c}_{i}e^{\lambda_{i}t}$. To
compare with the work of Huang & Bai (2022), we perform their tests B and C.
Both tests consider one gas and two dust species with constant stopping times
and homogeneous density and velocity backgrounds. The cases differ by a factor
of 1000 in damping timescale. In terms of the velocity, the analytical
solution corresponds to

$v_{\rm analytical}=v_{\rm CM}+c_{1}e^{\lambda_{1}t}+c_{2}e^{\lambda_{2}t}\,,$
(30)

where $v_{\rm CM}$ is the velocity of the center of mass. The values of the
eigenvalues $\lambda_{1}$, $\lambda_{2}$ and the coefficients $c_{1}$ and
$c_{2}$ can be obtained from table 1 of Huang & Bai (2022).

In Figure 1 we show a comparison between the numerical and the analytical
solutions, together with the numerical error. We include the first-order
solution obtained with the method of Benítez-Llambay et al. (2019) (green
circle curve). In addition, the relative error with respect to the analytic
solution as a function of the time step is displayed in the rightmost panel of
Figure 1. The error is obtained as

${\rm e}_{1}=\frac{1}{t_{\rm max}-t_{\rm min}}\sum_{k=0}^{2}\frac{\lvert
v_{{\rm numeric},k}(\Delta t)-v_{{\rm analytic},k}\rvert}{v_{{\rm
analytic},k}}\Delta t\,,$ (31)

where $\Delta t$ is the time step and $k=0,1,2$ is the index of gas, and dust
species, respectively.

The error defined by Eq. (31) corresponds to the L1 norm when $\Delta
t\lesssim t_{\rm damp}$, where $t_{\rm damp}$ is a characteristic damping time
to reach the equilibrium velocity. For the case of test B, $t_{\rm damp}\simeq
0.01$, whereas in test C, $t_{\rm damp}\simeq 10.0$. We calculate the
numerical solution for different time steps starting at $\Delta t=10^{-4}$
until a final time $t_{\rm final}>t_{\rm damp}$. In the regime of $\Delta
t<t_{\rm damp}$, we obtain ${\rm e}_{1}\sim\Delta t^{2}$ for both values of
$\gamma$. When $\gamma=1+1/\sqrt{2}$, the absolute error is significantly
larger in comparison with the solution obtained for $\gamma=1-1/\sqrt{2}$.
However, for $\gamma=1-1/\sqrt{2}$ the numerical solution may introduce small
oscillations about the correct equilibrium. Fortunately, since the method is
asymptotically stable, any initial oscillations are damped to the correct
equilibrium. This is shown in the left panel of Figure 1 where the gas
velocity initially undershoots the analytical solution (see Appendix D for
additional comments on the value of $\gamma$).

When $\Delta t\gg t_{\rm damp}$, from Equation (31) we can measure the
characteristic asymptotic convergence rate of the numerical method. For large
time steps, the error displayed in Fig. 1 converges to a constant value.
Therefore, the numerical solutions approaches the equilibrium as $\sim
1/\Delta t$. This convergence rate is in agreement with the asymptotic
behaviour derived in Appendix A.

Figure 2: Relative velocity between the gas and dust species for the damping
these with a constant external force (see Section 3.2). The numerical solution
is shown with orange and blue circles. The dashed line corresponds to the
analytical solution.

###  3.2 Damping with an external force

In this section, we extend the test of Section 3.1 by including an external
force. We illustrate the convergence properties of our method when the time-
explicit and time-implicit integrators are coupled together as described in
Section 2.2.

We again consider one gas and two dust species and solve the momentum equation
using Eq. (2.2), neglecting the continuity equation as well as the divergence
of fluxes. Under these assumptions, the momentum equation becomes

$\partial_{t}{\bf u}={\bf M}{\bf u}+{\bf G}\,$ (32)

with ${\bf G}=(G_{0},0,0)$. The analytical solution of Eq. (32) corresponds
to,

${\bf u}_{\rm analytical}={\bf
C}\,\,[t,e^{\lambda_{1}t},e^{\lambda_{2}t},1]^{\rm T}\,,$ (33)

where $\lambda_{1}$ and $\lambda_{2}$ are the non-zero eigenvalues of ${\bf
M}$ and ${\bf C}$ is a coefficient matrix that depends on the initial
condition, external force, and eigenvalues of ${\bf M}$. We set as initial
condition $\rho_{\rm g}=1.0$, $\rho_{\rm d1}=0.1$, $\rho_{\rm d2}=0.1$ for the
density and $v_{\rm g}=2.0$, $v_{\rm d1}=0.1$, $v_{\rm d2}=-0.5$ for the
velocity. The external force is set to $G_{0}=1$ whereas the drag force
coefficients are set to $\alpha_{\rm d,1}=1.0$, $\alpha_{\rm d,2}=0.75$,
respectively. With these values we obtain $\lambda_{1}=-1.125$ and
$\lambda_{2}=-0.8$ and the coefficient matrix ${\bf C}$ is

$\displaystyle{\bf C}$
$\displaystyle=\begin{pmatrix}0.833333&0.145014&0.059615&1.79537\\\
0.083333&-0.116011&0.029808&0.096204\\\
0.083333&-0.029003&-0.089423&0.068426\\\ \end{pmatrix}\,.$ (34)

In Figure 2 we show the relative velocity between the gas and the dust
species. The dashed line corresponds to the analytical solution, the agreement
between the analytical and numerical solutions is excellent. In the limit of
$\Delta t\rightarrow\infty$, the momentum of the $j$-th species (see Eq, (33))
converges to

$\displaystyle u_{j}$ $\displaystyle=P_{j0}\left(\left({\bf P}^{-1}{\bf
G}\right)_{0}\,t+\left({\bf P}^{-1}{\bf
u}\right)_{0}\right)-\sum^{N-1}_{k=1}P_{jk}\frac{\left({\bf P}^{-1}{\bf
G}\right)_{k}}{\lambda_{k}}\,$ (35)

where we assume $\lambda_{0}=0$ and $P_{jk}$ denotes the $jk$ entry of the
matrix ${\bf P}$, with ${\bf P}$ and ${\bf P}^{-1}$ the matrices such that
${\bf D}={\bf P}^{-1}{\bf M}{\bf P}$ is diagonal, that is, $\lambda_{k}={\bf
D}_{kk}$.

After dividing Eq. (35) by the corresponding density we can obtain the
relative velocity between species. First, note that
$P_{j0}/\rho_{j}-P_{i0}/\rho_{i}=0$, for any $i$ and $j$ components. This
follows from the properties of the matrix ${\bf M}$ described below Eq. (10)
and using that the first (zero) column of the matrix ${\bf P}$ is the
eigenvector $\tilde{\bf 1}^{T}$. Therefore, the relative velocity between two
species is obtained as

$\lim_{\Delta
t\rightarrow\infty}v_{j}-v_{i}=\frac{1}{\rho_{i}}\sum^{N-1}_{k=1}P_{ik}\frac{\left({\bf
P}^{-1}{\bf
G}\right)_{k}}{\lambda_{k}}-\frac{1}{\rho_{j}}\sum^{N-1}_{k=1}P_{jk}\frac{\left({\bf
P}^{-1}{\bf G}\right)_{k}}{\lambda_{k}}$ (36)

In this case we obtain $\lim_{\Delta t\rightarrow\infty}|v_{\rm g}-v_{\rm
d1}|=0.8333$ and $\lim_{\Delta t\rightarrow\infty}|v_{\rm g}-v_{\rm
d2}|=1.111$, which are the asymptotic values displayed in Figure 2. We
conclude that the numerical method converges to the correct equilibrium for
the case with external constant force. We study this convergence property
further in Appendix C.

Figure 3: Numerical (open circles) and analytical (solid lines) solution of
the five-fluid test described in Sec 3.3 in (Benítez-Llambay et al., 2019). We
plot the time evolution of the normalized density. The blue circles correspond
to the gas, while the other color corresponds to the dust species. All
solutions were obtained at $x=0$. The normalized density is defined as
$\delta\hat{\rho}/(A\rho^{0})$, with $A=10^{-4}$. We include in the plot the
$L_{1}$ error composed of the sum errors of density and velocity (not shown)
components of the five fluids. Individual errors are calculated over the whole
grid after $t=2.0$. Figure 4: Numerical (open circles) and analytical (solid
lines) solution of the five-fluid test described in Sec 3.2 in (Benítez-
Llambay et al., 2019). The numerical solution was obtained at time $t=500$,
starting from the initial jump condition. The plot shows the density of the
gas (blue) and dust (orange,green,red) species, respectively. In the large
panel, to allow the quality of the asymptotic behavior far away from the shock
to be assessed, the sampling rate for the open circles was reduced to 1:4 of
the original data. We additionally plot as an inset a zoomed region about the
shock showing the full sampling, where the unﬁlled circles correspond to the
actual grid points. The code resolves the shock within a single cell, even
when an increasing number of ﬂuids is considered. The overall agreement
between the numerical and analytical solutions is excellent.

###  3.3 Damping of a soundwave

In this section we reproduce the test described in Section 3.2 of Benítez-
Llambay et al. (2019). For this purpose, we have developed a simple code that
solves the Equations (1)-(4) in a periodic domain with the second-order Runge-
Kutta method of Eq. (2.2). We utilize a uniform mesh in the $x$-coordinate
with $N_{x}$ grid cells. After obtaining the numerical solution, we compare it
with the analytical one and show that our method described in Section 2.2
achieves the expected second-order accuracy.

In Figure 3, we show the analytic (solid lines) and numerical (open circles)
solutions for the normalized density, defined as $\delta\hat{\rho}/(\rho A)$.
The test includes five species represented by different colors. Numerical and
analytical solutions are in very good agreement. We study the convergence of
the numerical scheme computing the numerical error as

${\rm e}_{1}=\frac{1}{N_{x}}\sum_{k=1}^{5}\sum_{i=1}^{2}\left\lvert\delta
f_{i,{\rm numeric},k}-\delta f_{i,{\rm analytic},k}\right\rvert\,,$ (37)

with $k$ the index of the fluid and $i$ the index of the density and velocity
field, respectively. In the bottom panel of Figure 3 we show that the error
effectively scales as ${\rm e}_{1}\sim N_{x}^{-2}$.

###  3.4 Shock test

In this section we test our numerical method by reproducing the shock test
from Lehmann & Wardle (2018), extended to multiple fluids in (see Section 3.3
Benítez-Llambay et al., 2019). We solve Eqs. (1)-(4) in an evenly spaced mesh
along the $x$-coordinate with 400 grid cells. In this shock test, we assume
that Eq. (5) is constant, that is, $\rho_{\rm g}\alpha_{{\rm g},i}=\rho_{{\rm
d},i}\alpha_{{\rm d},i}\equiv K_{i}$. Since the drag-force solver has been
written in terms of the parameter $\epsilon_{i}$, the drag-force coefficient
$\alpha_{{\rm d},i}$ must be replaced by $\alpha_{{\rm d},i}=K_{i}/\rho_{{\rm
d},i}$ in the matrix ${\bf M}$. The details of the test together with initial
and boundary conditions are in section 3.3.2 from Benítez-Llambay et al.
(2019).

In Figure 4, we plot the exact (solid line) and numerical (open circles)
solutions for the density of the four fluid shock test. To compare the
numerical solution with the exact one, we shif the exact solution 2.97 in a
$x$-positive direction until the shock position as explained in (Benítez-
Llambay et al., 2019). The unfilled circles are a sub-sampling (1:4) of the
grid points. Inside, we plot a zoomed-in region containing the discontinuity
to show the quality of the numerical solution across the shock in the actual
grid (all the points are shown). The Riemann solver captures the shock within
one cell. For the dust fluid we utilize the solver described in Huang & Bai
(2022).

###  3.5 Jeans Instability

In order to expand our benchmarking to include codes with self-gravity, we
simulate a multifluid Jeans instability in an infinite uniform medium. To our
knowledge this is the first implementation of such a test, and so for
completeness we present the analytic solution in detail below.

We add the gravitational potential term into Equations (3) and (4), now
written as

$\displaystyle\rho_{\rm g}\frac{D{\bf v}_{\rm g}}{Dt}$
$\displaystyle=\rho_{\rm g}\nabla\Phi-\nabla P-\sum^{N-1}_{i=1}\rho_{{\rm
g}}\alpha_{{\rm g},i}\left({\bf v}_{\rm g}-{\bf v}_{{\rm d},i}\right),$ (38)
$\displaystyle\rho_{{\rm d},i}\frac{D{\bf v}_{{\rm d},i}}{Dt}$
$\displaystyle=\rho_{{\rm d},i}\nabla\Phi-\rho_{{\rm d},i}\alpha_{{\rm
d},i}\left({\bf v}_{{\rm d},i}-{\bf v}_{\rm g}\right),$ (39)

with $D/Dt=\partial/\partial_{t}+{\bf v}\cdot\nabla$ the material derivate.
The Poisson equation that governs the gravitational potential of the system is

$\nabla^{2}\Phi=4\pi G\left(\rho_{\rm g}+\sum_{i=1}^{N-1}\rho_{{\rm
d},i}\right)\,.$ (40)

In order to study the Jeans instability, we linearize Equations (1)-(2), and
(38)-(40) with respect to an equilibrium solution obtained after neglecting
the contribution of the background gravitational field in the Poisson equation
(Jeans Swindle) (Binney & Tremaine, 2008). The equilibrium gas density is
constant and denoted as $\rho^{0}_{\rm g}$. The equilibrium dust density of
the $i$-th species is also a constant, equal to $\rho^{0}_{{\rm
d},i}=\epsilon_{i}^{0}\rho_{\rm g}^{0}$. Finally, the equilibrium velocities
are set to zero. We define small density perturbations $\delta\rho_{\rm g}$
and $\delta\rho_{{\rm d},i}$, whereas the velocity perturbations are
${\delta\bf v}_{\rm g}=\delta v_{\rm g}\hat{x}$, $\delta{\bf v}_{{\rm
d},i}=\delta v_{{\rm d},i}\hat{x}$. We set $\rho^{0}_{g}=1.0$ and the total
dust-to-gas mass ratio $\epsilon^{0}=\sum^{N-1}_{i=1}\epsilon^{0}_{i}=0.01$,
with $\epsilon^{0}_{i}=\epsilon^{0}/N$. We consider $N=128$ fluids with a
uniform distribution of stopping times between $t_{\rm s}=10^{-4}-10$

The linearized continuity, momentum, and Poisson equation are

$\displaystyle\frac{\partial\delta\rho_{\rm g}}{\partial t}+\rho_{\rm
g}^{0}\frac{\partial\delta v_{\rm g}}{\partial x}$ $\displaystyle=0\,,$ (41)
$\displaystyle\frac{\partial\delta\rho_{{\rm d},i}}{\partial t}+\rho_{{\rm
d},i}^{0}\frac{\partial\delta v_{{\rm d},i}}{\partial x}$ $\displaystyle=0\,,$
(43) $\displaystyle\frac{\partial\delta v_{\rm g}}{\partial
t}+\sum_{m=1}^{N-1}\frac{\epsilon_{m}^{0}}{t_{{\rm s},m}}(\delta v_{\rm
g}-\delta v_{{\rm d},m})+\frac{c_{\rm
s}^{2}}{\rho_{g}^{0}}\frac{\partial\delta\rho_{\rm g}}{\partial x}$
$\displaystyle=-\frac{\partial\delta\phi}{\partial x}\,,$ (44)
$\displaystyle\frac{\partial\delta v_{{\rm d},i}}{\partial t}+\frac{1}{t_{{\rm
s},i}}(\delta v_{{\rm d},i}-\delta v_{\rm g})$
$\displaystyle=-\frac{\partial\delta\phi}{\partial x}\,,$ (45) $\displaystyle
4\pi G(\delta\rho_{\rm g}+\sum_{m=1}^{N-1}\delta\rho_{{\rm d},m})$
$\displaystyle=\nabla^{2}\delta\phi\,.$ (46)

For $N$ fluids, these $2N+1$ equations are coupled, due to the gravitational
field and the drag-force.

After assuming plane-wave perturbations $\delta f=\delta\hat{f}e^{i{\bf
k}\cdot{\bf x}-\omega t}$ with ${\bf k}$ a real wave-number, we obtain the
eigenvalues and eigenvectors of the system. We normalize the set of linearized
equations defining

$\displaystyle\tilde{\omega}=\omega/(c_{\rm s}k_{\rm
J})\,,\quad\quad\tilde{k}=k/k_{\rm J}$ (47) $\displaystyle\tilde{\rho}_{\rm
g}=\delta\rho_{\rm g}/\rho_{\rm g}^{0}\,,\quad\quad\tilde{\rho}_{{\rm
d},i}=\delta\rho_{{\rm d},i}/\rho_{\rm g}^{0}$ (48)
$\displaystyle\tilde{v}_{\rm g}=\delta v_{\rm g}/c_{\rm
s}\,,\quad\quad\tilde{v}_{{\rm d},i}=\delta v_{{\rm d},i}/c_{\rm s}$ (49)
$\displaystyle\beta_{i}=1/(t_{{\rm s},i}c_{\rm s}k_{\rm J})$ (50)

where

$\displaystyle k_{\rm J}=\sqrt{4\pi G\rho_{\rm g}^{0}(1+\epsilon^{0})}/c_{\rm
s},\,\rm{and}\quad\epsilon^{0}=\sum_{m=1}^{N-1}\epsilon_{m}^{0}.$ (51)

After normalization, (41)$-$(46) become

$\displaystyle-\tilde{\omega}\tilde{\rho}_{\rm g}+i\tilde{k}\tilde{v}_{\rm
g}=0\,,$ (52) $\displaystyle-\tilde{\omega}\tilde{\rho}_{{\rm
d},i}+i\epsilon^{0}_{i}\tilde{k}\tilde{v}_{{\rm d},i}=0\,,$ (53)
$\displaystyle-\tilde{\omega}\tilde{v}_{\rm
g}+\sum_{m=1}^{N-1}\epsilon_{m}^{0}\beta_{m}(\tilde{v}_{\rm g}-\tilde{v}_{{\rm
d},i})+i\tilde{\rho}_{\rm
g}\tilde{k}=-i\tilde{k}\frac{\delta\hat{\phi}}{c_{\rm s}^{2}}\,,$ (54)
$\displaystyle-\tilde{\omega}\tilde{v}_{{\rm d},i}+\beta_{i}(\tilde{v}_{{\rm
d},i}-\tilde{v}_{\rm g})=-i\tilde{k}\frac{\delta\hat{\phi}}{c_{\rm s}^{2}}\,,$
(55) $\displaystyle\frac{\delta\hat{\phi}}{c_{\rm
s}^{2}}=-\frac{1}{\tilde{k}^{2}(1+\epsilon^{0})}\left(\tilde{\rho}_{\rm
g}+\sum_{m}^{N-1}\tilde{\rho}_{{\rm d},m}\right)\,,$ (56)

and can be reduced to a $2N$ equation system, written as

$\displaystyle i\tilde{k}\tilde{v}_{\rm g}$
$\displaystyle=\tilde{\omega}\tilde{\rho}_{\rm g}\,,$ (57) $\displaystyle
i\tilde{k}\epsilon^{0}_{i}\tilde{v}_{{\rm d},i}$
$\displaystyle=\tilde{\omega}\tilde{\rho}_{{\rm d},i}\,,$ (58)
$\displaystyle\tilde{\omega}\tilde{v}_{\rm g}$
$\displaystyle=\left(i\tilde{k}-\frac{i}{\tilde{k}(1+\epsilon^{0})}\right)\tilde{\rho}_{\rm
g}-\frac{i}{\tilde{k}(1+\epsilon^{0})}\sum_{m=1}^{N-1}\tilde{\rho}_{{\rm
d},m}$ (59)
$\displaystyle+\sum_{m=1}^{N-1}\epsilon_{m}^{0}\beta_{m}\tilde{v}_{\rm
g}-\sum_{m=1}^{N-1}\epsilon_{m}^{0}\beta_{m}\tilde{v}_{{\rm d},m}\,$
$\displaystyle\tilde{\omega}\tilde{v}_{{\rm d},i}$
$\displaystyle=-\frac{i}{\tilde{k}(1+\epsilon^{0})}\tilde{\rho}_{\rm
g}-\frac{i}{\tilde{k}(1+\epsilon^{0})}\sum_{m=1}^{N-1}\tilde{\rho}_{{\rm
d},m}$ (60) $\displaystyle-\beta_{i}\tilde{v}_{\rm g}+\beta_{i}\tilde{v}_{{\rm
d},i}\,.$

The equivalent matrix equation is ${\bf A}{\bf V}=\tilde{\omega}{\bf V}$,
${\bf V}=(\tilde{\rho}_{\rm g},\tilde{\rho}_{{\rm
d},i},\dotsc,\tilde{\rho}_{{\rm d},N-1},\tilde{v}_{\rm g},\tilde{v}_{{\rm
d},i},\dotsc,\tilde{v}_{{\rm d},N-1})$. The eigenvalues and eigenvectors of
the problem, in our tests, have been obtained with the function eig() of the
Numpy Python package.

Figure 5 displays the grow rate of the dusty Jeans instability
$\sigma=-\min\\{\Re(\tilde{\omega})\\}\sqrt{1+\epsilon^{0}}$ as a function of
$k/k_{\rm J}$. Note that when $\epsilon_{0}=0$, $\sigma$ reduces to the
normalized growth rate of the gas-only problem. Thus, the inclusion of dust
enables higher growth rates as shown in Figure 5 (see solid line vs dashed
line (dust-free case). Comparing the gas only and dusty cases illustrates that
while both feature a growth rate that is larger for lower wave numbers, the
critical $k/k_{\rm J}$ for stability is shifted to larger values as dust
density increases. In Figure 5 we also the numerical solution obtained with
MDIRK described in Section 2.2. We calculate the growth rate varying the
number of cells from $N_{\rm cells}=16$ to $N_{\rm cells}=128$. The growth
rate converges to the analytic solution for all unstable modes considered in
this work222The largest wavenumber considered in this work corresponds to the
cutoff mode of the dust-free case, e.i., $k/k_{\rm
J}=1/\sqrt{1+\epsilon^{0}}$.

Figure 5: Test resolution convergence for Jeans instability including drag
forces. The dashed line denotes the analytic solution for only gas. The solid
line shows the analytic solution for 128 fluids including gas a dust species
with a linear distribution from $t_{\rm s}=10^{-4}$ to $10$ and uniform
density distribution with a total dust-to-gas ratio $\epsilon^{0}=0.01$. The
open circles indicate that the numerical solutions for different numbers of
cells converge to the solid line. The vertical dotted line indicates $k/k_{\rm
J}=1/\sqrt{1+\epsilon^{0}}$.

###  3.6 Axisymmetric steady-state drift solution

To demonstrate the scaling of our method with the number of species, we
recover the steady-state gas and dust drift solution in a shearing box, as in
section 3.5.2 of Benítez-Llambay et al. (2019). This test consists of a 1D
shearing box along the x direction, with size $L=1$ and with constant
frequency $\Omega_{0}=1$. The shear parameter is set to $q=3/2$. We solve the
continuity and $x$-$y$ momentum equations for each fluid. Gas and dust
momentum equations include the centrifugal force (${\bf f}_{\rm ce}=\rho
2q\Omega^{2}_{0}x\hat{x}$) and the Coriolis force ${\bf
f}_{c}=-2\rho\Omega_{0}\hat{z}\times(v_{x},v_{y})$. In addition, we add an
external constant force on the gas fluid that mimics the global pressure
gradient in a protoplanetary disk, ${\bf f}_{\rm g}=\chi_{0}\Omega_{0}\hat{x}$
with $\chi_{0}=0.005$. We add these external forces to the time-explicit
hydrodynamic solver in an unsplit fashion (see e.g. Gressel & Ziegler, 2007).
We apply periodic boundary conditions in $\rho,\rho v_{x}$. Whereas for
$v_{y}$ we add the relative shear velocity after applying periodic boundary
conditions, that is $\rho v_{y}(x)=\rho v_{y}(x\mp L)\pm q\Omega_{0}L\rho$
(Hawley et al., 1995). The initial condition is set to the exact equilibrium
following Eqs. (81-84) of (Benítez-Llambay et al., 2019). We integrate the
system until $t=20$ and vary the number of grid points in $x$ from $N_{x}=32$
to $N_{x}=1024$ . In all cases, our numerical method conserves the initial
equilibrium to machine precision.

In Figure 6, we show the normalized CPU time obtained for this test, where we
vary the number of dust species from 1 to 63. As expected, the computational
cost is proportional to the number of fluids (complexity $\mathcal{O}(N)$). We
moreover compare our result with the second-order multifluid method RK2Im of
(Huang & Bai, 2022), which is found to be close to $\mathcal{O}(N^{3})$) for
large number of fluids. It is interesting to note that the complexity of Huang
& Bai (2022) method can be reduced to $\mathcal{O}(N)$ by utilizing our
analytical solution in Appendix A.

##  4 Discussion and conclusions

In this work, we have presented MDIRK, a multifluid asymptotically stable
numerical method to solve the momentum exchange between gas and dust species.
We have shown that a hydrodynamical solver based on a Runge-Kutta implicit-
explicit (IMEX) algorithm converges to the correct equilibrium solution with
second-order accuracy. The time-implicit part of the IMEX utilizes a
diagonally implicit Runge–Kutta method (DIRK) (Alexander, 1977). The major
advantage of a diagonally implicit method for the drag force is that it admits
an analytical solution that reduces the complexity of the problem from
$\mathcal{O}(N^{3})$ (needed to invert the matrix numerically) to
$\mathcal{O}(N)$ (needed to evaluate the solution). Our second-order implicit-
explicit solver can be easily implemented in state-of-the-art hydrodynamics
codes (e.g., PLUTO (Mignone et al., 2007), ATHENA++ (Stone et al., 2020),
among others) and it allows for efficient parallelization by fluids through
the utilization of of weighted momenta and density (c.f. Krapp & Benítez-
Llambay, 2020). Extension to different second-order methods can be also
obtained utilizing the analytical solution of Appendix A.

Our results suggest that IMEX Runge-Kutta methods with diagonally implicit
drag-force schemes are an interesting alternative for developing high-order
multifluid codes. Examples of third-order and four-order IMEX methods can be
found in (Ascher et al., 1997; Pareschi & Russo, 2005). Note however that a
high-order integration will require monotonic preserving spatial
reconstruction techniques with adequate accuracy (Mignone, 2014; Felker &
Stone, 2018). In addition, we stress that only a subset of high-order Runge-
Kutta methods are consistent with a total-variation-diminishing property (TVD)
(Shu & Osher, 1988; Gottlieb & Shu, 1998) and therefore special care must be
taken with high-order methods to prevent the presence of spurious oscillation.

While we focus on 1D tests, our method can be utilized for 2D and 3D with
static and adaptive mesh refinement with the appropriate extension of the
hydrodynamics (Lebreuilly et al., 2019; Huang & Bai, 2022). Also, in this work
we did not take into account mass transfer between species, however, future
extensions could explore this feature by combining our method with a dust-
growth solution based on discretized dust bins method as described in Lombart
& Laibe (2021). Another important future extension of our method would be the
inclusion of energy exchange between dust and gas (Muley et al., 2023), since
collisional cooling could play an important role in the thermodynamics of
protoplanetary disks. In radiation hydrodynamics, multi-fluid dynamics can
also impact the opacity distribution and observational signatures of
accreting, embedded planets (Krapp et al., 2022). We emphasize that the
complexity $\mathcal{O}(N)$ combined with the stability properties will
expedite the study of particle-size distribution, especially in problems where
hundreds and/or thousands of bins are required to properly model the dust
component (see e.g., Krapp et al., 2019; Zhu & Yang, 2021).

Figure 6: Normalized CPU time vs the number of fluids for the test described
in Section 3.6. As expected, the complexity of MDIRK is $\mathcal{O}(N)$. For
comparison, we also include the results of the multifluid second-order method
(RK2Im) from Huang & Bai (2022) which has complexity $\mathcal{O}(N^{3})$.

L. K. and K. M. K acknowledge support from the Heising-Simons 51 Pegasi b
postdoctoral fellowship and TCAN grant 80NSSC19K0639. P. B. L. acknowledges
support from ANID, QUIMAL fund ASTRO21-0039 and FONDECYT project 1231205. J.G.
acknowledge support by ANID, – Millennium Science Initiative Program –
NCN19_171, CCJRF2105, FONDECYT Regular grant 1210425 and CATA-BASAL project
FB210003.

All codes and data utilized in this work are available in
bitbucket.org/secondorderriemannsolver/somrs.

## References

  * Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756, doi: 10.1143/PTP.56.1756
  * Alexander (1977) Alexander, R. 1977, SIAM Journal on Numerical Analysis, 14, 1006, doi: 10.1137/0714068
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Ascher et al. (1997) Ascher, U. M., Ruuth, S. J., & Spiteri, R. J. 1997, Applied Numerical Mathematics, 25, 151, doi: https://doi.org/10.1016/S0168-9274(97)00056-1
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
  * Balsara et al. (2009) Balsara, D. S., Tilley, D. A., Rettig, T., & Brittain, S. D. 2009, Monthly Notices of the Royal Astronomical Society, 397, 24, doi: 10.1111/j.1365-2966.2009.14606.x
  * Benítez-Llambay et al. (2019) Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25, doi: 10.3847/1538-4365/ab0a0e
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition, rev - revised, 2 edn. (Princeton University Press).  http://www.jstor.org/stable/j.ctvc778ff
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45, doi: 10.3847/2041-8213/aaf743
  * Booth et al. (2015) Booth, R. A., Sijacki, D., & Clarke, C. J. 2015, Monthly Notices of the Royal Astronomical Society, 452, 3932, doi: 10.1093/mnras/stv1486
  * Chen & Lin (2018) Chen, J.-W., & Lin, M.-K. 2018, Monthly Notices of the Royal Astronomical Society, 478, 2737, doi: 10.1093/mnras/sty1166
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Epstein (1924) Epstein, P. S. 1924, Phys. Rev., 23, 710, doi: 10.1103/PhysRev.23.710
  * Felker & Stone (2018) Felker, K. G., & Stone, J. M. 2018, Journal of Computational Physics, 375, 1365, doi: 10.1016/j.jcp.2018.08.025
  * Fung & Muley (2019) Fung, J., & Muley, D. 2019, ApJS, 244, 42, doi: 10.3847/1538-4365/ab45f3
  * Gottlieb & Shu (1998) Gottlieb, S., & Shu, C.-W. 1998, Mathematics of Computation, 67, 73.  http://www.jstor.org/stable/2584973
  * Gressel & Ziegler (2007) Gressel, O., & Ziegler, U. 2007, Computer Physics Communications, 176, 652, doi: 10.1016/j.cpc.2007.01.010
  * Hanasz et al. (2010) Hanasz, M., Kowalik, K., Wóltański, D., Pawłaszek, R., & Kornet, K. 2010, in EAS Publications Series, Vol. 42, EAS Publications Series, ed. K. Gożdziewski, A. Niedzielski, & J. Schneider, 281–285, doi: 10.1051/eas/1042030
  * Harten et al. (1983) Harten, A., Lax, P. D., & Leer, B. v. 1983, SIAM Review, 25, 35, doi: 10.1137/1025002
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742, doi: 10.1086/175311
  * Hopkins et al. (2020) Hopkins, P. F., Squire, J., & Seligman, D. 2020, MNRAS, 496, 2123, doi: 10.1093/mnras/staa1046
  * Huang & Bai (2022) Huang, P., & Bai, X.-N. 2022, ApJS, 262, 11, doi: 10.3847/1538-4365/ac76cb
  * Johansen et al. (2004) Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361, doi: 10.1051/0004-6361:20034417
  * Krapp & Benítez-Llambay (2020) Krapp, L., & Benítez-Llambay, P. 2020, Research Notes of the American Astronomical Society, 4, 198, doi: 10.3847/2515-5172/abc7be
  * Krapp et al. (2019) Krapp, L., Benítez-Llambay, P., Gressel, O., & Pessah, M. E. 2019, ApJ, 878, L30, doi: 10.3847/2041-8213/ab2596
  * Krapp et al. (2022) Krapp, L., Kratter, K. M., & Youdin, A. N. 2022, ApJ, 928, 156, doi: 10.3847/1538-4357/ac5899
  * Laibe & Price (2012) Laibe, G., & Price, D. J. 2012, MNRAS, 420, 2345, doi: 10.1111/j.1365-2966.2011.20202.x
  * Laibe & Price (2014) —. 2014, MNRAS, 444, 1940, doi: 10.1093/mnras/stu1367
  * Lebreuilly et al. (2019) Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96, doi: 10.1051/0004-6361/201834147
  * Lehmann & Wardle (2018) Lehmann, A., & Wardle, M. 2018, MNRAS, 476, 3185, doi: 10.1093/mnras/sty450
  * LeVeque (1992) LeVeque, R. J. 1992, Numerical methods for conservation laws (2. ed.)., Lectures in mathematics (Birkhäuser), 1–214
  * Leveque (2004) Leveque, R. J. 2004, Journal of Hyperbolic Differential Equations, 01, 315, doi: 10.1142/S0219891604000135
  * Lombart & Laibe (2021) Lombart, M., & Laibe, G. 2021, MNRAS, 501, 4298, doi: 10.1093/mnras/staa3682
  * Mignone (2014) Mignone, A. 2014, Journal of Computational Physics, 270, 784, doi: 10.1016/j.jcp.2014.04.001
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228, doi: 10.1086/513316
  * Mignone et al. (2019) Mignone, A., Flock, M., & Vaidya, B. 2019, ApJS, 244, 38, doi: 10.3847/1538-4365/ab4356
  * Miniati (2010) Miniati, F. 2010, Journal of Computational Physics, 229, 3916, doi: https://doi.org/10.1016/j.jcp.2010.01.034
  * Muley et al. (2023) Muley, D., Melon Fuksman, J. D., & Klahr, H. 2023, arXiv e-prints, arXiv:2308.03504, doi: 10.48550/arXiv.2308.03504
  * Paardekooper & Mellema (2006) Paardekooper, S. J., & Mellema, G. 2006, A&A, 453, 1129, doi: 10.1051/0004-6361:20054449
  * Pareschi & Russo (2005) Pareschi, L., & Russo, G. 2005, Journal of scientific computing, 25, 129–155, doi: 10.1007/s10915-004-4636-4
  * Pollack et al. (1994) Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615, doi: 10.1086/173677
  * Price et al. (2018) Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, Publications of the Astronomical Society of Australia, 35, e031, doi: 10.1017/pasa.2018.25
  * Rosenbrock (1963) Rosenbrock, H. H. 1963, The Computer Journal, 5, 329, doi: 10.1093/comjnl/5.4.329
  * Rusanov (1970) Rusanov, V. 1970, Journal of Computational Physics, 5, 507, doi: https://doi.org/10.1016/0021-9991(70)90077-X
  * Shu & Osher (1988) Shu, C.-W., & Osher, S. 1988, Journal of Computational Physics, 77, 439, doi: https://doi.org/10.1016/0021-9991(88)90177-5
  * Squire & Hopkins (2018) Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011, doi: 10.1093/mnras/sty854
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Stoyanovskaya et al. (2018) Stoyanovskaya, O. P., Glushko, T. A., Snytnikov, N. V., & Snytnikov, V. N. 2018, Astronomy and Computing, 25, 25, doi: 10.1016/j.ascom.2018.08.004
  * Strang (1968) Strang, G. 1968, SIAM Journal on Numerical Analysis, 5, 506, doi: 10.1137/0705041
  * Toro (2009) Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd edn. (Berlin, Heidelberg: Springer-Verlag), doi: 10.1007/b79761
  * van Leer (1977) van Leer, B. 1977, Journal of Computational Physics, 23, 276, doi: 10.1016/0021-9991(77)90095-X
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Yang & Johansen (2016) Yang, C.-C., & Johansen, A. 2016, ApJS, 224, 39, doi: 10.3847/0067-0049/224/2/39
  * Youdin & Johansen (2007) Youdin, A., & Johansen, A. 2007, ApJ, 662, 613, doi: 10.1086/516729
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459, doi: 10.1086/426895
  * Zhu & Yang (2021) Zhu, Z., & Yang, C.-C. 2021, MNRAS, 501, 467, doi: 10.1093/mnras/staa3628

##  Appendix A Asymptotic stability

In this section, we show that MDIRK is asymptotically stable. As in Section
2.1, we will only focus on the system of momentum equations. Moreover, we only
consider the implicit part of the method. In this case, MDIRK reduces to a
two-stage diagonally implicit method of the form

$\displaystyle\left({\bf I}-\gamma\Delta t{\bf M}\right){\bf k}_{1}$
$\displaystyle={\bf M}{\bf u}^{n}\,,$ (A1) $\displaystyle\left({\bf
I}-\gamma\Delta t{\bf M}\right){\bf k}_{2}$ $\displaystyle={\bf M}\left({\bf
u}^{n}+\beta\Delta t{\bf k}_{1}\right)$ (A2) $\displaystyle{\bf u}^{n+1}$
$\displaystyle={\bf u}^{n}+b_{1}\Delta t{\bf k}_{1}+b_{2}\Delta t{\bf
k}_{2}\,,$ (A3)

with $\beta=1-\gamma$, $b_{1}=1-\gamma$ and $b_{2}=\gamma$. The analytical
solution of Eq. (A3) can be obtained after solving Eqs. (A1)-(A2) for ${\bf
k}_{1}$ and ${\bf k}_{2}$ utilizing the following summations

$\displaystyle\mathcal{A}$ $\displaystyle=\sum^{N-1}_{k=1}\frac{\alpha_{{\rm
d},k}u_{{\rm d},k}}{(1+\gamma\Delta t\alpha_{{\rm
d},k})}\,,\quad\mathcal{B}=\sum^{N-1}_{k=1}\frac{\epsilon_{k}\alpha_{{\rm
d},k}}{(1+\gamma\Delta t\alpha_{{\rm d},k})}\,,$ $\displaystyle\mathcal{C}$
$\displaystyle=\sum^{N-1}_{k=1}\frac{\alpha^{2}_{{\rm d},k}u_{{\rm
d},k}}{(1+\gamma\Delta t\alpha_{{\rm
d},k})^{2}}\,,\quad\mathcal{D}=\sum^{N-1}_{k=1}\frac{\epsilon_{k}\alpha^{2}_{{\rm
d},k}}{(1+\gamma\Delta t\alpha_{{\rm d},k})^{2}}\,,$ (A4)

with

$\displaystyle k_{1{\rm g}}$ $\displaystyle=\frac{\mathcal{A}-u_{\rm
g}^{n}\mathcal{B}}{1+\gamma\Delta t\mathcal{B}}$ (A5) $\displaystyle k_{2{\rm
g}}$ $\displaystyle=k_{1{\rm g}}\left(\frac{1+(\gamma-\beta)\Delta
t\mathcal{B}+\gamma\beta\Delta t^{2}\mathcal{D}}{1+\gamma\Delta
t\mathcal{B}}\right)+\frac{\beta\Delta t}{(1+\gamma\Delta
t\mathcal{B})}\left(u_{\rm g}^{n}\mathcal{D}-\mathcal{C}\right)$ (A6)

and

$\displaystyle k_{1{\rm d},i}$ $\displaystyle=\frac{\alpha_{i}}{1+\gamma\Delta
t\alpha_{i}}\left(\epsilon_{i}u_{\rm g}-u_{{\rm d},i}+\gamma\Delta
t\epsilon_{i}k_{1{\rm g}}\right)$ (A7) $\displaystyle k_{2{\rm d},i}$
$\displaystyle=\frac{\left(\alpha_{i}+(\gamma-\beta)\Delta
t\alpha^{2}_{i}\right)}{\left(1+\gamma\Delta
t\alpha_{i}\right)^{2}}\left(\epsilon_{i}u_{\rm g}-u_{{\rm
d},i}\right)+\frac{\alpha_{i}\beta\Delta t\epsilon_{i}}{\left(1+\gamma\Delta
t\alpha_{i}\right)^{2}}k_{1{\rm g}}+\frac{\alpha_{i}\gamma\Delta
t\epsilon_{i}}{\left(1+\gamma\Delta t\alpha_{i}\right)}k_{2{\rm g}}\,.$ (A8)

MDIRK does not utilize the full analytical solutions given in Eqs. (A5)-(A7)
since the calculation of ${\bf k}_{1}$ and ${\bf k}_{2}$ are combined with the
explicit updates (see Section 2.2). However, this analytical solution could be
implemented in the method of (Huang & Bai, 2022) reducing the from
$\mathcal{O}(N^{3})$ complexity to $\mathcal{O}(N)$ without shifting to a
mixed implicit explicit scheme. Another alternative could be a method based on
Strang-Splitting (Strang, 1968), however, we anticipate that strang-splitting
may not converge to the correct equilibrium.

####  A.0.1 Asymptotic stability

We will show that this solution is asymptotically stable, that is,

$\lim_{\Delta t\rightarrow\infty}{\bf u}^{n+1}-{\bf c}=0\,,$ (A9)

where ${\bf c}$ is a constant vector. In fact, the terminal velocity of all
fluids converges to a steady state given by the velocity of the center of
mass, $v_{\rm cm}$, defined as

$v_{\rm cm}=\frac{u_{\rm g}+\sum_{k}u_{{\rm d},k}}{\rho_{\rm
g}+\sum_{k}\rho_{{\rm d},k}}=\frac{v_{\rm g}+\epsilon_{k}v_{{\rm
d},k}}{1+\sum_{k}\epsilon_{k}}\,.$ (A10)

To demonstrate the asymptotic convergence we define $u_{\rm d}=\sum_{k}u_{{\rm
d},k}$, $\epsilon=\sum_{k}\epsilon_{k}$, and $X=\left(u_{\rm d}-\epsilon
u_{\rm g}\right)$. With these definitions the summations from Eq. (A) can be
reduced to

$\displaystyle\lim_{\Delta t\rightarrow\infty}\Delta t\mathcal{A}$
$\displaystyle=\gamma^{-1}u_{\rm d}\,,\quad\lim_{\Delta
t\rightarrow\infty}\Delta t\mathcal{B}=\gamma^{-1}\epsilon\,,$
$\displaystyle\lim_{\Delta t\rightarrow\infty}\Delta t^{2}\mathcal{C}$
$\displaystyle=\gamma^{-2}u_{{\rm d}}\,,\quad\lim_{\Delta
t\rightarrow\infty}\Delta t^{2}\mathcal{D}=\gamma^{-2}\epsilon\,.$ (A11)

Using the limits from Eq. (A.0.1) and the coefficients $\beta=b_{1}=1-\gamma$
and $b_{2}=\gamma$ it can be shown that

$\displaystyle\lim_{\Delta t\rightarrow\infty}\Delta tk_{\rm 1g}$
$\displaystyle=\gamma^{-1}\frac{X}{1+\epsilon}$ $\displaystyle\lim_{\Delta
t\rightarrow\infty}\Delta tk_{\rm 2g}$
$\displaystyle=\gamma^{-1}\frac{X}{1+\epsilon}-\gamma^{-2}\frac{(1-\gamma)}{1+\epsilon}X$
$\displaystyle\lim_{\Delta t\rightarrow\infty}\Delta tk_{{\rm 1d},k}$
$\displaystyle=\gamma^{-1}\left(\epsilon_{k}u_{\rm g}-u_{{\rm
d},k}\right)+\gamma^{-1}\frac{\epsilon_{k}}{1+\epsilon}X$
$\displaystyle\lim_{\Delta t\rightarrow\infty}\Delta tk_{{\rm 2d},k}$
$\displaystyle=\gamma^{-2}(2\gamma-1)\left(\epsilon_{k}u_{\rm g}-u_{{\rm
d},k}\right)+\frac{\epsilon_{k}}{1+\epsilon}X\left(\gamma^{-1}-\gamma^{-2}(1-\gamma)\right)\,.$
(A12)

Taking the limit of $\Delta t\rightarrow\infty$ to Eq. (A3) using the limits
from Eq. (A.0.1) we obtain

$\displaystyle\lim_{\Delta t\rightarrow\infty}u^{n+1}_{\rm g}$
$\displaystyle=u_{\rm
g}^{n}+\left((1-\gamma)\gamma^{-1}+1-\gamma^{-1}(1-\gamma)\right)\frac{X}{(1+\epsilon)}$
$\displaystyle=\frac{\left(u_{\rm d}+u_{\rm g}\right)}{(1+\epsilon)}\,,$ (A13)
$\displaystyle\lim_{\Delta t\rightarrow\infty}u^{n+1}_{{\rm d},k}$
$\displaystyle=u_{{\rm
d},k}^{n}\left(1-\gamma^{-1}(1-\gamma)-\gamma^{-1}(2\gamma-1)\right)+u_{\rm
g}\epsilon_{k}\gamma^{-1}\left((1-\gamma)+(2\gamma-1)\right)$
$\displaystyle+\frac{\epsilon_{k}}{1+\epsilon}X\left((1-\gamma)\gamma^{-1}+\gamma\left(\gamma^{-1}-\gamma^{-2}(1-\gamma)\right)\right)$
$\displaystyle=u_{\rm g}\epsilon_{k}+\frac{\epsilon_{k}}{(1+\epsilon)}X$
$\displaystyle=\epsilon_{k}\frac{\left(u_{\rm d}+u_{\rm
g}\right)}{(1+\epsilon)}\,.$ (A14)

In this derivation, it seems that the result does not depend on the value of
$\gamma$. However, this is only due to the particular choice of the
coefficients $\beta$, $b_{1}$ and $b_{2}$. A different set of coefficients
will lead to the condition $\gamma=1\pm 1/\sqrt{2}$, together with constraint
for second-order accuracy.

##  Appendix B Second-order time accuracy

To demonstrate the second-order accuracy in time of the method described in
Section 2.1 we consider an equation of the form

$\frac{{\rm d}{\bf u}}{{\rm d}t}=\bf{G}({\bf u})+{\bf M}{\bf u}$ (B1)

where ${\bf G}({\bf u})$ is a vector of any external forces i.e. those that do
not involve coupling between different fluids, that is ${\bf G}({\bf
u})_{i}=G_{i}(u_{i})$. Forces due to coupling enter only through M. To proceed
we define ${\bf u}(t_{n})={\bf u}^{n}$ and emphasise that the notation ${\bf
q}^{n}$, $\left({\bf q}\right)^{n}$, and $\left[{\bf q}\right]^{n}$ indicates
the function $q(t)$ evaluated when $t=t_{n}$, for any arbitrary function
$q(t)$. The analytical solution of Eq. (B1) is obtained as

${\bf u}(t)={\bf u}^{n}+\int^{t}_{t_{n}}{\bf G}({\bf u})\,{\rm
d}t+\int^{t}_{t_{n}}{\bf M}{\bf u}\,{\rm d}t\,.$ (B2)

We now approximate the solution (B2) to second-order in time. First, the
conserved variable ${\bf u}(t)$ and the force, ${\bf G}({\bf u})$, can be
approximated (to first order in time) as

$\displaystyle{\bf u}(t)$ $\displaystyle\simeq{\bf
u}^{n}+(t-t_{n})\left(\frac{{\rm d}{\bf u}}{{\rm d}t}\right)^{n}$ (B3)
$\displaystyle{\bf G}({\bf u})$ $\displaystyle\simeq{\bf G}({\bf
u}^{n})+(t-t_{n})\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}\left(\frac{{\rm
d}{\bf u}}{{\rm d}t}\right)^{n}$

where $\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}$ is a diagonal operator
with diagonal elements $\left[\nabla_{u}{\bf G}({\bf
u})\right]^{n}_{ii}=\left({\rm d}G_{i}(u_{i})/{\rm d}u_{i}\right)^{n}$. We
emphasize again that $G_{i}$ is only a function of $u_{i}$. With this
consideration, the solution (B2) can be approximated as follows

$\displaystyle{\bf u}(t)$ $\displaystyle\simeq{\bf u}^{n}+{\bf G}({\bf
u}^{n})(t-t_{n})+\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}\left(\frac{{\rm
d}{\bf u}}{{\rm d}t}\right)^{n}\frac{\left(t-t_{n}\right)^{2}}{2}$
$\displaystyle+{\bf M}\left({\bf u}^{n}(t-t_{n})+\left(\frac{{\rm d}{\bf
u}}{{\rm d}t}\right)^{n}\frac{\left(t-t_{n}\right)^{2}}{2}\right)\,.$ (B4)

Eq. (B1) implies that

$\left(\frac{{\rm d}{\bf u}}{{\rm d}t}\right)^{n}={\bf G}({\bf u}^{n})+{\bf
M}{\bf u}^{n}\,,$ (B5)

we finally obtain

$\displaystyle{\bf u}(t)$ $\displaystyle\simeq{\bf u}^{n}+(t-t_{n}){\bf
G}({\bf u}^{n})+(t-t_{n}){\bf M}{\bf
u}^{n}+\frac{(t-t_{n})^{2}}{2}\left(\left[\nabla_{u}{\bf G}({\bf
u})\right]^{n}+{\bf M}\right){\bf G}({\bf
u}^{n})+\frac{(t-t_{n})^{2}}{2}\left(\left[\nabla_{u}{\bf G}({\bf
u})\right]^{n}{\bf M}+{\bf M}^{2}\right){\bf u}^{n}\,.$ (B6)

We now apply the numerical scheme from Section 2.2 with the coefficients
$\beta=b_{1}=1-\gamma$ and $b_{2}=\gamma$. A full update is obtained as
follows333Note that we already shown that this update is equivalent to the
scheme described in Section 2.2.

$\displaystyle\left({\bf I}-\gamma\Delta t{\bf M}\right){\bf k}_{1}$
$\displaystyle={\bf M}\left({\bf u}^{n}+\gamma\Delta t{\bf G}({\bf
u}^{n})\right)$ (B7) $\displaystyle{\bf u}^{*}$ $\displaystyle={\bf
u}^{n}+\gamma\Delta t{\bf k}_{1}+\Delta t\gamma{\bf G}({\bf u}^{n})$ (B8)
$\displaystyle\left({\bf I}-\gamma\Delta t{\bf M}\right){\bf k}_{2}$
$\displaystyle={\bf M}\left({\bf u}^{n}+(1-\gamma)\Delta t{\bf
k}_{1}+\delta\Delta t{\bf G}({\bf u}^{n})+(1-\delta)\Delta t{\bf G}({\bf
u}^{*})\right)$ $\displaystyle{\bf u}^{n+1}$ $\displaystyle={\bf
u}^{n}+(1-\gamma)\Delta t{\bf k}_{1}+\gamma\Delta t{\bf k}_{2}+\delta\Delta
t{\bf G}({\bf u}^{n})+(1-\delta)\Delta t{\bf G}({\bf u}^{*})$ (B9)

Defining

${\bf A}=\left({\bf I}-\gamma\Delta t{\bf M}\right)^{-1}{\bf M}\,,$ (B10)

the solution follows as

$\displaystyle{\bf u}^{n+1}$ $\displaystyle=\left({\bf I}+\Delta t{\bf
A}+\gamma(1-\gamma)\Delta t^{2}{\bf A}^{2}\right){\bf u}^{n}$ (B11)
$\displaystyle+\left((\gamma(1-\gamma)+\gamma\delta)\Delta t^{2}{\bf
A}+\gamma^{2}(1-\gamma)\Delta t^{3}{\bf A}^{2}+\delta\Delta t{\bf
I}\right){\bf G}({\bf u}^{n})$ $\displaystyle+\left(\gamma(1-\delta)\Delta
t^{2}{\bf A}+(1-\delta)\Delta t{\bf I}\right){\bf G}({\bf u}^{*})\,.$

To find the second order approximation of the numerical solution (B) we first
note that

$\displaystyle{\bf A}$ $\displaystyle\simeq{\bf M}+\Delta t{\bf M}^{2}\,,$
(B12) $\displaystyle{\rm G}({\bf u}^{*})$ $\displaystyle\simeq{\bf G}({\bf
u}^{n})+\gamma\Delta t\left(\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}{\bf
G}({\bf u}^{n})+\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}{\bf M}{\bf
u}^{n}\right)\,,$ (B13)

where we have used that ${\bf k}_{1}\simeq\Delta t{\bf M}{\bf u}^{n}$ to first
order.

Replacing these two expansions in Eq. (B), and neglecting terms with order
higher than $(\Delta t)^{2}$, we obtain

$\displaystyle{\bf u}^{n+1}$ $\displaystyle\simeq\left({\bf I}+\Delta t({\bf
M}+\Delta t{\bf M}^{2})+\gamma(1-\gamma)\Delta t^{2}{\bf M}^{2}\right){\bf
u}^{n}$ $\displaystyle+\left((\gamma(1-\gamma)+\gamma\delta)\Delta t^{2}{\bf
M}+\delta\Delta t{\bf I}\right){\bf G}({\bf u}^{n})$ (B14)
$\displaystyle+\left(\gamma(1-\delta)\Delta t^{2}{\bf M}+(1-\delta)\Delta
t{\bf I}\right)$ $\displaystyle\times\left[{\bf G}({\bf u}^{n})+\gamma\Delta
t\left(\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}{\bf G}({\bf
u}^{n})+\left[\nabla_{u}{\bf G}({\bf u})\right]^{n}{\bf M}{\bf
u}^{n}\right)\right]\,.$ (B15)

which can be cast as

$\displaystyle{\bf u}^{n+1}$ $\displaystyle\simeq{\bf u}^{n}+\Delta t({\bf
M}+\Delta t{\bf G}({\bf u}^{n})$ (B16) $\displaystyle+\Delta
t^{2}\left(\gamma(2-\gamma){\bf M}^{2}+\gamma(1-\delta)\left[\nabla_{u}{\bf
G}({\bf u})\right]^{n}{\bf M}\right){\bf u}^{n}$ (B17) $\displaystyle+\Delta
t^{2}\left(\gamma(1-\delta)\left[\nabla_{u}{\bf G}({\bf
u})\right]^{n}+\gamma(2-\gamma){\bf M}\right){\bf G}({\bf u}^{n})\,.$ (B18)

In order to match the Taylor series (B6) for $t=t_{n}+\Delta t$, we must have

$\gamma(2-\gamma)=\frac{1}{2}\,,\quad\quad\quad\gamma(1-\delta)=\frac{1}{2}$
(B19)

These two conditions give

$\gamma=1\pm\frac{1}{\sqrt{2}}\,,\quad{\rm
and,}\quad\delta=1-\frac{1}{2\gamma}\,,$ (B20)

which is the solution from Ascher et al. (1997).

##  Appendix C Convergence to the correct equilibrium

To study the convergence properties of the numerical method we consider the
case of Eq. (B1) with constant external forces ${\bf g}({\bf u})={\bf g}_{0}$.
In this case, the second-order Runge-Kutta IMEX method described in Section
2.1 reduces to

$\displaystyle{\bf u}^{n+1}$ $\displaystyle=\left({\bf I}+\Delta t{\bf
A}+\gamma(1-\gamma)\Delta t^{2}{\bf A}^{2}\right){\bf u}^{n}$
$\displaystyle+\left(\gamma(2-\gamma)\Delta t^{2}{\bf
A}+\gamma^{2}(1-\gamma)\Delta t^{3}{\bf A}^{2}+\Delta t{\bf I}\right){\bf
G}_{0}\,.$ (C1)

Since ${\bf A}$ and ${\bf A}^{2}$ are diagonalizable under the same
transformation (see Section C.1), we apply the transformation matrix ${\bf
P}^{-1}$ to Eq. (C) and obtain

$\displaystyle\hat{u}^{n+1}_{k}$ $\displaystyle=\hat{u}^{n}_{k}\left(1+\Delta
t\frac{\lambda_{k}}{1-\gamma\Delta t\lambda_{k}}+\gamma(1-\gamma)\Delta
t^{2}\left(\frac{\lambda_{k}}{1-\gamma\Delta t\lambda_{k}}\right)^{2}\right)$
$\displaystyle+\hat{G}_{k}\left[\gamma(2-\gamma)\Delta
t^{2}\frac{\lambda_{k}}{1-\gamma\Delta
t\lambda_{k}}+\gamma^{2}(1-\gamma)\Delta
t^{3}\left(\frac{\lambda_{k}}{1-\gamma\Delta t\lambda_{k}}\right)^{2}+\Delta
t\right]\,,$ (C2)

where $\hat{{\bf u}}={\bf P}^{-1}{\bf u}$ and $\hat{u}_{k}$ corresponds to the
entry of $\hat{{\bf u}}$ associated to the eigenvalue $\lambda_{k}$ of the
matrix ${\bf M}$ (equivalently for $\hat{G}_{k}$). As mentioned in Section 2,
all eigenvalues of the matrix ${\bf M}$ are negative, except for a unique
eigenvalue $\lambda_{j}=0$. We define $x=|\lambda_{k}|\Delta t$ (with
$|\lambda_{k}|=-\lambda_{k}$) and write Eq. (C) as

$\displaystyle\hat{u}^{n+1}_{k}=\hat{u}^{n}_{k}\left(\frac{x(2\gamma-1)+1}{(1+\gamma
x)^{2}}\right)+\frac{\hat{G}_{k}}{|\lambda_{k}|}\frac{x(\gamma^{2}x+1)}{(1+\gamma
x)^{2}}\,.$ (C3)

For the case of the eigenvalue $\lambda_{j}=0$,
$\hat{u}^{n+1}_{j}=\hat{u}^{n}_{j}+\Delta t\hat{G}_{j}$, which corresponds to
a uniformly accelerated center of mass of the system. For the case with
$\lambda_{k}\neq 0$, the limit of $x\rightarrow\infty$ gives

$\lim_{x\rightarrow\infty}\hat{u}^{n+1}_{k}=\lim_{x\rightarrow\infty}\left(\hat{u}^{n}_{k}\frac{2\gamma-1}{\gamma^{2}x}+\frac{\hat{G}_{k}}{|\lambda_{k}|}\left(1+\frac{1-2\gamma}{\gamma^{2}x}\right)\right)=\frac{\hat{G}_{k}}{|\lambda_{k}|}\,.$
(C4)

This limit is identical to that obtained after applying the transformation
${\bf P}^{-1}$ to the analytical solution (B2) with a constant external force
${\bf G}={\bf G}_{0}$ (see Eq. (35)). Thus, the solution corresponds to a
system with constant velocities relative to the center of mass. Note as well
that in cases with no external acceleration
$\lim_{x\rightarrow\infty}\hat{u}^{n+1}_{k}\propto\lim_{x\rightarrow\infty}1/x$,
which explains the asymptotic convergence of the error as $e\sim 1/\Delta t$
obtained in Section 3.1.

###  C.1  The matrix ${\bf A}$ is diagonalizable

First note that ${\bf M}$ is diagonalizable since it is similar to a real
symmetric matrix as shown in Benítez-Llambay et al. (2019). Now, we define
${\bf P}$ as the transformation matrix such that ${\bf P}^{-1}{\bf M}{\bf
P}={\bf\Lambda}$, with ${\bf\Lambda}$ a diagonal matrix. The transformation
defined by ${\bf P}$ also diagonalize ${\bf I}-\gamma\Delta t{\bf M}$,
together with the inverse $\left({\bf I}-\gamma\Delta t{\bf M}\right)^{-1}$.
Therefore,

$\displaystyle{\bf P}^{-1}{\bf A}{\bf P}$ $\displaystyle={\bf
P}^{-1}\left(\left({\bf I}-\gamma\Delta t{\bf M}\right)^{-1}{\bf M}\right){\bf
P}$ $\displaystyle=\left({\bf P}^{-1}\left({\bf I}-\gamma\Delta t{\bf
M}\right)^{-1}{\bf P}\right)\left({\bf P}^{-1}{\bf M}{\bf P}\right)={\bf
D}{\bf\Lambda}={\bf A}_{\rm D}$ (C5)

where ${\bf D}$ and ${\bf A}_{\rm D}$ are diagonal matrices.

##  Appendix D On the value of $\gamma$

We stress that $\gamma=1-1/\sqrt{2}$ (instead of $\gamma=1+1/\sqrt{2}$) will
provide the most accurate solution since the truncation error is minimized for
this value. However, depending on the problem at hand, the value of
$\gamma=1-1/\sqrt{2}$ may not be consistent with a TVD Runge-Kutta in the form
of (Shu & Osher, 1988). For example, the explicit part of the solver in
Section 2.2 can be written as

$\displaystyle{\bf u}^{(0)}$ $\displaystyle={\bf u}^{n}$ $\displaystyle{\bf
u}^{(1)}$ $\displaystyle=\alpha_{10}{\bf u}^{0}+\beta_{10}\Delta
t\mathcal{H}({\bf u}^{0})$ $\displaystyle{\bf u}^{(2)}$
$\displaystyle=\alpha_{20}{\bf u}^{(0)}+\alpha_{21}{\bf
u}^{(1)}+\beta_{20}\Delta t\mathcal{H}({\bf u}^{(0)})+\beta_{21}\Delta
t\mathcal{H}({\bf u}^{(1)})$ $\displaystyle{\bf u}^{n+1}$ $\displaystyle={\bf
u}^{(2)}$ (D1)

with $\alpha_{10}=1$, $\alpha_{20}=1$, $\alpha_{21}=0$, $\beta_{10}=\gamma$,
$\beta_{20}=\delta$, $\beta_{21}=1-\delta$. These values are consistent with
TVD second-order Runge-Kutta provided $\gamma=1+1/\sqrt{2}$, since all the
coefficients of the explicit update are non-negative (Gottlieb & Shu, 1998).
However, when $\gamma=1-1/\sqrt{2}$ the coefficient $\beta_{20}=\delta<0$.

There is yet another concern about the value $\gamma=1-1/\sqrt{2}$. To
showcase the problem, we consider a reduced system where only drag forces are
taken into account. In this case, the numerical solution of the system of
equations (3)-(4) reduces to the form

${\bf u}^{n+1}=\left({\bf I}+\Delta t{\bf A}+\gamma(1-\gamma)\Delta t^{2}{\bf
A}^{2}\right){\bf u}^{n}\equiv{\bf T}(\Delta t,\gamma){\bf u}^{n}\,$ (D2)

where ${\bf A}$ is defined in Eq. (B10). If for a given time step the
$\rm{det}\left({\bf T}(\Delta t,\gamma)\right)=0$, the solution ${\bf
u}^{n+1}$ may deviate from the correct instantaneous equilibrium. First, note
that the eigenvalues, $a_{k}$, of ${\bf A}=({\bf I}-\gamma\Delta t{\bf
M})^{-1}{\bf M}$ are

$a_{k}=\frac{\lambda_{k}}{1-\gamma\Delta t\lambda_{k}}\,,$ (D3)

with $\lambda_{k}$ eigenvalue of ${\bf M}$. Now, we define $\mu_{k}$ the
eigenvalues of ${\bf T}$ and note that

$\mu_{k}=1+\Delta t\frac{\lambda_{k}}{1-\gamma\Delta
t\lambda_{k}}+\gamma(1-\gamma)\Delta
t^{2}\left(\frac{\lambda_{k}}{1-\gamma\Delta t\lambda_{k}}\right)^{2}\,.$ (D4)

We introduce the variable $x=\Delta t|\lambda_{k}|$ as in Section C with
$|\lambda_{k}|=-\lambda_{k}$. Thus,

$\mu_{k}(x)=1-\frac{x}{1+x}+\gamma(1-\gamma)\frac{x^{2}}{(1+\gamma x)^{2}}\,.$
(D5)

which can be reduced to

$\mu_{k}(x)=\frac{x(2\gamma-1)+1}{(1+\gamma x)^{2}}\,.$ (D6)

We conclude that $\mu_{k}(x)=0$ for $x=1/(1-2\gamma)$ and therefore

${\rm det}({\bf T})=\prod_{k}\mu_{k}=0\iff x=1/(1-2\gamma)\,.$ (D7)

In terms of $\Delta t$, the condition (D7) implies

${\rm det}({\bf T})=0\iff\Delta t=\frac{|\lambda_{k}|^{-1}}{1-2\gamma}\,,$
(D8)

for any $k\in[0,N-1]$. In the case of $\gamma=1-1/\sqrt{2}$, $1-2\gamma>0$ and
therefore there is $\Delta t>0$ such that ${\rm det}({\bf T})=0$. In the other
hand, if $\gamma=1+1/\sqrt{2}$, $1-2\gamma<0$ and therefore ${\rm det}({\bf
T})\neq 0$ for all $\Delta t>0$. The same result also holds for the second-
order Runge-Kutta of Pareschi & Russo (2005). Finally, the same analysis can
be applied to a case with constant external force, however, in that case, the
condition that could affect the solution corresponds to $x=-1/\gamma^{2}$,
which is never satisfied for $\Delta t>0$.

##  Appendix E Flux calculation

To calculate the mass and momentum flux ($\mathcal{L}$ operator neglecting
external forces) we consider a 1D Cartesian evenly spaced mesh along the
$x$-direction. For each fluid, the operator $\mathcal{L}$ in Eq. (2.2) returns
the divergence of the flux as follows

$\mathcal{L}({\bf U}_{i}^{n})=-\frac{1}{\Delta x}\left({\bf
F}_{i+1/2}({\bf{U}}^{n})-{\bf F}_{i-1/2}({\bf{U}}^{n})\right)\,,$ (E1)

where the index $i$ denotes the spatial position in the discretized mesh and
$i\pm 1/2$ indicates the position of the interfaces where fluxes are
calculated. ${\bf U}=(\rho,\rho v)$ correspond to the vector of conserved
variables. The flux, ${\bf F}_{i+1/2}$, is obtained by solving the Riemann
problem with initial condition

${\bf U}(x,0)=\left\\{\begin{array}[]{ccc}{\bf U}_{{\rm
L},i+1/2}&\text{if}&x<x_{i+1/2}\,,\\\ {\bf U}_{{\rm
R},i+1/2}&\text{if}&x>x_{i+1/2}\,,\\\ \end{array}\right.$ (E2)

where ${\bf U}_{{\rm L},i+1/2}$ and ${\bf U}_{{\rm R},i+1/2}$ correspond to
the left and right states with respect to the interface at $x=x_{i+1/2}$. In
this work, we use a piece-wise linear reconstruction with a Van-Leer slope
limiter (van Leer, 1977) on the primitive variables ${\bf V}=(\rho,v)$. We
approximate the left and right states as follows

${\bf V}_{{\rm L},i+1/2}={\bf V}_{i}+\frac{1}{2}\Delta_{x}{\bf
V}_{i}\,,\quad{\bf V}_{{\rm R},i+1/2}={\bf V}_{i+1}-\frac{1}{2}\Delta_{x}{\bf
V}_{i+1}\,,$ (E3)

with

$\Delta_{x}{\bf V}_{i}=\frac{{\rm max}\left\\{({\bf V}_{i+1}-{\bf V}_{i})({\bf
V}_{i}-{\bf V}_{i-1}),0\right\\}}{{\bf V}_{i+1}-{\bf V}_{i-1}}\,.$ (E4)

After the reconstruction of the left and right states calculate the fluxes
utilizing the approximate solution from Harten, Lax and van Leer (HLL) (Harten
et al., 1983),

${\bf F}_{i+1/2}=\left\\{\begin{array}[]{ccc}{\bf F}_{{\rm L}}&\textrm{if
}&0<S_{\rm L}\,,\\\\[-5.0pt] \\\ \frac{\displaystyle S_{\rm R}{\bf F}_{\rm
L}-S_{\rm L}{\bf F}_{\rm R}+S_{\rm R}S_{\rm L}({\bf U}_{\rm R}-{\bf U}_{\rm
L})}{\displaystyle S_{\rm R}-S_{\rm L}}&\textrm{if}&S_{\rm L}\leq 0\leq S_{\rm
R}\,,\\\\[-5.0pt] \\\ {\bf F}_{{\rm R}}&\textrm{i}f&0>S_{\rm R}\,,\\\
\end{array}\right.$ (E5)

where ${\bf F}_{{\rm L}}$ and ${\bf F}_{{\rm R}}$ are fluxes computed with
${\bf V}_{{\rm L},i+1/2}$ and ${\bf V}_{{\rm R},i+1/2}$, respectively. Thus,
in terms of primitive variables, we define the vector fluxes as ${\bf F}_{\rm
g}=(\rho_{\rm g}{v}_{\rm g},\rho_{\rm g}{v}_{\rm g}^{2}+P)$, ${\bf F}_{\rm
d}=(\rho_{\rm d}{v}_{\rm d},\rho_{\rm d}{v}_{\rm d}^{2})$, for gas and dust,
respectively. For the gas species, the speeds $S_{\rm R}$ and $S_{\rm L}$ are
estimated as

$S_{\rm L}={\rm min}\\{v_{\rm L}-c_{{\rm s}},v_{\rm R}-c_{{\rm s}}\\}\,,\quad
S_{\rm R}={\rm max}\\{v_{\rm L}+c_{{\rm s}},v_{\rm R}+c_{{\rm s}}\\}\,,$ (E6)

where $v_{\rm L}$ and $v_{\rm R}$ are obtained after transforming the
conservative left and right states into primitive variables. For pressureless
fluids (dust), the sound speed is set to zero, $S_{\rm L}=-S_{\rm R}$, and
$S_{\rm R}={\rm max}\\{|v_{{\rm L}}|,|v_{{\rm R}}|\\}$. Therefore, for the
dust species the flux in Equation (E5) reduces to the Rusanov flux (Rusanov,
1970). We have also implemented the Riemann solver for the dust described in
(Huang & Bai, 2022). We utilize this method in the shock test and steady-state
gas and dust drift calculation. While in this work we adopted the (perhaps)
simpler flux estimation, there are different alternatives for the gas (see
e.g., Toro, 2009), and the dust as a pressureless fluid Leveque (2004);
Paardekooper & Mellema (2006).

##  Appendix F Gravitational potential

The gravitational potential is included as a source term in Eqs. (3)-(4). To
incorporate the potential into our numerical method we solve the discretized
form of the Poisson equation using a sparse linear solver with periodic
boundary conditions,

$\frac{\Phi_{i+1}-2\Phi_{i}-\Phi_{i-1}}{\Delta x^{2}}=4\pi G\rho_{i}^{n}\,,$
(F1)

and obtain $\Phi_{i}$ and include the potential into the $\mathcal{L}$
operator as

$\displaystyle\mathcal{L}({\bf U}^{n}_{i})$ $\displaystyle=-\frac{\Delta
t}{\Delta x}\left(\left({\bf F}_{i+1/2}({\bf U}^{n})-{\bf F}_{i-1/2}({\bf
U}^{n})\right)-\frac{\rho_{i}^{n}}{2}\left(\Phi_{i+1}-\Phi_{i-1}\right)\right)\,.$
(F2)