© The Authors 2025

1 Introduction

Dust, alongside gas, is one of the fundamental components of most astrophysical systems. From the cold interstellar medium to molecular clouds and protoplanetary disks, dust is observed ubiquitously across the Milky Way and other galaxies. Importantly, in many cases, dust is the only directly observable component, detected through its millimeter and infrared thermal emission. Gas properties are often inferred indirectly from dust observations, typically assuming a constant dust-to-gas ratio of 10⁻². Dust also plays a critical role in a variety of astrophysical processes. It contributes to the cooling of gas in star-forming regions (Vogelsberger et al. 2019), serves as a surface for chemical reactions that form complex molecules (Minissale et al. 2016), and provides the primary opacity source in many environments, impacting radiative transfer. Understanding the behavior of dust is therefore crucial not only for interpreting observations but also for modeling the physical conditions and evolution of astrophysical systems.

While dust is generally assumed to be perfectly coupled to the gas, this assumption breaks down either for larger dust-to-gas ratios or for larger grain sizes. Protoplanetary disks are the astrophysical environment where dust dynamics has been most thoroughly studied. Here, the uncoupling of dust and gas is essential to understand the evolution of the disk and the early stages of planet formation (Birnstiel 2024). The radial drift of dust grains inward, one of the most important dynamical mechanisms to understand the evolution of disks, is caused by a slight difference between the rotational velocity of gas and of dust.

But disks are not the only objects in which the dynamics of dust grains play a role. The outflows of AGB stars are also believed to be accelerated by radiation pressure on dust grains, which consequently accelerates gas by aerodynamic drag. The same mechanism is also proposed to fuel outflows from “dusty torus” regions around AGNs (Soliman & Hopkins 2023).

Given the importance of dust as a key component of many astrophysical systems, the inclusion of dust dynamics in modern hydrodynamical simulation codes has seen significant growth in recent years. One widely used approach involves post-processing gas-only simulations to study the behavior of dust. While this method improves upon the assumption of perfect coupling between dust and gas, it neglects the critical back-reaction of dust on gas motion. This omission can lead to inaccuracies, especially in scenarios with substantial dust-to-gas ratios or large grain sizes, where dust can significantly influence the gas dynamics. Self-consistent models of dust dynamics have been implemented in Eulerian codes, both as tracer particles, such as in AREPO (McKinnon et al. 2018), or as additional fluids, like in ATHENA++ (Huang & Bai 2022) and FARGO3D (Benítez-Llambay et al. 2019). In Smoothed Particle Hydrodynamics (SPH) frameworks, the model that has achieved the most success in simulating dust dynamics is the One-Fluid model developed by (Laibe & Price 2014a). This approach represents a significant departure from traditional multi-fluid methods by describing the motion of a gas-dust mixture self-consistently, rather than treating gas and dust as separate fluids coupled via aerodynamic drag. A notable implementation of the One-Fluid model, in the terminal velocity approximation limit, is in the SPH code PHANTOM (Price et al. 2018). While this model has been largely successful at simulating the dynamics of small dust grains, it still presents challenges when simulating large dust grains. Most notably, (Laibe & Price 2014b) observed that interpenetrating dust flows of weakly coupled grains can violate the assumptions underlying the fluid approximation for a pressureless dust model, causing the code to crash. More subtly, the inability of a pressureless fluid to capture the diffusive nature of a multi-valued velocity distribution leads to artificial clumping of dust, which can produce numerical failures or result in an overestimate of planetesimal formation timescales. This has led to the One-Fluid model to be used mostly in the terminal velocity approximation, where dust is assumed to couple efficiently to the gas on timescales smaller than the numerical integration timestep. One possible solution to this long-standing issue of simulating dust dynamics is the introduction of some form of dust pressure. For instance, Lynch et al. (2024) found that a treatment of dust pressure is needed in order to correctly simulate dust dynamics in debris disks.

In this paper, we will describe the implementation of the full (i.e., not in the terminal velocity approximation) One-Fluid model for dust dynamics in OpenGadget3. In particular, we will focus on improvements on the original implementation derived in (Laibe & Price 2014b), such as time-dependent artificial viscosity and conductivity, and a model of dust diffusion via a dust pressure-like term to treat clumping in simulations of weakly coupled dust. While the limitations of the full One-Fluid model are still present in the most extreme cases, we show how these improvements can make the One-Fluid model suitable for a larger set of astrophysical applications compared to the terminal velocity approximation.

The paper is organized as follows. In Section 2 we show how the One-Fluid model is derived from the equations describing the motion of dust and gas separately and we show the discretization of the continuum equations in the SPH framework, including the improved viscosity and conductivity terms and the dust-pressure diffusion model. In Section 3, we show the performance of the code on a suite of standard numerical tests, as well as on more complex simulations such as the Cold Keplerian Disk and the Kelvin–Helmholtz instability. Finally, in Section 4, we present our conclusions and outline potential avenues for future development.

2 Methods

2.1 The One-Fluid model

Our implementation of dust dynamics relies and improves on the One-Fluid model introduced in Laibe & Price (2014a). The equations governing gas dynamics in the presence of dust closely resemble the standard Euler equations, with additional terms in the momentum and energy conservation equations to account for the drag exerted by dust on the gas $${\partial }_t{\rho }_g+\nabla ·\left({\rho }_g{\mathit{v}}_g\right)=0,$$(1) $${\partial }_t{\mathit{v}}_g+({\mathit{v}}_g·\nabla ){\mathit{v}}_g=-\frac{\nabla P}{{\rho }_g}+\frac{K}{{\rho }_g}({\mathit{v}}_d-{\mathit{v}}_g)+\mathit{f},$$(2) $${\partial }_{t}u+({\mathit{v}}_g·\nabla )u=-\frac{P}{{\rho }_g}(\nabla ·{\mathit{v}}_g)+\frac{K}{{\rho }_g}{({\mathit{v}}_d-{\mathit{v}}_g)}^2+\Lambda ,$$(3) where u and P are the internal energy and the pressure of the gas, respectively, and K is the drag coefficient determining how strongly the gas and the dust are coupled. For simplicity, we will ignore any external force acting on the fluid, f, as well as any external heating, Λ.

The equations for dust are essentially describing a pressureless fluid, with an additional term in the momentum equation describing the drag of gas onto the dust, similar and opposite to the one in Eq. (2) $${\partial }_t{\rho }_d+\nabla ·\left({\rho }_d{\mathit{v}}_d\right)=0,$$(4) $${\partial }_t{\mathit{v}}_d+({\mathit{v}}_d·\nabla ){\mathit{v}}_d=-\frac{K}{{\rho }_d}({\mathit{v}}_d-{\mathit{v}}_g)+\mathit{f}.$$(5) As in the case of the gas momentum equation, we will ignore the term f as well as any effect of gas pressure on the dust.

Laibe & Price (2014a) show that, by introducing a new and convenient set of variables, the equations describing the motion of gas and dust can be combined into a set of equations describing the motion of the mixture of both fluids. In particular, the following variables are used $$\rho ={\rho }_g+{\rho }_d,$$(6) $$ϵ=\frac{{\rho }_d}{{\rho }_g+{\rho }_d},$$(7) $$\mathit{v}=\frac{{\rho }_g{\mathit{v}}_g+{\rho }_d{\mathit{v}}_d}{{\rho }_g+{\rho }_d},$$(8) $$\Delta v=v_d-v_g,$$(9) which described, respectively, the total density (ρ), the fraction of dust density compared to the total density (ε), the barycenter velocity of the mixture (v), and the velocity difference between the dust and the gas phases (Δ v). By using these variables, the Eqs. (1)–(5) can be rewritten as $$\frac{d\rho }{dt}=-\rho (\nabla ·\mathit{v})$$(10) $$\frac{dϵ}{dt}=-\frac{1}{\rho }\nabla ·\left[ϵ(1-ϵ)\rho \Delta \mathit{v}\right]$$(11) $$\frac{d\mathit{v}}{dt}=-\frac{\nabla P}{\rho }-\frac{1}{\rho }\nabla ·\left[ϵ(1-ϵ)\rho \Delta \mathit{v}\Delta \mathit{v}\right]$$(12) $$\frac{d\Delta \mathit{v}}{dt}=\frac{\nabla P}{\rho (1-ϵ)}-\frac{\Delta \mathit{v}}{t_s}-(\Delta \mathit{v}·\nabla )\mathit{v}+\Delta \mathit{v}(\Delta \mathit{v}·\nabla )ϵ+(2ϵ-1)(\Delta \mathit{v}·\nabla )\Delta \mathit{v}$$(13) $$\frac{du}{dt}=-\frac{P}{\rho (1-ϵ)}\nabla ·(\mathit{v}-ϵ\mathit{\Delta }\mathit{v})+ϵ(\mathit{\Delta }\mathit{v}·\nabla )u+ϵ\frac{\mathit{\Delta }{\mathit{v}}^2}{t_s},$$(14) where we already introduced the Lagrangian, or co-moving, temporal derivative $$\frac{d}{dt}=\frac{\partial }{\partial t}+(v·\nabla ),$$(15) and the stopping time t_s, which is related to the drag coefficient by $$t_s=\frac{\rho ϵ(1-ϵ)}{K}.$$(16) Notice that Eq. (13) is slightly different from the form presented in Laibe & Price (2014a) due to an error in the original derivation.

