© The Authors 2025

1 Introduction

Despite various efforts to decipher the nature of dark matter (DM) with laboratory experiments, all evidence for its existence still stems from astrophysical and cosmological observations, such as the rotation curves of galaxies (Rubin et al. 1980; Bosma 1981). Even though the collisionless cold dark matter (CDM) model is quite successful in explaining several observations, including the cosmological large-scale structure, it does not provide much insight into the particle nature of DM. Whether DM has any interaction other than the gravitational force remains an open question. Potentially, DM particles could interact with each other through novel physics of a dark sector or end up coupled to the particles of the standard model (SM).

Both scenarios have the potential to change the evolution of astrophysical objects such as galaxies and galaxy clusters, compared to the case of collisionless DM in cases where the interactions would be strong enough. This opens up a window onto probing the particle nature of DM via astronomy and eventually discovering new physics of the dark sector. It also paves an avenue for probing the cross-section for DM self-interactions and DM-baryon interactions.

In the first case of self-interacting dark matter (SIDM), various astrophysical systems and their observables have been used to constrain the strength of self-interactions and explain potential discrepancies between CDM predictions and observations (described in the review articles by Tulin & Yu 2018; Adhikari et al. 2022). The strongest upper bounds on the crosssection come from galaxy clusters, while low-mass systems such as dwarf galaxies offer hints that DM may have strong self-interactions at low velocities. These efforts have been supported by semi-analytical and numerical modelling of SIDM. In particular, for the latter, it is possible to run full physics cosmological simulations to obtain predictions for the formation of objects covering a large mass range from dwarf galaxies (e.g. Vogelsberger et al. 2014; Correa et al. 2025) through galaxy clusters (e.g. Robertson et al. 2019; Ragagnin et al. 2024; Despali et al. 2025).

In contrast, for the case of interacting dark matter (IDM), there are no cosmological simulations available to model the DM-baryon interactions in situ via the simulation codes. The astrophysical constraints are based on analytical calculations and simulations of the linear evolution of the matter power spectrum (e.g. Sigurdson et al. 2004; Dvorkin et al. 2014; Boddy & Gluscevic 2018; Ali-Haïmoud et al. 2024). They have, for example, been derived from cosmic microwave background (CMB) observations (Boddy et al. 2018; Gluscevic & Boddy 2018), the Lyman-alpha forest and large scale structure (Dvorkin et al. 2014; He et al. 2023, 2025), and the Milky Way satellite galaxy abundance (Maamari et al. 2021; Crumrine et al. 2025). A recent study (He et al. 2025) constrained velocity-dependent DM-baryon interactions employing CMB data together with observations of the large scale structure. Constraints on DM-proton and DM-electron scattering were also derived by Nguyen et al. (2021) and Buen-Abad et al. (2022) based on CMB anisotropies, baryon acoustic oscillations, the Lyman-α forest, and the abundance of Milky Way subhalos. The implications at low redshifts due to an altered matter power spectrum by DM-baryon scattering at high redshifts were recently studied using N-body simulations, assuming collisionless DM (e.g. Zhang et al. 2024; Nadler et al. 2025a,b; An et al. 2025). Moreover, galaxy clusters provide a promising probe for significant DM-baryon interactions at late times. By studying the heat exchange between the DM and the intra-cluster medium, constraints on IDM can be inferred (e.g. Shoji et al. 2024; Stuart & Pardo 2024). Furthermore, the ionisation of molecular clouds has been used to constrain DM-proton scattering (Prabhu & Blanco 2023; Blanco et al. 2024) and the orbital decay of pulsars has been employed in constraining IDM as well (Lucero et al. 2024).

More generally speaking, signatures of DM-baryon scattering at low redshift may have similarities to DM self-interactions, since IDM is also expected to alter the matter distribution on small scales. While SIDM is only capable of altering the properties of the baryons indirectly via changes in the DM gravitational potential, the effects of IDM go beyond that. In particular, DM-baryon interactions can lead to energy exchange between DM and baryons, which means that they could cool or heat the interstellar or intra-cluster medium, depending on the particle physics model. Moreover, these interactions can make relative motions between DM and baryons decay or affect the ionisation fraction of the gas. In turn, these effects can potentially impact other processes, such as star formation, and alter the evolution of galaxies. Overall, the phenomenology of IDM is rich and can vary substantially between models with interaction cross-sections that differ in their angular and velocity dependence as well as the mass ratio of the interacting particle species.

The extent to which astrophysical probes can be used to constrain DM-baryon interactions is limited by our ability to model these interactions. The use of N-body simulations is a common technique for studying various astrophysical systems at different scales, for example, the formation and evolution of galaxies within the cosmological context. The evolution of galaxies is shaped by various physical processes beyond gravity, such as gas dynamics, the evolution of the stellar component, and black holes (BH). The non-linear interplay of different physical processes often hinders a precise analytic description and requires numerical simulations. Unfortunately, however, it has not yet been possible to successfully include the effects of DM-baryon interactions in such simulations and thus exploit the full potential of astrophysical probes to learn about those interactions.

On the other hand, significant efforts have been undertaken to constrain DM interactions with direct detection experiments (see for example Cirelli et al. 2024, and the references therein). In particular, ground-based experiments have mainly been carried out, but studies are not limited to these (e.g. Emken et al. 2019; Du et al. 2024). In addition, DM interactions are also tested by searching for annihilation and decay products in cosmic and gamma rays. This implies a wealth of constraints that could be complemented by studies of the dynamical impact of DM physics.

While numerous constraints on DM-baryon interactions exist, leveraging astrophysical probes with the help of new techniques to model these interactions can provide novel constraints that may allow us to probe uncharted parts of the DM parameter space. Towards this end, we introduce a novel numerical scheme that allows us to simulate the interaction between SM and DM particles in N-body simulations in situ. We study its numerical behaviour, show its abilities, and apply it to the collapse of an overdensity and halo formation. Our aim is to develop a scheme that is suitable for application in cosmological simulations of galaxy formation involving various physical processes. However, the exploration of this possibility lies beyond the scope of this paper. Instead, we focus on the numerical foundation to describe the DM interactions.

This work is structured as follows. In Sect. 2 we explain the scheme for the DM-baryon interactions. It follows a set of test problems to study its numerical properties (Sect. 3). We study the effect of IDM on the halo formation by simulating the collapse of an overdensity (Sect. 4). In Sect. 5, we discuss the limitations of this work as well as directions for further research. Finally, we present our summary and conclusions in Sect. 6. Additional information is provided in the appendices.

Fig. 1 Illustration of the numerical scheme for the DM-baryon interactions. The numerical particles are shown together with the physical particles they represent. The velocities of the physical particles are indicated by small arrows. The different stages in treating the interaction between baryons and DM for a single pair of numerical particles are illustrated from the left to the right. As shown here, the case of an interaction between a numerical baryonic particle and a numerical DM particle involves a change in their bulk motion and heating of the baryons. Finally, we note that here we also illustrate the case where the baryons consist of two species, but only one interacts with DM.

2 Numerical methods

In this section, we introduce our novel formulation of the interactions between DM and baryons. A sketch of the idea can be found in Fig. 1. We explain the numerical scheme and describe the implementation.

2.1 Formulation of DM-baryon interactions

We aim to describe the DM-baryon interactions within N-body simulations where the baryons are modelled with smoothed particle hydrodynamics (SPH), meshless finite mass (MFM), or other schemes. Those codes are commonly used for cosmological simulations and studies of galaxy formation and evolution.

Overall, DM is represented by numerical particles characterised by the mass, m_DM, position, x_DM, and velocity, v_DM. This can be understood as a particle representing a phase-space patch with the physical DM particles having the same velocity. The numerical particles representing the baryons analogously have a mass, m_bary, and position, x_bary. In addition, the bulk velocity of the physical particles is described by v_bary and their random motions, assumed to follow a Maxwell-Boltzmann distribution, are characterised by the internal energy per mass, u_bary.

Similarly to the numerical schemes of SIDM, we formulated the interaction based on pairs of numerical particles. A pair consists always of one DM and one baryonic particle. The interaction only takes place if they are close enough to each other. The particles are assigned a kernel with a size, h, determined by the N_ngb next neighbours of the same particle specified (i.e. DM or baryons)¹. They only interact if their kernels are overlapping. During a time step, for a particle, all the pairs it forms with its neighbours from the other particle species (DM or baryons) are used to model the DM-baryon interactions. In the following, we describe the details of this interaction.

The physical particles represented by a numerical DM-baryon pair would scatter in different centre-of-mass systems and their post-scattering velocities would point in various directions. It is impossible to accommodate these distributions with the two numerical particles of the pair. Instead, we use a stochastic description and develop a Monte Carlo scheme, as commonly done in SIDM as well. However, in contrast to the DM particle, the baryonic one represents a whole distribution of velocities and not a single one. This is because for the baryons we can assume that the velocities locally follow a Maxwell-Boltzmann distribution, whereas the DM particles do not follow a specific velocity distribution, making it necessary to resolve the velocity space². To account for this, we take a random velocity from the distribution of baryonic velocities. This velocity is used to form an interaction partner for the numerical DM particle.

Specifically, we create a particle called a ‘virtual’ particle when we want to compute the interaction between a DM-baryon pair and we destroy the virtual particle when we completed the interaction of this pair (as illustrated in Fig. 1). This implies that we create many virtual particles from each considered baryonic particle per time step, but they never exist simultaneously. This is because all pairs that a baryonic particle forms are computed in a well-defined consecutive manner and a virtual particle exists only for the time we are considering a specific pair of a numerical baryonic and DM particle. This consecutive order is necessary to ensure energy conservation. If, in contrast, we were to execute the computations of two pairs that share a common particle at the same time, we would use the same initial properties for the computations. Thus, we might end up with two incompatible sets for the post-scattered properties. Despite being consecutive for every particle, the specific order in which the pairwise computations are executed does not matter and, thus, it should not impact the accuracy of our simulations. The only purpose of the virtual particles is to allow us to formulate the interaction between the DM and baryons in an manner that explicitly conserves the mass, energy, and momentum³.

The virtual particle sits at the same position as the baryonic one, x_virt = x_bary, and it has the same kernel size, h_virt = h_bary. Its velocity is given by the bulk motion of the baryons plus a random component, expressed as $${\mathbf{v}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}={\mathbf{v}}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}+{\mathbf{v}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}(u_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}).$$(1)

We note that the random component depends on the internal energy per mass of the baryonic particle, u_bary. It is drawn from a Maxwell-Boltzmann distribution, $$f(v)=\sqrt{\frac{2}{\pi }}\frac{v^2}{a^3}{e}^{-\frac{v^2}{2a^2}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}a=\sqrt{\frac{2}{3}{u}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}.$$(2)

To avoid very high velocities for the virtual particle, we cut the high-velocity end of the Maxwell-Boltzmann distribution. Velocities higher than v_cut = ζ a are reduced to v_cut. In practice, we used ζ = 5, which should be large enough such that its effect is negligible compared to all scatterings. More precisely, ζ = 5, implies a relative error for the energy represented by the velocity distribution of ≈10⁻⁵. However, this v² weight may underestimate the impact on the simulation results. A better estimate might be to weight the velocities by v⁵, as motivated by the effective or characteristic cross-sections for SIDM (Yang & Yu 2022; Yang et al. 2023). With such a weight, the relative error becomes ≈3 × 10⁻⁴. This implicitly assumes a velocity-independent crosssection. However, if the cross-section decreases with velocity, the error would be smaller.

It should be noted that with Eq. (2), we assume that the baryons consist only of one type of particle, namely, the one that interacts with the DM. However, we could also accommodate for more complicated situations, which may require a different value for the distribution parameter, a.

To model the scattering kinematics correctly, the numerical particles must have the same mass ratio, r, as the physical particles⁴. For the mass of the virtual particle, we obtain $$m_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}=r{m}_{\mathrm{D}\mathrm{M}}.$$(3)

Based on the DM and virtual particle, we can compute the scattering as done in SIDM codes (e.g. Koda & Shapiro 2011; Fry et al. 2015; Robertson et al. 2017; Yang & Yu 2022). In practice, we follow Fischer et al. (2021a). The scattering can alter the velocities of the two numerical particles. Hence we obtain the post-scattered velocity, v′_DM, for the DM particle and v′_virt for the virtual particle. The post-scattered velocities are obtained by rotating the momentum vectors in the centre-of-mass frame. In the case of large-angle scattering, we first compute a probability of determining whether the two particles interact or not (e.g. Burkert 2000; Rocha et al. 2013). In particular, we follow the scheme for rare interactions of Fischer et al. (2021a). The interaction strength for the particles i and j, depends on a geometrical factor, Λ_ij, based on the kernels, W(x, h), assigned to the numerical particles, as $${\mathrm{\Lambda }}_{ij}=\int W(|\mathbf{x}-{\mathbf{x}}_i|,h_i)W(|\mathbf{x}-{\mathbf{x}}_j|,h_j){\mathrm{d}}^3\mathbf{x}.$$(4)

