Issue 
A&A
Volume 666, October 2022



Article Number  A95  
Number of page(s)  15  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202243496  
Published online  20 October 2022 
Cosmological simulations of selfinteracting BoseEinstein condensate dark matter
Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern, 0315 Oslo, Norway
email: stian.hartman@gmail.com
Received:
8
March
2022
Accepted:
16
August
2022
Fully 3D cosmological simulations of scalar field dark matter with selfinteractions, also known as BoseEinstein condensate dark matter, are performed using a set of effective hydrodynamic equations. These are derived from the nonlinear Schrödinger equation by performing a smoothing operation over scales larger than the de Broglie wavelength, but smaller than the selfinteraction Jeans’ length. The dynamics on the de Broglie scale become an effective thermal energy in the hydrodynamic approximation, which is assumed to be subdominant in the initial conditions, but become important as structures collapse and the fluid is shockheated. The halos that form have NavarroFrenkWhite envelopes, while the centers are cored due to the fluid pressures (thermal + selfinteraction), confirming the features found by Dawoodbhoy et al. (2021, MNRAS, 506, 2418) using 1D simulations under the assumption of spherical symmetry. The core radii are largely determined by the selfinteraction Jeans’ length, even though the effective thermal energy eventually dominates over the selfinteraction energy everywhere, a result that is insensitive to the initial ratio of thermal energy to interaction energy, provided it is sufficiently small to not affect the linear and weakly nonlinear regimes. Scaling relations for the simulated population of halos are compared to Milky Way dwarf spheroidals and nearby galaxies, assuming a Burkert halo profile, and are found to not match, although they conform better with observations compared to fuzzy dark matteronly simulations.
Key words: cosmology: theory / dark matter
© S. T. H. Hartman et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Unveiling the fundamental nature of dark matter (DM) – the dominant matter component of our Universe – has been one of the ultimate goals of physics for many decades. Although its particle identity has remained elusive, the collisionless and cold DM (CDM) paradigm, along with a cosmological constant and inflationary initial conditions, has proven extremely successful at explaining a wide range of observables, and is known as the ΛCDM model (Davis et al. 1985; Percival et al. 2001; Tegmark et al. 2004; TrujilloGomez et al. 2011; Vogelsberger et al. 2014; Planck Collaboration XIII 2016; Riess et al. 2016; Cyburt et al. 2016). Nevertheless, the quest for a complete theory of DM has spawned numerous expanded models, many of which aim to reproduce the success of CDM while attempting to fix discrepancies between ΛCDM and observations, finding physics beyond the standard model of particle physics (SM), or both.
Some of the tensions in ΛCDM are related to the formation and shape of smallscale structures. CDM, due to its cold and collisionless nature, is very efficient at forming structure at all scales. For instance, Nbody simulations predict CDM halos to follow the universal NavarroFrenkWhite (NFW) profile (Navarro et al. 1996, 1997), which diverges near the center as ρ(r)∼r^{−1}, whereas measurements instead indicate the central regions of DM profiles of lowmass halos to be flatter, or rather, that lowmass halos are cored rather than cuspy. Furthermore, ΛCDM predicts a large number of lowmass halos, as well as massive subhalos that should be too big to fail at forming stars, yet these are largely absent in our galactic neighborhood. These issues of ΛCDM, known as the “cuspcore”, “missinghalo”, and “toobigtofail” problems, respectively, might be due to limitations in accurately modeling baryonic physics and processes that are too small to be resolved in largescale simulations (see Weinberg et al. 2015; Del Popolo & Le Delliou 2017; Bullock & BoylanKolchin 2017, and references therein for further discussion on these smallscale discrepancies in ΛCDM). Alternatively, the solution might lie in the dark sector.
A popular class of DM models beyond CDM involve ultralight scalar, or pseudoscalar, particles that have a sufficiently small mass to exhibit wavelike behaviour on astrophysical scales. The freefield case, termed Fuzzy DM (FDM; Dine & Fischler 1983; Preskill et al. 1983; Hu et al. 2000; Marsh 2016; Hui et al. 2017), with masses as low as 10^{−22} eV, have large de Broglie wavelengths λ_{dB} = 1/mv and produce solitonic cores of order 1 kpc at the centers of DM halos. Unlike CDM, which clusters on all scales, FDM suppresses smallscale structure because of an effective Jeans’ length due to the large de Broglie wavelength. Furthermore, FDM produces interference patterns in its density distribution, which is a very distinct feature of FDM compared to CDM and other DM models, such as warm DM. However, a recent bound on the FDM mass from the Lymanalpha forest constrains it to m > 2 × 10^{−20} eV, and therefore seems to rule out the canonical massrange generally considered to be needed to solve the smallscale problems of ΛCDM (Rogers & Peiris 2021).
Another kind of ultralight DM are scalar fields with interactions. These can vary in complexity, from singlefield DM with selfinteractions (Lee & Koh 1996; Peebles 2000; Goodman 2000; Arbey et al. 2002, 2003; Böhmer & Harko 2007; Chavanis 2011; RindlerDaller & Shapiro 2012), to multifield DM with nontrivial couplings that give rise to exotic properties (Matos et al. 2000; Bettoni et al. 2014; Berezhiani & Khoury 2015; Khoury 2016; Ferreira et al. 2019), though a common feature among these models is for the interactions to give rise to a fluid pressure that produces halo cores. The simplest scenario is singlefield DM with quartic selfinteractions and a negligible de Broglie wavelength, which, along with FDM, is often referred to as BoseEinstein condensed DM, because they can be regarded as the zerotemperature limit of a boson gas, which in meanfield theory is described by a classical field (Pitaevskii & Stringari 2016). The term superfluid DM is also frequently used due to the close relationship between BoseEinstein condensates and superfluidity, though in this work we will instead use “selfinteracting BoseEinstein condensed” (SIBEC) DM, in an effort to emphasize the importance of the selfinteractions on the dynamics of the scalar field.
The nonlinear structure of SIBECDM has largely been investigated near hydrostatic equilibrium, which finds SIBECDM halo cores to be independent of the total halo mass and core density. Fitting SIBECDM to nearby galaxies reproduces the observed rotation curves for core radii around 1 kpc and larger (Zhang et al. 2018; Crăciun & Harko 2020). However, to obtain more realistic SIBECDM halo profiles from cosmological structure formation, which are expected transition to into the NFW profile outside the cores, one needs to go beyond hydrostatic considerations and use numerical simulations, though some challenges arise in such an endeavor. For instance, the equation of motion for nonrelativistic scalar fields, the nonlinear Schrödinger equation (NLSE), is very computationally demanding to solve, even in the FDMlimit where the de Broglie wavelength is of astrophysical scale. For SIBECDM, where λ_{dB} is much smaller, but still needs to be resolved, the NLSE is even more computationally demanding. Simply setting the terms that describe the dynamics at λ_{dB}scales to zero, which works in the hydrostatic case, does not work in nonlinear simulations, as they play an important role in regularizing discontinuities and keeping the numerical solution wellbehaved. A scheme that somehow incorporates the de Broglie scale without actually needing to resolve it is therefore desired, and was recently employed by Dawoodbhoy et al. (2021) and Shapiro et al. (2021). They used a hydrodynamic formulation of the NLSE that incorporates the dynamics on the de Broglie scale as an effective thermal pressure, and performed 1D simulations of collapsing spherical overdensities of SIBECDM. They found that in a noncosmological setting, SIBECDM halos with cores of order 1 kpc provide a solution to the smallscale issues of the standard model. In fact, they found SIBECDM to provide a better solution to both the toobigtofail and cuspcore problems than FDM, which struggles to fix these two issues at the same time (Robles et al. 2019). Moving to a cosmological setting, on the other hand, results in an overly large suppression of the formation of smallscale halos for these same core radii. The cores therefore have to be smaller than what is needed to provide a solution to, for instance, the corecusp issue, in order to be compatible with the halo mass function. A similar conclusion was drawn by Hartman et al. (2022) using largescale observables, albeit to a much lesser degree. Exploring the formation of nonlinear structure in a universe with SIBECDM and investigating the properties of SIBECDM halos is nevertheless a worthwhile endeavor in the effort to further understand the rich phenomenology that is found in scalar field DM models. In this paper we therefore build on the work of Dawoodbhoy et al. (2021) and Shapiro et al. (2021) by implementing the de Brogliesmoothed hydrodynamic formulation of SIBECDM into a modified version of the RAMSES code (Teyssier 2002). This code is then used to perform fully 3D cosmological simulations to explore the formation of SIBECDM structures and their nonlinear dynamics. The paper is structured as follows: In Sect. 2 the basic theoretical framework for modelling scalar field DM is reviewed, along with the smoothing procedure used to derive the hydrodynamic approximation for SIBECDM. In Sect. 3 the numerical aspects of this paper are discussed, including a brief description of how the RAMSES code works, the initial setup of the SIBECDM fluid, and a discussion of the limitations of the simulations. Our main results are presented and discussed in Sect. 4, with final conclusions in Sect. 5. We use, for the most part, natural units.
2. Theory of selfinteracting BoseEinstein condensates dark matter
Ultralight scalar DM can be described by a classical field theory, and one of the simplest realizations is the Lagrangian for a complex scalar field Ψ with selfinteractions (Li et al. 2014),
$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{2m}{g}_{\mu \nu}{\partial}^{\mu}{\mathrm{\Psi}}^{\ast}{\partial}^{\nu}\mathrm{\Psi}\frac{1}{2}{m\mathrm{\Psi}}^{2}\frac{1}{2}g{\mathrm{\Psi}}^{4},\end{array}$$(1)
where m is the mass of the scalar field, and g is the coupling parameter for the selfinteraction. In the early universe, where DM is much denser and the selfinteraction of the scalar field dominates the energy density, the averaged pressure and energy of the field behave as radiation (Matos & Arturo UreñaLópez 2001; Li et al. 2014). As the universe expands and becomes more diluted, the interaction energy eventually becomes smaller than the rest energy of the field, and it ceases to be radiationlike. In this nonrelativistic phase we can define a new field, Ψ = ψe^{−imt}, whose dynamics are given by the NLSE (Ferreira 2021),
$$\begin{array}{c}\hfill i\frac{\partial \psi}{\partial t}=[\frac{{\mathrm{\nabla}}^{2}}{2m}+V]\psi ,\end{array}$$(2)
where V = gψ^{2} + mϕ includes twobody contact interactions and the gravitational potential ϕ. The NLSE can be reformulated in a hydrodynamical form by substituting for the wavefunction
$$\begin{array}{c}\hfill \psi =\sqrt{n}{e}^{\mathit{iS}}=\sqrt{\frac{\rho}{m}}{e}^{\mathit{iS}},\end{array}$$(3)
where n is the particle number density, ρ = mn the mass density, and the velocity field is defined by v = ∇S/m. This gives the Madelung equations (Madelung 1926),
$$\begin{array}{cc}& \frac{\partial \rho}{\partial t}+\mathbf{\nabla}\xb7(\rho \mathit{v})=0,\hfill \end{array}$$(4)
$$\begin{array}{cc}& \frac{\partial \mathit{v}}{\partial t}+(\mathit{v}\xb7\mathbf{\nabla})\mathit{v}+\mathbf{\nabla}(\frac{g\rho}{{m}^{2}}\frac{1}{2{m}^{2}}\frac{{\mathrm{\nabla}}^{2}\sqrt{\rho}}{\sqrt{\rho}}+\varphi )=\mathbf{0},\hfill \end{array}$$(5)
which are readily recognizable as an equation for mass conservation and a variant of the momentum equation. The largest difference from standard hydrodynamics is the absence of an energy equation, the presence of a selfinteraction pressure
$$\begin{array}{c}\hfill {P}_{\mathrm{SI}}=\frac{g{\rho}^{2}}{2{m}^{2}},\end{array}$$(6)
and a socalled quantum potential
$$\begin{array}{c}\hfill Q=\frac{1}{2{m}^{2}}\frac{{\mathrm{\nabla}}^{2}\sqrt{\rho}}{\sqrt{\rho}}\xb7\end{array}$$(7)
The phenomenology of ultralight DM depends crucially on whether the interaction pressure or the quantum potential dominates. This becomes apparent when considering halo properties of ultralight DM at equilibrium. In the FDMlimit, the de Broglie wavelength of the ultralight DM particles is much larger than the Jeans’ length of the selfinteraction and is of astrophysical size, such that smallscale structure is stabilized against collapse due to the formation of solitonic cores for masses m ≲ 10^{−21}. The core radii of FDM halos at hydrostatic equilibrium are approximately (Membrado et al. 1989)
$$\begin{array}{c}\hfill {R}_{\mathrm{FDM}}=\frac{10}{GM{m}^{2}},\end{array}$$(8)
where M is the total mass of the halo. The opposite limit Q ≪ P_{SI}/ρ, which is the case we refer to as SIBECDM, but is also known as the ThomasFermi approximation, results in the hydrostatic density profile (Goodman 2000)
$$\begin{array}{c}\hfill \rho (r)={\rho}_{0}\frac{sin(Ar)}{\mathit{Ar}},\end{array}$$(9)
where $A=\sqrt{4\pi G{m}^{2}/g}$, and goes to zero at
$$\begin{array}{c}\hfill r={R}_{\mathrm{c}}=\sqrt{\frac{g\pi}{4G{m}^{2}}}\xb7\end{array}$$(10)
The soliton cores of massive FDM halos are smaller than lowmass halos, a feature confirmed by simulations (Chan et al. 2022), whereas the core radius R_{c} of halos in the ThomasFermi limit is independent of the total halo mass and central density, depending only on the combination g/m^{2}.
Both the NLSE and its hydrodynamic Madelung formulation are widely used to study ultralight DM across a wide range of scales, from the properties of lowmass halos and galaxies (Membrado et al. 1989; Goodman 2000; Böhmer & Harko 2007; Chavanis 2011; Chavanis & Delfini 2011; RindlerDaller & Shapiro 2014; Hui et al. 2017; Zhang et al. 2018; Berezhiani et al. 2019; Crăciun & Harko 2020; Lancaster et al. 2020) to largescale simulations (Schive et al. 2014a,b; Schwabe et al. 2016; Mocz et al. 2017, 2018; Veltmaat et al. 2018; Nori & Baldi 2018; Mina et al. 2020, 2022; Nori & Baldi 2021; May & Springel 2021). However, only the FDMlimit has been thoroughly investigated using fully 3D cosmological simulations. No corresponding work has so far been carried out in the ThomasFermi limit, in part due to a number of challenges in simulating this kind of system. The hydrodynamic equations obtained by neglecting the quantum potential appear equivalent to regular hydrodynamics with an adiabatic index γ = 2, such that it seems standard numerical schemes should work, but as shock fronts and discontinuities arise, the ThomasFermi approximation breaks down and the full equations are needed. However, both the NLSE and Madelung equations are computationally demanding to solve due to the high resolution needed to resolve the scales on which the quantum potential operates. This is already a challenge in FDM simulations, but even more so for SIBECDM where the characteristic scale of Q is much smaller than the scales of interest given by the interaction pressure, λ_{dB} ≪ R_{c}.
2.1. Smoothed SIBEC hydrodynamics
An alternative hydrodynamic approximation of SIBECDM that includes the effect of the quantum potential without actually needing to resolve λ_{dB} was recently used by Dawoodbhoy et al. (2021) and Shapiro et al. (2021). It involves using a smoothed phase space representation of the wavefunction, known as the Husimi representation (Husimi 1940), which is obtained by essentially smoothing over ψ with a Gaussian window function of width η and Fourier transforming;
$$\begin{array}{c}\hfill \mathrm{\Psi}(\mathit{x},\mathit{p},t)={\displaystyle \int}\frac{{\mathrm{d}}^{3}y}{{(2\pi )}^{3/2}}\phantom{\rule{0.166667em}{0ex}}\frac{{e}^{{(\mathit{x}\mathit{y})}^{2}/2{\eta}^{2}}}{{(\eta \sqrt{\pi})}^{3/2}}\psi (\mathit{y},t){e}^{i\mathit{p}\xb7(\mathit{y}\mathit{x}/2)}.\end{array}$$(11)
Defining the phase space distribution function
$$\begin{array}{c}\hfill \mathcal{F}(\mathit{x},\mathit{p},t)={\mathrm{\Psi}(\mathit{x},\mathit{p},t)}^{2},\end{array}$$(12)
and computing the derivative ∂ℱ/∂t, using the NLSE i∂ψ/∂t = [ − ∇^{2}/2m + V]ψ, and smoothing over scales larger than the de Broglie wavelength λ_{dB} ≪ η, gives an approximate equation of motion that is of the same form as the collisionless Boltzmann equation (Skodje et al. 1989; Widrow & Kaiser 1993);
$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathcal{F}}{\mathrm{d}t}=\frac{\partial \mathcal{F}}{\partial t}+\frac{{p}_{i}}{m}\frac{\partial \mathcal{F}}{\partial {x}_{i}}\frac{\partial V}{\partial {x}_{i}}\frac{\partial \mathcal{F}}{\partial {p}_{i}}=0.\end{array}$$(13)
The gravitational potential and the meanfield potential due to selfinteractions in V are now given by the smoothed density field, which is obtained from ℱ in the same way as from a standard phase space distribution function,
$$\begin{array}{c}\hfill \rho (\mathit{x},t)=m{\displaystyle \int}{\mathrm{d}}^{3}p\phantom{\rule{0.166667em}{0ex}}\mathcal{F}(\mathit{x},\mathit{p},t)={\displaystyle \int}{\mathrm{d}}^{3}y\phantom{\rule{0.166667em}{0ex}}\frac{{e}^{{(\mathit{x}\mathit{y})}^{2}/{\eta}^{2}}}{{(\eta \sqrt{\pi})}^{3}}m{\psi (\mathit{y},t)}^{2}.\end{array}$$(14)
Following Dawoodbhoy et al. (2021), the hydrodynamic equations for SIBECDM in this approximation are obtained from the zeroth, first, and second moments of Eq. (13), and it is assumed that the phase space distribution function is isotropic around its mean momentum ⟨p⟩, that is, that ℱ(x, p, t) is skewless in p. This ansatz is supported by the results of Veltmaat et al. (2018), who performed 3D simulations of FDM halo mergers using the NLSE and found the velocity distribution of plane waves inside the virial radius to be Maxwellian. Before the growth of nonlinear structures, and outside of halos, such as in regions that have not undergone violent relaxation and shell crossing, the assumption of skewlessness in ℱ may not hold true. On the other hand, the fluid is supersonic and ballistic in these regions, as well as dominated by the interaction pressure, so it should not matter in these situations. During mergers, however, the skewless assumption might also break down and have significant implications if the merging halos have large planewave velocity dispersions. We nevertheless use this assumption throughout our simulations. The resulting equations are
$$\begin{array}{cc}& \frac{\partial \rho}{\partial t}+\mathbf{\nabla}\xb7\mathit{j}=0,\hfill \end{array}$$(15)
$$\begin{array}{cc}& \frac{\partial \mathit{j}}{\partial t}+\mathbf{\nabla}(P+{P}_{\mathrm{SI}})+\rho (\mathit{v}\xb7\mathbf{\nabla})\mathit{v}+\mathit{v}(\mathbf{\nabla}\xb7\rho \mathit{v})=\rho \mathbf{\nabla}\varphi ,\hfill \end{array}$$(16)
$$\begin{array}{cc}& \frac{\partial E}{\partial t}+\mathbf{\nabla}\xb7[(E+P+{P}_{\mathrm{SI}})\mathit{v}]=\mathit{j}\xb7\mathbf{\nabla}\varphi .\hfill \end{array}$$(17)
There are two major changes in the above hydrodynamic description of BECs compared to the Madelung equations; first, there is no quantum potential; and second, there is an energy equation in addition to the continuity and momentum equations, with the total energy given by
$$\begin{array}{c}\hfill E=\frac{1}{2}\rho {v}^{2}+U+{U}_{\mathrm{SI}},\end{array}$$(18)
where $P=\frac{2}{3}U$ and P_{SI} = U_{SI}. In fact, the smallscale dynamics of the Madelung equations driven by the quantum potential, after smoothing over scales larger than the de Broglie wavelength, behave effectively as an ideal gas with adiabatic index γ = 5/3, with an additional pressure and internal energy due to the selfinteractions. The fluid variables, such as the fluid density and momentum, and the effective thermal internal energy and pressure, are defined in the same way as in classical statistical mechanics (Huang 1987), but using the smoothed phase space representation of the wavefunction rather than a phase space distribution function of classical point particles;
$$\begin{array}{cc}& n(\mathit{x},t)={\displaystyle \int}{\mathrm{d}}^{3}p\phantom{\rule{0.166667em}{0ex}}\mathcal{F}(\mathit{x},\mathit{p},t),\hfill \end{array}$$(19)
$$\begin{array}{cc}& \mathit{j}(\mathit{x},t)=n\mathit{p}={\displaystyle \int}{\mathrm{d}}^{3}p\phantom{\rule{0.166667em}{0ex}}\mathit{p}\phantom{\rule{0.166667em}{0ex}}\mathcal{F}(\mathit{x},\mathit{p},t),\hfill \end{array}$$(20)
$$\begin{array}{cc}& P(\mathit{x},t)={\displaystyle \int}{\mathrm{d}}^{3}p\phantom{\rule{0.166667em}{0ex}}\frac{{\mathit{p}\mathit{p}}^{2}}{3m}\mathcal{F}(\mathit{x},\mathit{p},t),\hfill \end{array}$$(21)
$$\begin{array}{cc}& U(\mathit{x},t)={\displaystyle \int}{\mathrm{d}}^{3}p\phantom{\rule{0.166667em}{0ex}}\frac{{\mathit{p}\mathit{p}}^{2}}{2m}\mathcal{F}(\mathit{x},\mathit{p},t),\hfill \end{array}$$(22)
$$\begin{array}{cc}& E(\mathit{x},t)={U}_{\mathrm{SI}}+{\displaystyle \int}{\mathrm{d}}^{3}p\phantom{\rule{0.166667em}{0ex}}\frac{{p}^{2}}{2m}\mathcal{F}(\mathit{x},\mathit{p},t).\hfill \end{array}$$(23)
In the following, we will sometimes refer to the “effective thermal” energy and pressure as simply “thermal”.
2.2. Supercomoving hydrodynamics
Now that we have a set of fluid equations for SIBECDM, we would like to apply these to the growth of structure in a cosmological setting. A convenient set of coordinates for this purpose are the socalled supercomoving coordinates (Martel & Shapiro 1998), which uses the following change of variables in the case of a flat universe dominated by matter and a cosmological constant;
$$\begin{array}{cc}\hfill \stackrel{\sim}{\mathit{x}}& =\frac{1}{a}\frac{\mathit{x}}{L},\phantom{\rule{1em}{0ex}}\mathrm{d}\stackrel{\sim}{t}={H}_{0}\frac{\mathrm{d}t}{{a}^{2}},\phantom{\rule{1em}{0ex}}\stackrel{\sim}{\rho}={a}^{3}\frac{\rho}{{\mathrm{\Omega}}_{\mathrm{m}0}{\rho}_{\mathrm{c}0}},\hfill \\ \hfill \stackrel{\sim}{\mathit{u}}& =a\frac{\mathit{u}}{{H}_{0}L},\phantom{\rule{1em}{0ex}}\stackrel{\sim}{P}={a}^{5}\frac{P}{{\mathrm{\Omega}}_{\mathrm{m}0}{\rho}_{\mathrm{c}0}{H}_{0}^{2}{L}^{2}},\hfill \\ \hfill \stackrel{\sim}{m}& =m{H}_{0}L,\phantom{\rule{1em}{0ex}}\stackrel{\sim}{g}=g\frac{{\mathrm{\Omega}}_{\mathrm{m}0}{\rho}_{\mathrm{c}0}}{a},\hfill \end{array}$$(24)
where the peculiar motion u = v − Hr is introduced. L is a free length scale parameter, Ω_{m0} is the fraction of matter today, and ρ_{c0} the critical energy density today. In these coordinates, the fluid equations in an expanding universe are on a similar form as the standard equations on a static background;
$$\begin{array}{cc}& \frac{\partial \stackrel{\sim}{\rho}}{\partial \stackrel{\sim}{t}}+\stackrel{\sim}{\mathbf{\nabla}}\xb7\stackrel{\sim}{\mathit{j}}=0,\hfill \end{array}$$(25)
$$\begin{array}{cc}& \frac{\partial \stackrel{\sim}{\mathit{j}}}{\partial \stackrel{\sim}{t}}+\stackrel{\sim}{\mathbf{\nabla}}({\stackrel{\sim}{P}}_{\mathrm{tot}})+\stackrel{\sim}{\rho}(\stackrel{\sim}{\mathit{u}}\xb7\stackrel{\sim}{\mathbf{\nabla}})\stackrel{\sim}{\mathit{u}}+\stackrel{\sim}{\mathit{u}}(\stackrel{\sim}{\mathbf{\nabla}}\xb7\stackrel{\sim}{\rho}\stackrel{\sim}{\mathit{u}})=\stackrel{\sim}{\rho}\stackrel{\sim}{\mathbf{\nabla}}\stackrel{\sim}{\varphi},\hfill \end{array}$$(26)
$$\begin{array}{cc}& \frac{\partial \stackrel{\sim}{E}}{\partial \stackrel{\sim}{t}}+\stackrel{\sim}{\mathbf{\nabla}}\xb7[(\stackrel{\sim}{E}+{\stackrel{\sim}{P}}_{\mathrm{tot}})\stackrel{\sim}{\mathit{u}}]+\mathcal{H}(3{\stackrel{\sim}{P}}_{\mathrm{tot}}2{\stackrel{\sim}{U}}_{\mathrm{tot}})=\stackrel{\sim}{\mathit{j}}\xb7\stackrel{\sim}{\mathbf{\nabla}}\stackrel{\sim}{\varphi}.\hfill \end{array}$$(27)
where ${\stackrel{\sim}{P}}_{\mathrm{tot}}=\stackrel{\sim}{P}+{\stackrel{\sim}{P}}_{\mathrm{SI}}$ and ${\stackrel{\sim}{U}}_{\mathrm{tot}}=\stackrel{\sim}{U}+{\stackrel{\sim}{U}}_{\mathrm{SI}}$. The supercomoving Hubble parameter is given by $\mathcal{H}={a}^{1}\mathrm{d}a/\mathrm{d}\stackrel{\sim}{t}$, the gravitational potential by
$$\begin{array}{c}\hfill {\stackrel{\sim}{\mathrm{\nabla}}}^{2}\stackrel{\sim}{\varphi}=\frac{3}{2}a{\mathrm{\Omega}}_{\mathrm{m}}(\stackrel{\sim}{\rho}1),\end{array}$$(28)
and expressions for the energies, pressures, and so on, are the same as before, but with supercomoving quantities. A notable difference from regular hydrodynamics is the presence of the Hubble drag in the supercomoving energy equation, which vanishes for the ideal gas pressure with γ = 5/3, but not for the selfinteraction pressure.
3. Numerical implementation
A modified version of RAMSES (Teyssier 2002) is used to simulate the formation of structure of SIBECDM in a cosmological setting. The halo catalogs are obtained with Amiga’s Halo Finder (AHF; Gill et al. 2004; Knollmann & Knebe 2009), and plots of simulation snapshots are made using the python package YT (Turk et al. 2010). RAMSES is a grid based code that uses adaptive mesh refinement (AMR) to focus computational resources through a hierarchy of nested grids, increasing the local resolution in regions of interest. This is done by recursively splitting individual cells into 2^{N} smaller cells, where N is the dimensionality, until a set of refinement criteria are satisfied, doubling the spatial resolution at each new level of refinement. For a box of length L, the cell length at refinement level l is Δx^{l} = L/2^{l}.
RAMSES was originally developed for cosmological simulations of the formation and evolution of structure at both large and small scales, using an Nbody particlemesh code for CDM, with which the phase space distribution of DM is sampled through a computationally tractable number of collisionless macroparticles. The gravitational forces acting on these macroparticles are obtained by projecting the DM mass onto the AMR grid and solving the Poisson equation. RAMSES also includes code for studying gas dynamics, using the grid and a Godunov scheme to solve the hydrodynamic equations in supercomoving coordinates. It is this latter part of RAMSES that has been modified in order to solve the SIBECDM hydrodynamic equations, and we use a secondorder Godunov scheme with the Rusanov flux (an HLLE Riemann solver) and the van Leer slope limiter to compute the flux at the cell interfaces (Toro 2006). The Rusanov flux is computed using the conservative variables $\stackrel{\sim}{\rho}$, $\stackrel{\sim}{\mathit{j}}$, and $\stackrel{\sim}{E}$, and Eqs. (25)–(27), while the secondorder MUSCLHancock and slope limiting step is done using the primitive variables $\stackrel{\sim}{\rho}$, $\stackrel{\sim}{\mathit{u}}$, and $\stackrel{\sim}{P}$ (the thermal pressure).
A 3D test of the collapse of a spherically symmetric overdensity at two different grid resolutions (refinement levels 7 to 11, and 8 to 12) are included in Fig. 1. A periodic box of size L = 50 kpc was used, with the initial density ρ(r) = ρ_{0}[1 + e^{−(r/R)2}], where R = 5 kpc and ρ_{0} = 0.3ρ_{c0}. The initial thermal pressure was ζ = P/P_{SI} = 10^{−1}, and the simulations were run until 3 t_{dyn}, where ${t}_{\mathrm{dyn}}=\sqrt{3/8\pi G{\rho}_{0}}$ is the dynamical timescale of the initial overdensity. The energy density profiles, with
$$\begin{array}{c}\hfill W=\frac{GM(<r)\rho}{r}\xb7\end{array}$$(29)
Fig. 1. Selfinteraction energy density U_{SI}, thermal energy density U, and gravitational potential energy density W of the spherically symmetric test halo with R_{c} = 1 kpc, relative to the interaction energy of the initial background ${\overline{U}}_{\mathrm{SI}}$. The result of the simulation with refinement levels 8 to 12 are shown in solid, while 7 to 11 is shown in dotted. The 3D simulation reproduces the 1D results of Ahn & Shapiro (2005) and Dawoodbhoy et al. (2021), in particular the large thermal energy in the wake of the accretion shock, its sharp decrease to below the selfinteraction energy inside the core, and the fit U_{SI} ∝ ρ^{2} ∝ r^{−24/7} in the envelope. 
agrees with the results of the 1D simulations of Ahn & Shapiro (2005) and Dawoodbhoy et al. (2021), in particular the large thermal energy in the wake of the accretion shock, its sharp decrease to below the selfinteraction energy inside the core, and the fit ρ ∝ r^{−12/7} in the envelope, which gives U_{SI} ∝ ρ^{2} ∝ r^{−24/7}.
A number of simulations were run, with different SIBECDM interaction strengths (given by the core radius R_{c}), initial conditions, and box sizes to balance reaching sufficient spatial resolution to resolve the cores of the emerging SIBECDM halos, while also forming enough such halos to study their properties. Initial conditions were generated with the MUSIC code (Hahn & Abel 2011) using the standard CDM matter power spectrum, but with a cutoff in power below a certain scale to include the effect of a nonzero sound speed,
$$\begin{array}{c}\hfill {P}_{\mathrm{SIBEC}}(k)={f}_{\mathrm{cut}}(k){P}_{\mathrm{CDM}}(k).\end{array}$$(30)
Comparison with output from the Boltzmann code CLASS, modified to include SIBECDM (Hartman et al. 2022), shows that the function
$$\begin{array}{c}\hfill {f}_{\mathrm{cut}}(k)={e}^{{(k/{k}_{\mathrm{cut}})}^{3}},\end{array}$$(31)
with k_{cut} given by the effective sound horizon r_{s} of SIBECDM,
$$\begin{array}{c}\hfill {r}_{\mathrm{s}}={\displaystyle {\int}_{0}^{{\tau}_{\mathrm{init}}}\mathrm{d}\tau \phantom{\rule{0.166667em}{0ex}}{c}_{\mathrm{s}}={\int}_{{z}_{\mathrm{init}}}^{\infty}\mathrm{d}z\frac{{c}_{\mathrm{s}}}{H},}\end{array}$$(32)
provides a good fit to the suppression of power in a universe with SIBECDM relative to CDM. The SIBECDM sound speed c_{s}, which includes the early radiationlike epoch during which the interaction energy dominates, is approximately given by (Hartman et al. 2022)
$$\begin{array}{c}\hfill {c}_{\mathrm{s}}^{2}=\frac{\partial \overline{P}}{\partial \overline{\rho}}=\frac{1}{3}\frac{1}{1+{a}^{3}/6{\omega}_{0}},\end{array}$$(33)
where ω_{0} is the SIBECDM equation of state today,
$$\begin{array}{c}\hfill {\omega}_{0}=\frac{{\overline{P}}_{0}}{{\overline{\rho}}_{0}}=\frac{g{\overline{\rho}}_{0}}{2{m}^{2}}\approx {10}^{15}{\left(\frac{{R}_{\mathrm{c}}}{1\phantom{\rule{0.166667em}{0ex}}\mathrm{kpc}}\right)}^{2}.\end{array}$$(34)
The resulting cutoff scale is
$$\begin{array}{c}\hfill {k}_{\mathrm{cut}}=\frac{1}{2.2}{k}_{\mathrm{s}}=\frac{1}{2.2}\frac{2\pi}{{r}_{\mathrm{s}}}\approx 0.3\phantom{\rule{0.166667em}{0ex}}h\phantom{\rule{0.166667em}{0ex}}{\mathrm{Mpc}}^{1}{\left(\frac{{R}_{\mathrm{c}}}{1\phantom{\rule{0.166667em}{0ex}}\mathrm{kpc}}\right)}^{2/3},\end{array}$$(35)
which is in close agreement with previous results on the SIBECDM transfer function and cutoff scale (Shapiro et al. 2021), apart from a slightly different prefactor and slope (0.3 → 0.2 and −2/3 → −1/2), obtained by considering the first kmode that enters the horizon while it is subJeans, instead of the effective sound horizon. Although this cutoff provides initial conditions appropriate for SIBECDM as described by the Lagrangian in Eq. (1), it poses a computational challenge that unfortunately could not be overcome at the present time. The large difference between the cutoff scale λ_{cut} = π/k_{cut} and the Jeans’ scale R_{c}, λ_{cut} ≫ R_{c}, means the simulations require both a large box and very high spatial resolution to simultaneously investigate the formation of SIBECDM structure and to resolve the SIBECDM halo cores. For instance, R_{c} = 1 kpc gives a cutoff at λ_{cut} ≈ 15 Mpc. Including an extra order of magnitude at each end to get sufficient resolution to resolve halo cores, as well as to have a large enough simulated volume to form halos, ends up being at least 6 orders of magnitude of spatial resolution that has to be simulated. We therefore chose to forgo the realistic initial conditions, and to instead use a simple Heaviside function f_{cut}(k) = Θ(k_{cut} − k), with k_{cut} equal to the Jeans’ length k_{J} = π/R_{c} of SIBECDM at the time the simulation starts at z = 50. In this way there is enough power in the initial conditions to give rise to a population of halos that can be analysed, while removing all perturbations on subJeans’ scales, although at the expense of a realistic cosmological history for SIBECDM.
The strong cutoff in Eq. (35) also has important consequences for the canonical parameter space of SIBECDM. By considering the largest halo masses that are affected by the cutoff,
$$\begin{array}{c}\hfill {M}_{\mathrm{cut}}\approx \frac{4\pi}{3}{\rho}_{\mathrm{dm},0}{\left(\frac{\pi}{{k}_{\mathrm{cut}}}\right)}^{3}\approx 5\times {10}^{14}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}\phantom{\rule{0.166667em}{0ex}}{\left(\frac{{R}_{\mathrm{c}}}{1\phantom{\rule{0.166667em}{0ex}}\mathrm{kpc}}\right)}^{2},\end{array}$$(36)
it is clear that SIBECDM with core radii of R_{c} ∼ 1 kpc results in a suppression in structure formation that is clearly at odds with observations. This was first pointed out and discussed in detail by Shapiro et al. (2021). In their paper, Shapiro et al. (2021) computed the halo mass function (HMF) for SIBECDM from linear theory, and found that a SIBECDM Jeans’ length as low as R_{c} ∼ 10 pc is needed to be consistent with observations, which poses a serious challenge for SIBECDM as a candidate to solve, for instance, the cuspcore problem. In the present work we nevertheless use R_{c} ∼ 1 kpc in our effort to explore the basic properties of SIBECDM halos. Our results are therefore in a sense more directly relatable to previous studies in which SIBECDM with halo cores of order 1 kpc is assumed to have a matter power spectrum very close to CDM (Harko 2011; Harko & Mocanu 2012; Velten & Wamba 2012; Freitas & Gonçalves 2013; Bettoni et al. 2014; de Freitas & Velten 2015; Berezhiani & Khoury 2015; Khoury 2016; Hartman et al. 2022), usually justified by invoking a phase transition from an initial CDMlike state into a latetime SIBECDM fluid, or an additional temperature dependence in the selfcoupling in order to make it weaker at high redshifts, although no concrete realization of any such mechanism or transition has yet been suggested.
Initial conditions for the effective thermal energy of SIBECDM present in Eqs. (25)–(27) must also be provided and is in principle given by the underlying wavefunction, although in the hydrodynamic approximation it has been integrated away with the smoothing kernel in Eq. (11). We therefore make the ansatz that the effective thermal pressure P is small compared to the interaction pressure P_{SI} at the start of the simulation, and set $\stackrel{\sim}{P}({z}_{\mathrm{init}})=\zeta {\stackrel{\sim}{P}}_{\mathrm{SI}}({z}_{\mathrm{init}})$ with ζ ≪ 1. The supercomoving effective thermal pressure at the background level is constant in time, while the interaction pressure decreases as ${\stackrel{\sim}{P}}_{\mathrm{SI}}\propto {a}^{1}$, so ζ should be sufficiently small that $\stackrel{\sim}{P}$ is largely negligible until the first structures start to form, at which point the effective thermal energy increases rapidly as matter is accreted onto halos and is heated.
The simulations performed in this work are summarized in Table 1, and an overview of the conservation of energy, resolution, and runtime for a few of these are shown in Fig. 2 at two different initial refinement levels. The conservation of energy is satisfied to within 1%, and the Jean’s length is resolved in overdense regions with the use of AMR at all times, more so at later times as the initial perturbations collapse to form halos. On the other hand, the initial resolution of the grid, which is at the lowest level of refinement, places a lower bound on the halo masses that are accurately evolved. An estimate of this bound is the total mass enclosed by a minimum number of cells N_{min} on the initial grid,
$$\begin{array}{c}\hfill {M}_{\mathrm{min}}={\rho}_{\mathrm{dm},0}\phantom{\rule{0.166667em}{0ex}}{N}_{\mathrm{min}}{\left(\frac{L}{{2}^{{l}_{\mathrm{init}}}}\right)}^{3}.\end{array}$$(37)
Fig. 2. Energy conservation (top), the smallest resolved physical cell length (middle), and the runtime in CPU hours (bottom) for the main suite of simulations. Runs with a minimum refinement level of 8 are shown in solid lines, and 7 in dotted lines. 
Overview of the simulation runs.
We find that N_{min} = 300, the mass enclosed by an approximately 7 × 7 × 7 cube, provides a reasonable estimate, as shown in Fig. 3. In this case M_{min} ≈ 10^{7} M_{⊙} for l_{init} = 8, whereas at l_{init} = 7 the mass limit is M_{min} ≈ 10^{8} M_{⊙}. The profiles are largely the same except in the lowest mass bin, which is below M_{min} for l_{init} = 7. Additionally, the HMF is suppressed below M_{min}, hence halos with M_{200} < M_{min} are discarded in the following analysis. It should be noted that the higher resolution yields slightly denser halos even for the masses that are most resolved, and that the HMF above M_{min} is not exactly identical at the two refinement levels, as one would expect. This indicates our simulations would benefit from further grid refinement, in particular the Rc1 run, which has the highest ratio of both M_{min}/M_{cut} and Δx_{min}/R_{c}. Unfortunately, the computational cost of further increasing the resolution was prohibitive, since at l_{init} = 8 the simulation used up to a hundred thousand CPU hours to reach z = 0.5, while an increase to the next refinement level l_{init} = 9 is expected to use a few million CPU hours.
Fig. 3. Binned SIBECDM halo profiles and the HMF for R_{c} = 1 kpc at z = 1.5, with minimum refinement level 7 (dotted) and 8 (solid). The HMF of CDM for the same box size and initial resolution is included in black. The halo mass limit M_{min} are indicated by dots. The standard deviation from the binned halo profile mean for level 8, as well as the Poisson error in the HMF, are shown in shaded. 
4. Results
In the cosmological simulations, the SIBECDM collapses and forms halos with a cored inner structure and a CDMlike envelope, which we find to be wellfitted by a cored NFW (NFWc) profile,
$$\begin{array}{c}\hfill {\delta}_{\mathrm{NFWc}}(r)={[{\delta}_{\mathrm{c}}^{1}+{\delta}_{\mathrm{NFW}}^{1}]}^{1}.\end{array}$$(38)
The profile is given in terms of overdensity $\delta =\rho /\overline{\rho}$, where δ_{c} is the central overdensity, and δ_{NFW} is the NFW profile,
$$\begin{array}{c}\hfill {\delta}_{\mathrm{NFW}}(r)=\frac{{\delta}_{\mathrm{s}}}{\frac{r}{{r}_{\mathrm{s}}}{(1+\frac{r}{{r}_{\mathrm{s}}})}^{2}}\xb7\end{array}$$(39)
We define the core radius r_{c} as where the density has dropped to 50% of its central value, δ(r_{c}) = 0.5δ(0). At hydrostatic equilibrium we would have r_{c} ≈ 0.6 R_{c}. A number of other cored profiles could also be used to fit the SIBECDM halos. In Li et al. (2020) several DM profiles were fitted to the Spitzer Photometry & Accurate Rotation Curves (SPARC) dataset of rotation curves of nearby galaxies (Lelli et al. 2016), such as the Burkert profile (Burkert 1995),
$$\begin{array}{c}\hfill {\delta}_{\mathrm{Burkert}}(r)=\frac{{\delta}_{\mathrm{c}}}{(1+\frac{r}{{r}_{\mathrm{s}}})[1+{\left(\frac{r}{{r}_{\mathrm{s}}}\right)}^{2}]},\end{array}$$(40)
the Einasto profile (Einasto 1965),
$$\begin{array}{c}\hfill {\delta}_{\mathrm{Einasto}}(r)={\delta}_{\mathrm{s}}exp\{\frac{2}{{\alpha}_{\u03f5}}[{\left(\frac{r}{{r}_{\mathrm{s}}}\right)}^{{\alpha}_{\u03f5}}1]\}\end{array}$$(41)
and the socalled Lucky13 profile (Li et al. 2020),
$$\begin{array}{c}\hfill {\delta}_{13}(r)=\frac{{\delta}_{\mathrm{c}}}{[1+{\left(\frac{r}{{r}_{\mathrm{s}}}\right)}^{3}]}\xb7\end{array}$$(42)
These are fitted to the SIBECDM halos by minimizing the cost function
$$\begin{array}{c}\hfill {\chi}_{\nu}^{2}=\frac{1}{{N}_{\mathrm{d}}{N}_{\mathrm{f}}}{\displaystyle \sum _{i=1}^{{N}_{\mathrm{d}}}}\frac{{[{\delta}_{i}\delta ({r}_{i})]}^{2}}{{\delta}_{i}^{2}},\end{array}$$(43)
with respect to the N_{f} free parameters of the density function δ(r), where N_{d} is the number of radial points δ_{i} from the simulation that we are fitting the density profile to. Each halo in the simulation gives a value for ${\chi}_{\nu}^{2}$, and the corresponding cumulative distribution function (CDF) for the set of these tells us what fraction of halos has ${\chi}_{\nu}^{2}$ smaller than a given value, and therefore how well the fitting function generally describes the halos. In Fig. 4 the CDF for the different types of profiles are shown, from which we find that the NFWc profile provides the best fit overall. Of the cored profiles used by Li et al. (2020) to fit the DM halos in the SPARC dataset, the Burkert and Einasto profiles are best suited and perform about the same, but as Burkert is the simpler of the two, we will use it in the following.
Fig. 4. Cumulative distribution functions for different cored halo profiles from SIBECDM simulations with R_{c} = 1 kpc (upper), R_{c} = 3 kpc (middle), and R_{c} = 10 kpc (lower) at z = 0.5. This shows the fraction of the halos that has ${\chi}_{v}^{2}$ smaller than a given value, and therefore how well the fitting function generally describes the halos in our simulation. 
Both NFWc and Burkert fitted to binned SIBECDM halos with R_{c} = 1 kpc at z = 0.5 are shown in Fig. 5. The two profiles yield slightly different results for the core radius. Burkert is a bit steeper than NFWc near the core, and therefore generally gives a core radius r_{c} of around half compared to NFWc (according to our definition δ(r_{c}) = 0.5δ(0)). The core density δ_{c} is also about twice as high for Burkert compared to NFWc. In the following we use the Burkert profile, even though the fit using NFWc is better, for two reasons; it is a simpler function, and we can readily compare the fitting parameters of the SIBECDM halos to observations, because the Burkert profile has previously been used for the DM component of nearby galaxies in the SPARC dataset (Li et al. 2020) and the dwarf spheroidal (dSph) satellites of the Milky Way (Salucci et al. 2012).
Fig. 5. Mean (solid) and standard deviation (shaded) of binned SIBECDM halo profiles with R_{c} = 1 kpc at z = 0.5. The binned profiles are fitted to NFWc (dashed) and Burkert (dashdotted). The fits of the NFW profile to the halo envelopes are also shown (dotted). 
In order to investigate the general trends in the SIBECDM halos in our simulations, and compare these to observations, we introduce scaling functions for the core radius r_{c}, the central density δ_{c}, and the mass enclosed inside the core radius, M_{c};
$$\begin{array}{cc}& {r}_{\mathrm{c}}({M}_{200})={r}_{\mathrm{c},10}{\left(\frac{{M}_{200}}{{10}^{10}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}}\right)}^{\alpha},\hfill \end{array}$$(44)
$$\begin{array}{cc}& {\delta}_{\mathrm{c}}({M}_{200})={\delta}_{\mathrm{c},10}{\left(\frac{{M}_{200}}{{10}^{10}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}}\right)}^{\beta},\hfill \end{array}$$(45)
$$\begin{array}{cc}& {M}_{\mathrm{c}}({M}_{200})={M}_{\mathrm{c},10}{\left(\frac{{M}_{200}}{{10}^{10}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}}\right)}^{\gamma}.\hfill \end{array}$$(46)
These are fitted using TheilSen regression (Theil 1950; Sen 1968). This method finds the slope m and the yintercept b of the linear function y(x) = b + mx (which Eqs. (44)–(46) are in logspace) by first computing the median of slopes between all pairs of points (y_{i}, x_{i}), which gives m, and then finding the median of y_{i} − mx_{i} for all points, which gives b. This fitting procedure is robust against outliers, and provides a simple measure of the variation in m and b present in the data through, for example, the first and third quartiles of the pairwise slopes and individual yintercepts. Tables listing the scaling parameters fitted using this procedure are included in Appendix A for several redshifts.
The properties of the Burkert profile fitted to SIBECDM halos are shown in Fig. 6 for R_{c} = 1 kpc, 3 kpc, and 10 kpc, along with the above scaling relations. The fit to the nearby galaxy rotation curves in the SPARC dataset (Li et al. 2020) and Milky Way dSphs (Salucci et al. 2012) are also included in Fig. 6.
Fig. 6. Core radii r_{c} (left), core densities δ_{c} (middle), and core masses M_{c} (right) versus M_{200} for halos in cosmological simulations of SIBECDM with R_{c} = 1 kpc, 3 kpc, and 10 kpc, using the Burkert profile. The upper plots show the scatter of halos at z = 0.5, with the fitted scaling Eqs. (44)–(46) shown in dashed lines, compared to the SPARC dataset (Li et al. 2020) and Milky Way dSphs (Salucci et al. 2012). The lower plots show the fitted median (solid) and the first and third quartiles (shaded) at several redshifts using TheilSen regression. The results from SPARC and the dSphs are also shown. In the scatter plot for r_{c} the core radii R_{c} are indicated in dotted colored lines. 
While the halo catalogues from the simulation snapshots only span 2−3 orders in magnitude in mass and number 200−400 halos due to the limited volume of the simulated box, and the fitted parameters of the scaling functions show some timedependence (the prefactors more so than the exponents), there are a number of interesting features to note. The halo core radii r_{c} are scattered near R_{c}, and the dependence of r_{c} on M_{200} is nearly constant, largely as expected from hydrostatic considerations, although there is a slight positive trend, ${r}_{\mathrm{c}}\sim {M}_{200}^{\alpha}$ with α ≈ 0.05 − 0.1. Also, the SIBECDM cores are generally larger compared to what we find from Eq. (9). Both of these observations are consistent with the addition of thermal pressure in the halos centers, which is also slightly increasing for more massive halos, as seen in Fig. 7. In fact, the thermallike energy due to the hidden dynamics on the de Brogliescale dominates the central regions of the SIBECDM halos. This is illustrated in Fig. 8, where the selfinteraction, thermal, and gravitational potential energy densities of a sample halo are shown. As structure forms, the SIBECDM fluid is heated as it is accreted onto the growing halos, the kinetic energy of the infalling matter being converted into effective thermal energy. The result are halos with a thermal energy much larger than the interaction energy throughout the halo, even in the core where it dominates the overall dynamics. Additionally, the thermal energy is nearly isothermal. At the edge of the core the potential energy becomes larger than the internal energies, and the halo transitions to an NFWlike envelope.
Fig. 7. Ratio of average thermal and selfinteraction energy, U and U_{SI}, inside the halo cores r < r_{c}/2 in simulations with R_{c} = 3 kpc at z = 1. 
We illustrate with the sample halo in Fig. 8, which is relatively isolated in the simulated box, that the SIBECDM halos in the simulations achieve approximate virial equilibrium. The virial 𝒱 satisfies, under the assumption of spherical symmetry (Dawoodbhoy et al. 2021),
$$\begin{array}{c}\hfill \mathcal{V}=2\mathcal{T}+\mathcal{W}+2\mathcal{U}+3{\mathcal{U}}_{\mathrm{SI}}+\mathcal{S}+{\mathcal{S}}_{\mathrm{SI}}=0,\end{array}$$(47)
Fig. 8. Top panel: selfinteraction energy density U_{SI}, thermal energy density U, and gravitational potential energy density W of a sample halo with R_{c} = 1 kpc and M_{200} = 3 × 10^{9} M_{⊙}, relative to the interaction energy at the background level ${\overline{U}}_{\mathrm{SI}}$. The middle panel shows instead the specific energy, and the bottom panel the cumulative virial. The core radius r_{c} and halo radius r_{200} are indicated with the inner and outer dotted vertical lines, respectively. 
where the various terms are cumulative energy contributions and surface terms,
$$\begin{array}{cc}& \mathcal{T}=\frac{1}{2}\int {v}^{2}\mathrm{d}M2\pi {r}^{2}(\rho r{v}^{2}+\frac{1}{2}{r}^{2}\frac{\partial (\rho v)}{\partial t}),\hfill \end{array}$$(48)
$$\begin{array}{cc}& \mathcal{W}={\displaystyle \int}\frac{\mathit{GM}}{r}\mathrm{d}M,\hfill \end{array}$$(49)
$$\begin{array}{cc}& \mathcal{U}={\displaystyle \int}\frac{U}{\rho}\mathrm{d}M,\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\mathcal{U}}_{\mathrm{SI}}=\int \frac{{U}_{\mathrm{SI}}}{\rho}\mathrm{d}M,\hfill \end{array}$$(50)
$$\begin{array}{cc}& \mathcal{S}=4\pi {r}^{3}P,\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\mathcal{S}}_{\mathrm{SI}}=4\pi {r}^{3}{P}_{\mathrm{SI}},\hfill \end{array}$$(51)
and 𝒰_{tot} = 𝒰 + 𝒰_{SI}.
Our 3D simulations confirm many of the features of SIBECDM halos that were reported in Dawoodbhoy et al. (2021), such as the transition from a cored interior to a CDMlike envelope near R_{c}, and the domination of the thermal energy over selfinteraction energy outside the core due to shock heating. The only significant difference is the domination of thermal energy inside the cores of halos in 3D, whereas in 1D simulations the cores are not heated and are instead dominated by selfinteractions. In fact, the halos from our simulations look more like the doublepolytrope solutions of Dawoodbhoy et al. (2021) for SIBECDM with both the selfinteraction pressure and an isothermal thermal pressure at hydrostatic equilibrium. Dawoodbhoy et al. (2021) anticipated that after collapse, relaxation, and finally virial equilibrium, the thermal profile should be nearly isothermal, with core radii set by the selfinteraction pressure and thus on the order of R_{c}, which is what we find in our simulations.
We attribute the difference between our simulations and the ones in Dawoodbhoy et al. (2021) to the absence of mixing in spherical symmetry. In 1D, halos form by the accretion of spherically symmetric shells that are heated as they are slowed down, resulting in an accretion shock expanding outwards as shells collide, but do not cross. The cores themselves experience little heating as the fluid is decelerated smoothly by the selfinteraction pressure and without a shock. In fully 3D simulations the halos form as clumps – not perfect shells – of matter merge, which perturb the halos and causes the outer shock heated layers to mix with the interior. We demonstrate this difference between symmetrical and asymmetrical collapse using 2D simulations; one in which an initial stationary overdensity is centered in the simulated box and collapses symmetrically; and another with a second smaller overdensity offset at x_{2} that merges with the first in an asymmetrical manner. The initial density field in the two scenarios are ρ(x) = ρ_{0} + ρ_{0}Δ_{1}e^{−(x/R1)2}, and ρ(x) = ρ_{0} + ρ_{0}Δ_{1}e^{−(x/R1)2} + ρ_{0}Δ_{2}e^{−(x − x2/R2)2}, and we use a periodic box of size L = 40 kpc with R_{c} = 1 kpc, R_{1} = 5 kpc, R_{2} = 2 kpc, Δ_{1} = Δ_{2} = 99, x_{2} = 10 kpc, and ρ_{0} = Ω_{m0}ρ_{c0}, with the initial ratio ζ = P/P_{SI} = 10^{−2}. The “symmetrical” simulation was run for 20 t_{dyn}, while the “asymmetrical” one was run for 200 t_{dyn}, to give it enough time to virialize, where ${t}_{\mathrm{dyn}}=1/\sqrt{2\pi G{\rho}_{0}{\mathrm{\Delta}}_{1}}$ is the dynamical timescale of the central overdensity. The final distributions of U/U_{SI} are shown in Fig. 9. In the symmetrical case, the thermal energy is dominant everywhere except in the core, whereas in the asymmetrical case the merging with the second overdensity results in a strong mixing of the fluid layers, causing shockheated fluid in the region outside the core to penetrate into the core.
Fig. 9. Final distributions of U/U_{SI} in 2D simulations with symmetrical (upper) and asymmetrical (lower) collapse, centered on the density peak. In the asymmetrical case, the second smaller overdensity was initially located to the right. The minimum in the symmetrical case is around U/U_{SI} ≈ 2 × 10^{−2}, while in the asymmetrical case its is U/U_{SI} ≈ 2. 
In Fig. 10 the internal energies of SIBECDM with R_{c} = 1 kpc at z = 0.5 are shown, where the extended envelopes of high thermal energies compared to the interaction energy are clearly seen around the collapsed structures. In the same figure the density field is compared to standard CDM, which illustrates the smoothing of the DM density field, especially in the halo cores. The corresponding matter power spectrum is shown in Fig. 11 at several redshifts, as well as CDM with the same initial conditions, but without the cutoff, with
$$\begin{array}{c}\hfill {\mathrm{\Delta}}^{2}(k)=\frac{{k}^{3}P(k)}{2{\pi}^{2}},\end{array}$$(52)
Fig. 10. Projection plots of the DM density in cosmological simulations with L = 2 Mpc h^{−1} at z = 0.5 SIBECDM with R_{c} = 1 kpc (upper right) and CDM with the same initial conditions, but without a cutoff (upper left). The SIBCEDM internal energies for the same snapshot showing the effective thermal energy U (lower left) and the selfinteraction energy U_{SI} (lower right). 
Fig. 11. Matter power spectrum for cosmological simulations L = 2 Mpc h^{−1} at z = 0.5 for SIBECDM with R = 1 kpc (solid) and CDM (dashed). The comoving cutoff k_{cut} is indicated with a dotted vertical line, and the spectrum at z = 50 is multiplied by 50 in the figure. 
where
$$\begin{array}{c}\hfill {(2\pi )}^{3}P(k){\delta}^{3}(\mathit{k}{\mathit{k}}^{\prime})=\widehat{\delta}(\mathit{k})\widehat{\delta}({\mathit{k}}^{\prime}),\end{array}$$(53)
and
$$\begin{array}{c}\hfill \widehat{\delta}(\mathit{k})={\displaystyle \int}{\mathrm{d}}^{3}x\phantom{\rule{0.166667em}{0ex}}\delta (\mathit{x}){e}^{i\mathit{k}\xb7\mathit{x}},\end{array}$$(54)
with δ^{3} being the 3D Dirac delta function. At the largest scales the SIBECDM modes grow like CDM, while at smaller scales SIBECDM is suppressed compared to CDM due to the DM fluid pressure.
The simulations are insensitive to the initial ratio P/P_{SI} as long as it is much smaller than unity, as shown in Figs. 12 and 13. Below ζ = P/P_{SI} ≲ 0.1, there is little change in the final halos. As soon as structure begins to form, the thermal energy increases due to the conversion of kinetic energy to internal energy, hence, a small amount of initial thermal energy matters little. Nevertheless, despite the effective thermal energy dominating the DMhalo cores, it is the strength of the selfinteraction that determines the characteristic size of the cores. This is because it is at the Jeans’ scale due to the selfinteractions that the pressure forces begin opposing gravity, creating the shocks that cause the heating of the fluid.
Fig. 12. Evolution of the ratio of the total thermal energy and total selfinteraction energy for the simulations run. 
Fig. 13. Scaling functions at several redshifts for the core radii r_{c} (left), the core density δ_{c} (middle), and the core mass M_{c} (right) for SIBECDM halos with R_{c} = 3 kpc and different initial ζ = P/P_{SI} (upper), and different initial power cutoffs (lower). The shaded areas give the first and third quartiles for the scaling parameters obtained from the TheilSen regression. 
Another interesting feature is the dependence of δ_{c} on the halo mass. The Burkert profile fitted to observational data prefers a negative slope, meaning massive halos are less dense. In Nbody simulations of CDM, considering instead the characteristic overdensity δ_{s} of the halo in Eq. (39), such mass dependence arises, in very simple terms, because δ_{s} is proportional to the mean density at the time the halo forms, with lowmass halos generally collapsing at a higher redshift, when the universe was denser (Navarro et al. 1996). We might therefore expect processes that transform the CDM halo cusps into cores, for instance baryonic feedback, to produce core densities that inherit the same halo mass dependence (Ogiya et al. 2014). SIBECDM halos prefer instead a positive slope β ≈ 0.5 at z = 0.5, such that more massive halos are also denser. Furthermore, the SIBECDM halo core masses M_{c} scales approximately as ${M}_{\mathrm{c}}\sim {M}_{200}^{\gamma}$, with γ ≈ 0.75 at z = 0.5. The scaling of both M_{c} and δ_{c} can be qualitatively understood under the assumption of velocity dispersion tracing outside the core, meaning that the circular orbital velocity is nearly constant (Chavanis 2019a,b; Padilla et al. 2019; Dawoodbhoy et al. 2021),
$$\begin{array}{c}\hfill {v}_{\mathrm{c}}^{2}=\frac{G{M}_{\mathrm{c}}}{{r}_{\mathrm{c}}}\approx {v}_{200}^{2}=\frac{G{M}_{200}}{{R}_{200}},\end{array}$$(55)
which we indeed see to be approximately true in Fig. 8, since v^{2} = −W/ρ. Inserting that r_{c} is constant, and ${M}_{200}\sim {r}_{200}^{3}$, gives
$$\begin{array}{c}\hfill {M}_{\mathrm{c}}\sim {M}_{200}^{2/3},\end{array}$$(56)
and
$$\begin{array}{c}\hfill {\delta}_{\mathrm{c}}\propto \frac{{M}_{\mathrm{c}}}{{r}_{\mathrm{c}}^{3}}\sim {M}_{200}^{2/3},\end{array}$$(57)
which are in the neighborhood of what is obtained in the simulations. We can improve this scaling argument by including a small mass dependence in ${r}_{\mathrm{c}}\sim {M}_{200}^{\alpha}$, which gives
$$\begin{array}{cc}& {M}_{\mathrm{c}}\sim {M}_{200}^{\frac{2}{3}+\alpha},\hfill \end{array}$$(58)
$$\begin{array}{cc}& {\delta}_{\mathrm{c}}\sim {M}_{200}^{\frac{2}{3}2\alpha}.\hfill \end{array}$$(59)
Inserting α = 0.05 gives ${M}_{\mathrm{c}}\sim {M}_{200}^{0.72}$ and ${\delta}_{\mathrm{c}}\sim {M}_{200}^{0.57}$, while α = 0.1 gives ${M}_{\mathrm{c}}\sim {M}_{200}^{0.77}$ and ${\delta}_{\mathrm{c}}\sim {M}_{200}^{0.47}$. For comparison, independent FDM simulations, using the full NLSE or the Madelung equations, find varying exponents, with ${M}_{\mathrm{c}}\sim {M}_{200}^{1/3}$ (Schive et al. 2014a,b; Veltmaat et al. 2018), ${M}_{\mathrm{c}}\sim {M}_{200}^{5/9}$ (Mocz et al. 2017; Mina et al. 2022), or ${M}_{\mathrm{c}}\sim {M}_{200}^{0.6}$ (Nori & Baldi 2021), all of which are less steep than what we find for SIBECDM. Furthermore, as expected from hydrostatic equilibrium, the size of FDM halo cores have been found to be decreasing with the virial mass (Chan et al. 2022), as opposed to slightly increasing for SIBECDM.
Generally, we find the trends in r_{c}, δ_{c}, and M_{c} for the values of R_{c} and initial conditions tested in this paper to be in conflict with the SPARC dataset and Milky Way dSphs, although their slopes are in less tension with observations compared to FDM. The trends are mostly insensitive to the initial ratio P/P_{SI} ≪ 1, but they are dependent on the cutoff scale in the initial matter power spectrum (when it is not scaled along with R_{c} to match the Jeans’ length), as shown in Fig. 13. By moving the initial cutoff to larger scales, such that there is less initial power at small scales, the scaling exponents shift towards those inferred from the SPARC dataset and the Milky Way dSphs. On the flipside, the prefactors shift away from the observations. For example, the total ratio P/P_{SI} increases as k_{cut} is lowered, meaning there is more thermal pressure, causing halos to have larger cores and lower central densities. However, as we saw in Fig. 6, decreasing the interaction strength, or equivalently, the core radius R_{c}, hardly affects the scaling exponents, but moves the prefactors in the desired direction. Based on this observation, one could imagine that a sufficiently small R_{c} and strong cutoff might create a population of halos that alleviates to some degree the current tension with the observed trends. According to Eqs. (58) and (59), a stronger mass dependence in the core radius ${r}_{\mathrm{c}}\sim {M}_{200}^{\alpha}$ will bring the SIBEC trends towards the observed ones, and can come about by having a larger thermal Jeans’ length for massive halos, such that more massive halos have a larger ratio of thermal energy to selfinteraction energy, U/U_{SI}, due to increased heating, as observed in our simulations.
However, a good deal of caution is called for in making such a statement. For instance, the range of halo masses represented in the simulations is rather small, and baryonic physics are not included at all, but are known to affect the distribution of DM in halos. Furthermore, the scaling function parameters have some timedependence, primarily the prefactors, and our simulations were only run until z = 0.5, whereas the observed galaxies and dSphs are at z ≈ 0. This is particularly true for the Rc3b run with k_{cut} = k_{J}/2, which has about a fourth of the halos compared to k_{cut} = k_{J}, and is therefore more sensitive to disruptions in parts of the simulated halo population due to events such as mergers. We see this in Fig. 13 at z = 0.5, where the scaling parameters suddenly change, and the first and third quartileregion becomes very large, signaling that the scaling relations generally do not provide a good fit to the halo population at that point. Finally, but not least, the SIBECDM parameter space and initial conditions used in our simulations are not appropriate for a realistic realization of scalar field DM with selfinteractions as given by the Lagrangian in Eq. (1), for which R_{c} ≳ 1 kpc are ruled out and the suppression in the initial power spectrum is much stronger (Shapiro et al. 2021; Hartman et al. 2022).
5. Conclusion
In this paper, a hydrodynamic approximation of the NLSE that includes the de Brogliescale dynamics as an effective thermal energy was used to simulate the formation of structure of SIBECDM in a fully 3D cosmological setting. The advantage of such an approach is that the dynamics on the de Brogliescale need not be resolved on the numerical grid, making the task of simulating SIBECDM simpler and computationally tractable. This is of particular importance when the de Broglie wavelength is much smaller than the selfinteraction Jeans’ scale that we are interested in. On the other hand, this approach introduces an additional equation for the effective thermal energy that must be evolved through time, information about the underlying wavefunction is lost, and some simplifying assumptions were made, such as the skewlessness of the phase space distribution function.
Our simulations reproduce many of the features reported in Dawoodbhoy et al. (2021). The SIBECDM halos have cores with radii of order R_{c}, which is only weakly dependent on the halo mass. Outside the cores the halos transition to a CDMlike profile, which in our 3D simulations is the NFW profile. The SIBECDM halos are therefore generally welldescribed by a cored NFWprofile such as NFWc, Eq. (38), or the Burkert profile, Eq. (40). Despite the selfinteractions determining the scale of the cores, the effective thermal energy ends up dominating over the selfinteractions throughout the halos due to significant heating during collapse, even in the central regions. This is in contrast to the 1D simulations of Dawoodbhoy et al. (2021), who found the SIBECDM halo cores to not experience much heating, and the thermal energy to fall sharply below the selfinteraction energy inside the core radius. However, Dawoodbhoy et al. (2021) anticipated that after collapse and relaxation, the SIBECDM halos should be largely isothermal with the central thermal energy on the order of the selfinteraction energy. Spherically symmetric 1D simulations do not include the mixing of fluid layers that allows for such heating of the core, while in 3D, matter collapses onto the halos in clumps rather than shells, causing the outer shockheated layers to be mixed with the central fluid, resulting in an isothermal profile. Our simulations therefore confirm this anticipated result of Dawoodbhoy et al. (2021).
Scaling relations for the core radii r_{c}, core densities δ_{c}, and core masses M_{c} as functions of the total halo mass M_{200} were fitted to the simulated halo populations, which largely agree with hydrostatic considerations of the halo cores where r_{c} is nearly constant, as well as velocity dispersion tracing in the halo envelope, ${v}_{\text{c}}^{2}\approx {v}_{200}^{2}$. However, these trends do not agree with those obtained by fitting the Burkert profile to nearby galaxies in the SPARC dataset and the classical Milky Way dSphs. This poses an issue for SIBECDM with R_{c} ≳ 1 kpc and a largely CDMlike matter power spectrum at late times (Harko 2011; Harko & Mocanu 2012; Velten & Wamba 2012; Freitas & Gonçalves 2013; Bettoni et al. 2014; de Freitas & Velten 2015; Hartman et al. 2022), as was used in our simulations, although these scenarios are not wellmotivated. For SIBECDM given by the field Lagrangian in Eq. (1), the selfinteraction is constrained to R_{c} < 1 kpc, otherwise an early radiationlike period and a large comoving Jeans’ length washes out too much structure to be consistent with observations (Shapiro et al. 2021; Hartman et al. 2022). In fact, Shapiro et al. (2021) found by using constraints on FDM as a proxy for SIBECDM, and matching their transfer function cutoffs and HMFs, that the SIBECDM selfinteraction should be as low as R_{c} ∼ 10 pc to not be in conflict with observations. We were unable to probe SIBECDM with initial conditions and parameters consistent with the Lagrangian in Eq. (1), since the large gap between the halo cores and the cutoff scale requires both a large simulation box and very high spatial resolution. It should be noted that our SIBECDMonly simulations do provide a better agreement with the slopes in observed scaling relations than FDM. In particular, FDM simulations generally find ${M}_{\mathrm{c}}\sim {M}_{200}^{\gamma}$ with 1/3 < γ < 0.6, while we find γ ≈ 0.75, which is closer to the observed γ ≈ 1.1. Additionally, FDM halos have core radii that generally decrease with the halo mass, while we find a slightly increasing trend due to larger halos experiencing more thermal heating, although not as steep as in the SPARC dataset and the Milky Way dSphs.
The main results of the present work is summarized as follows:

