A Fast Second-order Solver for Stiff Multifluid Dust and Gas Hydrodynamics

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 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 O(N) operations, instead of standard matrix inversion, which is 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 multifluid 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.


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 solids 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 Corresponding author: Leonardo Krapp krapp@arizona.edumulti-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;Bai & Stone 2010;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 Semianalytic 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 (and time scales) 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 firstorder time-implicit scheme for the drag force with an explicit first-order hydrodynamics solver.This multispecies numerical solver is unconditionally and asymptotically stable and conserves momentum to machine precision.Moreover, the complexity of the scheme is order 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 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 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.

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.Each dust species is characterized by a density, velocity, and a characteristic stopping time1 , t s , or equivalent a rate of momentum exchange with the gas, α d .Gas density and velocity are defined as ρ g and v g , respectively.Equivalently, for the i-th dust species we define the density and velocity as ρ d,i and v 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 s ρ g , with c s as the sound speed.Under these assumptions, the continuity and momentum equations for gas and dust species are with i = 1, . . ., N − 1.The right-hand sides of Equations (3) and (4) containt the drag force terms and any external forces G g and G d,i , respectively.The rates of momentum exchange are α g,i and α d,i .Momentum conservation implies therefore α g,i = ϵ i α d,i , with ϵ i = ρ d,i /ρ 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 now, Eqs. ( 3)-( 4) can be cast into the form where M is the N × N matrix.
Note that in our case the matrix 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 M BKP19 , defined in equation ( 11) of Benítez-Llambay et al. (2019) (see Appendix G).In this work the matrix M = −RM BKP19 R −1 , where R a diagonal matrix such that R ii = ρ 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 λ k is an eigenvalue of M, then −λ k is an eigenvalue of M BKP19 .Note that the eigenvalues of M 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 x is an eigenvector of M BKP19 , then x = −Rx is an eigenvector of the matrix M. Since

Implicit solutions of the drag-force term
Before introducing a numerical method to solve Eqs. ( 1)-( 4), we will discuss three different alternatives to solve the momentum equation ( 9) neglecting all terms, except for the drag force.The three methods correspond to the first-order fully implicit method of Benítez-Llambay et al. (2019), a second-order fully implicit method (see e.g, Huang & Bai 2022), and the diagonally implicit Runge-Kutta method (DIRK) utilized in this work.Our goal is to highlight the major differences between the three methods when solving the momentum equation We define the numerical solution at t = t n +∆t of Eq.( 5) as u n+1 .A first-order fully implicit solution can be obtained after solving whereas a second-order fully implicit method can be cast as These two fully implicit methods converge to the correct equilibrium solution and are asymptotically stable, that is, In particular, utilizing the properties of the matrix M it can be shown that the second-order method asymptotically converges as lim ∆t→∞ |u n+1 − u(∆t)| ∼ 1/∆t 2 , whereas the first-order method converges as lim ∆t→∞ |u n+1 − u(∆t)| ∼ 1/∆t.Therefore, the fully implicit second-order method provides a more accurate solution for all ∆t.
The advantage of the first-order method is that no matrix inversion is required to obtain the solution (Krapp & Benítez-Llambay 2020).The matrix inversion can add a complexity O(N 3 ) to the problem, making the second-order method ∼ ×N 2 more expensive than the first-order method.The extension of the analytical solution of Krapp & Benítez-Llambay (2020) to the second order method is challenging since I − ∆tM + ∆t 2 2 M 2 is not a sparse matrix.
Ideally, we would like a numerical solution that takes advantage of the reduced complexity of the first-order method and preserves the accuracy of the fully implicit second-order method.Thus, in this work, we study and test the properties of a second-order diagonally implicit Runge-Kutta (DIRK) method of the form where γ is known as the singular value of the DIRK method (Alexander 1977).This two-stage method is second-order accurate (see Appendix B).Moreover, the solution of each stage requires the inversion of the matrix (I − γ∆tM) and therefore an extension of the analytical solution of (Krapp & Benítez-Llambay 2020) can be utilized to reduce the complexity to order O(N ).We will show the analytical solution of the DIRK method in Section 2.2.We also emphasize that coupling the solution (17) with other numerical methods for hydrodynamics is not a trivial task, as standard splitting may fail to converge to the correct equilibrium.In Section 2.2 we describe the Implicit-Explicit Runge-Kutta method from Ascher et al. (1997) to solve the full system of momentum equations ( 3)-( 4) utilizing the solution (17).Continuity equations are then included in the MDIRK method described in Section 2.3.