Although Eq. (13) is the most compact form for the evolution of Δ v, here we present two alternative formulations. The first formulation is the one derived in Lebreuilly et al. (2019) $$\begin{array}{rl}\frac{d\mathbf{\Delta }\mathit{v}}{dt}=& \frac{\mathrm{\nabla }P}{\rho (1-ϵ)}-\frac{\mathbf{\Delta }\mathit{v}}{t_s}-(\mathbf{\Delta }\mathit{v}\cdot \mathrm{\nabla })\mathit{v}+\frac{1}{2}\mathrm{\nabla }[(2ϵ-1)\mathbf{\Delta }\mathit{v}\cdot \mathbf{\Delta }\mathit{v}]\\ & +(1-ϵ)\mathbf{\Delta }\mathit{v}\times [\mathrm{\nabla }\times (1-ϵ)\mathbf{\Delta }\mathit{v}]-ϵ\mathbf{\Delta }\mathit{v}\times [\mathrm{\nabla }\times ϵ\mathbf{\Delta }\mathit{v}].\end{array}$$(17) The second version will be the one that we will use for the later SPH discretization $$\begin{array}{rl}& \frac{d\mathrm{\Delta }v}{dt}=\frac{\mathrm{\nabla }P}{\rho (1-ϵ)}-\frac{\mathrm{\Delta }v}{t_s}-(\mathrm{\Delta }v\cdot \mathrm{\nabla })v+\frac{1}{2}\mathrm{\nabla }[(2ϵ-1)\mathrm{\Delta }v\cdot \mathrm{\Delta }v]\\ & +\mathrm{\Delta }v(\mathrm{\Delta }v\cdot \mathrm{\nabla })ϵ+(2ϵ-1)(\mathrm{\Delta }v\cdot \mathrm{\nabla })\mathrm{\Delta }v-|\mathrm{\Delta }v{|}^2\mathrm{\nabla }ϵ-\frac{1}{2}(2ϵ-1)\mathrm{\nabla }|\mathrm{\Delta }v{|}^2\end{array}$$(18) While all these three formulations seem very different, they are actually equivalent to one another, and each may be useful under particular circumstances.

We note here that Eqs. (17) and (18) share a common property: under scalar product with Δv itself, the second row of both equations cancels out. It is easier to see in Eq. (17), where all terms in the second row are the result of cross-products with Δv, and thus perpendicular to it, but the same holds for Eq. (18) as well. This property will be crucial in deriving a suitable SPH discretization of Eq. (18).

2.2 SPH implementation of the One-Fluid model

In order to computationally solve the One-Fluid model presented above, we need to find a suitable discretization that conserves mass, momentum, and energy. We will closely follow the derivation of Laibe & Price (2014b) (hereafter LP14b) for the well-known Smoothed Particle Hydrodynamics (SPH) framework, but corrected for the new terms in Equation (18). We will also add time-dependent versions of the artificial viscosity and conductivity, a new dust diffusion term, and the explicit form of the equations for a parametrization of the dust fraction.

2.2.1 Density estimate