We compute Λ_ij as described in Appendix A by Fischer et al. (2021a).

The probability that the DM particle and the virtual particle interact is given as $$P_{ij}=\frac{\sigma (v)}{m_{\chi }}\frac{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\mu }{m}_{\mathrm{D}\mathrm{M}}v\mathrm{\Delta }t{\mathrm{\Lambda }}_{ij}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}v=|{\mathbf{v}}_i-{\mathbf{v}}_j|.$$(5)

The total cross-section is given by σ(ν), and the physical DM particle mass by m_χ. We employ the mass ratio μ = m_virt/m_bary. We note that here we use the velocity of the virtual particle, not the baryonic one. This is a consequence of the fact that the baryonic particle represent a distribution of velocities and the velocity of the virtual particle is a random velocity drawn from this distribution. Using the velocity from the virtual particle gives the correct interaction probability and allows us to account for arbitrary velocity dependencies of the interaction cross-section. In addition, we introduced the parameter f_bary to specify the mass fraction of the baryonic particle taking part in the interaction. It is worth mentioning that here we assume a physical particle scatters only once with a particle represented by the other numerical particle of the pair per time step. This still allows for a particle to scatter multiple times, but only with partners from different pairs. The probability that a physical particle scatters twice is $P_{ij}^2$. Hence, we have to choose the time step, Δt, that would be small enough to make $P_{ij}^2$ negligible. This constraint is a consequence of modelling the interactions between two numerical particles analogously to a single physical scattering event; namely we produce the post-scattered velocity distribution by assuming a single scattering event. For multiple scattering events per particle, this distribution would look differently and, thus, the contribution of those must be kept small to accurately model the interactions.

Similarly to the interaction probability, we can formulate a drag force term for small-angle scattering analogously to the scheme for frequent self-interactions by Fischer et al. (2021a). To characterise the strength of the interaction, we use the momentum transfer cross-section, $${\sigma }_{\mathrm{T}}=2\pi {\int }_{-1}^1\frac{\mathrm{d}\sigma }{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}}(1-\cos {\theta }_{\mathrm{c}\mathrm{m}\mathrm{s}})\mathrm{d}\cos {\theta }_{\mathrm{c}\mathrm{m}\mathrm{s}}.$$(6)

The drag force for our DM-baryon interactions is given as $$F_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}=\frac{{\sigma }_{\mathrm{T}}(v)}{m_{\chi }}\frac{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\mu (1+r)}{m}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}{m}_{\mathrm{D}\mathrm{M}}{v}^2{\mathrm{\Lambda }}_{ij}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}v=|{\mathbf{v}}_i-{\mathbf{v}}_j|.$$(7)

A derivation of the interaction probability, P_ij (Eq. (5)), and the drag force, F_drag (Eq. (7)), can be found in Appendix A. The final step of the interaction is to destroy the virtual particle or in other words, thermalise it back into the baryonic particle. The underlying idea is that the interactions between the physical baryonic particles are strong enough to quickly distribute the exchanged momentum and energy from the DM-baryon scattering over the physical baryonic particles (represented by the single numerical particle under consideration). However, this thermalisation timescale must be small enough (relative to the numerical time step, which is roughly speaking set by the minimum of the local dynamical time and the inverse scattering rate), so that our previous assumption of a Maxwell-Boltzmann distribution is valid as well.

Next, we can derive the post-scattered properties of the baryonic particle. To do so, we start with momentum conservation to obtain the new bulk motion; namely, the momentum change that the virtual particle has experienced is then applied to the baryonic particle. $${\mathbf{v}}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^{\prime }={\mathbf{v}}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}+\mu ({\mathbf{v}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}^{\prime }-{\mathbf{v}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}).$$(8)

We also update the internal energy per mass fulfilling energy conservation, namely, the baryonic particle must account for energy change that the virtual particle experiences. $$u_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^{\prime }=u_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}+\frac{1}{2}[{\mathbf{v}}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^2-{\mathbf{v}}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^{\prime 2}+\frac{m_{\mathrm{D}\mathrm{M}}}{m_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}({\mathbf{v}}_{\mathrm{D}\mathrm{M}}^2-{\mathbf{v}}_{\mathrm{D}\mathrm{M}}^{\prime 2})].$$(9)

This enables the internal energy of the baryons to increase or decrease, effectively allowing for heat flow between the two components in each direction.

We expect the numerical scheme we present here to converge to the true physical solution in the simultaneous limit of N → ∞, N_ngb → ∞, N/N_ngb → ∞, m_bary/m_DM → ∞ and Δt → 0⁵. We note that for very anisotropic cross-sections in the frequent limit, we no longer require m_bary/m_DM → ∞. This limit is relevant to reduce the effect a single numerical particle interaction has on the baryonic particle, but in the frequent interaction scheme for small-angle scattering, this is already achieved by Δt → 0 or N_ngb → ∞. In the frequent interaction scheme, an effective scattering angle determined by a drag force description is used (Fischer et al. 2021a). It also depends on the size of the time step and the neighbour number. In this case, the effect on the baryonic particle per numerical particle interaction can simply be reduced by choosing smaller time steps or increasing the neighbour number. This implies, for the limit of very anisotropic cross-sections, that we do not require the numerical mass of the baryonic particle to be much higher than the one of the DM particle.

2.2 Negative internal energy problem

The energy required to generate the numerical virtual particle from the baryonic one can potentially be greater than the internal energy of the baryonic particle. Precisely speaking, in the rest frame of the baryonic particle before creating the virtual particle, the sum of the kinetic energy of the virtual particle and the remaining baryonic one (minus the mass of the virtual particle) could be greater than its initial internal energy. This can cause a problem as in the case when the following inequality is violated, $$\underset{=E_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l},\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\underbrace{u_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{m}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}}>\underset{=E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}}{\underbrace{\frac{1}{2}{m}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}{\mathbf{v}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}^2}}+\underset{=E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}}{\underbrace{\frac{1}{2}\frac{m_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}}{m_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}-m_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}}{\mathbf{v}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}^2}}.$$(10)

If this condition is not fulfilled, the virtual particle could lose enough energy via scattering with a DM particle to the extent that even after the virtual and baryonic particle are rejoined, the internal energy of the baryonic particle is non-positive. This would cause a severe problem for the numerical scheme to model the hydrodynamics.

To ensure that the specific internal energy, u, does not become zero or negative, we find that the velocity of the virtual particle must fulfil $$|{\mathbf{v}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}|<a\sqrt{3(\frac{m_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{r{m}_{\mathrm{D}\mathrm{M}}}-1)}=a\sqrt{3(\frac{m_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{m_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}}-1)}.$$(11)

This can be seen by requiring that the remaining internal energy of the baryonic particle must be positive after creating the virtual particle.

When replacing the left-hand side of Eq. (11) with v_cut = a ζ, because |v_rand| is chosen not to be larger than v_cut, we can derive a bound on the mass ratio of the baryonic and DM particles, $$\frac{m_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{m_{\mathrm{D}\mathrm{M}}}>r(\frac{{\zeta }^2}{3}+1).$$(12)

This constrains the mass ratio of the numerical particles required to prevent negative internal energies. We have to note that Eq. (12) only applies to large-angle scattering as the maximal change in internal energy that we use here only arises from a large scattering angle.

In contrast, small scattering angles typically do not lead to enough energy exchange for creating negative values of the internal energy. But if the numerical particle mass of the baryons is much smaller than the one of the DM particles even smallangle scattering could run into the problem of negative internal energies given that the time step is not chosen small enough⁶. However, this also depends on the baryon temperature and DM velocity dispersion involved. Also, there is no equation as simple as Eq. (12) for such a case.

For some models, Eq. (12) might be to restrictive. Those featuring large-angle scattering and having a sizeable value of r would require numerical particles for the DM being much less massive than the numerical baryonic particles. This makes simulating models with a large r, fairly expensive or even unfeasible. However, at the same time, scatterings in such a model are much less likely to lead to negative internal energy as they have a smaller impact on the baryonic scattering partner. Consequently, it might be more useful to choose a set-up that does not satisfy Eq. (12). To ensure at the same time that internal energies stay positive, we explicitly check the internal energy after each interaction and reject the post-scattered values if the internal energy has a non-positive value. We recompute the scattering event from the pre-scattered values (for large-angle scattering this includes the decision if the particles interact) and repeat this until the internal energy is positive. In practice, for large-angle scattering, this turns most of the problematic scattering events into nonscatters. If too many scatterings are rejected this would lead to an underestimate of the effect of DM-baryon interactions. We note that the scatterings with a high relative velocity are predominately rejected, namely, the ones that lead to the largest energy and momentum exchange between DM and baryons.

We note that when the simulation parameters are chosen carefully, this rejection scheme has a negligible impact on the simulation results, as only a small number of scatterings end up being rejected. With this in mind, we set up the simulations for this work to ensure that rejections either do not occur at all, which applies for most of our simulations, or occur very rarely. In the latter case, the number of rejections is many orders of magnitude smaller than the total number of interactions. In principle, the number of rejections depends on the choice of ζ ; for a larger value, more rejections could occur.

2.3 The role of viscosity and heat conduction

To model DM-baryon interactions accurately with the scheme described above, viscosity plays a crucial role. This can be illustrated by considering a situation where the DM and baryon distributions are at rest, but have a different velocity dispersion (temperature)⁷. Physically speaking, the interactions would cool the DM (reduce their velocity dispersion) and heat the baryons or vice versa, while no net momentum would be exchanged and both components would remain at rest. However, in our numerical scheme, we also modify the bulk motion, i.e. the velocity of the numerical baryonic particles. This bulk motion is random and not coherent among the numerical particles. It might be viewed as an artificial small-scale turbulence. The strength of this turbulence could potentially grow over the course of the simulation. How strong it is depends on the numerical set-up, for example, the value of μ.

To ensure accurate results, the artificial turbulent motion of the baryons (random bulk motion) must be kept small. If it reaches a significant strength though, the temperature of the baryons would appear to be cooler compared to the exact physical solution. If the baryons are subject to viscous forces, turbulent motion, including artificial turbulence, is dampened. Viscosity effectively transfers the kinetic energy induced by the DM-baryon interactions into the internal energy of the baryons and thus increases the temperature of the baryons. Overall, viscosity must be strong enough to ensure that the heat exchanged between baryons and DM does not significantly end up as kinetic energy of the baryons, but is transferred to internal energy. Otherwise, the temperature of the baryons is expected to be inaccurate.

We note that even if the internal energy of the baryons is significantly off this must not imply that the heat exchange between the DM and baryonic component is off as well. This is because the artificial small scale turbulence, i.e. the velocities of the numerical particles representing the baryons is present in calculations for the DM-baryon interactions. However, it is no longer possible to clearly distinguish between turbulence and temperature of the baryons on small scales.

For the test cases in Sect. 3, we explicitly checked the kinetic energy of the baryonic particles to understand how strong the numerically induced turbulence is. In general, the strength of the required viscosity to dampen the turbulence depends on the strength of the DM-baryon interactions. If artificial viscosity is used to dampen the numerically induced turbulence, it might need to only act on the smallest resolved flows, as larger scales are not likely to be affected by this numerical artefact.

Another relevant aspect is that the interactions artificially increase the temperature variation among the numerical baryonic particles. Thus, it could be favourable to have heat conduction that reduces this variation and prevents the internal energy of particles from approaching zero. Furthermore, it helps prevent negative internal energies (see also Sect. 2.2) and thus makes the numerical scheme more robust.

The problems described above limit the physical cases to which our numerical scheme can be applied. Astrophysical gases may not in general have a physical viscosity and heat conduction being large enough to mitigate the described issues when needed. Artificial viscosity and artificial heat conduction can be applied instead if they allow for approximating the underlying physical system well enough, for example, ensuring they do not destroy relevant physical turbulence.

2.4 Implementation

We implemented the scheme for IDM in the cosmological N-body code OpenGadget3, a successor of Gadget-2 described by Springel (2005). The SIDM implementation in OpenGadget3 was written by Fischer et al. (2021a), and is the basis for our implementation of IDM. Further improvements to the SIDM module are described by Fischer et al. (2021b, 2022, 2024b). We generalised the module and re-used it for the DM-baryon scattering. In consequence, we employ the same scattering routine which has the advantage of being capable of simulating the limit of very anisotropic cross-sections (fSIDM, see Kahlhoefer et al. 2014; Fischer et al. 2021a). However, also we took advantage of the explicit energy conservation in parallel computations for the scattering of the particles. The details of the parallelisation scheme, realised with the message passing interface (MPI), are explained in Fischer et al. (2024b).

The domain decomposition and the neighbour search in OpenGadget3 have been described by Ragagnin et al. (2016). Furthermore, the simulation code offers an SPH implementation (Beck et al. 2016) and an MFM scheme (Groth et al. 2023). We have implemented the interaction for both fluid schemes. This allowed us to study the differences between SPH and MFM and evaluate which is better suited for our purpose. In practice, the code does not save the internal energy but the entropy, this implies that we have to convert between entropy and internal energy as required.

For the scattering, we used a kernel assigned to each particle. The kernel size of the DM particles is determined by searching for the next N_ngb,DM neighbours. We search only for DM particles and ignore all other particle types. In the case of the baryonic particles, the scheme to model hydrodynamics already employs a kernel size, h_hydro, based on a neighbour number, N_ngb,hydro. However, h_hydro might be larger or smaller than what we preferably use for the DM-baryon interactions. Instead of performing another neighbour search, we rescale the kernel sizes for the interactions as we describe in the subsequent section (Sect. 2.4.1). These rescaled values are employed in Eq. (4) to compute the kernel overlap integral. In our simulations, we used N_ngb,DM = 64 for the DM particles, N_ngb,SPH = 230 for the SPH particles and N_ngb,MFM = 32 for the MFM particles. The first number was motivated by an aim to obtain a sufficiently accurate estimate of the local DM density, so that the kernel size rescaling will work as expected. It is likely that a somewhat smaller number is sufficient as well. The numbers for the hydrodynamic schemes are chosen in the range of commonly employed values, which are known to provide accurate results.

The adaptive time-stepping scheme in OpenGadget3 leads to active and passive particles, with the status of a particle being dependent on the time step that is computed. A detailed description can be found in the Gadget- 2 paper (Springel 2005). Given that not only active particles interact with each other, but also active-passive particle pairs are possible, the search for scattering partners must be performed from both directions. This means that we take all active DM particles and search for baryonic scattering partners, including active and passive particles. As we might have missed pairs consisting of an active baryonic particle and a passive DM particle, we perform another search around all active baryonic pairs to find DM particles. We note that the simulations presented in this paper employ a fixed time step for all particles implying that all particles are always active. The bidirectional search described above implies that we find every pair consisting of two active particles twice. However, for the interactions, we consider them only once per time step. This is realised with a criterion based on a unique identification number assigned to every particle.

2.4.1 Rescaling of the local kernel size

In our IDM simulations, we have two species that are interacting with each other, this makes controlling the number of interactions between the numerical particles more complicated than in the single species case usually being present in SIDM studies. The number of interactions does not simply depend on the kernel sizes chosen for baryons (h_bary) and DM (h_DM), but also on the local number density of the numerical particles of each species, namely, n_bary and n_DM. The expected maximum number of interactions that a particle could undergo per time step can be expressed as, $$N=\frac{4}{3}\pi (h_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}+h_{\mathrm{D}\mathrm{M}}{)}^{3}max(n_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}},n_{\mathrm{D}\mathrm{M}}).$$(13)