On the value of γ and stiff condition
Alternative versions of the method of Alexander (1977) can be obtained by varying the coefficients of the Runge-Kutta.Depending on the coefficients, the value of γ is obtained either by the order conditions (e.g., second-order accuracy as ∆t → 0 and/or as ∆t → ∞) or by the asymptotic stability condition lim ∆t→∞ |u n+1 − u(∆t)| → 0. Since the form of Eqs. ( 15)-( 17) provides an asymptotically stable solution for any value of γ (see Appendix C), we will use the order conditions to set its value.
Instead of utilizing the Taylor series at ∆t → 0, one could also set the value of γ by matching the series expansion of u n+1 and u(∆t) at ∆t → ∞.It can be shown that u n+1 asymptotically converges to the correct solution as , for some λ k (non-zero) eigenvalue of the matrix M. In this case, second-order accuracy is obtained with γ = 1/2, otherwise, the method will converge as a first-order method.
We note that the two order conditions can not be simultaneously satisfied by the same value of γ.Therefore, there must be a criterion condition that sets the value of γ for a given time step ∆t.We define this criterion as, where max(α −1 d,i ) is the largest stopping time among the dust species.The choice of ∆t > max(α −1 d,i ) is only approximate since formally one should look into the non-zero eigenvalues of the matrix M. For instance, the timescale for the condition (18) could be better defined in terms of the largest eigenvalue (nonzero and smallest in absolute value) of the matrix M, e.g., ∆t > |max(λ k )| −1 .However, it can be shown that |max(λ k )| −1 is bound between the two largest stopping times of the system (see Appendix D Benítez-Llambay et al. 2019).This justify our choice of ∆t > max(α −1 d,i ).As we describe below, when ∆t ∼ max(α −1 d,i ) the truncation error is in the asymptotic regime.We therefore define the system of Eqs ( 1)-( 4) as a "stiff system" when ∆t ≥ max(α −1 d,i ).To highlight the definition of the stiff condition and further compare the methods described in this section, we solve Eq. 11 for a system with N = 17 fluids (16 dust species).For each dust fluid, the dust-to-gas mass ratios are obtained as s,max ), and ϵ = 1.0.The stopping time t s,k = α −1 d,k is uniformly spaced in logscale between t s,min = 10 −3 and t s,max = 10.Note that t s,max = max(α −1 d,i ) in Eq. ( 18).In Figure 1 we show the truncation error, e(∆t), for the three different methods described in this section.The error is obtained as As ∆t → 0, the difference between the Taylor series implies that e(∆t) ∼ ∆t 2 and e(∆t) ∼ ∆t 3 for the first and second-order methods, respectively.These scaling are shown in Figure 1 for ∆t ≲ 0.001.As already mentioned in this section, the DIRK method is more accurate for In the regime where t s,min ≲ ∆t ≲ t s,max , the DIRK method shows an error comparable to that obtained with the firstorder implicit method.Similarly, the second-order fully implicit method is also more accurate in this regime.
Comparison of the (local) truncation error for different methods described in Section 2.1.The blue solid line corresponds to the first-order implicit method of Eq. ( 12) whereas the red solid line shows the error for the second-order implicit method of Eq. ( 13).The solutions obtained with the DIRK method are shown in blue and red dashed lines.The black solid line indicates the transition to a stiff regime where γ is set to γ = 1/2 for the DIRK method.For ∆t ≳ 0.002 (vertical dotted line) time explicit methods will be unstable.The analytical solution, u k (∆t), is obtained by solving the eigenvalue problem of the linear equation ( 11).
As the time-step increases the system reaches the stiff regime where ∆t ≥ t s,max = 10.In this regime, the numerical solutions show the expected asymptotic convergences e(∆t) ∼ ∆t −1 and e(∆t) ∼ ∆t −2 for the first and second-order methods, respectively.If condition ( 18) is not applied to the DIRK method, i.e. γ = 1 ± 1/ √ 2 for all ∆t, the asymptotic convergence will be that of the first-order method.Therefore, without condition (18) for γ there is not a clear advantage of utilizing a DIRK method in the stiff regime.

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.3.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 sta-ble (L-stable), and converges to the correct equilibrium.We can write Eq. ( 9) as, with H(u) = −∂ x F + G.The brackets [] exp and [] 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 O(N 3 ) to O(N ), making it suitable for a system with a large number of species.Consider the momenta vector u(t n ) = u n .The full update to a new time t = t n + ∆t is obtained as follows: where the coefficients γ, δ and β are defined in Eq. ( 30).
The vectors w 1 and w 2 are the solution of the following equations: Alternatively, the solution can be also obtained in terms of k 1 = Mw 1 and k 2 = Mw 2 , by applying the M matrix to Eqs .(22)-( 23).
The major advantage of the diagonally implicit method is that one first solves for k 1 and then for k 2 , which simplifies the calculations.Note that this is not the case in fully implicit methods since Eq. ( 24) will also depend on k 2 .Moreover, we can find the analytical solution of the linear system, which has a complexity O(N ) (c.f.Krapp & Benítez-Llambay 2020).For instance, consider the arbitrary vector q n = (q n g , q n d,1 , . . ., q n d,N −1 ) and the equation for i = 1, . . ., N − 1, with It can be shown that γ = 1 ± 1/ √ 2 are the only two possible values for γ that provide an asymptotically stable and second-order accurate solution (for ∆t → 0) with two-stages (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) An alternative version of the second-order Implicit-Explicit Runge-Kutta can be found in Pareschi & Russo (2005) 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 γ = 1 + 1/ √ 2, provided a Courant number C ≲ 0.58.