Using a hydrodynamic approximation and fully 3D cosmological simulations, SIBECDM halos are found to have NFWlike envelopes with central core radii of order R_{c}. This confirms the results reported in Shapiro et al. (2021), based on 1D spherically symmetric simulations in the same hydrodynamic approximation, when the initial perturbation is chosen to yield an infall rate that matches the mass assembly history of 3D Nbody simulations of CDM. The density profiles are therefore wellfitted by cored NFW profiles such as NFWc, Eq. (38), and the Burkert profile, Eq. (40).

The SIBECDM core radii are only weakly dependent on the halo mass, largely as expected from hydrostatic equilibrium. The slight increase is due to larger halos generally experiencing more heating, and hence an increase in their thermal Jeans’ length.

Despite the selfinteraction energy initially dominating the internal energy and fluid pressure of the SIBECDM, as well as determining the general scale of the SIBECDM cores, the effective thermal energy ends up dominating throughout the halo. This result is insensitive to the initial ratio of the thermal energy relative to the selfinteraction energy, as long as it is small, since the final thermal energy comes from heating as SIBECDM collapses and forms structure. Below ζ = P/P_{SI} ≲ 0.1, there are hardly any changes in the final halos.

The SIBECDM halo core densities δ_{c} and core masses M_{c} are found to scale with the virial mass M_{200} as ${\delta}_{\mathrm{c}}\sim {M}_{200}^{0.5}$ and ${M}_{\mathrm{c}}\sim {M}_{200}^{0.75}$. This result is insensitive to changes in R_{c} that is matched by a corresponding change in the cutoff scale, k_{cut} ∼ 1/R_{c}. Velocity tracing in the halo envelope and a nearly constant core radius predict these relations to be ${\delta}_{\mathrm{c}}\sim {M}_{\mathrm{c}}\sim {M}_{200}^{2/3}$. Including a small mass dependence in r_{c} improves the agreement with the simulations, ${r}_{\mathrm{c}}\sim {M}_{200}^{0.05}$, gives ${M}_{\mathrm{c}}\sim {M}_{200}^{0.72}$ and ${\delta}_{\mathrm{c}}\sim {M}_{200}^{0.57}$.