To effectively control the performance of the simulations it is important to control N. To do so, we introduce a factor, ξ, to obtain a rescaled kernel size of h^* = h ξ. We note that every numerical particle has its individual factor ξ. Altering the kernel sizes via ξ allows us to control N, our target value is N_idm. Initially we assume ξ = 1 and update ξ every time step, $${\xi }^{\prime }=\frac{1}{2h}\sqrt[3]{\frac{N_{\mathrm{i}\mathrm{d}\mathrm{m}}}{n_{\mathrm{m}\mathrm{a}\mathrm{x}}}},n_{\mathrm{m}\mathrm{a}\mathrm{x}}=max(\frac{N_{\mathrm{n}\mathrm{g}\mathrm{b},\mathrm{o}\mathrm{w}\mathrm{n}}}{h^3},\frac{N_{\mathrm{n}\mathrm{g}\mathrm{b},\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}}}{{h_{\mathrm{m}\mathrm{i}\mathrm{n}}}^3}).$$(14)

Here, we employ the minimum kernel size, h_min, that a particle has seen over the last time step. Moreover, N_ngb,own, refers to the number of neighbours that were used to set the kernel size, h, of the species under consideration, i.e. the one to which the particle for which we want to compute ξ′ belongs. Similarly, N_ngb,other, is the number of neighbours that were used for the other species and is relevant for h_min. We note that we derived Eq. (14) by assuming that h^* does only depend on the location but not on the species, namely, particles of different species should have the same rescaled kernel size if they are at the same location.

In practice, we employ the kernel size rescaling based on the DM and SPH/MFM kernel sizes. It allows us to significantly speed up the simulations, depending on the chosen value for N_idm. Importantly, the rescaled kernel sizes are also employed in the search for scattering partners.

A reasonable choice for N_idm might be somewhat larger than typical neighbour numbers used in SIDM simulations. We note, from Eq. (14) it follows for the single species case that N = 8 N_ngb. This may give an idea of how to compare N_idm to the neighbour numbers employed in SIDM simulations. For our simulations we use N_idm = 384 if not stated otherwise. This number is motivated by typical choices for SIDM. In addition, we tested to check that a larger number does not lead to a significant improvement, as described in Appendix C.

2.4.2 Time step criterion

Analogously to time step criteria for SIDM (Fischer et al. 2021b, 2024b), we derive the IDM time step criterion based on Eqs. (5) and (7) with the aim to limit the interaction probability or respectively the fractional velocity change. In the case of a velocity-independent cross-section, the interaction probability and the fractional velocity change increase monotonically with relative velocity. Consequently, we estimate the maximal scattering velocity that each particle could encounter. Therefore we compute for every particle, v_max = max(v + v_cut), by taking the maximum over all pairs. In the case of a velocity-dependent cross-section, there might eventually exist a finite velocity v_e for which the interaction probability and fractional velocity change become maximal. This is for example the case for the velocitydependent cross-section employed by Fischer et al. (2024b)⁸. In this case we can combine v_max and v_e to obtain a critical velocity v_c = min(v_max, v_e) to formulate the time step criteria. For a velocity-independent cross-section, we can simply set v_c = v_max.

The time step criterion for the frequent interactions, that is small-angle scattering, depends on the momentum transfer cross-section (see Eq. (6)) and is given by $$\mathrm{\Delta }{t}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{q}}=\frac{{\tau }_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{q}}}{m_{\mathrm{D}\mathrm{M}}{v}_c{\mathrm{\Lambda }}_{ix}}\frac{\mu }{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\left(\frac{{\sigma }_{\mathrm{T}}(v_c)}{m_{\chi }}\right)}^{-1}.$$(15)

Analogously, we can give the time step criterion for rare interactions, this is large-angle scattering, employing the total cross-section, $$\mathrm{\Delta }{t}_{\mathrm{r}\mathrm{a}\mathrm{r}\mathrm{e}}=\frac{{\tau }_{\mathrm{r}\mathrm{a}\mathrm{r}\mathrm{e}}}{m_{\mathrm{D}\mathrm{M}}{v}_c{\mathrm{\Lambda }}_{ix}}\frac{\mu }{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\left(\frac{\sigma (v_c)}{m_{\chi }}\right)}^{-1}.$$(16)

The factors τ_freq and τ_rare, allow the size of the time step to be controlled, namely, by limiting the fractional velocity change and the interaction probability, respectively. In contrast to the SIDM time step criterion employed by Fischer et al. (2024b), we do not use the self-overlap Λ_ii. In the case of two separate species when the kernel size is determined using particles of the same species only, it can happen that the kernel sizes differ between the species vastly, leading to a dramatic over or underestimate of the time step when using Λ_ii. Instead, we employ Λ_ix when computing the time step of particle i, where x refers to the particle with the smallest kernel size that particle i has seen over the last time step. For the computation of Λ_ix, however, we assume that both particles are located at the same position; this is the same as for the self-overlap but with two different kernel sizes⁹.

2.4.3 Considerations of the energy conservation

Although the DM-baryon interactions are formulated in a strict energy-conserving manner, this does not necessarily imply that they do not harm energy conservation. While SPH and MFM can be formulated in a strict conserving manner as well, the combination of an explicit conservative hydrodynamics scheme with IDM might result in a loss of conservation of total energy. The IDM scheme changes the velocity and in particular the internal energy of the baryonic particles. Within pure hydro schemes, those quantities are not expected to be altered at the stage where the IDM kicks happen. This can invalidate previously computed hydrodynamical accelerations and cause energy non-conservation.

The IDM case is much more complicated than the one of SIDM. In the latter, we only have to deal with the DM and its accelerations are computed from gravity only. In consequence, the DM acceleration does not contain any velocity-dependent term¹⁰, which would be affected by the SIDM kicks, but is purely based on the positions of the DM particles. In contrast, for IDM we might have velocity-dependent terms such as viscosity in the hydrodynamics scheme; importantly, however, the hydrodynamical accelerations depend on the internal energy, which is modified by the DM-baryon interactions.

As a consequence, it matters where in the time integration scheme the IDM interactions are executed. In OpenGadget3, we use a leapfrog scheme in the kick-drift-kick (KDK) formulation. For the simulations in Sect. 3 and Sect. 4, we implemented the IDM interactions between the two half-step kicks. But given that we alter the internal energy of the baryons we invalidate the acceleration used in the second half step-kick. As a consequence, energy is not explicitly conserved any more. We note, that the issue of energy conservation might be more complicated than what we describe here; for example additional complications may arise from using variable time steps and a wake-up scheme employed as part of the solver for hydrodynamics (e.g. Saitoh & Makino 2009; Pakmor et al. 2012). We leave a detailed investigation for future work. Nevertheless, it is in principle possible to improve on energy conservation by modifying the time integration while accepting higher computational costs as we discuss and show in Appendix B. Though for the test problems we simulated (see Sect. 3) this has not been of much relevance as the energy error is small enough to hardly affect the agreement between the numerical and exact solution.

3 Test problems

In this section, we study multiple test problems, including heat conduction and momentum transfer between baryons and DM. Firstly we consider a set-up where heat flows from the DM to the baryons and one with the opposite case where heat flows from baryons to DM (Sect. 3.1). Secondly, we study a problem in which DM and baryons initially move relative to each other and exchange momentum (Sect. 3.2). While we first limit those tests to equal mass ratios and small-angle scattering, we also test the implementation for an unequal mass ratio (Sect. 3.3) and for large scattering angles, namely isotropic scattering (Sect. 3.4).

In this work, we consider two differential cross-sections, a forward-dominated model and an isotropic one. The two crosssections are velocity-independent and the interactions are elastic. The forward-dominated cross-section is given by the limit where the transfer cross-section, σ_τ (Eq. (6)), is held constant while the scattering angles approach zero. This model could also be expressed as $${\frac{\mathrm{d}\sigma }{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}}|}_{\mathrm{f}\mathrm{w}\mathrm{d}}=\underset{ϵ\to 0}{lim}\frac{{\sigma }_{\mathrm{T}}}{8\pi \ln ({ϵ}^{-2})}\frac{1}{{({ϵ}^2+{\sin }^2{\theta }_{\mathrm{c}\mathrm{m}\mathrm{s}}/2)}^2}.$$(17)

A similar formulation but for identical particles has been used in several studies of SIDM, usually referred to as frequent scattering (e.g. Kahlhoefer et al. 2014; Kummer et al. 2019; Sabarish et al. 2024; Arido et al. 2025). The isotopic model can be formulated with the total cross-section, σ, $${\frac{\mathrm{d}\sigma }{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}}|}_{\mathrm{i}\mathrm{s}\mathrm{o}}=\frac{\sigma }{4\pi },$$(18)

and would for example follow from hard-sphere scattering.

3.1 Heat conduction

We simulate DM and baryons together, while both species interact with each other. However, at the same time we neglect gravitational forces. The test set-up for our heat conduction problem consists of a cube with a side length of 10 kpc and periodic boundary conditions. Both components initially have the same density and are at rest. The velocities of the DM component follow a Maxwell-Boltzmann distribution. The baryonic particles have no bulk velocity but a non-zero temperature. The one-dimensional velocity dispersion of the DM particles is ν = 2km s⁻¹ and the baryons have an internal energy that corresponds to 10% of the kinetic energy of the DM. We have 10⁵ DM particles and 46 656 baryonic particles arranged in a lattice with 36 particles per dimension. The particle numbers are chosen such that the baryonic particles are roughly twice as massive as the DM particles. Each component accounts for a mass of 10¹⁰M_⊙. This implies a DM density, ρ_DM = 10⁷ M_⊙ kpc⁻³ and a baryonic density, ρ_bary = 10⁷ M_⊙ kpc⁻³, while the energy densities are w_DM = 6 × 10⁷ M_⊙ km² s⁻² kpc⁻³ and w_bary = 6 × 10⁶ M_⊙ km² s⁻² kp c⁻³. To compute the DM-baryon interaction we employ a cubic spline kernel (Monaghan & Lattanzio 1985). It is used in Eq. (4) to compute the kernel overlap relevant for the interaction probability (Eq. (5)) and strength of the drag force (Eq. (7)). We simulate the set-up with an extremely anisotropic cross-section (Eq. (17)) and use a velocity-independent momentum transfer cross-section of σ_τ/m_χ = 10 cm² g⁻¹. In addition, we use the MPI parallelisation scheme for all tests of this Section.