MDIRK
To solve the fluid equations we have reduced the Runge-Kutta described in Section 2.2 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 U n g = (ρ n g , ρ n g v n g ) = (U n g0 , U n g1 ) and K n g = (0, k n g ) for the gas, and ) for the i-th dust-species, with k n g and k n d,i defined in Eq. ( 27) and Eq. ( 28), respectively.Therefore, K g and K 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 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 fluid2 , a full update of the conserved variables from time step t n to time step t n+1 = t n +∆t is obtained as follows with The coefficients γ, δ, β, b 1 and b 2 are given in Eq. ( 30).
The vectors K * and K p are obtained from the analytical solution given in Eqs. ( 27)-( 28).When utilizing these analytical solutions, the densities and momenta must be taken from the partial updates U * and U p , respectively.For example, with Note that when the vector K is set to zero (no drag forces), second-order accuracy can be obtained by setting γ = 1 and δ = 1/2.

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 where for any arbitrary q (which for example can be taken as u * or u p ), the vector k satisfies From left to right we plot the evolution of the velocity for the gas and two dust species.The dashed line corresponds to the analytical solution, given by Eq (39).The circles indicate the numerical solutions obtained with different integrators.
Green color shows the first-order method, orange and dark lines the second-order method with γ = γ √ 2, respectively.The upper and lower panels show the results for the tests 1 and 2, respectively.In the top and bottom panels, the time step is ∆t1 = 0.005, ∆t2 = 0.05, respectively (see Section 3.1).In the rightmost panel, we plot, the error by Eq. ( 40) as a function of the time step ∆t, for the first-order method (green color) and the DIRK method with γ = 1 ± 1/ √ 2. The vertical solid line corresponds to ∆t = ts,max, where condition ( 18) is applied for the DIRK method.

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 multifluid shock, an idealized, multi-fluid Jeans instability, and steady-state drift in a 2D shearing box.

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 for a given initial condition u 0 .After solving the eigenvalue problem Mx = λx, the solution of Eq. ( 38) can be expressed as u(t) = N −1 k=0 ck e λ k 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 where v CM is the velocity of the center of mass.The values of the eigenvalues λ 1 , λ 2 and the coefficients c 1 and c 2 can be obtained from table 1 of Huang & Bai (2022).
In Figure 2 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 2. The error is obtained as where ∆t i is the time step and M is the number of time steps during the integration.The index k = 0, 1, 2 is the index of gas, and dust species, respectively.We calculate the numerical solution for different time steps starting at ∆t = 10 −4 until a final time t final > t damp , where t damp ≃ |max(λ k ) −1 | is a characteristic damping time to reach the equilibrium velocity.When ∆t ≤ t damp , the number of times steps M is given by the ratio t damp /∆t.When ∆t ≥ t damp we only have one timestep, i.e., M = 1.For the case of test B, t final = 2 and t damp ≃ 0.007, whereas in test C, t final = 2 × 10 3 and t damp ≃ 1.9.
In Figure 2 we show that the error as a function of ∆t.In the regime of ∆t < t damp , we obtain e 1 ∼ ∆t 2 for both values of γ of the DIRK method.When γ = 1 + 1/ √ 2, the absolute error is significantly larger in comparison with the solution obtained for γ = 1 − 1/ √ 2. However, for γ = 1 − 1/ √ 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 2 where the gas velocity initially undershoots the analytical solution (see Appendix D for additional comments on the value of γ).
When ∆t ≫ t damp , the error obtained from Eq. ( 40) should agree with the characteristic asymptotic convergence rate of the numerical method.For large time steps, the error displayed in Figure 2 converges as ∼ 1/∆t for the first-order method, whereas for the DIRK method, the error goes as ∼ 1/∆t 2 .Therefore, the DIRK method provides a second-order accurate solution in the stiff regime.We significantly improve the convergence rate of the DIRK method utilizing condition for γ given in Eq. ( 18).The condition sets γ = 1/2 if ∆t ≥ t s,max .In the top panel, t s,max = 0.01 whereas in the bottom panel t s,max = 2.Note that these stopping times are upper bounds for t damp as expected from the discussion in Section 2.3 and the results on Appendix D of Benítez-Llambay et al. (2019).

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.3.
We again consider one gas and two dust species and solve the momentum equation using Eq. ( 31), neglecting the continuity equation as well as the divergence of fluxes.Under these assumptions, the momentum equation becomes  with G = (G 0 , 0, 0).The analytical solution of Eq. ( 41) corresponds to, where λ 1 and λ 2 are the non-zero eigenvalues of M and C is a coefficient matrix that depends on the initial condition, external force, and eigenvalues of M. We set as initial condition ρ g = 1.0, ρ d1 = 0.1, ρ d2 = 0.1 for the density and v g = 2.0, v d1 = 0.1, v d2 = −0.(43) In Figure 3 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 ∆t → ∞, the momentum of the j-th species (see Eq, ( 42)) converges to where we assume λ 0 = 0 and P jk denotes the jk entry of the matrix P, with P and P −1 the matrices such that  (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 δ ρ/(Aρ 0 ), with A = 10 −4 .We include in the plot the L1 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.
After dividing Eq. ( 44) by the corresponding density we can obtain the relative velocity between species.First, note that P j0 /ρ j − P i0 /ρ i = 0, for any i and j components.This follows from the properties of the matrix M described below Eq. (G59) and using that the first (zero) column of the matrix P is the eigenvector 1T .Therefore, the relative velocity between two species is obtained as 45) In this case we obtain lim ∆t→∞ |v g − v d1 | = 0.8333 and lim ∆t→∞ |v g − v d2 | = 1.111, which are the asymptotic values displayed in Figure 3.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.

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 secondorder Runge-Kutta method of Eq. ( 31).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 de-  (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 unfilled circles correspond to the actual grid points.The code resolves the shock within a single cell, even when an increasing number of fluids is considered.The overall agreement between the numerical and analytical solutions is excellent.
scribed in Section 2.3 achieves the expected second-order accuracy.
In Figure 4, we show the analytic (solid lines) and numerical (open circles) solutions for the normalized density, defined as δ ρ/(ρ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 with k the index of the fluid and i the index of the density and velocity field, respectively.In the bottom panel of Figure 4 we show that the error effectively scales as e 1 ∼ N −2 x .

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, ρ g α g,i = ρ d,i α d,i ≡ K i .Since the dragforce solver has been written in terms of the parameter ϵ i , the drag-force coefficient α d,i must be replaced by α d,i = K i /ρ d,i in the matrix M. The details of the test together with initial and boundary conditions are in section 3.3.2from Benítez-Llambay et al. (2019).
In Figure 5, 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 xpositive 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).

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 with D/Dt = ∂/∂ t + v • ∇ the material derivate.The Poisson equation that governs the gravitational potential of the system is In order to study the Jeans instability, we linearize Equations ( 1)-(2), and ( 47)-( 49) with respect to an equilibrium solution obtained after neglecting the contribution of the background gravitational field in the Poisson equation (Jeans Swindle) (Jeans 1902;Binney & Tremaine 2008).The equilibrium gas density is constant and denoted as ρ 0 g .The equilibrium dust density of the i-th species is also a constant, equal to ρ 0 d,i = ϵ 0 i ρ 0 g .Finally, the equilibrium velocities are set to zero.We define small density perturbations δρ g and δρ d,i , whereas the velocity perturbations are δv g = δv g x, δv d,i = δv d,i x.We set ρ 0 g = 1.0 and the total dust-to-gas mass ratio ϵ 0 = N −1 i=1 ϵ 0 i = 0.01, with ϵ 0 i = ϵ 0 /N .We consider N = 128 fluids with a uniform distribution of stopping times between t s = 10 −4 − 10 The linearized continuity, momentum, and Poisson equation are For N fluids, these 2N + 1 equations are coupled, due to the gravitational field and the drag-force.
Figure 6 displays the grow rate of the dusty Jeans instability σ = − min{ℜ(ω)} √ 1 + ϵ 0 as a function of k/k J .Note that when ϵ 0 = 0, σ reduces to the normalized growth rate of the gas-only problem.Thus, the inclusion of dust enables higher growth rates as shown in Figure 6 (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 J for stability is shifted to larger values as dust density increases.In Figure 6 we also the numerical solution obtained with MDIRK described in Section 2.3.We calculate the growth rate varying the number of cells from N cells = 16 to N cells = 128.The growth rate converges to the analytic solution for all unstable modes considered in this work 3 .

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 Ω 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 (f ce = ρ2qΩ 2 0 xx) and the Coriolis force f c = −2ρΩ 0 ẑ × (v x , v y ).In addition, we add an external constant force on the gas 3 The largest wavenumber considered in this work corresponds to the cutoff mode of the dust-free case, i.e., k/k J = 1/ √ 1 + ϵ 0 fluid that mimics the global pressure gradient in a protoplanetary disk, f g = χ 0 Ω 0 x with χ 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 ρ, ρv x .Whereas for v y we add the relative shear velocity after applying periodic boundary conditions, that is ρv y (x) = ρv y (x ∓ L) ± qΩ 0 Lρ (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 7, 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 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 O(N 3 )) for a large number of fluids (N > 10).It is interesting to note that the complexity of Huang & Bai (2022) method can be reduced to O(N ) by utilizing our analytical solution in Appendix A.

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 secondorder 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 O(N 3 ) (needed to invert the matrix numerically) to 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 dustgrowth 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 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 7. Normalized CPU time vs the number of fluids for the test described in Section 3.6.As expected, the complexity of MDIRK is O(N ).For comparison, we also include the results of the multifluid second-order fully implicit method (RK2Im) from Huang & Bai (2022) which has complexity O(N 3 ).
In order to match the Taylor series (B20) for t = t n + ∆t, we must have These two conditions give which is the solution from Ascher et al. (1997).

C. CONVERGENCE TO THE CORRECT EQUILIBRIUM
To study the convergence properties of the numerical method we consider the case of Eq. (B15) with constant external forces g(u) = g 0 .In this case, the second-order Runge-Kutta IMEX method described in Section 2.2 reduces to Since A and A 2 are diagonalizable under the same transformation (see Section C.2), we apply the transformation matrix P −1 to Eq. ( C39) and obtain where û = P −1 u and ûk corresponds to the entry of û associated to the eigenvalue λ k of the matrix M (equivalently for Ĝk ).As mentioned in Section 2, all eigenvalues of the matrix M are negative, except for a unique eigenvalue λ j = 0. We define x = |λ k |∆t (with |λ k | = −λ k ) and write Eq. (C36) as For the case of the eigenvalue λ j = 0, ûn+1 j = ûn j + ∆t Ĝj , which corresponds to a uniformly accelerated center of mass of the system.For the case with λ k ̸ = 0, the limit of x → ∞ gives ) This limit is identical to that obtained after applying the transformation P −1 to the analytical solution (B16) with a constant external force G = G 0 (see Eq. ( 44)).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→∞ ûn+1 k ∝ lim x→∞ 1/x, which explains the asymptotic convergence of the error as e ∼ 1/∆t obtained in Section 3.1.

C.1. Strang-splitting
We now show that a strang-splitting method does not converge to the correct equilibrium solution.Consider a strangsplitting method where the drag force is integrated on ∆t/2 before and after a full integration of the hydrodynamics.
To study the convergence of the strang-splitting method we assume u n = 0 since the initial condition should not affect the convergence to the correct equilibrium.Thus, we skip the first implicit update on the half step.The solution of the strang-splitting method is obtained as follows Following the diagonalization steps and normalization described for Eq.(C36) we obtain for some eigenvalue λ k ̸ = 0. Now we take the limit of x → ∞ and obtain The convergence to the correct equilibrium implies 4γ − 2/γ 2 = 1 which gives γ = 2 ± √ 2. However, this value is not consistent with the second-order accuracy requirement discussed in Section 2.3.Therefore, the strang-splitting fashion will not converge to the correct steady-state solution.

C.2. The matrix A is diagonalizable
First note that M is diagonalizable since it is similar to a real symmetric matrix as shown in Benítez-Llambay et al. (2019).Now, we define P as the transformation matrix such that P −1 MP = Λ, with Λ a diagonal matrix.The transformation defined by P also diagonalize I − γ∆tM, together with the inverse (I − γ∆tM) −1 .Therefore, where D and A D are diagonal matrices.

D. ON THE VALUE OF γ
We stress that γ = 1 − 1/ √ 2 (instead of γ = 1 + 1/ √ 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 γ = 1 − 1/ √ 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.3 can be written as u (0) = u n u (1) = α 10 u 0 + β 10 ∆tH(u 0 ) u (2) = α 20 u (0) + α 21 u (1) + β 20 ∆tH(u (0) ) + β 21 ∆tH(u (1) ) with α 10 = 1, α 20 = 1, α 21 = 0, β 10 = γ, β 20 = δ, β 21 = 1 − δ.These values are consistent with TVD second-order Runge-Kutta provided γ = 1 + 1/ √ 2, since all the coefficients of the explicit update are non-negative (Gottlieb & Shu 1998).However, when γ = 1 − 1/ √ 2 the coefficient β 20 = δ < 0. There is yet another concern about the value γ = 1 − 1/ √ 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 u n+1 = I + ∆tA + γ(1 − γ)∆t 2 A 2 u n ≡ T(∆t, γ)u n (D44) where A is defined in Eq. (B24).If for a given time step the det (T(∆t, γ)) = 0, the solution u n+1 may deviate from the correct instantaneous equilibrium.First, note that the eigenvalues, a k , of A = (I − γ∆tM) −1 M are with λ k eigenvalue of M. Now, we define µ k the eigenvalues of T and note that where F L and F R are fluxes computed with V L,i+1/2 and V R,i+1/2 , respectively.Thus, in terms of primitive variables, we define the vector fluxes as F g = (ρ g v g , ρ g v 2 g + P ), F d = (ρ d v d , ρ d v 2 d ), for gas and dust, respectively.For the gas species, the speeds S R and S L are estimated as where v L and v 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 L = −S R , and S R = max{|v L |, |v R |}.Therefore, for the dust species the flux in Equation (E55) 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).

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, and obtain Φ i and include the potential into the L operator as Figure2.From left to right we plot the evolution of the velocity for the gas and two dust species.The dashed line corresponds to the analytical solution, given by Eq (39).The circles indicate the numerical solutions obtained with different integrators.Green color shows the first-order method, orange and dark lines the second-order method with γ = γ+ = 1 + 1/ √ 2 and γ = γ − = 1 − 1/√ 2, respectively.The upper and lower panels show the results for the tests 1 and 2, respectively.In the top and bottom panels, the time step is ∆t1 = 0.005, ∆t2 = 0.05, respectively (see Section 3.1).In the rightmost panel, we plot, the error by Eq. (40) as a function of the time step ∆t, for the first-order method (green color) and the DIRK method with γ = 1 ± 1/ √ 2. The vertical solid line corresponds to ∆t = ts,max, where condition (18) is applied for the DIRK method.

Figure 3 .
Figure3.Relative velocity between the gas and dust species for the damping test 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.
5 for the velocity.The external force is set to G 0 = 1 whereas the drag force coefficients are set to α d,1 = 1.0, α d,2 = 0.75, respectively.With these values we obtain λ 1 = −1.125 and λ 2 = −0.8 and the coefficient matrix C is

Figure 4 .
Figure 4. 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 δ ρ/(Aρ 0 ), with A = 10 −4 .We include in the plot the L1 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 5 .
Figure5.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 unfilled circles correspond to the actual grid points.The code resolves the shock within a single cell, even when an increasing number of fluids is considered.The overall agreement between the numerical and analytical solutions is excellent.

Figure 6 .
Figure 6.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 ts = 10 −4 to 10 and uniform density distribution with a total dust-to-gas ratio ϵ 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/kJ = 1/ √ 1 + ϵ 0 .