The scaling relations for r_{c}(M_{200}), δ_{c}(M_{200}), and M_{c}(M_{200}), Eqs. (44)–(46), assuming the Burkert profile for the SIBECDM halos, generally do not agree with trends of the classical Milky Way dSphs and nearby galaxies in the SPARC dataset for the SIBECDM parameters tested in this work. Nevertheless, the slopes of the scaling relations in our SIBECDMonly simulation are in better agreement with observations compared to FDMonly simulations.
In future work, larger highresolution simulations should be carried out to further investigate SIBECDM, as our current simulations have a limited box size and halo population. This is particularly relevant for probing SIBECDM with realistic initial conditions and selfinteractions that are consistent with the HMF (Shapiro et al. 2021) and largescale observables (Hartman et al. 2022). It will be important to test the validity and accuracy of the smoothing procedure and the resulting hydrodynamic approximation used in this work, especially if it continues to be used to test SIBECDM models. Dawoodbhoy et al. (2021) showed analytically for a few idealized 1D cases that smoothing the exact solutions of the NLSE gives the same largescale fluid properties, such as pressure and density, as the smoothed phase space distribution function. These tests should be extended to more complicated scenarios, such as selfgravitation and the formation of shock fronts, to check, for instance, if the heating of the effective thermal energy is accurately captured by the hydrodynamic approximation, or if the skewlessness of ℱ is a valid assumption, both inside virialized halos and during the infall and collapse. In particular, the assumption of skewlessness should be tested in mergers, since the merging SIBECDM halos can have large planewave velocity dispersions, which would have significant implications for the merger event, although this will require using large and costly numerical simulations.
Acknowledgments
We thank the Research Council of Norway for their support, and the anonymous referee for their helpful comments and suggestions. We are also grateful to Matteo Nori and Marco Baldi for our discussions on simulations of ultralight DM, and to Taha Dawoodbhoy and Paul R. Shapiro for their feedback on this work. Computations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway.
References
 Ahn, K., & Shapiro, P. R. 2005, MNRAS, 363, 1092 [NASA ADS] [CrossRef] [Google Scholar]
 Arbey, A., Lesgourgues, J., & Salati, P. 2002, Phys. Rev. D, 65, 083514 [NASA ADS] [CrossRef] [Google Scholar]
 Arbey, A., Lesgourgues, J., & Salati, P. 2003, Phys. Rev. D, 68, 023511 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhiani, L., & Khoury, J. 2015, Phys. Rev. D, 92, 103510 [Google Scholar]
 Berezhiani, L., Elder, B., & Khoury, J. 2019, JCAP, 2019, 074 [Google Scholar]
 Bettoni, D., Colombo, M., & Liberati, S. 2014, JCAP, 2014, 004 [CrossRef] [Google Scholar]
 Bullock, J. S., & BoylanKolchin, M. 2017, ARA&A, 55, 343 [Google Scholar]
 Burkert, A. 1995, ApJ, 447, L25 [NASA ADS] [Google Scholar]
 Böhmer, C. G., & Harko, T. 2007, JCAP, 2007, 025 [CrossRef] [Google Scholar]
 Chan, H. Y. J., Ferreira, E. G. M., May, S., Hayashi, K., & Chiba, M. 2022, MNRAS, 511, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Chavanis, P.H. 2011, Phys. Rev. D, 84, 043531 [NASA ADS] [CrossRef] [Google Scholar]
 Chavanis, P.H. 2019a, Phys. Rev. D, 100, 123506 [NASA ADS] [CrossRef] [Google Scholar]
 Chavanis, P.H. 2019b, Phys. Rev. D, 100, 083022 [NASA ADS] [CrossRef] [Google Scholar]
 Chavanis, P.H., & Delfini, L. 2011, Phys. Rev. D, 84, 043532 [NASA ADS] [CrossRef] [Google Scholar]
 Crăciun, M., & Harko, T. 2020, Eur. Phys. J. C, 80, 735 [CrossRef] [Google Scholar]
 Cyburt, R. H., Fields, B. D., Olive, K. A., & Yeh, T.H. 2016, Rev. Mod. Phys., 88, 015004 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]
 Dawoodbhoy, T., Shapiro, P. R., & RindlerDaller, T. 2021, MNRAS, 506, 2418 [NASA ADS] [CrossRef] [Google Scholar]
 de Freitas, R. C., & Velten, H. 2015, Eur. Phys. J. C, 75, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Del Popolo, A., & Le Delliou, M. 2017, Galaxies, 5, 17 [Google Scholar]
 Dine, M., & Fischler, W. 1983, Phys. Lett. B, 120, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Einasto, J. 1965, Trudy Astrofizicheskogo Instituta AlmaAta, 5, 87 [NASA ADS] [Google Scholar]
 Ferreira, E. G. M. 2021, A&ARv, 29, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Ferreira, E. G. M., Franzmann, G., Khoury, J., & Brandenberger, R. 2019, JCAP, 2019, 027 [CrossRef] [Google Scholar]
 Freitas, R. C., & Gonçalves, S. V. B. 2013, JCAP, 2013, 049 [CrossRef] [Google Scholar]
 Gill, S. P. D., Knebe, A., & Gibson, B. K. 2004, MNRAS, 351, 399 [Google Scholar]
 Goodman, J. 2000, New Astron., 5, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Hahn, O., & Abel, T. 2011, MNRAS, 415, 2101 [Google Scholar]
 Harko, T. 2011, Phys. Rev. D, 83, 123515 [NASA ADS] [CrossRef] [Google Scholar]
 Harko, T., & Mocanu, G. 2012, Phys. Rev. D, 85, 084012 [Google Scholar]
 Hartman, S. T. H., Winther, H. A., & Mota, D. F. 2022, JCAP, 2022, 005 [CrossRef] [Google Scholar]
 Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, K. 1987, Statistical Mechanics, 2nd edn. (John Wiley& Sons) [Google Scholar]
 Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541 [NASA ADS] [CrossRef] [Google Scholar]
 Husimi, K. 1940, Proc. PhysicoMath. Soc. Jpn., 22, 264 [Google Scholar]
 Khoury, J. 2016, Phys. Rev. D, 93, 103533 [Google Scholar]
 Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608 [Google Scholar]
 Lancaster, L., Giovanetti, C., Mocz, P., et al. 2020, JCAP, 2020, 001 [Google Scholar]
 Lee, J.W., & Koh, I.G. 1996, Phys. Rev. D, 53, 2236 [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157 [Google Scholar]
 Li, B., RindlerDaller, T., & Shapiro, P. R. 2014, Phys. Rev. D, 89, 083536 [NASA ADS] [CrossRef] [Google Scholar]
 Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2020, ApJS, 247, 31 [Google Scholar]
 Madelung, E. 1926, Naturwissenschaften, 14, 1004 [Google Scholar]
 Marsh, D. J. E. 2016, Phys. Rep., 643, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Martel, H., & Shapiro, P. R. 1998, MNRAS, 297, 467 [Google Scholar]
 Matos, T., & Arturo UreñaLópez, L. 2001, Phys. Rev. D, 63, 063506 [NASA ADS] [CrossRef] [Google Scholar]
 Matos, T., Guzmán, F. S., & UreñaLópez, L. A. 2000, Class. Quant. Grav., 17, 1707 [NASA ADS] [CrossRef] [Google Scholar]
 May, S., & Springel, V. 2021, MNRAS, 506, 2603 [NASA ADS] [CrossRef] [Google Scholar]
 Membrado, M., Pacheco, A. F., & Sañudo, J. 1989, Phys. Rev. A, 39, 4207 [NASA ADS] [CrossRef] [Google Scholar]
 Mina, M., Mota, D. F., & Winther, H. A. 2020, A&A, 641, A107 [EDP Sciences] [Google Scholar]
 Mina, M., Mota, D. F., & Winther, H. A. 2022, A&A, 662, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mocz, P., Vogelsberger, M., Robles, V. H., et al. 2017, MNRAS, 471, 4559 [Google Scholar]
 Mocz, P., Lancaster, L., Fialkov, A., Becerra, F., & Chavanis, P.H. 2018, Phys. Rev. D, 97, 083519 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
 Nori, M., & Baldi, M. 2018, MNRAS, 478, 3935 [Google Scholar]
 Nori, M., & Baldi, M. 2021, MNRAS, 501, 1539 [Google Scholar]
 Ogiya, G., Mori, M., Ishiyama, T., & Burkert, A. 2014, MNRAS, 440, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Padilla, L. E., RindlerDaller, T., Shapiro, P. R., Matos, T., & Vázquez, J. A. 2019, Phys. Rev. D, 103, 063012 [Google Scholar]
 Peebles, P. J. E. 2000, ApJ, 534, L127 [NASA ADS] [CrossRef] [Google Scholar]
 Percival, W. J., Baugh, C. M., BlandHawthorn, J., et al. 2001, MNRAS, 327, 1297 [Google Scholar]
 Pitaevskii, L. P., & Stringari, S. 2016, BoseEinstein Condensation and Superfluidity (Oxford: Oxford University Press) [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Preskill, J., Wise, M. B., & Wilczek, F. 1983, Phys. Lett. B, 120, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56 [Google Scholar]
 RindlerDaller, T., & Shapiro, P. R. 2012, MNRAS, 422, 135 [NASA ADS] [CrossRef] [Google Scholar]
 RindlerDaller, T., & Shapiro, P. R. 2014, Mod. Phys. Lett. A, 29, 1430002 [NASA ADS] [CrossRef] [Google Scholar]
 Robles, V. H., Bullock, J. S., & BoylanKolchin, M. 2019, MNRAS, 483, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, K. K., & Peiris, H. V. 2021, Phys. Rev. Lett., 126, 071302 [NASA ADS] [CrossRef] [Google Scholar]
 Salucci, P., Wilkinson, M. I., Walker, M. G., et al. 2012, MNRAS, 420, 2034 [NASA ADS] [CrossRef] [Google Scholar]
 Schive, H.Y., Chiueh, T., & Broadhurst, T. 2014a, Nat. Phys., 10, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Schive, H.Y., Liao, M.H., Woo, T.P., et al. 2014b, Phys. Rev. Lett., 113, 261302 [NASA ADS] [CrossRef] [Google Scholar]
 Schwabe, B., Niemeyer, J. C., & Engels, J. F. 2016, Phys. Rev. D, 94, 043513 [Google Scholar]
 Sen, P. K. 1968, J. Am. Stat. Assoc., 63, 1379 [CrossRef] [Google Scholar]
 Shapiro, P. R., Dawoodbhoy, T., & RindlerDaller, T. 2021, MNRAS, 509, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Skodje, R. T., Rohrs, H. W., & VanBuskirk, J. 1989, Phys. Rev. A, 40, 2894 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, ApJ, 606, 702 [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
 Theil, H. 1950, Nederl. Akad. Wetensch. Proc., 53, 386 [Google Scholar]
 Toro, E. 2006, Appl. Numer. Math., 56, 1464 [Google Scholar]
 TrujilloGomez, S., Klypin, A., Primack, J., & Romanowsky, A. J. 2011, ApJ, 742, 16 [Google Scholar]
 Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2010, ApJS, 192, 9 [Google Scholar]
 Velten, H., & Wamba, E. 2012, Phys. Lett. B, 709, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Veltmaat, J., Niemeyer, J. C., & Schwabe, B. 2018, Phys. Rev. D, 98, 043509 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [Google Scholar]
 Weinberg, D. H., Bullock, J. S., Governato, F., Naray, R. K. D., & Peter, A. H. G. 2015, Proc. Natl. Acad. Sci., 112, 12249 [Google Scholar]
 Widrow, L. M., & Kaiser, N. 1993, ApJ, 416, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, X., Chan, M. H., Harko, T., Liang, S.D., & Leung, C. S. 2018, Eur. Phys. J. C, 78, 346 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Fitted scaling relation parameters
The fitted parameters for the scaling relations of r_{c}(M_{200}), δ_{c}(M_{200}), and M_{c}(M_{200}), Eqs. (44), (45), and (46), are shown in Tables A.1, A.2, and A.3. These values correspond to the fits shown in Figs. 6 and 13 obtained using TheilSen regression, including the results for the SPARC dataset (Li et al. 2020) and the Milky Way dSphs (Salucci et al. 2012).
Fitted scaling parameters for r_{c}(M_{200}), Eq. (44), at various redshifts. The results from the SPARC dataset and the Milky Way dSphs are included under the label “Data”.
Fitted scaling parameters for δ_{c}(M_{200}), Eq. (45), at various redshifts. The results from the SPARC dataset and the Milky Way dSphs are included under the label “Data”.
Fitted scaling parameters for M_{c}(M_{200}), Eq. (46), at various redshifts. The results from the SPARC dataset and the Milky Way dSphs are included under the label “Data”.
All Tables
Fitted scaling parameters for r_{c}(M_{200}), Eq. (44), at various redshifts. The results from the SPARC dataset and the Milky Way dSphs are included under the label “Data”.
Fitted scaling parameters for δ_{c}(M_{200}), Eq. (45), at various redshifts. The results from the SPARC dataset and the Milky Way dSphs are included under the label “Data”.
Fitted scaling parameters for M_{c}(M_{200}), Eq. (46), at various redshifts. The results from the SPARC dataset and the Milky Way dSphs are included under the label “Data”.
All Figures
Fig. 1. Selfinteraction energy density U_{SI}, thermal energy density U, and gravitational potential energy density W of the spherically symmetric test halo with R_{c} = 1 kpc, relative to the interaction energy of the initial background ${\overline{U}}_{\mathrm{SI}}$. The result of the simulation with refinement levels 8 to 12 are shown in solid, while 7 to 11 is shown in dotted. The 3D simulation reproduces the 1D results of Ahn & Shapiro (2005) and Dawoodbhoy et al. (2021), in particular the large thermal energy in the wake of the accretion shock, its sharp decrease to below the selfinteraction energy inside the core, and the fit U_{SI} ∝ ρ^{2} ∝ r^{−24/7} in the envelope. 

In the text 
Fig. 2. Energy conservation (top), the smallest resolved physical cell length (middle), and the runtime in CPU hours (bottom) for the main suite of simulations. Runs with a minimum refinement level of 8 are shown in solid lines, and 7 in dotted lines. 

In the text 
Fig. 3. Binned SIBECDM halo profiles and the HMF for R_{c} = 1 kpc at z = 1.5, with minimum refinement level 7 (dotted) and 8 (solid). The HMF of CDM for the same box size and initial resolution is included in black. The halo mass limit M_{min} are indicated by dots. The standard deviation from the binned halo profile mean for level 8, as well as the Poisson error in the HMF, are shown in shaded. 

In the text 
Fig. 4. Cumulative distribution functions for different cored halo profiles from SIBECDM simulations with R_{c} = 1 kpc (upper), R_{c} = 3 kpc (middle), and R_{c} = 10 kpc (lower) at z = 0.5. This shows the fraction of the halos that has ${\chi}_{v}^{2}$ smaller than a given value, and therefore how well the fitting function generally describes the halos in our simulation. 

In the text 
Fig. 5. Mean (solid) and standard deviation (shaded) of binned SIBECDM halo profiles with R_{c} = 1 kpc at z = 0.5. The binned profiles are fitted to NFWc (dashed) and Burkert (dashdotted). The fits of the NFW profile to the halo envelopes are also shown (dotted). 

In the text 
Fig. 6. Core radii r_{c} (left), core densities δ_{c} (middle), and core masses M_{c} (right) versus M_{200} for halos in cosmological simulations of SIBECDM with R_{c} = 1 kpc, 3 kpc, and 10 kpc, using the Burkert profile. The upper plots show the scatter of halos at z = 0.5, with the fitted scaling Eqs. (44)–(46) shown in dashed lines, compared to the SPARC dataset (Li et al. 2020) and Milky Way dSphs (Salucci et al. 2012). The lower plots show the fitted median (solid) and the first and third quartiles (shaded) at several redshifts using TheilSen regression. The results from SPARC and the dSphs are also shown. In the scatter plot for r_{c} the core radii R_{c} are indicated in dotted colored lines. 

In the text 
Fig. 7. Ratio of average thermal and selfinteraction energy, U and U_{SI}, inside the halo cores r < r_{c}/2 in simulations with R_{c} = 3 kpc at z = 1. 

In the text 
Fig. 8. Top panel: selfinteraction energy density U_{SI}, thermal energy density U, and gravitational potential energy density W of a sample halo with R_{c} = 1 kpc and M_{200} = 3 × 10^{9} M_{⊙}, relative to the interaction energy at the background level ${\overline{U}}_{\mathrm{SI}}$. The middle panel shows instead the specific energy, and the bottom panel the cumulative virial. The core radius r_{c} and halo radius r_{200} are indicated with the inner and outer dotted vertical lines, respectively. 

In the text 
Fig. 9. Final distributions of U/U_{SI} in 2D simulations with symmetrical (upper) and asymmetrical (lower) collapse, centered on the density peak. In the asymmetrical case, the second smaller overdensity was initially located to the right. The minimum in the symmetrical case is around U/U_{SI} ≈ 2 × 10^{−2}, while in the asymmetrical case its is U/U_{SI} ≈ 2. 

In the text 
Fig. 10. Projection plots of the DM density in cosmological simulations with L = 2 Mpc h^{−1} at z = 0.5 SIBECDM with R_{c} = 1 kpc (upper right) and CDM with the same initial conditions, but without a cutoff (upper left). The SIBCEDM internal energies for the same snapshot showing the effective thermal energy U (lower left) and the selfinteraction energy U_{SI} (lower right). 

In the text 
Fig. 11. Matter power spectrum for cosmological simulations L = 2 Mpc h^{−1} at z = 0.5 for SIBECDM with R = 1 kpc (solid) and CDM (dashed). The comoving cutoff k_{cut} is indicated with a dotted vertical line, and the spectrum at z = 50 is multiplied by 50 in the figure. 

In the text 
Fig. 12. Evolution of the ratio of the total thermal energy and total selfinteraction energy for the simulations run. 

In the text 
Fig. 13. Scaling functions at several redshifts for the core radii r_{c} (left), the core density δ_{c} (middle), and the core mass M_{c} (right) for SIBECDM halos with R_{c} = 3 kpc and different initial ζ = P/P_{SI} (upper), and different initial power cutoffs (lower). The shaded areas give the first and third quartiles for the scaling parameters obtained from the TheilSen regression. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.