In Fig. 2, we show the results for the heat conduction problem when using SPH with 230 neighbours and a Wendland C⁶ kernel (Dehnen & Aly 2012). Neighbour numbers of this size are a typical choice for our SPH implementation and have proven to give reliable results. The DM neighbour number is N_ngb,DM = 64 and the interaction number is N_idm = 384. The first one is chosen large enough to give a sufficiently accurate estimate for the local DM density and is used for the kernel size rescaling only. The second one controls the number of particles that a given particle could interact with per time step. In Appendix C, we show that increasing the number any further does not significantly improve the results. Moreover, we employed artificial viscosity and artificial head conduction (Price 2012; Beck et al. 2016). The artificial viscosity is formulated to act only against highvelocity divergence, aiming to leave rotating or shearing flows unchanged. The settings were chosen to be the same as those for the Magneticum simulations (Dolag et al. 2025). Hence, we do not expect that they would harm the formation of galaxies, for example, affecting the star formation rate. Furthermore, we note that the same fixed time step of Δt_SPH = 0.024 Gyr is employed for all particles, this is true for all SPH simulations of this section. For the size of the time step, the value implied by the time step criterion for the hydrodynamics is taken as an orientation (Beck et al. 2016). A smaller choice for the time step could help to improve the results, as we show in Appendix C. From the figure, we can see that the total energy (black) stays constant, but the energy of the DM particles (violet) is decreasing over time and the energy of the baryons (orange) is increasing. We note, that we separately show the kinetic energy from the bulk motion of the baryons, i.e. the kinetic energy of the numerical particles. It is supposed to be zero and as we can see, it increases only slightly over the course of the simulation.

The set-up for the heat flow from baryons to DM is the same as for the heat conduction from the DM to the baryons, but here we interchange the energies, i.e. the baryons contain 10 times the energy of the DM. We simulate the set-up with the same cross-section as before (σ_T/m_χ = 10 cm² g⁻¹ ) and thus expect the thermalisation to take place at the same speed.

In Fig. 3, we show the results for the problem with heat conduction from the baryons to the DM using the SPH implementation. Here the energy densities of DM and baryons are swapped compared to the previous set-up. We can see that the thermalisation takes place at roughly the same speed as for the set-up with heat conduction in the other direction. Any difference in the thermalisation speed is not physical but due to numerical errors, a comparison to the exact solution follows next with Fig. 4. Moreover, the kinetic energy of the baryonic particles stays small over the course of the simulation, a detailed investigation follows later with Fig. 5.

In the following, we compare the test simulations in larger detail. We did not only use SPH but also ran the same simulations using MFM with 32 neighbours. With this, we follow a common choice for MFM simulations (e.g. Gaburov & Nitadori 2011). Moreover, we employ the cubic spline kernel (Monaghan & Lattanzio 1985) for MFM and use Eq. (14) to rescale the kernel size for the DM-baryon interactions as we do for SPH too, while keeping the DM neighbour number at 64 as before. We note, that for MFM we do not employ artificial viscosity and artificial heat conduction, the MFM scheme itself already gives rise to numerical viscosity. In addition, we employ a fixed time step of Δt_MFM = 0.006 Gyr for all particles. This is a quarter of the value for the SPH simulations and used for all MFM simulations of this section. Again, we took the corresponding time step criterion as an orientation (Groth et al. 2023). Decreasing the time step can help to improve the results as we show in Appendix C.

In Fig. 4, we show the kinetic energy of the DM particles. Depending on the set-up it is increasing or decreasing over time. We find that a significant difference between SPH and MFM is present which is growing over time. The energy of the DM particles in the MFM runs is always lower than in the corresponding SPH runs, no matter in which direction the heat is flowing. This is related to energy non-conservation in the MFM runs as we see later.

Figure 4 also allows us to compare the simulation results to the exact solution indicated by the red lines. The analytic description follows from Dvorkin et al. (2014) (see also Munoz et al. 2015) and simplifies in our case with a velocity-independent cross-section to $$E_{\mathrm{D}\mathrm{M}}(t)=E_{\mathrm{D}\mathrm{M},\mathrm{e}\mathrm{q}}+(E_{\mathrm{D}\mathrm{M},\mathrm{i}\mathrm{n}\mathrm{i}.}-E_{\mathrm{D}\mathrm{M},\mathrm{e}\mathrm{q}})e^{-\kappa t}.$$(19)

Here, E_DM,ini is the initial energy of the DM and E_DM,eq = (E_tot r)/(1 + r), with E_tot being the total energy (i.e. the sum of DM and baryons). The speed at which the heat transfer happens is set by $$\kappa =\frac{8}{\sqrt{\pi }}{\left(\frac{2}{3}\right)}^{3/2}\frac{\sqrt{w_{\mathrm{t}\mathrm{o}\mathrm{t}}{\rho }_{\mathrm{t}\mathrm{o}\mathrm{t}}}}{1+r}\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }}.$$(20)

The total energy density w_tot = w_DM + w_bary is given by the sum of the energy density (i.e. energy per volume) of DM and baryons. Analogously the total matter density is ρ_tot = ρ_DM + ρ_bary. We note that the equations above are only valid when ρ_DM = ρ_bary, given that we made this assumption for obtaining a simpler expression (a more general formulation can be found in Appendix D).

Overall, the simulation results agree well with the exact solution, in particular the SPH runs (see Fig. 4). The MFM scheme agrees very well at the beginning of the simulation and gives abit too low energies at the late stages of the evolution. The coupling of IDM to the hydro schemes can lead to different numerical artefacts affecting how well the simulation results agree with the exact solution. In the following, we look closer into this.

Next, we study the kinetic energy of the baryonic particles. It is supposed to stay zero as DM and baryons are at rest and the interactions only heat up or cool down the baryons. However, in practice, our numerical scheme must allow for an exchange of momentum in each DM-baryon interaction. In the convergence limit, the kinetic energy of the baryons would stay zero for our heat conduction test problems. How much it deviates from zero can only (in a limited sense) be considered a measure of how accurate the simulations are (see Sect. 2.3). In Fig. 5, the kinetic energy of the baryonic particles is shown. It is visible that initially, the energy is increasing steeply for all runs. For the simulation with heat flow from the baryons to DM (green) a plateau is reached after the sharp increase and the SPH and MFM runs exhibit a very similar behaviour. In general, the difference between the simulated set-ups is larger then between the numerical schemes for modelling hydrodynamics. The simulations for the set-up with heat flow from DM to baryons reaches a peak about twice as high as for the test with the heat flow in the opposite direction. After reaching the peak, the kinetic energy of the baryonic particles declines and reaches almost the value of the other simulations. Here, viscous forces reduce the artificial small-scale turbulence.

The last aspect we studied with our simulations for the heat conduction problem is energy conservation. In Fig. 6, we show how well total energy is conserved as a function of time. Ideally, we expect energy to be perfectly conserved as indicated by the red line. However, in practice, it is subject to numerical error. We can see that the energy errors are monotonically increasing over the simulation time. Interestingly, the errors are much smaller for SPH compared to MFM, where for the later one they rise up to a few percent of the total energy. As a consequence, the SPH implementation appears to be preferable in terms of energy conservation, with energy errors below one percent. In Sect. 2.4.3, we discuss reasons for non-conservation, and in Appendix B, we demonstrate that the conservation of total energy can be improved for the case of SPH. Moreover, we ran tests varying the time step size and the interaction number N_idm, they are presented in Appendix C.

Lastly, we want to mention that we performed a test with a similar set-up in an expanding space. In Appendix D, we show the results and demonstrate that we can also accurately model IDM when using comoving integration.

Fig. 2 Heat conduction problem where energy flows from the dark matter to the baryons. Different types and components of the energy are shown as a function of time. The total energy of the system is illustrated in black. For DM the energy is shown in violet, it corresponds to the kinetic energy of the particles as other contributions are zero. In orange, we illustrate the energy of the baryons. It consists of the internal or thermal energy and the kinetic energy (bulk motion). The kinetic part is also displayed separately (dashed line). The simulation was conducted employing a forward dominated cross-section and SPH to describe the gas (baryons).

Fig. 3 Same as in Fig. 2, but for the heat conduction problem with energy flowing from the baryons to the DM.

Fig. 4 Kinetic energy of the DM as a function of time for both heat conduction set-ups. For heat flowing from the DM to baryons (DM→Gas) the results are given in blue and for heat flowing from baryons to DM (Gas→DM) in green. We display the results for both methods SPH (dotted) and MFM (dashed). The red lines give the exact solution according to Eqs. (19) and (20). The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution.

Fig. 5 Kinetic energy of the baryonic particles in units of the total energy is shown as a function of time. Labels, colours, and line types are the same as in Fig. 4. We note that this figure is basically a zoom-in on the orange dashed curves displayed in Figs. 2 and 3.

Fig. 6 Total energy in units of the initial energy as a function of time to study energy conservation. The results are displayed for the two heat conduction set-ups and both fluid methods, SPH, and MFM. Labels, colours and line types are chosen as in Fig. 4.

3.2 Momentum transfer

We considered a problem similar to the previous one, but now the DM and baryonic components are moving relative to each other and contain the same energy. Their matter densities are the same as in Sect. 3.1 (ρ_DM = ρ_bary = 10⁷ M_⊙ kpc⁻³). The onedimensional velocity dispersion of the DM is ν = 0.5 km s⁻¹ and the relative velocity between the two components is 4km s⁻¹. The corresponding energy densities (excluding the bulk motion) are w_DM = w_bary = 3.75 × 10⁶ M_⊙ km² s⁻² kpc⁻³. Thus, the two components have the same energy in their common centre-of-mass frame. As a result of having an equal mass ratio (r = 1), we do not expect any heat transfer to occur between the two components, but they should decelerate relative to each other and in consequence heat up, while the total energy of each component stays constant. Again we consider velocity-independent forward dominated scattering with a momentum transfer cross-section per DM particle mass of σ_τ/m_χ = 10 cm² g⁻¹.

Based on the work by Dvorkin et al. (2014), we can analytically estimate the solution for this test problem assuming a velocity-independent cross-section. In the limit where the relative velocity is much larger than the velocity dispersion of the two components, the evolution can be well described by a drag force decelerating the DM and baryons. The drag deceleration acting on the DM component is given by $$a_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}=-\frac{{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{v}_{\mathrm{r}\mathrm{e}\mathrm{l}}^2}{1+r}\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }}.$$(21)

In the regime where the velocity dispersion dominates over the relative velocity, the deceleration can be expressed as $$a_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p}}=-\frac{8\sqrt{2}}{3\sqrt{\pi }}\frac{{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{v}_{\mathrm{r}\mathrm{e}\mathrm{l}}}{1+r}\sqrt{{\nu }_{\mathrm{D}\mathrm{M}}^2+{\nu }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^2}\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }}{G}_0\left(\frac{v_{\mathrm{r}\mathrm{e}\mathrm{l}}^2}{{\nu }_{\chi }^2+{\nu }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^2}\right).$$(22)

Here, ${\nu }_{\mathrm{D}\mathrm{M}}^2$ and ${\nu }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^2$ are the one-dimensional velocity dispersion of the DM and baryons respectively. The function G₀ can be expressed as a series, $$G_0(X)=1+\frac{X}{10}-\frac{X^2}{280}+\frac{X^3}{5040}+...,$$(23)

where the variable $X=v_{\mathrm{r}\mathrm{e}\mathrm{l}}^2/({\nu }_{\chi }^2+{\nu }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}^2)$. To estimate the deceleration we compute G₀ up to the third order and interpolate between the two regimes, $$\frac{\mathrm{d}{v}_{\mathrm{D}\mathrm{M}}}{\mathrm{d}t}=a_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}f(X)+a_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p}}(1-f(X)).$$(24)

For the interpolation, we used a logistic weighting function, $$f(X)=\frac{1}{1-e^{a(1-bX)}}.$$(25)

For the comparison, we employ the following parameters for the weight, a = 3.0 and b = 0.1. Moreover, we note that in our set-up ρ_DM = ρ_bary and thus the two components experience the same deceleration. To obtain the time evolution of the relative velocity between the two components we integrate Eq. (24) numerically.