In SPH, the fluid is modeled as a collection of particles, each representing a portion of its mass. Any hydrodynamical quantity at a point in space r is then represented by a weighted average of the contribution of all neighboring particles, weighted by a kernel 𝒲(r_a, h_a), depending on the distance r_a between the particle a and the point r and the smoothing length h_a of the particle. The kernel should satisfy several properties: compact support, continuity, radial symmetry, and the following normalization constraint $${\int }_V\mathcal{W}(|{\mathit{r}}^{\text{'}}-{\mathit{r}}_a|,h_a)d{V}^{\text{'}}=1.$$(19) Some common choices for SPH kernels include the cubic (Monaghan & Lattanzio 1985) and quintic spline (Morris 1996) and the Wedland C2/C4 and C6 kernels (Wendland 1995; Dehnen & Aly 2012).

The density of a given particle is then determined by a weighted sum over the particle’s neighbors $$\rho \left({\mathit{r}}_a\right)={\sum }_{b\in NGB}{m}_b\mathcal{W}(|{\mathit{r}}_a-{\mathit{r}}_b|,h_a),$$(20) where m_b is the mass of the b particle. Eq. (20) can be shown to be a solution to the continuity equation Eq. (10). The variable smoothing length h_a of each particle is computed to keep the mass inside the kernel volume of each particle fixed. In three dimensions, h_a is computed to fulfill the equation $$\frac{4}{3}\pi h_a^3{\rho }_a=N_{\mathrm{ngb}}{m}_a$$(21) where N_ngb is the chosen number of neighbors.

2.2.2 Differential operators

To obtain a suitable discretization of the One-Fluid model equations we will employ the following SPH differential operators (Price 2012) $$\nabla ·\mathit{A}\underset{\mathrm{SPH}}{=}{\rho }_a{\sum }_b{m}_b[\frac{{\mathit{A}}_a·\nabla {\mathcal{W}}_{ab}\left(h_a\right)}{{\Omega }_a{\rho }_a^2}+\frac{{\mathit{A}}_b·\nabla {\mathcal{W}}_{ab}\left(h_b\right)}{{\Omega }_b{\rho }_b^2}],$$(22) $$(\mathit{A}·\nabla )f\underset{\mathrm{SPH}}{=}\frac{1}{{\rho }_a{\Omega }_a}{\sum }_b{m}_b(f_b-f_a){\mathit{A}}_a·\nabla {\mathcal{W}}_{ab}\left(h_a\right),$$(23) $$(\mathit{A}·\nabla )\mathit{B}\underset{\mathrm{SPH}}{=}\frac{1}{{\rho }_a{\Omega }_a}{\sum }_b{m}_b({\mathit{B}}_b-{\mathit{B}}_a){\mathit{A}}_a·\nabla {\mathcal{W}}_{ab}\left(h_a\right),$$(24) where ${\mathcal{W}}_{ab}\left(h_a\right)=\mathcal{W}(|{\mathit{r}}_a-{\mathit{r}}_b|,h_a),{\Omega }_a=(1+\frac{h_a}{{\rho }_{a}d}\frac{\partial {\rho }_a}{\partial h_a})$, with d the number of dimensions, and A, B and f are arbitrary vector and scalar fields.

We will apply these operators and, step by step, ensure that mass, total momentum, and total energy are conserved.

2.2.3 Mass conservation

The first requirement that the SPH discretization should meet is the conservation of mass. The total mass of the mixture is conserved by construction by the SPH framework, and since there is no dust-mass evolution mechanism taken into account in the One-Fluid model (e.g., sputtering, erosion, or astration), the individual dust and gas masses should also be conserved. Dust mass conservation requires that $$\frac{d{M}_d}{dt}\stackrel{=}{\mathrm{SPH}}{\sum }_a{m}_a\frac{d{ϵ}_a}{dt}=0.$$(25) By applying the divergence operator Eqs. (22)–(11) we obtain the following SPH discretization $$\begin{array}{rl}\frac{d{ϵ}_a}{dt}=-\sum _b{m}_b[\frac{{ϵ}_a(1-{ϵ}_a)}{{\mathrm{\Omega }}_a{\rho }_a}& \mathrm{\Delta }{\mathit{v}}_a\cdot {\mathrm{\nabla }}_a{\mathcal{W}}_{ab}\left(h_a\right)\\ & +\frac{{ϵ}_b(1-{ϵ}_b)}{{\mathrm{\Omega }}_b{\rho }_b}\mathrm{\Delta }{\mathit{v}}_b\cdot {\mathrm{\nabla }}_a{\mathcal{W}}_{ab}\left(h_b\right)]\end{array}$$(26) Using the anti-symmetry of the kernel gradient ∇ 𝒲_ab(h_a) = −∇ 𝒲_ab(h_b), it is straightforward to show that the total mass of dust is conserved.

2.2.4 Momentum conservation

While dust and gas can exchange momentum, the total momentum of the mixture should be conserved. A suitable SPH discretization of Eq. (12) is constrained by $${\sum }_a{m}_a\frac{d{\mathit{v}}_a}{dt}=0.$$(27) Compared to the gas-only momentum conservation equation, the One-Fluid model version, Eq. (12), has a first term that represents the usual pressure gradient acceleration and a new term specific to the One-Fluid model. The pressure gradient term can be solved as usual in SPH $${\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{\nabla P}=-{\sum }_b{m}_b[\frac{P_a}{{\Omega }_a{\rho }_a^2}\nabla {\mathcal{W}}_{ab}\left(h_a\right)+\frac{P_b}{{\Omega }_b{\rho }_b^2}\nabla {\mathcal{W}}_{ab}\left(h_b\right)].$$(28) The Eq. (18) also presents a similar term, which we discretize in the same manner $$\begin{array}{rl}{\left(\frac{d\mathbf{\Delta }{\mathit{v}}_a}{dt}\right)}_{\mathrm{\nabla }P}& -\frac{1}{1-{ϵ}_a}{\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{\mathrm{\nabla }P}\\ & =\frac{1}{1-{ϵ}_a}\sum _b{m}_b[\frac{P_a}{{\mathrm{\Omega }}_a{\rho }_a^2}\mathrm{\nabla }{\mathcal{W}}_{ab}\left(h_a\right)+\frac{P_b}{{\mathrm{\Omega }}_b{\rho }_b^2}\mathrm{\nabla }{\mathcal{W}}_{ab}\left(h_b\right)]\end{array}$$(29) Back to the momentum conservation, the second term of Eq. (12) can be discretized by applying the differential operator Eq. (22) $$\frac{-b\pm \sqrt{b^2-4ac}}{2a}$$(30) Both Eqs. (28) and (30) conserve momentum as required by Eq. (27) and can thus be used as SPH discretizations of the corresponding continuum equations.

2.2.5 Energy conservation

The mixture of dust and gas should also conserve its total energy, and from this requirement, we will derive the SPH discretizations for both the internal energy u and the remaining terms in the Δ v evolution.

The total energy of the mixture is the sum of the total energy of the gas phase (kinetic and internal energy) and the total energy of the dust phase (which only has a kinetic component, as it is modeled as a pressureless fluid) $$\frac{-b\pm \sqrt{b^2-4ac}}{2a}$$(31) where in the second row we introduced the One-Fluid model variables and in the final row we discretized the integral as a sum over all particles.

By computing the time derivative of Eq. (31) we obtain the constraint that the SPH discretizations for Eqs. (18) and (14) should satisfy $$\begin{array}{rl}\frac{dE}{dt}=& \sum _a{m}_a[{\mathit{v}}_a\cdot \frac{d{\mathit{v}}_a}{dt}+{ϵ}_a(1-{ϵ}_a)\mathrm{\Delta }{\mathit{v}}_a\cdot \frac{d\mathbf{\Delta }{\mathit{v}}_a}{dt}\\ & +((1-2{ϵ}_a)\frac{\mathbf{\Delta }{\mathit{v}}_a^2}{2}-u_a)\frac{d{ϵ}_a}{dt}+(1-{ϵ}_a)\frac{d{u}_a}{dt}]=0\end{array}$$(32) The first term we want to discretize is the PdV term in Eq. (14). By combining all terms in which pressure is present, we obtain the following constraint from energy conservation $$\begin{array}{rl}\sum _a{m}_a(1& -{ϵ}_a){\left(\frac{d{u}_a}{dt}\right)}_{PdV}\\ =& -\sum _a{m}_a[{\mathit{v}}_a\cdot {\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{\mathrm{\nabla }P}+{ϵ}_a(1-{ϵ}_a)\mathbf{\Delta }{\mathit{v}}_a\cdot {\left(\frac{d\mathbf{\Delta }{\mathit{v}}_a}{dt}\right)}_{\mathrm{\nabla }P}]\\ =& -\sum _a{m}_a[({\mathit{v}}_a-{ϵ}_a\mathbf{\Delta }{\mathit{v}}_a)\cdot {\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{\mathrm{\nabla }P}]\\ =& \sum _a{m}_a\sum _b{m}_b\frac{P_a}{{\mathrm{\Omega }}_a{\rho }_a^2}({\mathit{v}}_a-{ϵ}_a\mathbf{\Delta }{\mathit{v}}_a)\cdot \mathrm{\nabla }{\mathcal{W}}_{ab}\left(h_a\right)\\ & +\sum _a{m}_a\sum _b{m}_b\frac{P_b}{{\mathrm{\Omega }}_b{\rho }_b^2}({\mathit{v}}_a-{ϵ}_a\mathbf{\Delta }{\mathit{v}}_a)\cdot \mathrm{\nabla }{\mathcal{W}}_{ab}\left(h_b\right)\end{array}$$(33) where in the last step we have substituted ${\left(\frac{d{v}_a}{dt}\right)}_{\nabla P}$ from Eq. (28). Swapping the summation indices of the second term and using the anti-symmetry of the kernel, we obtain $${\left(\frac{d{u}_a}{dt}\right)}_{PdV}=\frac{P_a}{{\Omega }_a{\rho }_a^2(1-{ϵ}_a)}{\sum }_b{m}_b[({\mathit{v}}_a-{ϵ}_a\mathit{\Delta }{\mathit{v}}_a-{\mathit{v}}_b+{ϵ}_b\mathit{\Delta }{\mathit{v}}_b)·\nabla W_{ab}\left(h_a\right)].$$(34) The second term in Eq. (14) is also constrained from energy conservation. By comparing the terms in Eq. (32) involving the internal energy $${\sum }_a{m}_a(1-{ϵ}_a){\left(\frac{d{u}_a}{dt}\right)}_{(\Delta v·\nabla )u}={\sum }_a{m}_a{u}_a\frac{d{ϵ}_a}{dt}.$$(35) By substituting Eq. (26) and swapping indices in the double summation, we obtain $${\left(\frac{d{u}_a}{dt}\right)}_{(\mathit{\Delta }v·\nabla )u}=\frac{{ϵ}_a}{{\Omega }_a{\rho }_a}{\sum }_b{m}_b(u_b-u_a)\mathit{\Delta }{\mathit{v}}_a·\nabla {\mathcal{W}}_{ab}\left(h_a\right).$$(36) We will now deal with Eq. (18), describing the evolution of Δ v. The first term, similar to the pressure gradient force term, was already dealt with and its SPH discretization is described in Eq. (29). The second term in Equation (18) relating to the drag will be integrated implicitly and addressed later, together with the last term in Equation (14). The SPH discretization for the third term in Eq. (18) can be derived from energy conservation, by equating the first two terms in Eq. (32): $${\sum }_a{m}_a{ϵ}_a(1-{ϵ}_a)\mathit{\Delta }{\mathit{v}}_a·{\left(\frac{d\mathit{\Delta }{\mathit{v}}_a}{dt}\right)}_{-(\mathit{\Delta }\mathit{v}·\mathit{\nabla })\mathit{v}}=-{\sum }_a{m}_a{\mathit{v}}_a·{\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{OFM}.$$(37) Also in this case, we simply substitute Eq. (30) and swap summation indices to obtain $${\left(\frac{d\mathit{\Delta }{\mathit{v}}_a}{dt}\right)}_{-(\mathit{\Delta }\mathit{v}·\nabla )\mathit{v}}=-\frac{1}{{\Omega }_a{\rho }_a}{\sum }_b{m}_b[({\mathit{v}}_b-{\mathit{v}}_a)\mathit{\Delta }{\mathit{v}}_a·\nabla {\mathit{W}}_{ab}\left(h_a\right)].$$(38) The fourth term in Eq. (18) is also constrained by energy conservation, but comparing now the second and third terms in Eq. (32) $${\sum }_a{m}_a{ϵ}_a(1-{ϵ}_a)\mathit{\Delta }{\mathit{v}}_a·{\left(\frac{d\mathit{\Delta }{\mathit{v}}_a}{dt}\right)}_4=-{\sum }_a{m}_a(1-2{ϵ}_a)\frac{\mathit{\Delta }{\mathit{v}}_a^2}{2}\frac{d{ϵ}_a}{dt}.$$(39) By substituting Eq. (26) and swapping summation indices we obtain $$\begin{array}{rl}{\left(\frac{d\mathbf{\Delta }{\mathit{v}}_a}{dt}\right)}_4=& \frac{1}{2{\mathrm{\Omega }}_a{\rho }_a}\sum _b{m}_b[(1-2{ϵ}_a)\mathbf{\Delta }{\mathit{v}}_a^2-\\ & -(1-2{ϵ}_b)\mathbf{\Delta }{\mathit{v}}_b^2]\mathrm{\nabla }{\mathit{W}}_{ab}\left(h_a\right)\end{array}$$(40) It is crucial to note here that the remaining terms in Eq. (18) cannot be constrained by energy conservation, as they cancel out when multiplied scalarly by Δv, as noted when they were first introduced (see Section 2.1). Hence, their SPH discretization will be obtained directly by use of SPH differential operators (Eqs. (22)–(24)), and we will only later check that energy is conserved. We propose the following SPH discretizations $$\mathit{\Delta }\mathit{v}(\mathit{\Delta }\mathit{v}·\nabla )ϵ\underset{\mathrm{SPH}}{=}\frac{\mathit{\Delta }{\mathit{v}}_a}{{\Omega }_a{\rho }_a}{\sum }_b{m}_b({ϵ}_b-{ϵ}_a)\mathit{\Delta }{\mathit{v}}_a·\nabla {\mathcal{W}}_{ab}\left(h_a\right),$$(41) $$(2ϵ-1)(\Delta \mathit{v}·\nabla )\Delta \mathit{v}\underset{\mathrm{SPH}}{=}\frac{2{ϵ}_a-1}{{\Omega }_a{\rho }_a}{\sum }_b{m}_b(\Delta {\mathit{v}}_b-\Delta {\mathit{v}}_a)\Delta {\mathit{v}}_a·\nabla {\mathcal{W}}_{ab}\left(h_a\right),$$(42) $$-{\left|\mathit{\Delta }v\right|}^2\nabla {ϵ}_{\mathrm{SPH}}^{=}-\frac{{\left|\mathit{\Delta }{\mathit{v}}_a\right|}^2}{{\Omega }_a{\rho }_a}{\sum }_b{m}_b({ϵ}_b-{ϵ}_a)\nabla {\mathcal{W}}_{ab}\left(h_a\right)$$(43) $$-\frac{1}{2}(2ϵ-1)\nabla {\left|\mathit{\Delta }v\right|}^2\underset{\mathrm{SPH}}{=}-\frac{2{ϵ}_a-1}{{\Omega }_a{\rho }_a}{\sum }_b{m}_b(\mathit{\Delta }{\mathit{v}}_b-\mathit{\Delta }{\mathit{v}}_a)·\mathit{\Delta }{\mathit{v}}_a\nabla {\mathcal{W}}_{ab}\left(h_a\right)$$(44) It can be shown that these SPH discretizations conserve energy if inserted into Eq. (32).

One quantity that this set of SPH discretizations fails to conserve is the dust momentum in the absence of drag forces. However, deviations from conservation are only of order O(Δv²), and thus negligible for strong drag. In the weak drag case, we also observed little deviation from conservation in the benchmark tests we show in Section 3. In general, the conservation of dust momentum for zero drag is a property that has still not been enforced by previous SPH discretizations of the One-Fluid model, and could be an area of future improvement.

2.3 Timestepping

OpenGadget3 employs, for integrating in time all hydrodynamical variables, a leapfrog scheme in a Kick-Drift-Kick (KDK) form. This scheme ensures second-order accuracy for all variables and is implemented as in Springel (2005).

We indicate with W=(ε, v, Δ v, u) the vector of primitive quantities of the One-Fluid model, except for density, which is computed by the kernel density estimate in Eq. (20). The W vector and the positions of the particles, x, are evolved using predicted quantities to compute the increments of the second kick $$\begin{array}{rl}& {\mathit{W}}_{n+1/2}={\mathit{W}}_n+\frac{1}{2}\mathrm{\Delta }t{\left(\frac{d\mathit{W}}{dt}\right)}_n,\\ & {\mathit{W}}_{n+1}^{\ast }={\mathit{W}}_n+\mathrm{\Delta }t{\left(\frac{d\mathit{W}}{dt}\right)}_n,\\ & {\mathit{x}}_{n+1}={\mathit{x}}_n+\mathrm{\Delta }t{\mathit{v}}_{n+1/2},\\ & {\mathit{W}}_{n+1}={\mathit{W}}_{n+1/2}+\frac{1}{2}\mathrm{\Delta }t{\left(\frac{d\mathit{W}}{dt}\right)}_{n+1}^{\ast }.\end{array}$$(45)

2.4 Implicit integration of the drag term

The last two terms in Eqs. (18) and (14) can be numerically integrated both explicitly and implicitly¹. As shown in LP14b, explicit integration is subject to the timestep requirement Δ t<t_s, which becomes prohibitive for strong drag. Implicit integration, instead, does not require such a timestep requirement and can be implemented via operator splitting.

We show here the procedure for the first kick in the KDK scheme, but the same is applied also to compute the predicted quantities during the drift and for the second kick. By separating the drag term from all the other contributions, we can rewrite Eqs. (18) and (14) as $$\{\begin{array}{c}\frac{d\mathit{\Delta }\mathit{v}}{dt}={\left(\frac{d\mathit{\Delta }\mathit{v}}{dt}\right)}_0-\frac{\mathit{\Delta }\mathit{v}}{t_s},\hfill \\ \frac{du}{dt}={\left(\frac{du}{dt}\right)}_0+ϵ\frac{\mathit{\Delta }{\mathit{v}}^2}{t_s}.\hfill \end{array}$$(46) To use operator splitting, we first solve the corresponding homogeneous system by ignoring the drag terms and using the SPH discretization derived above $$\{\begin{array}{c}{\overline{\mathit{\Delta }v}}_{n+1/2}=\mathit{\Delta }{\mathit{v}}_n+\frac{1}{2}\Delta t{\left(\frac{d\mathit{\Delta }\mathit{v}}{dt}\right)}_{n,0},\hfill \\ {\overline{u}}_{n+1/2}=u_n+\frac{1}{2}\Delta t{\left(\frac{du}{dt}\right)}_{n,0}.\hfill \end{array}$$(47) The Ordinary Differential Equation corresponding to the system in Eq. (46) can be solved analytically as $$\{\begin{array}{c}\frac{d\mathit{\Delta }\mathit{v}}{dt}=-\frac{\mathit{\Delta }\mathit{v}}{t_s},\hfill \\ \frac{du}{dt}=ϵ\frac{\mathit{\Delta }{\mathit{v}}^2}{t_s},\hfill \end{array}\to \{\begin{array}{c}\mathit{\Delta }\mathit{v}\left(t\right)=\mathit{\Delta }{\mathit{v}}_0{e}^{-\frac{t}{t_s}},\hfill \\ u\left(t\right)=u_0-\frac{1}{2}ϵ\mathit{\Delta }{\mathit{v}}_0^2[e^{-\frac{2t}{t_s}}-1].\hfill \end{array}$$ By discretizing this solution to half a timestep as in the first kick, and by using the solutions to the homogeneous system as initial conditions, we obtain the final result $$\{\begin{array}{c}\mathit{\Delta }{\mathit{v}}_{n+1/2}=[\mathit{\Delta }{\mathit{v}}_n+\frac{1}{2}\Delta t{\left(\frac{d\mathit{\Delta }\mathit{v}}{dt}\right)}_{n,0}]e^{-\frac{\Delta t}{2t_s}}\hfill \\ u_{n+1/2}=u_n+\frac{1}{2}\Delta t{\left(\frac{du}{dt}\right)}_{n,0}-\frac{1}{2}{ϵ}_{n+1/2}{\overline{\mathit{\Delta }v}}_{n+1/2}^2[e^{-\frac{\Delta t}{t_s}}-1]\hfill \end{array}$$(48)

2.5 Time-dependent artificial viscosity

In order to dampen post-shock oscillations, prevent particle interpenetration, and treat contact discontinuities, SPH frameworks often employ artificial viscosity mechanisms (Price 2008). Our artificial viscosity implementation follows the steps of LP14b, with some modifications and improvements that make the implementation time-dependent. We show both the conservative (C) and the non-conservative (NC) approaches to artificial viscosity, as shown in LP14b, an artificial dissipation term for Δv, and a simple unifying scheme to combine both conservative and non-conservative viscosities.

2.5.1 Conservative artificial viscosity

A simple momentum-conserving formulation for artificial viscosity can be written as $${\left(\frac{d\mathit{v}}{dt}\right)}_{AV}^C={\sum }_b{m}_b(1-{\overline{ϵ}}_{ab})v_{ab}\overline{\nabla }{\overline{\mathcal{W}}}_{ab},$$(49) where the overline stands for the mean value between particle a and particle b of a given quantity. For instance, ${\overline{ϵ}}_{ab}=\frac{1}{2}({ϵ}_a+{ϵ}_b)$. A possible choice for the viscosity v_ab is the one employed by LP14b $$v_{ab}=\frac{v_{ab}^{sig}}{{\overline{\rho }}_{ab}}{\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab},$$(50) where ${\mathit{v}}_{ab}^g={\mathit{v}}_a^g-{\mathit{v}}_b^g$ is the difference between the gas velocities of the particles, ${\hat{\mathit{r}}}_{ab}$ is the unit vector pointing from particle a to particle b and $v_{ab}^{\text{sig}}$ is the signal velocity.

The signal velocity is used to switch on and off the artificial viscosity depending on whether the gas in the particles is approaching or distancing $$v_{ab}^{sig}=\{\begin{array}{cc}\frac{1}{2}(c_{s,a}+c_{s,b}-\beta {\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab})\hfill & {\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab}\le 0,\hfill \\ 0\hfill & {\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab}>0,\hfill \end{array}$$(51) where c_s is the sound speed of the particle and β=3 following Beck et al. (2016).

In our implementation, we adopt two additional terms in the formulation of viscosity $$v_{ab}=\frac{v_{ab}^{sig}}{{\overline{\rho }}_{ab}}{\alpha }_{ab}{\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab},$$(52) where α_ab is the symmetrized viscosity coefficient, and is used to reduce artificial viscosity when it is not needed, i.e., away from shocks. This framework was already implemented (Dolag et al. 2005; Beck et al. 2016) and extensively tested in OpenGadget3 by Marin-Gilabert et al. (2022), we only added the precaution of using the gas velocity instead of the mixture velocity when computing α_ab.

In order to make this formulation conservative, we need to add matching terms in the equations of Δv and u. In particular, for Δv we simply have $${\left(\frac{d\mathbf{\Delta }{\mathit{v}}_a}{dt}\right)}_{AV}^C=\frac{1}{1-{ϵ}_a}\sum _b{m}_b(1-{\overline{ϵ}}_{ab})v_{ab}{\overline{\mathrm{\nabla }\mathcal{W}}}_{ab},$$(53) while the term in the internal energy equation is constrained from the equivalence already shown in Eq. (33) $${\sum }_a{m}_a(1-{ϵ}_a){\left(\frac{d{u}_a}{dt}\right)}_{AV}^C=-{\sum }_a{m}_a[{\mathit{v}}_a^g·{\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{AV}^C],$$(54) which gives, after inserting Eq. (49), expanding $\overline{\nabla {\mathcal{W}}_{ab}}$ and swapping summation indices $${\left(\frac{d{u}_a}{dt}\right)}_{AV}^C=-\frac{1}{2(1-{ϵ}_a)}{\sum }_b{m}_b(1-{\overline{ϵ}}_{ab})v_{ab}{\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab}{\overline{F}}_{ab},$$(55) where the kernel gradient has been expanded into direction and magnitude: ∇𝒲_ab(h_a)=r̂_ab F_ab(h_a).

As described by LP14b, only using Eqs. (49), (53) and (55) leaves some oscillations in Δv in post-shock regions, which can be treated by using an additional dissipative term in the Δ v equation. In particular, a momentum and energy-conserving formulation is $${\left(\frac{d\mathit{\Delta }{\mathit{v}}_a}{dt}\right)}_{AD}=\frac{1}{{ϵ}_a(1-{ϵ}_a)}{\sum }_b\frac{v_{sig}^{\mathit{\Delta }\mathit{v}}}{{\overline{\rho }}_{ab}}{\overline{ϵ}}_{ab}(1-{\overline{ϵ}}_{ab})\mathit{\Delta }{\mathit{v}}_{ab}·{\widehat{\mathit{r}}}_{ab}\overline{\nabla {\mathcal{W}}_{ab}}.$$(56) In order to conserve energy, a matching term in the internal energy equation must be added $${\left(\frac{d{u}_a}{dt}\right)}_{AD}=-\frac{1}{2(1-{ϵ}_a)}{\sum }_b{m}_b\frac{v_{sig}^{\Delta v}}{{\overline{\rho }}_{ab}}{\overline{ϵ}}_{ab}(1-{\overline{ϵ}}_{ab}){(\Delta {\mathit{v}}_{ab}·{\widehat{\mathit{r}}}_{ab})}^2{\overline{F}}_{ab}.$$(57) The signal velocity used here is different from the one used in the artificial viscosity. In particular, LP14b propose the following definition $$v_{ab}^{sig,\Delta v}=\{\begin{array}{cc}-{\alpha }_{\Delta v}\mathit{\Delta }{\mathit{v}}_{ab}·{\widehat{\mathit{r}}}_{ab}\hfill & \mathit{\Delta }{\mathit{v}}_{ab}·{\widehat{\mathit{r}}}_{ab}\le 0,\hfill \\ 0\hfill & \mathit{\Delta }{\mathit{v}}_{ab}·{\widehat{\mathit{r}}}_{ab}>0,\hfill \end{array}$$(58) with α_{Δ v}=1. This concludes the conservative formulation for artificial viscosity. As already found by LP14b, this formulation works well for low Δv, but fails in shocks with large Δv, as can be seen in the DUSTYSHOCK benchmark test, at low drag. For this reason, LP14b introduced an alternative, nonconservative, formulation of artificial viscosity that shows a much better performance on the DUSTYSHOCK test.

2.5.2 Nonconservative artificial viscosity

The nonconservative formulation of artificial viscosity we employed is very similar to the one derived by LP14b, with the only improvement of making it adaptive and time-dependent by using the viscosity in Eq. (52) $${\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{AV}^{NC}={\sum }_b{m}_b{v}_{ab}\overline{\nabla }{\overline{\mathit{W}}}_{ab},$$(59) $${\left(\frac{d\mathit{\Delta }{\mathit{v}}_a}{dt}\right)}_{AV}^{NC}=\frac{1}{1-{ϵ}_a}{\sum }_b\frac{v_{ab}^{sig}}{{\overline{\rho }}_{ab}}\mathit{\Delta }{\mathit{v}}_{ab}·{\widehat{\mathit{r}}}_{ab}\overline{\nabla W_{ab}},$$(60) where the signal velocity employed in the second equation is the one in Eq. (51).

The viscous heating term is also added as in Eq. (55), slightly modified to adjust to the different term in the $\frac{dv}{dt}$ equation $${\left(\frac{d{u}_a}{dt}\right)}_{AV}^{NC}=-\frac{1}{2(1-{ϵ}_a)}{\sum }_b{m}_b{v}_{ab}{\mathit{v}}_{ab}^g·{\widehat{\mathit{r}}}_{ab}{\overline{F}}_{ab}\left(h_a\right).$$(61)

2.5.3 Unified framework

In order to derive a unified framework that could be used under all circumstances, instead of relying on two options to be left to the user, we use the shock indicator R_a to interpolate between the two formulations. The shock indicator is defined as in Cullen & Dehnen (2010) $$R_a=\frac{1}{{\rho }_a}{\sum }_{b}sign(\nabla ·{\mathit{v}}_b)m_b{\mathit{W}}_{ab}.$$(62) For non-shocked regions, R_a ≃ 0, while for shocked regions R_a ≃ ± 1. This indicator is already used to compute the symmetrized viscosity coefficient α_ab, to determine the strength of the shock.

Starting from the shock indicator, we define a new variable, σ_ab, to interpolate between the conservative and the non-conservative formulations of artificial viscosity. σ_ab is defined as a sigmoid function $${\sigma }_{ab}=\frac{1}{1+e^{\alpha (|\overline{R}{|}_{ab}-\beta )}},$$(63) where appropriate values for α and β where found to be α=20 and β=0.25. This makes the sigmoid function steep enough to quickly transition from one formulation to the other. In particular, for shocked regions σ_ab ≃ 0, while for non-shocked regions σ_ab ≃ 1.

The interpolation between conservative and non-conservative formulations is performed as follows: $${\left(\frac{d\mathit{v}}{dt}\right)}_{AV}={\left(\frac{d\mathit{v}}{dt}\right)}_{AV}^C{\sigma }_{ab}+{\left(\frac{d\mathit{v}}{dt}\right)}_{AV}^{NC}(1-{\sigma }_{ab}),$$(64) $${\left(\frac{d\Delta \mathit{v}}{dt}\right)}_{AV}=[{\left(\frac{d\Delta \mathit{v}}{dt}\right)}_{AV}^C+{\left(\frac{d\Delta \mathit{v}}{dt}\right)}_{AD}]{\sigma }_{ab}+{\left(\frac{d\Delta \mathit{v}}{dt}\right)}_{AV}^{NC}(1-{\sigma }_{ab}),$$(65) $${\left(\frac{du}{dt}\right)}_{AV}=[{\left(\frac{du}{dt}\right)}_{AV}^C+{\left(\frac{du}{dt}\right)}_{AD}]{\sigma }_{ab}+{\left(\frac{du}{dt}\right)}_{AV}^{NC}(1-{\sigma }_{ab}).$$(66) This is the framework that will be used throughout all the benchmark tests and simulations.

2.6 Time-dependent artificial conductivity

Artificial conductivity is introduced to treat discontinuities in internal energy (Price 2008). Here we employ its time-dependent formulation introduced in Beck et al. (2016), which is also used in OpenGadget3 (see Marin-Gilabert et al. 2022; Groth et al. 2023), adapted for the One-Fluid model., $${\left(\frac{d{u}_a}{dt}\right)}_{AC}=\frac{1}{1-{ϵ}_a}{\sum }_b\frac{m_b}{{\overline{\rho }}_{ab}}(u_b-u_a){\alpha }_{ab}^C{v}_{ab}^{sig,c}{\overline{F}}_{ab},$$(67) where the signal velocity for the artificial conductivity is defined as $$v_{ab}^{sig,c}=\sqrt{\frac{|P_a-P_b|}{\overline{\rho {(1-ϵ)}_{ab}}}},$$(68) and the conductivity coefficient ${\alpha }_a^C$ is made time-dependent by computing it from the internal energy and its gradient: $${\alpha }_a^C=\frac{h_a}{2}\frac{\left|\nabla u_a\right|}{\left|u_a\right|}.$$(69)

2.7 Dust fraction parametrization

In order to ensure that the dust fraction, ε, remains bounded between 0 and 1, a suitable parametrization is needed. Here we follow (Ballabio et al. 2018) and implement the following parametrization $$ϵ=\frac{s^2}{1+s^2}\to s=\sqrt{\frac{ϵ}{1-ϵ}}.$$(70) The corresponding continuum equation for the variable s can be obtained by substituting the parametrization in Eq. (11) and by using that $$\frac{ds}{dt}=\frac{{(1+s^2)}^2}{2s}\frac{dϵ}{dt},$$(71) which results in $$\frac{ds}{dt}=-\frac{{(1+s^2)}^2}{2\rho s}\nabla ·\left[\frac{s^2}{{(1+s^2)}^2}\rho \mathit{\Delta }\mathit{v}\right].$$(72) Since the s at the denominator could cause numerical issues if a given SPH particle has very little or no dust, we split the s²=s · s inside the divergence operator to obtain the following $$\frac{ds}{dt}=-\frac{{(1+s^2)}^2}{2\rho }\nabla ·\left[\frac{s}{{(1+s^2)}^2}\rho \Delta v\right]-\frac{1}{2}(\Delta v·\nabla )s.$$(73) This last formulation is thus discretized using the SPH operators Eqs. (22)–(23) and rearranged to obtain $$\begin{array}{rl}\frac{d{s}_a}{dt}=-\frac{1}{2}\sum _b{m}_b{s}_b& [\frac{1}{{\mathrm{\Omega }}_a{\rho }_a}\mathbf{\Delta }{\mathit{v}}_a\cdot \mathrm{\nabla }{\mathcal{W}}_{ab}\left(h_a\right)\\ & +\frac{{(1+s_a^2)}^2}{{\mathrm{\Omega }}_b{\rho }_b{(1+s_b^2)}^2}\mathbf{\Delta }{\mathit{v}}_b\cdot \mathrm{\nabla }{\mathcal{W}}_{ab}\left(h_b\right)].\end{array}$$(74) The other equations and corresponding SPH discretizations do not need to be recomputed; we simply replace ε when needed with the parametrization in Eq. (70), in terms of the new variable s.

2.8 Artificial diffusion of dust

One of the consequences of modeling dust as a pressureless fluid is its tendency, when the drag coefficient is sufficiently low, to form clumps of dust that cannot self-regulate and dissipate. This leads to large gradients in ε, which will lead, eventually, to code crashes. In order to alleviate this tendency, we tried including some form of dust diffusion in the One Fluid model. Our first attempt was to tackle directly the evolution of ρ_d by adding a Gradient Diffusion flux to its continuity equation. While this may certainly be a valid approach in Eulerian simulations, we found this approach not suitable for our SPH framework, as it would require adding a term in the continuity equation for the total mixture density, and hence a mass transfer term between SPH particles.

We thus resorted to a different approach, the inclusion of a “pressure-like” term in the dust momentum conservation equation, as described by Klahr & Schreiber (2021) and Binkert (2023) $$\frac{\partial }{\partial t}\left({\rho }_d{\mathit{v}}_d\right)+\nabla ·\left({\rho }_d{\mathit{v}}_d{\mathit{v}}_d\right)=-\nabla \left(\frac{1}{3}\frac{D}{t_s}{\rho }_d\right)-\frac{{\rho }_d}{t_s}({\mathit{v}}_d-{\mathit{v}}_g).$$(75) Translating this approach in the One Fluid model equations, we obtain the following continuum equations $${\left(\frac{d\mathit{v}}{dt}\right)}_{diff}=-\frac{DK}{3\rho {(1-ϵ)}^2}\nabla ϵ=-\frac{2DKs}{3\rho }\nabla s,$$(76) $${\left(\frac{d\mathit{\Delta }\mathit{v}}{dt}\right)}_{diff}=-\frac{DK}{3\rho ϵ{(1-ϵ)}^2}\nabla ϵ=-\frac{2sDK}{3\rho ϵ}\nabla s,$$(77) where we already provide the formulation for the dust fraction parametrization (Eq. (70)) that we will be using throughout the paper.

We propose the following suitable (momentum- and energyconserving) SPH discretization for the previous equations $${\left(\frac{d{\mathit{v}}_a}{dt}\right)}_{diff}=-{\sum }_b{m}_b\frac{2{\overline{D}}_{ab}{\overline{s}}_{ab}K}{3{\rho }_a{\rho }_b}(s_b-s_a)\overline{\nabla {\mathcal{W}}_{ab}},$$(78) $${\left(\frac{d\mathit{\Delta }{\mathit{v}}_a}{dt}\right)}_{diff}=-{\sum }_b{m}_b\frac{2{\overline{D}}_{ab}{\overline{s}}_{ab}K}{3{\rho }_a{\rho }_b{\overline{ϵ}}_{ab}}(s_b-s_a)\overline{\nabla {\mathcal{W}}_{ab}},$$(79) $${\left(\frac{d{u}_a}{dt}\right)}_{diff}={\sum }_b{m}_b\frac{2{\overline{D}}_{ab}{\overline{s}}_{ab}K}{3{\rho }_a{\rho }_b(1-{\overline{ϵ}}_{ab})}(s_b-s_a)({\mathit{v}}_a^d-{\mathit{v}}_b^d)·{\widehat{r}}_{ab}{\overline{F}}_{ab},$$(80) where the additional heating term was added to conserve energy, similarly to how Eq. (55) is introduced.

We only need to specify what diffusion coefficient D to use. While some implementations of dust diffusion leave this choice to the user (as in Huang & Bai (2022)), we want our implementation to be adaptive and only trigger when needed in order to make the dust density field smoother in case of large gradients. We tried different approaches, but what was found to have the best results in terms of stability was the following. We define a new variable, χ_a, as the ratio between the dust fraction of particle a and the mean of the dust fractions of all neighboring particles $${\chi }_a=\frac{{ϵ}_a}{\frac{1}{N_{ngb}}{\sum }_b{ϵ}_b}.$$(81) We want diffusion to act on particles for which χ_a departs from unity. A suitable definition for D was found to be $$D_a=\frac{{({\chi }_a-1)}^2}{{\chi }_a}.$$(82) In the SPH discretizations we use the symmetrical $\overline{D_{ab}}$ to ensure momentum and energy conservation.

This implementation could also be extended to describe a physical (or non-artificial) diffusion. This would require modeling the diffusion coefficient D based on gas or dust properties, rather than on the numerical prescription Eq. (82).

In principle, defining a spatially varying diffusion coefficient could allow shocks to develop, which would then require a dedicated artificial-viscosity term. To test whether this actually occurs, we examined the ratio between the compression and advection timescales as a shock indicator. These timescales are defined as $${\tau }_{\mathrm{compr}}=\frac{1}{|\nabla ·{\mathit{v}}_{\mathrm{d}}|},{\tau }_{\mathrm{adv}}=\frac{h}{\left|{\mathit{v}}_{\mathrm{d}}\right|}.$$(83) In a shock, compression is much faster than advection alone would allow, and thus τ_compr<τ_adv. We thus define the shock indicator for dust $$C=\frac{h·|\nabla ·{\mathit{v}}_{\mathrm{d}}|}{\left|{\mathit{v}}_{\mathrm{d}}\right|},$$(84) so that a shock happens in the dust when C > 1.

In Figure 1 we show the histograms of the values for C for the Cold Keplerian Disk and the Kelvin–Helmholtz instability simulations discussed, respectively, in Section 3.5 and Section 3.7. For all the drag coefficients K, and dust-to-gas ratios δ, in simulations the vast majority of particles show values of C < 1, indicating the absence of dust shocks.

Fig. 1 Distribution of the shock indicator defined in Eq. (84) for the Cold Keplerian Disk (top panel) and the Kelvin–Helmholtz instability (bottom panel) simulations.

3 Tests

In this section, we aim at showing both the accuracy and the stability of the implementation of the One-Fluid model in OpenGadget3 shown in the previous section. We first show that the code is able to reproduce the expected results from standard dust dynamics tests: the DUSTYBOX, DUSTYWAVE, and DUSTYSHOCK tests, as well as a dusty Sedov-Taylor blast wave. For all these tests, we will not use our new dust diffusion model. Afterward, we test the stability of the code under more complex flows for various dust-to-gas ratios and drag coefficients. One of the major challenges in simulating dust dynamics is the tendency in low drag and high dust-to-gas ratio environments for dust to clump and create large gradients in dust fraction. To study these behaviors and find numerical solutions to make the code stable under these extreme circumstances, we simulate a Cold Keplerian Disk, a dusty protoplanetary disk, and a Kelvin–Helmholtz instability. These more complex flows will be simulated with dust diffusion turned on. All the values reported in the following section are in internal code units.

Fig. 2 Result of the DUSTYBOX benchmark test. The continuous lines represent the analytical solution, while the dots are the simulation outputs. The results are shown for four different values of the drag coefficient K.

3.1 DUSTYBOX

The first benchmark test we perform is the so-called DUSTYBOX test. It consists of a mixture of gas and dust moving with respect to each other, but globally at rest. This is primarily used to test the integration of the drag term, since all other terms in the evolution equations are zero. The setup is extremely simple: 100 particles are equally spaced in one dimension, the density is set to ρ₀=1, and the dust fraction is set to ε₀=1/2, meaning that there is an equal amount of dust and gas. The barycenter velocity is zero, v=0, so the mixture is at rest and the particles do not move in space. Additionally, the pressure is constant and set to P₀=1.

The only dynamics arises from the dust-gas velocity difference, which is initialized to Δv₀=1. The expected behavior is for this initial velocity difference to decay exponentially, following the simple analytical solution $$\mathit{\Delta }\mathit{v}\left(t\right)=\mathit{\Delta }{\mathit{v}}_0{e}^{-\frac{t}{t_s}}.$$(85) The stronger the drag coefficient K, the smaller the stopping timescale t_s and the faster the exponential decay.

In Figure 2, we show the results of our simulations for four different values of the drag coefficient. The simulation outputs show very good agreement with the expected behaviors. More quantitatively, we found the L1 norm to be always smaller than <0.1%. In addition, energy is also very well conserved: as the simulation evolves, the kinetic energy stored in the second term of the expression for the total energy in Eq. (31) is converted into internal energy of the gas, via drag heating. The resulting energy is conserved up to differences of 10⁻⁷.

3.2 DUSTYWAVE

The next benchmark test is a simple linear wave of a mixture of dust and gas. Since dust is a pressureless fluid, the motion is primarily driven by gas pressure, which makes the gas wave move across the domain. The dust is dragged by the gas according to the strength of the drag coefficient.

The wave is initialized by perturbing an equilibrium particle distribution. The particles span a periodic domain x ∈ (0, 1). The initial conditions for all variables are $$\begin{array}{cc}& \rho (x,0)={\rho }_0+Asin\left(2\pi x\right),\hfill \\ & \mathit{v}(x,0)=Asin\left(2\pi x\right),\hfill \\ & ϵ(x,0)=1/2,\hfill \\ & \mathit{\Delta }\mathit{v}(x,0)=0,\hfill \end{array}$$ where A=10⁻⁴ is the amplitude of the perturbation and ρ₀=1 is the background density value. The pressure is initialized with an isothermal equation of state, $P=\rho c_s^2$, with sound speed c_s=1. What makes the DUSTYWAVE a useful benchmark test is the availability of an analytical solution, derived by Laibe & Price (2011).

The results of our simulations are shown in Figure 3, for three values of drag coefficient. Depending on the strength of the drag coefficient, the dynamics of dust and gas can differ widely. We find very good agreement between simulations and analytical solutions, with L₁ norms of the order of 10⁻⁷. Also in this case energy is very well conserved. This is also achieved by artificial viscosity being active only in its conservative formulation, since the shock indicator here is very close to 0.

Fig. 3 Results from the DUSTYWAVE benchmark test, for three different values of drag coefficient K. The continuous lines represent the analytical solution and the dots the output of the simulations. Gas and dust velocities are plotted separately.

3.3 DUSTYSHOCK

The last suite of benchmark tests consists of a set of shocks. This is useful to test fully the implementation of the One-Fluid model, as shock environments require all of our artificial viscosity and conductivity terms. We explore two opposite scenarios: first, a weak drag regime, in which dust is expected to remain still as the shock propagates in the gas; secondly, a strong drag regime, where dust and gas move tightly coupled, and the shock propagates as a gas-only shock with modified sound speed. For both these regimes, we can compute a numerical solution using a Riemann Solver (here we chose an Exact Riemann Solver).

The DUSTYSHOCK test was already used in LP14b to test the One-Fluid model implementation in one dimension, but here we extend the test also to two dimensions.

Fig. 4 DUSTYSHOCK benchmark test for zero drag. The blue and red dots represent the output of the simulation for gas and dust, respectively, while the continuous line shows the solution for the gas evolution. As expected, while the gas undergoes the typical Sod shock tube evolution, the dust remains unperturbed.

Fig. 5 DUSTYSHOCK benchmark test for strong drag (K=106). The blue and red dots represent the output of the simulation for gas and dust, respectively, while the continuous line shows the solution for the gas evolution. The evolution of the dust now closely follows the gas, and both evolve at a lower speed compared to a gas-only simulation, with the modified sound speed shown in Eq. (86).

3.3.1 Weak coupling

The gas is initialized as in a standard Sod shock tube test (Sod 1978). The left and right states are $$\begin{array}{ll}{\rho }_L^g=1& {\rho }_R^g=0.125,\\ {\mathit{v}}_L^g=0& {\mathit{v}}_R^g=0,\\ P_L=1& P_R=0.1.\end{array}$$ The dust, instead, is initialized as a constant background density at rest $$\begin{array}{rl}& {\rho }_L^d=0.125{\rho }_R^d=0.125,\\ & {\mathit{v}}_L^d=0{\mathit{v}}_R^d=0.\end{array}$$ The drag coefficient is set to K=0. The purpose of this test is to verify that the One-Fluid model implementation is able to evolve the gas shock without perturbing the dust background uniform density, which should remain unperturbed given the zero drag it feels from the gas.

The simulation is initialized with a total of N=1375 particles in a domain x ∈ (−0.5, 0.5). Boundary particles are placed outside each end of the domain to keep the particles confined. The internal energy is initialized and evolved following an ideal gas law, with constant γ=5/3.

The results are shown in Figure 4. As expected, the gas undergoes the standard shock as if the dust were not present, and the dust, conversely, stays unperturbed. This excellent performance in keeping the dust still is to be attributed to the Non-Conservative artificial viscosity scheme, which is here fully active in the shocked regions of the fluid.

3.3.2 Strong coupling

We also explore the opposite scenario, a Sod shock tube test, where dust is strongly coupled to the gas. While the gas initial left and right states are identical to the ones used for the weak coupling case, the dust is now initialized using the same density distribution as the gas, as opposed to the constant background used previously $$\begin{array}{ll}{\rho }_L^d=1& {\rho }_R^d=0.125,\\ {\mathit{v}}_L^d=0& {\mathit{v}}_R^d=0.\end{array}$$ The drag coefficient is set to K=10⁶ and a total of N=1700 particles are employed.

The expected behavior for this shock is to evolve as a standard Sod shock tube, but with a modified sound speed of $$\widetilde{c_s}=\frac{c_s}{\sqrt{1-ϵ}}.$$(86) As can be shown in Fig. 5, the result from the simulation follows closely the expected solution, and dust and gas remain tightly coupled.

Fig. 6 Same as in Figure 4, but now in two dimensions. The solution was obtained with a one-dimensional Exact Riemann Solver with an additional geometrical source term as described in Eq. (87). The dust now shows some additional wiggles compared to the one-dimensional case, which can be attributed to the lower effective resolution of this higher-dimensional test.

Fig. 7 Same as in Figure 5, but now in two dimensions. The solution was obtained with a one-dimensional Exact Riemann Solver with an additional geometrical source term as described in Eq. (87).

3.3.3 Two-dimensional shocks

In order to test the One-Fluid model implementation in two dimensions, we employed a strategy discussed in Toro (2009) to turn one-dimensional shocks into two-dimensional radially symmetric shocks. In particular, Toro (2009) shows that the radial evolution of a shock in multiple dimensions can be replicated in one-dimensional simulations by adding a geometrical source term to the Euler equations $${\partial }_t\mathit{U}+{\partial }_r\mathit{F}\left(\mathit{U}\right)=\mathit{S}\left(\mathit{U}\right),$$(87) where U is the vector of conserved variables, F(U) are the corresponding fluxes and the geometrical source term S(U) is defined as $$\mathit{S}\left(\mathit{U}\right)=-\frac{\alpha }{r}\left[\begin{array}{c}\rho v\\ \rho v^2\\ v(E+p)\end{array}\right].$$(88) By varying the α parameter we obtain a cylindrical (α=1) or a spherical (α=2) symmetry. With this procedure, we can obtain the expected result for two-dimensional simulations using a onedimensional Riemann solver.

The initial conditions for the shocks in two dimensions are essentially the same as in the one-dimensional shock. The particles are arranged in a hexagonal closed pack lattice in a periodic domain x ∈ (−1, 1), y ∈ (−1, 1). Particles inside a sphere centered around 0 and of radius R=0.5 have the “left” state initial conditions, while the particles outside have the “right” state initial conditions. The results for the two-dimensional shocks for both the weak and strong drag cases are shown in Fig. 6 and Fig. 7. The performances are slightly worse compared to the one-dimensional case, which can be explained by a lower effective resolution of the multidimensional cases, but the shocks all evolve as expected.

3.4 Dusty Sedov-Taylor blast wave

The Sedov–Taylor blast wave is a strong, radially symmetric shock wave, first introduced by Sedov (1946, 1959). It is widely used in astrophysics to model the evolution of supernova blast waves (Steinwandel et al. 2020) and to benchmark the shock-capturing capabilities and limitations of hydrodynamical simulation codes (Springel & Hernquist 2002; Rosswog & Price 2007; Groth et al. 2023).

For gas-only simulations, a self-similar solution exists and can be used to verify the performance of numerical schemes. Although no such solution is known for dust-gas mixtures, in the limit of strong drag and high dust-to-gas ratios, the evolution is expected to follow the gas-only solution but at a modified sound speed, as in the DUSTYSHOCK tests.

In the context of dust-gas mixtures, Laibe & Price (2012a) simulated a Sedov-Taylor blast wave using a two-fluid model for a moderate drag coefficient of K=1, since a high drag simulation would have been too computationally expensive, and a zero-drag simulation would have been trivial in the two-fluid model. Here we will simulate the Sedov-Taylor blast wave for a mixture of dust and gas in two cases, mirroring the ones used in the DUSTSHOCK tests: (i) a zero-drag simulation with dust-to-gas ratio of 0.01 and (ii) a high-drag (K=10⁶) simulation with dust-to-gas ratio of 1.

We initialize N=64³ particles in a regular grid x, y, z ∈ (−0.5, 0.5) with periodic boundary conditions, uniform gas density ρ_g=1 and vanishingly low pressure P=10⁻⁶, except for the particles closest to the center which are initialized with a large internal energy U=10. Consequently, a shock with Mach number M ∼ 2 · 10⁴ propagates through the domain. The simulation is evolved until t_f=0.02

In Figure 8, we show the radial density profiles obtained from the simulations for both regimes, zero and high drag, together with the self-similar solution. For K=0 (left), the gas undergoes the Sedov-Taylor blast wave and approaches the expected self-similar solution, while the dust remains largely at rest. For high drag (right), the dust and gas stay tightly coupled, and the mixture can be described by the self-similar solution evaluated at an earlier time, consistent with the modified sound speed given in Eq. (86).

Fig. 8 Dusty Sedov-Taylor blast wave for a dust and gas mixture in two regimes, K=0 (left) and K=106 (right). The radial density profiles are shown as blue dots for gas and red dots for dust, and the self-similar solution at t=tf is plotted as a continuous black line. The simulation correctly keeps dust at rest while the gas evolves following the self-similar solution, for uncoupled dust. For strongly coupled gas and dust, the mixture follows the self-similar solution evaluated at a previous time, consistent with the modified sound speed in Eq. (86), shown as a dashed black line.

3.5 Dusty Cold Keplerian Disk

The Cold Keplerian Disk is a standard hydrodynamical test used to study the stability of numerical frameworks (Cartwright et al. 2009; Hopkins 2015; Groth et al. 2023). SPH has already been shown to struggle in keeping a Cold Keplerian Disk stable for various orbits; thus, we are not interested in studying its stability. Instead, we will focus on how dust clumping behaves under various dust-to-gas ratios and drag coefficients and how the dust pressure model introduced in Section 2.7 can help keep the code stable.

The initial condition for the Cold Keplerian Disk is taken, for instance, from Groth et al. (2023). The gas surface density is given by $${\Sigma }_g=0.01+\{\begin{array}{cc}{(R/0.5)}^3\hfill & \text{if}R\le 0.5,\hfill \\ 1\hfill & \text{if}0.5\le R\le 2,\hfill \\ {(1+(R-2)/0.1)}^{-3}\hfill & \text{if}R>2.\hfill \end{array}$$(89) Particles are initially distributed in a hexagonal close-packed grid in a domain x ∈ (−2, 2), y ∈ (−2, 2), and then removed from R<0.5 and R>2 until the desired density profile is reached. The effective resolution of the initial grid is 256². The velocity is set to strictly Keplerian, and the pressure is set at a vanishingly low value of P=10⁻⁶. In adding the dust we simply increase the density uniformly by a factor (1+δ), where δ is the dust-to-gas ratio. The dust fraction is set accordingly as a uniform ε=δ /(1+δ). Dust is also initialized at a Keplerian velocity, hence, the initial difference velocity vector is Δv=0.

Such initial conditions should remain at equilibrium indefinitely, but in SPH the artificial viscosity eventually causes the disk to break up after roughly 5−10 orbits (Hopkins 2015). We conduct a series of simulations varying the initial dust-to-gas ratio and the drag coefficient K.

In Figure 10, we show two simulations, both with δ=0.01 and K=0.01, the only difference being the use of dust pressure as a diffusive term, as described in Section 2.8. The left half of the figure, corresponding to a simulation without any diffusion, shows severe clumping in dust density, which eventually causes a crash in the code. The snapshot shown here is the last available from the simulation. On the right half is shown a simulation where diffusion is switched on, at the same snapshot. As shown in the figure, the dust density field is less clumped, resulting in smoother gradients and a more stable simulation.

Figure 9 presents the results of three distinct simulation configurations at time t=120, corresponding to roughly 20 orbits at R=1. As anticipated, the disk becomes unstable over time in all cases due to the artificial viscosity inherent in SPH. The left panels depict a simulation with weak drag (K=0.01) and a low dust-to-gas ratio (δ=0.01). In this scenario, the dust gradually drifts inward as it loses angular momentum, but its influence on the gas dynamics is negligible due to the small dust-to-gas ratio. The middle panels illustrate a simulation with strong drag (K=100) and the same low dust-to-gas ratio (δ=0.01). Here, the dust closely follows the motion of the gas, as expected under strong coupling conditions. Finally, the right panels show a simulation with strong drag (K=100) and a high dust-to-gas ratio (δ=1). In this case, the disk undergoes significant disruption, as the pressure gradient force of the gas is insufficient to counteract the influences from the dust dynamics. The dust concentrates in spiral structures that actively drag the gas along, resulting in a highly coupled system where the gas motion is strongly influenced by the dust.

While the presence of dust diffusion helped stabilize the code under various conditions, especially for weak drag regimes, the code could not remain stable for the most extreme conditions of weak and intermediate drag (K=0.01 and K=1) for a dust-to-gas ratio of unity (δ=1). Although our code aims at being as general-purpose as possible, such extreme conditions are rarely encountered in common astrophysical applications and would require specialized techniques to be simulated.

Fig. 9 Dusty Cold Keplerian Disk simulations for three different sets of parameters. Gas density is shown in the top panels while dust density is shown in the bottom panels.

Fig. 10 Zoom-in in the central region of the Cold Keplerian Disk simulation with weak drag (K=0.01) and low dust-to-gas ratio (δ=0.01). The left half shows the last snapshot of a run without diffusion. The lack of a restoring force causes strong clumping of dust density, which leads to a code crash very early into the simulation. The addition of diffusion, as shown in the right half of the figure, dissipates these clumps and makes the code more stable.

3.6 Dusty protoplanetary disk

The Cold Keplerian Disk is a useful and challenging hydrodynamical test to investigate the stability of our implementation. However, a more standard and physically motivated setup is a model of an actual protoplanetary disk. Our configuration follows previous SPH studies of protoplanetary disks (Laibe & Price 2012b; Vericel et al. 2021), but restricted here to two dimensions.

The surface density distribution follows a power law with a smooth inner edge $$\Sigma \left(R\right)={\Sigma }_0{R}^{-p}(1-\sqrt{\frac{R_{\mathrm{in}}}{R}}),$$(90) with p=3/2 and R_in=10 au, and Σ₀ chosen so that the total mass of the disk is M_disk=0.01 M_⊙. To generate a particle distribution with this density profile, we first use inverse sampling of the cumulative distribution function to generate values for the radii, and random values for the ϕ coordinate. To reduce the statistical noise of the resulting point distribution, we apply the Centroidal Voronoi Tessellation (CVT) technique described in Du et al. (1999). The initial particle set is used to construct the corresponding Voronoi tessellation with the qvoronoi command from the Qhull library (Barber et al. 1996). We then compute the mass centroid of each Voronoi cell as $${\mathit{x}}_{j,i+1}=\frac{{\int }_{V_{j,i}}\mathit{x}\Sigma \left(\mathit{x}\right)dV}{{\int }_{V_{j,i}}\Sigma \left(\mathit{x}\right)dV},$$(91) where Σ is the target density profile in Eq. (90), the j index spans over the points in the distribution and the i denotes the iteration step. The new point distribution x_{j, i+1} is then used to draw a new Voronoi tessellation, and the procedure is iterated until convergence. The resulting point distribution is much smoother than the initial one obtained from inverse sampling.

The CVT was used to generate the positions of N=10⁶ points between R_in=10 au and R_out=400 au. The dust-to-gas ratio is set to 0.01, and both gas and dust follow the same initial density profile. The disk is locally isothermal, with the sound speed following the power-law $$c_s=c_{s,0}{R}^{-q/2},$$(92) with the index q=3/4 and the sound speed at radius R=1 au equal to the aspect ratio: c_{s, 0}=h=0.0281, as in Vericel et al. (2021).

The azimuthal velocity is set to exactly Keplerian for dust, and slightly sub-Keplerian for the gas $$v_{g,\varphi }=v_K\sqrt{1-h^2(p+q)},$$(93) where v_K is the Keplerian azimuthal velocity.

The dust drag coefficient is set through its Stokes number, and two values are used: St=0.01 and St=1, to study the different evolution of tightly and weakly coupled dust grains.

These initial conditions are evolved until t_max=20 kyr under the influence of the gravitational potential of the central star of mass M_*=1 M_⊙. During the evolution, particles that cross the inner radius R_in are removed from the simulation.

The gas and dust density maps for both St=0.01 and St=1 at t=t_max are shown in Figure 11. For tightly coupled grains, both the gas and dust density profiles remain stable for long timescales. The pressure support of the gas and the drag force acting on the dust make this test much more stable than the Cold Keplerian Disk counterpart. For St=1, the low dust content is not enough to modify the gas dynamics. However, dust undergoes more significant clumping and radial drift. In Figure 12 we compare the radial dust density profiles for the two runs, and it is evident that the St=1 run favors radial drift toward the pressure maximum created by the inner radius tapering in Eq. (90).

Fig. 11 Gas (left panels) and dust (right panels) density for a dusty protoplanetary disk for two values of Stokes number at t=20 kyr. For tightly coupled grains (St=0.01, top row), the gas and dust density remains at equilibrium, while for weakly coupled grains (St=1, bottom row), the dust undergoes both clumping and radial drift toward the pressure maximum.

Fig. 12 Dust radial density profile for both runs with St=0.01 and St=1. The coupled grains keep their original density profile, driven by the strong drag force, which forces them to follow the dynamics of the gas. The weakly coupled grains, on the other hand, undergo a much stronger radial drift toward the position of the pressure maximum.

Fig. 13 Dusty Kelvin–Helmholtz instability simulations for three different sets of parameters. Gas density is shown in the top panels while dust density is shown in the bottom panels.

Fig. 14 Time evolution of the Kelvin–Helmholtz instability mode amplitude for three dust-to-gas ratios, δ=(0.01, 0.1, 1), and three drag coefficients, K=(0.01, 1, 100). Time is shown in units of the growth timescale τK H.

3.7 Dusty Kelvin–Helmholtz instability

To explore the capability of the One Fluid model in simulating the dynamics of dust in turbulent motion, we simulate a Kelvin–Helmholtz instability for various choices of dust-to-gas ratios and drag coefficients. As before, we add the dust to a gas-particle setup by uniformly multiplying the density by (1+δ). We use the same initial conditions as, for instance, McNally et al. (2012), Hopkins (2015) or Groth et al. (2023): a periodic box in two dimensions is filled with two fluids of density ρ₁=1 and ρ₂=2 with opposite velocities v₁=0.5 x̂ and v₂=−0.5 x̂ and a small perturbation in the ŷ direction to trigger the instability. We expect the instability to grow with a timescale ${\tau }_{KH}=\frac{\lambda }{v_{x,1}-v_{x,2}}\frac{{\rho }_1+{\rho }_2}{\sqrt{{\rho }_1{\rho }_2}}$. In Figure 13 we show both the gas and dust densities at a time t ≃ 1.2 τ_{K H} for the same dust-to-gas ratios and drag coefficients that were shown in Figure 9 for the Cold Keplerian Disk.