For the simulations, we use the same numerical parameters as for the heat conduction problem and run the test with SPH and MFM. In particular, we employ the same resolution (N_DM = 10⁵, N_bary = 46 656), the same neighbour numbers (N_ngb,DM = 64, N_ngb,SPH = 230, N_ngb,MFM = 32, N_idm = 384), the same kernel functions and the same fixed time step (Δt_SPH = 0.024 Gyr, Δt_MFM = 0.006 Gyr). Again, the SPH run is conducted using artificial viscosity and artificial heat conduction.

In Fig. 7, we display the relative velocity between the DM and the baryons as a function of time and compare it to our analytic estimate. As expected the relative velocity decreases as the baryon-DM interactions lead to a momentum transfer between the two components. The simulation results agree well with the analytic estimate and the two hydro schemes SPH and MFM behave very similarly. We note that the accuracy of the analytic estimate, given that it is an interpolation of two regimes, might account for some of the deviation between numerical and analytic results.

We also show the kinetic energy of DM in Fig. 8. It is expected to remain constant over time as the energy of the bulk motion is converted into random motion. However, in the beginning, we see an increase in the kinetic energy, which is more pronounced for SPH. Later, the DM significantly loses energy in the MFM run, which is related to energy non-conservation as we see next. In contrast, the kinetic energy of the DM for the SPH run stays rather constant or slightly decreases at this stage.

Finally, we studied the energy conservation for this test problem. In Fig. 9, we can see that the energy is better conserved for SPH than for MFM. This is in line with our previous findings for the heat conduction problem illustrated in Fig. 6.

Fig. 7 Relative velocity between DM and baryons for the momentum transfer problem. The time evolution for our SPH and MFM simulations (black) as well as the analytic estimate (red) of the relative velocity are displayed. The latter is computed according to Eq. (24). The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution. We note, that the SPH and MFM results are very similar.

Fig. 8 Kinetic energy of the DM particle as a function of time for the momentum transfer problem. The dotted lines give the results using SPH and the dashed ones are for MFM.

Fig. 9 Energy conservation for the momentum transfer problem as a function of time. The total energy divided by the initial energy is given for SPH (dotted line) and MFM (dashed line).

3.3 Unequal mass ratio

In addition to the tests where particles with an equal mass ratio scatter, we also consider a problem where the physical DM particle has only 1/1000th of the mass of its baryonic scattering partner, i.e. r = 1000. Further we choose, f_bary = 1 and each mass component makes up 10¹⁰ M_⊙. We have an unequal number of baryonic and DM particles, N_bary = 46 656 and N_DM = 10⁵. As for the other test problems, we have a periodic box with a side length of 10 kpc containing a constant density (ρ_DM = ρbary = 10⁷ M_⊙ kpc⁻³). The system is evolved at a fixed time step of Δt = 0.015 Gyr. Small-angle scattering is simulated and a momentum transfer cross-section per DM particle mass of σ_T/m_χ = 1000 cm² g⁻¹ is employed. Initially, the two components contain the same energy, but due to the scattering with an unequal mass ratio energy is transferred from the baryons to the DM. In the equilibrium state, which is asymptomatically approached by the system, the DM has an energy that is r times, namely, 1000 times, the energy of the baryons. The initial energy of the DM is given by a one-dimensional velocity dispersion of ν = 2km s⁻¹. This implies that the initial energy densities are w_DM = w_bary = 6 × 10⁷ M_⊙ km² s⁻² kpc⁻³; whereas in the equilibrium state, they become w_DM = 11.988 × 10⁷ M_⊙ km² s⁻²kpc⁻³ and w_bary = 0.012 × 10⁷ M_⊙ km² s⁻² kpc⁻³.

In Fig. 10, the evolution of the different energy components is shown. As expected the energy of the DM particles increases over the course of the simulation while the internal energy of the baryons is decreasing. A detailed comparison of the simulation results to the exact solution given by Eqs. (19) and (20) follows in Fig. 11. Here, we can see that the simulation results agree well with the analytic description of the test problem. Thus we can conclude that we are able to simulate fairly unequal mass ratios.

However, although these results are promising, we have to note that due to the stochastic nature of the simulation scheme, some of the baryonic particles could reach internal energies close to zero way before the average internal energy of the baryons has reached similar values. This poses a limitation on what can be simulated, it makes the set-up prone to the negative internal energy problem discussed in Sect. 2.2 (see also Sect. 2.3). In practice, this problem is prevented by the baryons being subject to a strong enough heat conduction (among the baryons). If we would run the test without artificial heat conduction, it would not have been possible to evolve the system as far as shown in Figs. 10 and 11.

Fig. 10 Same as in Fig. 2, but for a heat conduction problem with an unequal mass ratio evolved using SPH.

Fig. 11 Kinetic energy of the DM particles as a function of time. The simulated (black) time evolution for the heat conduction problem with a mass ratio of r = 1000 is compared to the exact solution (red) from Eqs. (19) and (20). The upper panel gives the absolute values and the lower panel displays the ratios to the exact solution.

Fig. 12 Energy evolution for isotropic scattering with equal particle masses. The different energy components for DM and baryons are shown as a function of time. This is the same test problem as shown in Fig. 2 for small-angle scattering.

3.4 Isotropic scattering

We also run a simulation with an isotropic cross-section (Eq. (18)) and equal physical particle masses (r = 1). Following the explanations in Sect. 2.2, we choose the numerical baryonic particle mass to be larger than the numerical DM particle mass to avoid non-positive internal energies. For this test, we set the baryonic particle mass to about ten times the DM mass. We simulate a total cross-section per DM particle mass of σ/m_χ = 20 cm² g⁻¹, this corresponds to a momentum transfer cross-section of σ_τ/m_χ = 10cm² g⁻¹, which is relevant for the heat exchange between the DM and baryons. Furthermore, we assume f_bary = 1. The number of DM particles is N_DM = 10⁵ and the number of baryonic particles is N_bary = 9261, their masses are m_DM = 10⁵ M_⊙ and m_bary = 1.0798 × 10⁶ M_⊙. The corresponding matter densities are ρ_DM = ρ_bary = 10⁷ M_⊙ kpc⁻³ and the initial energy densities are the same as in Sect. 3.1 (w_DM = 6 × 10⁷M_⊙ km² s⁻² kpc⁻³ and w_bary = 6 × 10⁶ M_⊙ km² s⁻²kpc⁻³). This implies a one-dimensional DM velocity dispersion of ν = 2kms⁻¹. We employ the same neighbour numbers as before and evolve the system at a fixed time step Δt = 0.024 Gyr.

The simulation results are shown in Fig. 12. Overall they look promising as we do not find a drastic increase in total energy, it is conserved up to 0.1%. Moreover, the kinetic energy of baryons stays also rather low and the baryonic and DM components evolve smoothly towards the equilibrium state.

In Fig. 13, we display the kinetic energy of the DM as a function of time and compare it to the exact solution from Eqs. (19) and (20). As we can see, the simulation results agree well with the analytic description of test problem. Hence, we can conclude that we are not only able to simulate small-angle scattering but large-angle scattering as well, in this case, an isotropic cross-section.

Fig. 13 Time evolution of the kinetic energy of the DM. The simulation results (black) for the heat conduction problem with isotropic scattering are compared to the exact solution (red) from Eqs. (19) and (20). The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution.

4 Halo formation

In this section, we study the collapse of an overdensity and the formation of the halo under the influence of DM-baryon interactions. We begin with a description of our simulation set-up and subsequently show and discuss our results. This time we only use SPH as it gave better results for the test problems of the previous section.

4.1 Simulation set-up

To simulate the formation of a M ≈ 10¹² M_⊙ halo from the collapse of an overdensity we employ an idealized set-up. For the ICs, we sample only a single spherically symmetric overdensity at z_i = 1000. We compute the density distribution following Lu et al. (2006) (see also Sect. 3.3 by Nadler et al. 2017). Specifically we employed the universal form for halo mass accretion histories as found by Wechsler et al. (2002) from a fit to simulations, $$M(z)=M_0\exp \left[\frac{-S}{1+z_c}(\frac{1+z}{1+z_f}-1)\right].$$(26)

At redshift z the virial mass is given by M₀, while redshift z_c characterizes the point in time at which the mass accretion rate d(log M)/d(log a) falls below a critical value S = 2. The linear overdensity δ_i at redshift z_i for a perturbation of mass M that collapses at z is $${\delta }_i(M)=1.686\frac{D(z_i)}{D(z(M))}.$$(27)

For the linear growth factor D(z) we use the fitting formula by Carroll et al. (1992), $$D(z)=\frac{g(z)}{1+z},$$(28) $$g(z)\approx \frac{5}{2}\frac{{\mathrm{\Omega }}_{\mathrm{M}}(z)}{{\mathrm{\Omega }}_{\mathrm{M}}^{4/7}(z)-{\mathrm{\Omega }}_{\mathrm{\Lambda }}(z)+[1+\frac{{\mathrm{\Omega }}_{\mathrm{M}}(z)}{2}][1+\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}(z)}{70}]}.$$(29)

The cosmological parameters Ω_Μ(z) and Ω_Λ(z) denote the density parameter of non-relativistic matter and of the cosmological constant respectively at redshift z. For the mass, M, enclosed within a radius r_i at redshift z_i we use $$r_i(M)={\left\{\frac{3M}{4\pi \overline{\rho }(z_i)[1+{\delta }_i(M)]}\right\}}^{1/3},$$(30)

where $\overline{\rho }(z_i)={\rho }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}{\mathrm{\Omega }}_{\mathrm{M}}(0)(1+z_i{)}^3$, with the critical density ρ_crit,0 at z = 0. We sample the spherical overdensity by creating equal mass bins with radii according to Eq. (30). Next, we approximate the density within the bins by a power law and enforce continuity. The positions are sampled in a Monte Carlo fashion using direct sampling and are rearranged on spherical shells to reduce density fluctuations. For the simulations in this section we use z_i = 200, z_c = 3.0, z_f = 0.0, and M₀ = 10¹² M_⊙ to generate the ICs.

We embed the overdensity in a cubic box with a comoving side length of l_box = 3698.4 ckpc. The volume around the overdensity is filled with a constant density ρ(r_i(M₀)). The positions of the corresponding particles are sampled employing a Monte Carlo approach with rejection sampling. In total the ICs contain N_bary = 2.5 × 10⁵ and N_DM = 2.5 × 10⁶ particles. Initially, the velocity of all the numerical particles is set to zero, which implies that DM is cold. The temperature of the SPH particles is set to T = 547.73 K assuming a mean molecular weight of μ = 1.2 m_p for a neutral gas with primordial abundances of hydrogen and helium. Adopting the Planck 2018 results (Planck Collaboration I 2020), we employ the following cosmological parameters: H₀ = 67.66 kms⁻¹ Mpc⁻¹, Ω_M0 = 0.3106, Ω_Λ⊙ = 0.6894, Ω_B0 = 0.0489.

A comoving softening length of ε = 2.96 ckpc is in place and we use a fixed time step for all particles to achieve high accuracy. We note that the time stepping is done in η = log(a), with a being the scale factor. This implies that the time steps in physical units are not equal, but Δη is constant. However, given that the required size of Δη decreases substantially over the course of the simulation, we decrease the time step a few times over the course of the simulation for all particles when required. Additionally, we use artificial viscosity and artificial heat conduction as in the previous section for the SPH simulations. We do not expect that these terms have a large global effect, as they are implemented time and spatial dependent to act against local discontinuities only, e.g. the artificial heat conduction does not lead to a coherent heat flow as the gravitational forces are taken into account (the implementation has been described by Beck et al. 2016). Furthermore, we employ our default neighbour numbers: N_ngb,SPH = 230, N_ngb,DM = 64, and N_idm = 384. The simulations are executed using the MPI parallelisation of the IDM scheme.

We simulate the collapse of the overdensity with collisionless DM and a simple IDM model with a velocity-independent forward-dominated cross-section (Eq. (17)) and assume r = 1 and f_bary = 1. We chose the cross-section relative to the CMB constraints by Boddy & Gluscevic (2018). It corresponds to a cross-section of σ_τ/m_χ ≈ 0.1 cm² g⁻¹ for velocity-independent equal-mass scattering. The cross-sections we simulate are ${\sigma }_{\mathrm{T}}/m_{\chi }\in \{0.1,1.0,10.0\}({\mathrm{c}\mathrm{m}}^2{\mathrm{g}}^{-1})$.

In contrast to the tests of the previous section, we use here for our simulations comoving integration in an expanding space, which we tested (described in Appendix D). Moreover, we note that the local kernel size rescaling as described in Sect. 2.4.1 is very helpful to speed up these simulations.

Fig. 14 Density profile for a collapsing overdensity in the case of collisionless DM. The density is shown as a function of radius at several redshifts for DM (solid) and baryons (dashed). The density is not displayed for small radii where the particle number is too low to obtain a reliable value.

4.2 Results

In the following, we present the results of our simulations of a collapsing overdensity. To analyse the simulations, we make use of the peak finding algorithm based on the gravitational potential presented by Fischer et al. (2021b). It gives us the centre with respect to which we compute other quantities such as the density profiles.

Figure 14 gives the comoving matter density for the DM and the baryons as a function of comoving radius. The collisionless case is shown. Overall the comoving density is increasing for the DM and baryons with time as the overdensity collapses under gravity. The two components evolve very similarly, except that the density of the baryons is not growing as much as for the DM at the halos centre, because of the thermal pressure. We note that we did not include gas cooling or star formation. This allows us to compare CDM and IDM runs better and directly infer the impact of the DM-baryon scattering without being affected by additional physical processes as we do next. However, this limits the physical fidelity at the same time.

At high redshifts, the collisionless and interacting DM models behave very similarly and it is only at later times that the DM-baryon interactions lead to lower DM densities, compared to the collisionless case. As shown in the upper panel of Fig. 15 for z = 0, this can be understood on the basis of two mechanisms. Firstly, at small radii, where the baryons in the absence of interactions are hotter than the DM, the scattering leads to a flow of heat from the baryons to the DM. Secondly, the velocity dispersion gradient of the DM and the temperature gradient of the baryons are positive at small radii, i.e. velocity dispersion and temperature increase with radius (lower panel)¹¹. Similarly to SIDM halos the interactions can give rise to heat transport following the velocity dispersion/temperature gradient. Given that the gradient is positive, heat is transported inward, which contributes to the formation of a density core.

Interestingly, we find for the cross-section of σ_τ/m_χ = 1.0 cm² g⁻¹ that the baryon density at small radii increases (upper panel of Fig. 15). This can be understood with the heat exchange between DM and baryons at these radii. The interactions effectively cool the baryons (middle panel of Fig. 15), which makes them contract and leads to a higher baryon density in the centre. For the larger cross-section (σ_τ/m_χ = 10.0 cm² g⁻¹) we do not find this higher baryon density, because much more heat from larger radii is transported inwards leading to a higher baryon temperature compared to the σ_τ/m_χ = 1.0 cm² g⁻¹ case. The smallest cross-section of σ_τ/m_χ = 0.1 cm² g⁻¹ has only a small impact. It decreases the DM density and increases the baryonic density due to effectively cooling the baryons (Fig. 15).

In summary, we find that the interplay of local heat exchange between the two components and heat inflow from larger radii leads to a non-monotonic behaviour of the central baryon density and temperature with cross-section. In contrast, the central DM density and its velocity dispersion behave monotonically with cross-section for the range we are studying here.

We want to note that we additionally simulated the set-up with σ_τ/m_χ = 10.0 cm² g⁻¹, using the scheme for improved energy conservation (Appendix B). As visible in Fig. 15, the default scheme does not deviate significantly from the improved scheme for our set-up of the collapsing overdensity. In conclusion, the default scheme is accurate enough for these simulations.

In addition to density and velocity dispersion, we also computed the entropic function, $$A=(\gamma -1)u{\rho }^{1-\gamma },$$(31)

which is closely related to entropy. It is shown in the lower panel of Fig. 15. Here, we also find the non-monotonic impact of the DM-baryon interactions with cross-section on the baryon density.

Besides the comoving density profile, we also measure the average density and specific energy within a comoving radius of r = 2.7 ckpc, and compute the corresponding interaction rate of the DM particles. The results for our simulations are shown in Fig. 16. We note, that the shown specific energy and densities are not comoving but physical. For the specific energy, u, we compute the sum of the kinetic and internal energy divided by mass. This specific energy is then used to compute the interaction rate, R_DM, of the DM particles shown in the lower panel. For a velocity-independent cross-section, it is given by, $$R_{\mathrm{D}\mathrm{M}}=4\frac{\sigma }{m_{\chi }}{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}\sqrt{\frac{u_{\mathrm{D}\mathrm{M}}+u_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{3\pi }}.$$(32)

We note for the simulated forward-dominated model (Eq. (17)) the total cross-section σ does not have a finite value. Instead, we assume that the impact of the DM-baryon interactions could be matched to isotropic scattering by employing the momentum transfer cross-section (Eq. (6)). For isotropic scattering the two cross-sections are related to each other through σ|_iso = 2 σ_T|_iso.

In the upper panel of Fig. 16 we can see that the central density for DM and baryons decreases with time. At the same time, the specific energy (middle panel) increases and reaches very similar values for DM and baryons. We note that we only show the late part of the evolution where the DM-baryon interactions have a significant impact on the density profile. At the earlier times, which are not shown, CDM and IDM behave almost the same. In principle, the DM interaction rate is fairly high at early times but decreases very quickly and stays for almost all of the cosmic evolution much lower. Given that R_DM decreases very quickly, the high values at early times have very little impact. In the lower panel of Fig. 16, we show only the interaction rate at later times but also give the expected number of scattering events a particle has undergone since z = 200. For the collisionless case, the interaction rate is always zero.

Overall, our simulations demonstrate that we can model the impact of DM-baryon interaction on physically interesting systems. The most noticeable effects result from cross-sections that are large compared to the CMB bounds by, for instance, a factor of 100 (σ_τ/m_χ = 10.0 cm² g⁻¹). While such a large crosssection can lead to a sizeable density core, a cross-section of σ_T/m_χ = 0.1 cm² g⁻¹, which roughly corresponds to the CMB bound, leads to a much more subtle effect, but could still increase the central baryon density by a factor of two. To more fully unpack the effects of DM-baryon scattering at low redshifts, several steps beyond the scope of our work are necessary. First, we only studied one system with a mass of 10¹² M_⊙, but this could be extended to other systems as well. Second, we neglected various baryonic processes. The inclusion of gas cooling would lead to higher densities in the centre and thus increase the interaction rate (at least at a fixed relative DM-baryon velocity), whereas, for example, supernova feedback could have the opposite effect. However, when taking the velocity dependence into account, the picture could become more complicated. This leads us to the third point, we only considered a velocity-independent crosssection with scattering partners having equal masses. This could be extended to a velocity-dependent cross-section with unequalmass scattering. Fourth, we assumed that there is no relative motion between the DM and the baryons for the initial conditions of our simulations (at z = 200). However, there could be a significant relative velocity (e.g. Tseliakhovich & Hirata 2010; Dvorkin et al. 2014). If this is the case, the impact of the DM-baryon interactions might be stronger than what one would expect based on our simulations. IDM simulations like ours can be a valuable tool for deriving bounds on the interaction cross-section, particularly for the impact of DM-baryon scattering at low redshift. Finally, we assumed ICs (in the form of an initial density perturbation for halo collapse) consistent with CDM, whereas IDM can significantly suppress the linear matter power spectrum and thus delay (or even prevent) the formation of halos below a critical mass that depends on the IDM mass and cross-section (Nadler et al. 2019; Nadler et al. 2025a). We leave a self-consistent treatment of halo formation in IDM to future work.

Fig. 15 Density, velocity dispersion, and entropy profile of DM and baryons for CDM and IDM. We show the density as a function of radius (upper panel) at a redshift of z = 0. The velocity dispersion, including kinetic and internal energy but excluding radial bulk motion, is displayed in the middle panel. The results for collisionless DM are shown in black. In the lower panel, the entropic function (Eq. (31)) is displayed. The solid lines are for the DM and the dashed ones are for the baryons. Furthermore, the coloured curves are for IDM simulations employing different cross-sections. In addition, the results for a simulation that makes use of the scheme for improved energy conservation as explained in Appendix B are given too.

Fig. 16 Central density, specific energy and DM interaction rate as a function of scale factor. For the simulated DM models we show the evolution of different quantities. The central density (upper panel) and specific energy (middle panel) within a radius of r = 2.7 ckpc are shown for the DM and baryonic components. The lower panel gives the average DM interaction rate within the same radius according to Eq. (32) assuming an isotropic cross-section matched to the simulated forward dominated one with the momentum transfer cross-section (σ|iso = 2σT|iso).

5 Discussion

In this section, we discuss various aspects of our numerical scheme for IDM and the simulations we ran. In particular, we highlight the remaining challenges and discuss potential next steps for investigating the numerics and physics of DM-baryon interactions.

In Sect. 3, we tested our IDM implementation with two schemes for hydrodynamics, namely SPH and MFM. We found that those schemes show significant differences for the tests we conducted. These differences stem from how the fluid equations are solved. For example, to conduct an SPH simulation we need artificial viscosity while for MFM this is not the case. Moreover, we found differences in the capability of conserving the total energy, which is related to the time integration, as we discuss in Sect. 2.4.3 and Appendix B. There exist further numerical schemes to solve the fluid equations, for example, there are a variety of advanced SPH schemes (for a review see Zhang et al. 2022). It is left for future work to understand how those schemes can be best combined with IDM and to investigate which of them might be best suited for simulations incorporating DM-baryon interactions.

We have found that viscosity and heat conduction can play a crucial role in making the numerical scheme robust. In particular, we discuss the fact that a negative internal energy can arise if the interactions are strong compared to viscosity and heat conduction (see Sect. 2.2). Depending on the numerical scheme for hydrodynamics this can make the use of artificial viscosity and artificial heat conduction crucial. Ideally, this is done in a fashion that acts against local discontinuities only and does not alter the properties on larger scales of the simulated system. In this line, we used the implementation described by Beck et al. (2016). Moreover, we found that the numerical formulation of the DM-baryon interactions introduces small-scale turbulence, which can be suppressed by viscosity, but must not harm the overall evolution of the system as turbulent pressure compensates for the missing thermal pressure. However, at the same time, this limits the interpretation of the hydrodynamical properties as we might not be able to clearly distinguish between the temperature and small-scale motion of the baryons.

Another aspect that can be crucial in making the scheme numerically robust is the choice of the numerical particle masses. In general, it is preferable to choose the masses for the baryons larger than for DM. This becomes particularly important for large scattering angles (see Sect. 2.2). Unfortunately, this stands in contrast to the way particle masses are typically chosen. For cosmological simulations of galaxy formation, it is common to choose the mass of the baryonic particles smaller than the one of the DM particles. To achieve the mass resolution for the baryons, many more DM particles would be needed making the simulation more expensive. The baryonic particle mass should also not be chosen too large because gravitational interactions can then lead to artificial mass segregation (e.g. Ludlow et al. 2019; Ludlow et al. 2023).

We found that numerically ensuring energy conservation is more challenging for IDM than for SIDM simulations. These issues are related to the fact that for IDM we also change the internal energies of the baryonic particles and not only the velocities, potentially interfering with the time integration. At least for SPH, we have been able to demonstrate that it is possible to improve energy conservation by modifying the time integration scheme with a moderate increase in computational costs (see Appendix B). However, in the case of MFM, it is left for future work to investigate how improvements in energy conservation could be achieved. An important aspect in this context is the use of variable time steps. So far, we have restricted ourselves to fixed time steps because we found that variable time steps can lead to a somewhat lower accuracy. It remains an open question how the time integration scheme can be improved when using variable time steps. This is left for future work, but it constitutes an important aspect given that it would allow for drastically reduced computational costs for typical astrophysical systems while maintaining high accuracy.

Overall, the aspects we mention in this paper indicate that there are possibilities and parameters that can be tuned to achieve a high accuracy of the IDM simulations. However, at the same time, they come with an increase in computational costs. One may have to find the best trade-off in respect of the specific problem at hand. There might not be simple guidance given the non-trivial interplay of the different aspects such as the size of the time steps, the numerical particle mass, the interaction number N_idm, viscosity, heat conduction, and the size of the scattering angles.

We note that the inclusion of additional processes could interfere with the numerical scheme and degrade its numerical properties, leading to larger numerical errors. We expect that simulating gravity along with the DM-baryon interactions and hydrodynamics is unproblematic as long as the gravitational accelerations are computed based on the positions only, in contrast to higher-order schemes that employ the particle velocities as well. This is in line with the finding by Fischer et al. (2024a) that the inclusion of their scheme for DM self-interactions does not hurt the vanishing drift in total energy for a symplectic and time-symmetric leap-frog integrator. However, problems that are even present without IDM are not expected to vanish when it is included. This can be seen analogously to the difficulties found in simulations of gravothermally collapsing SIDM halos. The main problems that show up are already present in a gravity-only simulation, but are amplified in the collapsing halos due to the high concentration of these objects (Fischer et al. 2024a).

On a similar note, the simulations of a collapsing overdensity presented in Sect. 4 are challenging. They are resolved by 2.75 × 10⁶ particles. Although one may naively assume that this resolution is sufficiently high, we want to note that the identified minimum in the gravitational potential, which we use as the centre, fluctuates. Hence, inferred quantities suffer from this noise.

However, simply using the centre as specified when generating the ICs does not help. This problem stems from artificial fragmentation and improves when increasing the resolution. That artificial fragmentation can be problematic in N-body simulations is well known and can in principle be resolved by using a different numerical scheme (e.g. Hahn et al. 2013; Hahn & Angulo 2015).

In this paper, we consider only the case where DM particles scatter off one species of baryonic particles. However, an interesting case is scatterings of DM with the nuclei of Hydrogen and Helium, the two most abundant elements in the Universe. It would be straightforward to extend our numerical scheme to this case. For each species, the interactions can be computed following our description and employing the correct parameters, such as the mass ratio r and the mass fraction of the baryonic species f_bary. The scattering routine that we implemented would then simply be executed twice per time step, for each species once using the corresponding parameters. Besides, the velocity distribution (see Eq. (2)) from which the velocity of the virtual particle is drawn needs to be adjusted in the case of baryonic species with unequal masses such as Hydrogen and Helium. Moreover, each species would have its own time step criterion.

Our numerical method could be applied to various astrophysical systems, as we know it from SIDM, to constrain the DM-baryon interactions. The following systems could prove particularly interesting for IDM:

• a)

  have large difference in temperature between the DM and baryonic component:

• b)

  show a significant relative motion between the DM and the baryons.

Systems falling into category a) can give rise to a significant heat flux between DM and baryons due to the DM-baryon scattering. Studies in this regard have been conducted for galaxy clusters (e.g. Shoji et al. 2024; Stuart & Pardo 2024). Potentially also less massive systems could be of interest. An example for a system that falls into category b) would a dissociative galaxy cluster merger. Here the DM can move with a relatively high velocity relative to the intra-cluster medium. IDM could eventually produce small offsets between the distribution of galaxies and the DM component as it is known from SIDM (e.g. Robertson et al. 2017; Fischer et al. 2023; Valdarnini 2024).

When extending the application of the IDM scheme to other physical systems and higher fidelity of modelling them, it could become relevant to include gas cooling as well as stars and black holes in the simulations and their associated processes. Those processes when happening on scales that cannot be resolved are termed ‘subgrid physics’ and are described by effective models, for example, star formation. Depending on how these models are formulated, they could plausibly interfere with the DM-baryon interactions. As we discuss in Sect. 2.3, our scheme could lead to artificial small-scale turbulence and artificially large variations in the internal energy, potentially affecting subgrid models. Investigating the compatibility with various subgrid models is left to a future work.

6 Conclusions

In this work, we present a novel method for simulating interactions between DM and baryons within the framework of the N-body method. It allows us to simulate a wide range of models with different mass ratios of the interacting particles, comprising a cross-section of varying angular and velocity dependencies. Our scheme allows for the simulation of typical astrophysical set-ups and we have demonstrated its accuracy with several test problems introduced for this purpose. In addition, we simulated the collapse of an overdensity to study the effect of IDM on halo formation. Furthermore, we have discussed the limitations of the scheme and directions for future development. Our main conclusions are as follows:

1. It is possible to simulate various models of IDM self-consistently within N-body simulations. In particular, we are able to simulate unequal-mass scatterings more easily than in state-of-the-art SIDM schemes;

2. Artificial viscosity and heat conduction can play an important role in ensuring numerical stability for simulations involving a large cross-section compared to the viscosity and heat conduction of the baryons;

3. For simulations of low-viscosity systems, it might not always be possible to clearly distinguish between the turbulence and temperature of the baryons on small scales;

4. With our simulations of a collapsing overdensity, we have demonstrated that it is possible to simulate a physically relevant system and found that a velocity-independent crosssection at the level of the CMB bounds could exert only a small dynamical impact on the collapse of a Milky Way-mass halo; Since we made several highly simplifying assumptions, such as neglecting the gas cooling, a more sophisticated study is needed to make a definitive assessment;

5. Moreover, we found in those simulations that the central baryon density and temperature do not respond monotonically to an increase in the cross-section, which is in contrast to the behaviour of the DM (Fig. 15).

This paper constitutes only a first step towards exploring the full astrophysical phenomenology of DM-baryon interactions by providing a numerical method. Future studies that use complementary IDM simulation methods (Wen et al., in preparation) and explore the effects of these interactions on diverse astrophysical systems could help constrain the cross-section of DM-baryon interactions and shed new light on the nature of DM.

Acknowledgements

The authors thank Ludwig Daniel Schmidt for his feedback on the analytic solution to the heat conduction problem. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 “Origins” - 390783311. KD, MSF and FG acknowledge support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. Software: NumPy (Harris et al. 2020), Matplotlib (Hunter 2007).

Appendix A Derivation of interaction probability and drag force

In this appendix, we derive the interaction probability (Eq. 5) for the rare scattering scheme as well as the drag force (Eq. 7) for the frequent scatterings in the case of DM-baryon interactions. Our derivations follow closely the ones presented by Fischer et al. (2021a) for the case of SIDM.

A.1 Interaction probability

First, we derived the interaction probability by starting from the probability of aDM particle scattering about a baryonic particle. The probability, P_scatter, depends on the total cross-section, σ(ν), which can be a function of the scattering velocity, v. The DM particle is assumed to travel through a constant baryonic density ρ_bary for a time t. In case not all baryonic particles are interacting with the DM, f_bary gives the mass fraction of particles taking part in the interaction. The scattering probability is, $$P_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}}=\sigma (v)\frac{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{r{m}_{\chi }}vt.$$(A.1)

We note that m_χ denotes the mass of the DM particles and r is the mass ratio between the physical baryonic particle, for instance, a proton, and the physical DM particle.

The scattering velocity is ν = |v_i - v_j|, where one of the particles is the numerical DM particle and the other one is the virtual particle. We can compute the expected number of scatterings 〈N〉 between DM and baryons by integrating over the DM, ρ_DM, and baryon, ρ_bary, density. $$⟨N⟩=\int \frac{{\rho }_{\mathrm{D}\mathrm{M}}{f}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{r{m}_{\chi }^2}\sigma (v)v\mathrm{\Delta }t\mathrm{d}V.$$(A.2)

The densities can be interpreted as the mass represented by the numerical particles, i.e. ρ_i = m_i W(|x − x_i|, h_i). Using Eq. (4), we can write, $$⟨N⟩=N_{\mathrm{D}\mathrm{M}}{N}_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}\sigma (v)v\mathrm{\Delta }t{\mathrm{\Lambda }}_{ij}.$$(A.3)

Here, N_DM denotes the number of physical particles represented by the numerical particles. But N_bary = f_bary m_bary/(rm_χ) is only the number of baryonic particles that could scatter. As each scattering event involves one physical DM particle and one baryonic particle, N_DM = N_virt = m_virt/(rm_χ) must hold. With the expected number of scattering events, we can easily express the probability, P_ij = 〈N〉/N_DM, for two numerical particles to interact. $$P_{ij}=\frac{\sigma (v)}{m_{\chi }}\frac{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\mu }{m}_{\mathrm{D}\mathrm{M}}v\mathrm{\Delta }t{\mathrm{\Lambda }}_{ij}.$$(A.4)

A.2 Drag force

Second, we provide a derivation of the drag force for the DM-baryon interactions. As done by Kahlhoefer et al. (2014), we start with a physical DM particle travelling through a background density ρ. When the DM particle interacts with a baryonic particle, it scatters in the centre-of-mass frame about an angle θ_cms and changes its momentum. We are interested in the momentum change parallel to the direction of motion, $$\mathrm{\Delta }{p}_{\parallel }=p_{\mathrm{D}\mathrm{M}}(1-\cos {\theta }_{\mathrm{c}\mathrm{m}\mathrm{s}}),$$(A.5)

as it gives rise to the drag force. The idea behind this is that the momentum changes perpendicular to the direction of motion average out over many scattering events, but the ones parallel to the direction of motion accumulate and effectively decelerate the particle.

Per time dt, we have dC interactions. $$\mathrm{d}C=\frac{\rho }{r{m}_{\chi }}\frac{\mathrm{d}\sigma }{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}}v\mathrm{d}t\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}.$$(A.6)

The number of interactions depends on the differential crosssection, dσ/dΩ_cms and the relative velocity v. We can express the change of parallel momentum that a physical particle experiences as $$\mathrm{d}{p}_{\parallel }=\frac{\rho }{r{m}_{\chi }}\frac{\mathrm{d}\sigma (v)}{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}}v{p}_{\mathrm{D}\mathrm{M}}\mathrm{d}t\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}.$$(A.7)

Here, we use the momentum of the DM particle in the centre-of-mass frame. It is given by $$p_{\mathrm{D}\mathrm{M}}=m_{\chi }v\frac{r}{1+r}.$$(A.8)

After integrating over the differential cross-section one obtains the momentum transfer cross-section, σ_τ (see Appendix A by Kahlhoefer et al. 2014). This allows for the change in the parallel momentum of a physical DM particle travelling through a baryonic background density to be expressed as $$\frac{\mathrm{d}{p}_{\parallel }}{\mathrm{d}t}={\sigma }_{\mathrm{T}}(v)\frac{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}{v}^2}{1+r}.$$(A.9)

We note that f_bary ρ_bary gives the density of particles that the DM particle could interact with. The drag force acting on the numerical DM particle follows by integrating over the number density of DM particles it represents. $$F_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}=\int \frac{\mathrm{d}{p}_{\parallel }}{\mathrm{d}t}\frac{{\rho }_{\mathrm{D}\mathrm{M}}}{m_{\chi }}\mathrm{d}V.$$(A.10)

Using Eq. (4) and μ = m_virt/m_bary, we can derive the final equation for the drag force. For the implementation we use, $$F_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}=\frac{{\sigma }_{\mathrm{T}}(v)}{m_{\chi }}\frac{f_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{\mu (1+r)}{m}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}}{m}_{\mathrm{D}\mathrm{M}}{v}^2{\mathrm{\Lambda }}_{ij}.$$(A.11)

Appendix B Improving energy conservation

To understand how the coupling between the hydrodynamics scheme and the IDM kicks affect the energy conservation, we performed a few tests. Here, we describe how energy conservation can be improved in the case of SPH. We note, however, that based on what we found in our first tests, this does not apply to MFM in a straightforward way. In terms of energy conservation, our implementation for IDM with MFM works quite well and better compared to SPH as long as no modification of the internal energy is involved even if velocity-dependent terms such as viscosity are in place. However, when the internal energy of the gas particles is modified, SPH performs better and we are able to find a way to further improve energy conservation. Clearly, more work is needed for both MFM and SPH. To give a bit of guidance, we now explain how an improvement in the case of SPH is possible.

For improving the energy conservation, we alter the time integration and split the drift step of the KDK scheme into two half-step drifts and do the IDM kicks with the update of the internal energy in between those two half-step drifts. Previously, the DM-baryon interactions were computed between the two half-step kicks. In addition, we also recompute the hydrodynamic accelerations directly after the DM-baryon interactions and before the second half-step drift. Thus, the hydrodynamic accelerations are computed twice per time step, before the second half-step drift and before the second half-step kick, where they are usually computed only (see e.g. Price et al. 2018; Groth et al. 2023). These modifications allow us to significantly improve energy conservation as we show below.

Fig. B.1 Kinetic energy evolution of DM for the heat conduction problem. We show the results for the test problem with heat flow from the baryons to DM for default (blue) and improved (green) implementations. In addition, the red line displays the exact solution given by Eqs. (19) and (20).

Fig. B.2 Energy conservation for the heat conduction problem with the default and improved implementation. From the same simulations as used for Fig. B.1, we show the energy conservation as a function of time for the default (blue) and improved (green) implementations. The expectation, i.e. perfect energy conservation, is indicated by the red line.

In Fig. B.1, we display the evolution of the kinetic energy of DM for the heat transfer problem where heat flows from the baryons to the DM. This is the same set-up used in Sect. 3.1. We show the results from our two implementations, the default one and the improved one. Both show a very similar agreement with the exact solution. The improvement in energy conservation due to the modified time integration scheme is visible from Fig. B.2. We managed to largely improve compared to the default implementation. We note that for this test we only consider the case where all particles reside on the same fixed time step. An extension to variable time steps as well as improving the energy conservation when using MFM is left for future work.

Appendix C Convergence test

In addition to the test shown in Sect. 3, we also test if the results of the heat conduction problem with heat flowing from the baryons to the DM changes when varying the interaction number, N_idm, and the size of the time step.

Compared to our fiducial set-up, we run a simulation using half the time step and one using a four times larger interaction number. All SPH simulations were conducted using artificial viscosity and artificial heat conduction, in contrast to the MFM simulations. The set-up is exactly the same as in Sect. 3.1 except for the mentioned parameters.

In Fig. C.1, we compare the SPH simulation to the exact solution and find that all runs agree well with the expectation. Moreover, we also compare the conservation of total energy among the simulations in Fig. C.2. Here, we find that increasing N_idm does only very slightly improve energy conservation. In contrast, decreasing the time step to half its value allows fora significant reduction in energy errors. It is visible that the remaining error is less than half the error in the simulation with the larger time step. In addition, we also show the results for the MFM simulations in Figs. C.3 and C.4.

From this, we can conclude that increasing N_idm beyond 384 does not help much to reduce numerical errors. In contrast, reducing the time step turned out to be more helpful.

Fig. C.1 Test of varying the time step and the interaction number for the heat conduction problem for SPH. We show the kinetic energy of the DM particles for the problem with heat flow from the baryons to DM as a function of time. The red line indicates the exact solution given by Eqs. (19) and (20). The simulation results are marked by the dotted lines and the varied parameters are displayed in the legend.

Fig. C.2 Energy conservation as a function of time of the heat conduction problem for SPH. The results for the same simulations as in Fig. C.1 are shown. The expectation for perfect energy conservation is indicated by the red line and the numerical results as displayed by the dotted lines with the corresponding simulation parameters being specified by the legend.

Fig. C.3 Varying the time step and the interaction number for the heat conduction test problem. The kinetic energy of the DM is shown as a function of time. Simulation results are indicated by the dotted lines and the analytic solution is given by the solid red line. We note this is the same as in Fig. C.1 but for MFM.

Fig. C.4 Energy conservation for the heat conduction test problem as a function of time with varying time steps and interaction numbers. The results of the same simulations as used for Fig. C.3 are displayed by the dotted lines and the expectation for perfect energy conservation is indicated by the red line. We note this is the same as in Fig. C.2 but for MFM.

Appendix D Comoving integration test

In order to run cosmological simulations we also implemented the comoving integration for the DM-baryon interactions. To test the implementation we simulate a heat conduction problem similar to the tests in Sect. 3.1, but now in expanding space.

The set-up consists of two homogenous components, DM and baryons, which have the same mass. Together they account for a total mass of M_tot = 22.8465 × 10¹⁰M_⊙ h⁻¹. We generate the ICs fora scale factor of a = 0.5 with the mass residing within a cube with a comoving side length of 1400 kpc h⁻¹, implying a comoving density of 83.26 M_⊙ kpc⁻² h². The ICs contain, N_DM = 10⁵, numerical DM particle and, N_bary = 46656, SPH particles. The velocities of the DM initially follow a Maxwell-Boltzmann distribution, their kinetic energy sums up to E_DM = 2.74969 × 10¹¹ M_⊙ h⁻¹ km² s⁻². The bulk velocity of the SPH particle is set to zero but their internal energy makes up E_bary = 2.74158 × 10¹² M_⊙ h⁻¹ km² s⁻². For the simulation, we switch off gravity but use periodic boundary conditions and an expanding space (h = 0.7, Ω_M,0 = 0.3, Ω_Λ = 0.7). Furthermore, we employ artificial viscosity and artificial heat conduction (Price 2012; Beck et al. 2016). For the DM-baryon interactions, the default neighbour numbers of, N_hydro = 230, N_ngb = 64, and N_idm = 384 are used. We simulate a cross-section of σ_τ/m_χ = 2 × 10⁶ cm² g⁻¹ employing a forward dominated cross-section (Eq. 17) with equal masses of the scattering particles (r = 1).

In Fig. D.1, we show the results for our test set-up. In detail, we display the kinetic energy of the DM as a function of the scale factor, a. Here, the kinetic energy is computed as $$E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{D}\mathrm{M}}=\frac{1}{2}\sum _i{m}_i{\mathbf{w}}^2,$$(D.1)

Fig. D.1 Heat conduction problem for comoving integration. The kinetic energy in terms of the canonical momentum is shown as a function of the scale factor. The simulation results are indicated by the black dashed line and the exact solution (Eq. D.3) is shown in red. The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution.

with the canonical momentum w = av_pec. In contrast to the peculiar velocity, v_pec, the canonical momentum is conserved while the scale factor is increasing. This implies that also the kinetic energy when expressed by Eq. (D.1) is conserved. As a consequence, all of the change in energy shown in Fig. D.1 can be attributed to the DM-baryon interactions and is not affected by the expansion.

To compare the simulation results with the exact solution, we compute the solution building up on the work by Dvorkin et al. (2014). In the non-expanding case the power P_DM, arising from the heat conduction is given by $$\begin{gathered} P_{\mathrm{D}\mathrm{M}}=\frac{\mathrm{d}{E}_{\mathrm{D}\mathrm{M}}}{\mathrm{d}t}=-\frac{32}{\sqrt{27\pi }}\frac{{\rho }_{\mathrm{D}\mathrm{M}}{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{(1+r{)}^2}\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }}{(\frac{E_{\mathrm{D}\mathrm{M}}}{M_{\mathrm{D}\mathrm{M}}}+\frac{E_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}{M_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}})}^{1/2} \\ \times [E_{\mathrm{D}\mathrm{M}}(\frac{1}{{\rho }_{\mathrm{D}\mathrm{M}}}+\frac{r}{{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}})-\frac{r{E}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{{\rho }_{\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}}}]. \end{gathered}$$(D.2)

For Sect. 3.1, we computed the full solution analytically (Eqs. 19 and 20). For the comoving integration test the differential equation becomes, $$\frac{\mathrm{d}{E}_{\mathrm{D}\mathrm{M}}}{\mathrm{d}a}=\frac{P_{\mathrm{D}\mathrm{M}}(a)}{aH(a)}-\frac{2E_{\mathrm{D}\mathrm{M}}}{a}.$$(D.3)

Here we express the energy change as a derivative of the scale factor and we use the Hubble parameter H(a). The first term on the right-hand side of Eq. (D.3) corresponds to the energy change caused by heat conduction, whereas the second one corresponds to the expansion of space. Moreover, we note that the quantities in Eq. (D.2) and (D.3) are in terms of their peculiar values, i.e. ρ ∝ a⁻³ and E ∝ a⁻²¹². To obtain the exact solution, we solve the differential equation given by Eq. (D.3) numerically.

In Fig. D.1, we can see that the simulation result (black) agrees well with the exact solution (red). Hence, we can conclude that our implementation of DM-baryon interactions works for comoving integration.

References

All Figures

Fig. 1 Illustration of the numerical scheme for the DM-baryon interactions. The numerical particles are shown together with the physical particles they represent. The velocities of the physical particles are indicated by small arrows. The different stages in treating the interaction between baryons and DM for a single pair of numerical particles are illustrated from the left to the right. As shown here, the case of an interaction between a numerical baryonic particle and a numerical DM particle involves a change in their bulk motion and heating of the baryons. Finally, we note that here we also illustrate the case where the baryons consist of two species, but only one interacts with DM.

In the text

Fig. 2 Heat conduction problem where energy flows from the dark matter to the baryons. Different types and components of the energy are shown as a function of time. The total energy of the system is illustrated in black. For DM the energy is shown in violet, it corresponds to the kinetic energy of the particles as other contributions are zero. In orange, we illustrate the energy of the baryons. It consists of the internal or thermal energy and the kinetic energy (bulk motion). The kinetic part is also displayed separately (dashed line). The simulation was conducted employing a forward dominated cross-section and SPH to describe the gas (baryons).

In the text

	Fig. 3 Same as in Fig. 2, but for the heat conduction problem with energy flowing from the baryons to the DM.
In the text

Fig. 4 Kinetic energy of the DM as a function of time for both heat conduction set-ups. For heat flowing from the DM to baryons (DM→Gas) the results are given in blue and for heat flowing from baryons to DM (Gas→DM) in green. We display the results for both methods SPH (dotted) and MFM (dashed). The red lines give the exact solution according to Eqs. (19) and (20). The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution.

In the text

	Fig. 5 Kinetic energy of the baryonic particles in units of the total energy is shown as a function of time. Labels, colours, and line types are the same as in Fig. 4. We note that this figure is basically a zoom-in on the orange dashed curves displayed in Figs. 2 and 3.
In the text

	Fig. 6 Total energy in units of the initial energy as a function of time to study energy conservation. The results are displayed for the two heat conduction set-ups and both fluid methods, SPH, and MFM. Labels, colours and line types are chosen as in Fig. 4.
In the text

Fig. 7 Relative velocity between DM and baryons for the momentum transfer problem. The time evolution for our SPH and MFM simulations (black) as well as the analytic estimate (red) of the relative velocity are displayed. The latter is computed according to Eq. (24). The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution. We note, that the SPH and MFM results are very similar.

In the text

	Fig. 8 Kinetic energy of the DM particle as a function of time for the momentum transfer problem. The dotted lines give the results using SPH and the dashed ones are for MFM.
In the text

	Fig. 9 Energy conservation for the momentum transfer problem as a function of time. The total energy divided by the initial energy is given for SPH (dotted line) and MFM (dashed line).
In the text

	Fig. 10 Same as in Fig. 2, but for a heat conduction problem with an unequal mass ratio evolved using SPH.
In the text

	Fig. 11 Kinetic energy of the DM particles as a function of time. The simulated (black) time evolution for the heat conduction problem with a mass ratio of r = 1000 is compared to the exact solution (red) from Eqs. (19) and (20). The upper panel gives the absolute values and the lower panel displays the ratios to the exact solution.
In the text

	Fig. 12 Energy evolution for isotropic scattering with equal particle masses. The different energy components for DM and baryons are shown as a function of time. This is the same test problem as shown in Fig. 2 for small-angle scattering.
In the text

	Fig. 13 Time evolution of the kinetic energy of the DM. The simulation results (black) for the heat conduction problem with isotropic scattering are compared to the exact solution (red) from Eqs. (19) and (20). The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution.
In the text

	Fig. 14 Density profile for a collapsing overdensity in the case of collisionless DM. The density is shown as a function of radius at several redshifts for DM (solid) and baryons (dashed). The density is not displayed for small radii where the particle number is too low to obtain a reliable value.
In the text

Fig. 15 Density, velocity dispersion, and entropy profile of DM and baryons for CDM and IDM. We show the density as a function of radius (upper panel) at a redshift of z = 0. The velocity dispersion, including kinetic and internal energy but excluding radial bulk motion, is displayed in the middle panel. The results for collisionless DM are shown in black. In the lower panel, the entropic function (Eq. (31)) is displayed. The solid lines are for the DM and the dashed ones are for the baryons. Furthermore, the coloured curves are for IDM simulations employing different cross-sections. In addition, the results for a simulation that makes use of the scheme for improved energy conservation as explained in Appendix B are given too.

In the text

Fig. 16 Central density, specific energy and DM interaction rate as a function of scale factor. For the simulated DM models we show the evolution of different quantities. The central density (upper panel) and specific energy (middle panel) within a radius of r = 2.7 ckpc are shown for the DM and baryonic components. The lower panel gives the average DM interaction rate within the same radius according to Eq. (32) assuming an isotropic cross-section matched to the simulated forward dominated one with the momentum transfer cross-section (σ|iso = 2σT|iso).

In the text

	Fig. B.1 Kinetic energy evolution of DM for the heat conduction problem. We show the results for the test problem with heat flow from the baryons to DM for default (blue) and improved (green) implementations. In addition, the red line displays the exact solution given by Eqs. (19) and (20).
In the text

	Fig. B.2 Energy conservation for the heat conduction problem with the default and improved implementation. From the same simulations as used for Fig. B.1, we show the energy conservation as a function of time for the default (blue) and improved (green) implementations. The expectation, i.e. perfect energy conservation, is indicated by the red line.
In the text

Fig. C.1 Test of varying the time step and the interaction number for the heat conduction problem for SPH. We show the kinetic energy of the DM particles for the problem with heat flow from the baryons to DM as a function of time. The red line indicates the exact solution given by Eqs. (19) and (20). The simulation results are marked by the dotted lines and the varied parameters are displayed in the legend.

In the text

	Fig. C.2 Energy conservation as a function of time of the heat conduction problem for SPH. The results for the same simulations as in Fig. C.1 are shown. The expectation for perfect energy conservation is indicated by the red line and the numerical results as displayed by the dotted lines with the corresponding simulation parameters being specified by the legend.
In the text

	Fig. C.3 Varying the time step and the interaction number for the heat conduction test problem. The kinetic energy of the DM is shown as a function of time. Simulation results are indicated by the dotted lines and the analytic solution is given by the solid red line. We note this is the same as in Fig. C.1 but for MFM.
In the text

	Fig. C.4 Energy conservation for the heat conduction test problem as a function of time with varying time steps and interaction numbers. The results of the same simulations as used for Fig. C.3 are displayed by the dotted lines and the expectation for perfect energy conservation is indicated by the red line. We note this is the same as in Fig. C.2 but for MFM.
In the text

	Fig. D.1 Heat conduction problem for comoving integration. The kinetic energy in terms of the canonical momentum is shown as a function of the scale factor. The simulation results are indicated by the black dashed line and the exact solution (Eq. D.3) is shown in red. The upper panel gives the absolute values, and the lower panel displays the ratios to the exact solution.
In the text