Varying primordial state fractions in exo-and endothermic SIDM simulations of Milky Way-mass haloes

Self-interacting dark matter (SIDM) is increasingly studied as a potential solution to small-scale discrepancies between simulations of cold dark matter (CDM) and observations. We examine a physically motivated two-state SIDM model with both elastic and inelastic scatterings. In particular, endothermic, exothermic, and elastic scattering occur with equal probability at high relative velocities ( v rel ≳ 400 km / s ) . In a suite of cosmological zoom-in simulation of Milky Way-size haloes, we vary the primordial state fractions to understand the impact of inelastic dark matter self-interactions on halo structure and evolution. In particular, we test how the initial conditions impact the present-day properties of dark matter haloes. Depending on the primordial state fraction, scattering reactions will be dominated by either exothermic or endothermic effects for high and low initial excited state fractions respectively. We find that increasing the initial excited fraction reduces the mass of the main halo, as well as the number of subhaloes on all mass scales. The main haloes are cored, with lower inner densities and higher outer densities compared with CDM. Additionally, we find that the shape of the main halo becomes more spherical the higher the initial excited state fraction is. Finally, we show that the number of satellites steadily decreases with initial excited state fraction across all satellite masses.


INTRODUCTION
Dark matter (DM) forms about 85% of the mass budget of the Universe, but it has never been observed directly due to its minimal interaction with baryonic matter.Models for DM are often explored through cosmological simulations, which are then compared to observational constraints -for example, from rotation curves, strong or weak gravitational lensing, and the luminosity function of galaxies.The prevailing model is the cold DM plus cosmological constant (ΛCDM) model, where DM is assumed collisionless and has non-relativistic velocities since the early universe.Rather than committing to a description of DM on the particle physics scale, CDM is flexible and uses macroscopic properties that apply to a broad class of DM candidates.Weakly interacting massive particles (WIMPs) have been a particularly promising candidate for CDM since the freeze-out density naturally falls within the energy and satellite (Zavala et al. 2019;Hayashi et al. 2020) populations while simulated satellites are much more uniformly described by the NFW profile.
The missing satellites problem refers to the historical discrepancy where many more satellites were predicted in simulations than satellite galaxies in observations (Moore et al. 1999;Springel et al. 2008;Klypin et al. 2011;Garrison-Kimmel et al. 2014a).More recently, however, more faint dwarf galaxies have been observed that, combined with accounting for the completeness of surveys, bring the number counts for observed galaxies and simulated haloes into agreement (Kim et al. 2018;Drlica-Wagner et al. 2020).Finally, the too-big-to-fail problem refers to the scenario where there are more large DM haloes predicted than large galaxies seen, but the predicted haloes have a large enough gravitational potential that they should readily form galaxies (Boylan-Kolchin et al. 2011;Boyer et al. 2012;Garrison-Kimmel et al. 2014b).The potential should also be strong enough to resist losing the majority of the stars to stripping.
These discrepancies were primarily observed in DM-only simulations, where baryons are ignored.Because of the high dark-tobaryonic matter ratio, galaxies are expected to form alongside DM overdensities.However, despite making up only a fraction of the total mass, baryonic matter has a significant influence on galaxy formation.Feedback from supernovae or stellar winds, for example pushes, DM outwards (Navarro et al. 1996;Governato et al. 2012;Di Cintio et al. 2014a,b), which has produced cores in simulations across a range of galaxy sizes (Pontzen & Governato 2012;Keller et al. 2016;Burger et al. 2022).
Many prescriptions for baryonic processes rely on subgrid models, where the relevant physics occurs at scales less than the resolution of the simulation.The effects of these subgrid models are sensitive to model parameters and can vary significantly and impact key comparisons between simulations and observations (Vogelsberger et al. 2020).The hydrodynamics solver for gas, e.g. using the moving-mesh scheme AREPO or the smoothed particle hydrodynamics scheme GADGET, can change the size of disks that are produced (Torrey et al. 2012).The star formation models are also not well understood.Altering the density threshold for star formation in gas influences whether cores are formed in dwarf galaxies Benítez-Llambay et al. (2019), altering the primary factor that determines the star formation rate, e.g. a density threshold or an H 2based model, changes the number of satellites by a factor of two (Munshi et al. 2019), and altering the star formation efficiency affects the morphology of galaxies (Agertz & Kravtsov 2015).These parameters can thus be the difference between matching observations or producing unrealistic galaxies.
Although baryonic models may mitigate some of the smallscale discrepancies, it is difficult to resolve them all simultaneously.Since DM is also not well understood, there is significant interest in modifications to CDM that alleviate its small-scale problems while maintaining its large-scale accuracy (Cyr-Racine et al. 2016;Vogelsberger et al. 2016).In addition, it is also difficult to predict how effects from baryonic physics and alternative DM combine.Alternative DM models can alter star formation and gas distribution (Shen et al. 2024a) and the stellar density profile (Vogelsberger et al. 2014).
Early proposed DM candidates based on known physics are heavily constrained.Warm DM (WDM), for instance, that may be made up of particles like neutrinos (Dodelson & Widrow 1994;Shi & Fuller 1999), are strongly constrained by the Lyman-α forest (Viel et al. 2013;Schneider et al. 2014;Iršič et al. 2017;Hooper et al. 2022).Weakly interacting massive particles (WIMPs), while promising due to their interaction scale, have become increasingly limited in their allowed cross sections as both direct and indirect search efforts have yet to detect them (Marrodán Undagoitia & Rauch 2016;Strigari 2013).Like with baryonic models, it is difficult to simultaneously solve all the small-scale problems (Colín et al. 2008;Polisensky & Ricotti 2015;Lovell et al. 2017).
This has led to the increasing interest in self-interacting DM (SIDM) (Spergel & Steinhardt 2000), which produces distinct direct detection signals (Vogelsberger & Zavala 2013).This is a large class of particles that allow for self-interactions between DM particles but not between DM and baryons beyond gravity.With small interaction cross sections, SIDM does not alter the large-scale structure of the universe (Rocha et al. 2013;Vogelsberger et al. 2016;Sameie et al. 2019;Andrade et al. 2022) but can alleviate small-scale problems by transferring orbital energy between particles (Tulin & Yu 2018, and references therein), undergoing corecollapse (Kaplinghat et al. 2019), causing variations in galactic rotation curves (Creasey et al. 2017), and combining with tidal forces (Dooley et al. 2016;Sameie et al. 2020).
Additionally, SIDM scattering events can be elastic or inelastic.While many models initially explored elastic SIDM (Vogelsberger & Zavala 2013;Rocha et al. 2013;Zavala et al. 2013;Elbert et al. 2015;Andrade et al. 2022), models including inelastic collisions are rising in interest due to their ability to alter the energy within the DM particles (Medvedev 2010b(Medvedev ,a, 2014a,b;,b;Schutz & Slatyer 2015;Blennow et al. 2017;Zhang 2017;Vogelsberger et al. 2019;Shen et al. 2021Shen et al. , 2024b;;O'Neil et al. 2023) and have important cosmological consequences (e.g.Fan et al. 2013;Foot 2014;Vogelsberger et al. 2019;Xiao et al. 2021;Roy et al. 2023).One way to provide the mechanism for energy transfer during inelastic reactions is for the DM particles to exist in multiple states (e.g.Schutz & Slatyer 2015;Vogelsberger et al. 2019;O'Neil et al. 2023).Other mechanisms, like decay (e.g.Wang et al. 2014), could have similar effects, but the multi-state model provides a large amount of flexibility and is what we focus on in this paper.The energy difference between the states can then be exchanged with the kinetic energy of the particles.
The parameter space in the multistate model is then defined by the ground state mass (m χ 1 ), the energy difference between the two states (δ ), the dark force mediator mass (m φ ), and the coupling between the mediator and the DM particle (α) (Schutz & Slatyer 2015).By adjusting these parameters, a wide range of SIDM behaviours can be achieved, including elastic, exothermic, and endothermic reactions (Vogelsberger et al. 2019;O'Neil et al. 2023).It is not currently feasible to fully explore this parameter space with simulations, so we take representative models that emphasise different reactions.
In addition to the parameters of the SIDM model, it is also necessary to determine the initial abundance of states in the early universe, which depends on the cosmological history of the DM model.Previous work (e.g.Vogelsberger et al. 2019;O'Neil et al. 2023) has been restricted to initialising the DM particles entirely in one state or split evenly between the two states.The initial abundance of each state has marked effects on the resulting SIDM behaviour.If the particles are initially all in the ground state, for example, a model that is primarily exothermic will be indistinguishable from CDM. Generically, it is expected that the ground state is energetically favoured by the high-density environment of the early universe (Medvedev 2001(Medvedev , 2010b(Medvedev , 2014b)).Depending on the processes that determine the initial abundances, for example freeze-in or resonant annihilation, it is possible to obtain varying amounts of DM in the ground and excited states (Brahma et al. 2024).
In this paper, we expand on the work done in O'Neil et al. (2023).Using the described multi-state SIDM framework, they introduced a model where elastic, exothermic, and endothermic scattering events had similar cross sections and therefore all contributed to the SIDM effects in the halo.They then ran this model, along with several comparison models, in DM-only, Milky Way-sized zoom-in simulations.In isolation, elastic scattering will produce cores in a halo, exothermic scattering will push particles to larger radii and make it more difficult for satellites to form, and endothermic scattering will lead to high central densities in both Milky Waymass and satellite haloes.Thus, in isolation, endothermic scattering is expected to exacerbate small-scale DM discrepancies, but the non-linear interplay between the different reaction types makes it difficult to predict how DM haloes will behave when all reactions are able to occur.They compared this model to the model in Vogelsberger et al. (2019), which suppressed the endothermic reactions, and showed that viable DM haloes could be produced even with a significant endothermic cross section.
Here, we explore the effects of altering the initial abundance of DM particle states on the introduced endothermic SIDM model.O'Neil et al. (2023) explored only the case in which all particles started in the ground state.We study the effects of initial fractional abundance of each state between 0 and 1 in a suite of eleven Milky Way-mass zoom-in simulations using the same initial conditions as in O'Neil et al. (2023).Although not each of these initial abundances is equally likely, it is necessary to explore the intermediate abundances to understand how the halo properties transition.
The paper is organised as follows.In Section 2, we describe our simulations (2.1) and the SIDM model (2.2).In Section 3, we show the results of varying the initial abundances on properties of DM haloes: density profiles (3.1), abundances of each state at z = 0 (3.2), halo shape (3.3), velocity dispersion and anisotropy (3.4), and substructure (3.5).We also show the redshift evolution of the haloes (3.6).Finally, we summarise our conclusions in Section 4.

Simulations
The analysis in this paper is based on a suite of 11 dark matter-only (DMO) cosmological zoom-in simulations of Milky Way-mass haloes.The simulations are performed using the multi-physics, moving-mesh, hydrodynamic code AREPO (Springel 2010;Weinberger et al. 2020), which employs the tree-particle-mesh (Tree-PM) algorithm for gravity with periodic boundary conditions.The 10 6 10 7 10 8 10 9 10 10 10 11 10 12 10 13 10 14 Halo mass [M ] 10 0 10 1 10 2 10 3 The transfer cross sections as a function of relative particle velocity for each scattering reaction in the endothermic SIDM model.The top axis shows the halo M 200,mean with circular velocity v at R 200,mean .Elastic scattering between same state particles is very unlikely, while the elastic scattering between different state particles is similar in likelihood to downand up-scattering.Down-scattering is most likely at lower relative velocities, and there is a velocity threshold for up-scattering to ensure that particles have enough kinetic energy to transform to the energy difference between the ground and excited states.
main target haloes in the zoom-in region have M 200,mean ≈ 2 × 10 12 M ⊙ at z = 0, where M 200,mean is the mass within the radius where the average density is 200 times the mean density of the universe.The simulations use the following cosmological parameters: matter density fraction Ω m = 0.302, dark energy fraction Ω λ = 0.698, and Hubble constant H 0 = 100h kms −1 Mpc −1 where h = 0.691.The main halo and its surroundings are made up of 55,717,200 high-resolution DM particles.The high-resolution region is surrounded by particles of mass 1.67 × 10 6 M ⊙ and 1.13 × 10 8 M ⊙ farther away from the main halo.For the high-resolution region, the effective Plummer equivalent gravitational softening length of DM is ε = 0.153 kpc and the DM particle mass resolution is m zoom dm = 2.2 × 10 5 M ⊙ .The background particles get less resolved the farther they are from the main halo in particle masses in the range m back dm = 1.76 × 10 6 − 9.03 × 10 8 M ⊙ .The gravitational softening lengths for these particles are adjusted by a factor of (m back dm /m zoom dm ) 1/3 .The SIDM model is not implemented for the background particles, which behave as CDM.
The periodic simulation box has a side length of 100 Mpc/h.The initial conditions are generated at redshift z = 127 using MU-SIC (Hahn & Abel 2011).All the simulations in our suite were run using the same initial conditions.
Haloes are identified using a Friends-of-Friends (FoF) algorithm (Davis et al. 1985).This spatially identifies haloes by linking particles and substructures when they lie with a linking length of b (V /N) 1/3 .We use a value of b = 0.2, V is the volume of the box and N is the number of particles in the box.Gravitationally bound subhaloes are identified using SUBFIND (Springel et al. 2001;Dolag et al. 2009).The position of a subhalo is determined by the position of its most bound member particle, and the most massive subhalo in a FoF group is identified as the main halo.

SIDM model
In this paper, we explore a DM model where the DM particles are able to self-interact with small cross-sections in addition to gravitational interactions.The particles scatter only a few times during the formation of the halo, so the DM remains largely collisionless.We employ the endothermic model introduced in O'Neil et al. ( 2023) based on the theory described in Section 3.1 ofSchutz & Slatyer (2015).The code is the same as used in Vogelsberger et al. (2019), which modelled a similar two-state model with different cross sections.In this framework, DM particles have two states, and the energy difference between the two states is much less than the groundstate energy.The model is fully described with four parameters: the ground state mass m χ 1 , the mass splitting between the ground and excited states δ , the dark force mediator mass m φ , and the dark coupling constant α.The particles can scatter elastically, where there is no state change, or inelastically, where the particles change states.There are five ways for particles to scatter in our model: (i) Elastic scattering of two ground state particles: In this reaction, two ground-state particles scatter and remain in the ground state.The total kinetic energy remains the same, but energy can be transferred between the particles and therefore through the halo.Elastic scatterings tend to produce a thermalised core in DM haloes.
(ii) Inelastic up-scattering of two ground state particles: In this reaction, two ground state particles scatter to produce two excited state particles.In this case, kinetic energy from the two particles is transformed to the energy of excited states, which results in "dark cooling flow" towards the centre of the halo.
(iii) Elastic scattering of two excited state particles: This reaction is similar to the elastic scattering between ground state particles but occurs between excited state particles.In this model, the cross-section for these two elastic scattering reactions is the same.This model also contributes energy transfer between particles and core creation.
(iv) Inelastic down-scattering of two excited state particles: In this reaction, two excited state particles scatter to produce two ground state particles.The energy difference between the states is converted to the kinetic energy of DM particle.This can cause the particles to move farther out in the DM halo, or to become unbound entirely.
(v) Elastic scattering of a ground and an excited state particle: Similar to the other elastic scattering reactions, there is no state change effectively.The cross-section for this reaction is calculated using the analytic cross-sections derived in Schutz & Slatyer (2015) in the Born regime (αm χ 1 ≪ m φ ).To account for the heavy mediator in this model, the calculation is supplemented with expressions from Feng et al. (2010).
The values for the model we adopt in this work are m χ 1 = 2.3 MeV, δ = 0.48 eV, m φ = 7.3 MeV, and α = 0.17.These values lead to a significant probability of both up-and down-scattering in addition to elastic scattering between ground and excited state particles.Elastic scattering between either two ground state or two excited state particles is suppressed though does happen in small amounts.While there have been some constraints placed on elastic scattering (see Adhikari et al. 2022, for review), inelastic models are less constrained.Additionally, these cross sections are high enough to be on the extreme end of what might be seen, both making the influences clear and testing the limits of the model.
Figure 1 shows the transfer cross section per unit mass (σ T /m χ ) for each type of scattering.The transfer cross sections are a function of relative velocity between the two particles.The conversion between halo mass and relative velocity is done by assuming an isotropic Maxwell-Boltzmann velocity distribution of DM particles with one-dimensional velocity dispersion σ 1d ≃ 0.6V vir (Łokas & Mamon 2001), where V vir ≡ GM 200,mean /R 200,mean is the virial velocity of the halo.The elastic scattering between ground state particles (purple line) coincides with the cross section for elastic scattering between excited state particles (green dashed line) at a value of 0.21 cm 2 g −1 .The upscattering reaction (orange line) is only possible when the relative velocity is high enough such that there is sufficient kinetic energy to convert to the energy difference between the ground and excited states.At lower velocities, the down-scattering reaction (red line) becomes more likely, scaling with a power law of -1, which can give particles more kinetic energy.Elastic scattering between a ground and an excited state particle is approximately constant at a value of 15 cm 2 g −1 .
Within the simulation, individual DM particles are represented with simulation macroscopic particles with a mass as given by the resolution (i.e.2.2 × 10 5 M ⊙ for zoom particles).Each DM particle represents an ensemble of fundamental DM particles that are entirely in either the ground or excited state and all of which scatter simultaneously.These particles scatter less frequently than fundamental DM particles such that over the course of the simulation, the expected amount of DM mass experiences each scattering reaction.
We follow the algorithm described in Schutz & Slatyer (2015) and used in Vogelsberger et al. (2019) to implement the SIDM model in the simulations, which we briefly summarise here.Whether a particle scatters or not, and which scattering occurs, depends on both the scattering transfer cross sections and the local density of dark matter.The local density is determined by the 10 to 38 nearest neighbours smoothed by a cubic spline kernel function around its nearest neighbours.The cross-section at the relevant velocity between two particles is multiplied by the DM (simulation macroscopic) particle mass and the distance each particle would move in one time step to get the probability of scattering.The probabilities are appropriately adjusted by the kernel factor as determined by the distance between the two potential scattering particles.Whether a particle scatters with a nearby particle and which reaction occurs is determined by drawing a random number.
In the event that a scattering reaction occurs, an appropriate second particle is selected to scatter with.The timesteps are set to be smaller than the expected scattering frequency so that multiple scatters per particle per time step are rare.In the case that there are multiple scatters in a single time step, the second reaction is rejected and does not take place.This is because a particle's position and velocity can be updated only once per time step.Rejected scatters make up less than 1% of the total expected scatters, and further decreasing the time step to avoid these does not impact the results.We have tested this by running the same initial conditions with larger cross sections and varying the time step and rejection rates under 5 percent do not noticeable alter the particle distribution in the halo.The particles scatter isotropically, and any changes to the particles' energies are done in the centre of momentum frame of the two particles.
We run a suite of 11 simulations with varying initial excited and ground state fractions between 0% and 100%.

RESULTS
Our suite is made up of 11 simulations with the DM particles initially 0 − 100%, increasing in intervals of 10%, in the excited state.We refer to the initial excited state fraction as χ 2 init and the initial ground state fraction as χ 1 init .In Figure 2, we show the DM projection for our fiducial CDM simulation.In Figure 3, we show DM projections for all particles several simulations spanning the range of χ 2 init .We show the DM projections for excited and ground state particles separately in Figure 4.Each simulation has a similar virial radius, ranging from 370.58 kpc for χ 2 init = 100% to 401.27 kpc for CDM.

Effects on halo density profiles
The mass distribution within the haloes varies significantly.In Figure 5, we plot density profiles for the main halo at z = 0 in the top panel.The SIDM model gives a cored inner density profile for all initial excited fractions, compared to the cuspy profile of CDM.This is similar to the effects of classical elastic SIDM models previously explored in literature (e.g.Zavala et al. 2013;Rocha et al. 2013;Elbert et al. 2015) and is due to thermalisation of DM cores from heat conduction.
The two-state configuration of the SIDM model studied in this paper enables additional channels of heating or cooling in the dark sector.As shown in Figure 5, when χ 2 init is low, we find that the central density of the simulated Milky Way-mass halo is enhanced compared to CDM.This is similar to the gravothermal core-collapse studied in lower halo mass regimes and often in SIDM models with strong velocity-dependence of cross sections (e.g.Nadler et al. 2020;Sameie et al. 2020;Correa et al. 2022;Zeng et al. 2022).The up-scatterings of DM particles dissipate away their kinetic energy and can accelerate the gravothermal collapse of the halo (e.g.Essig et al. 2019).This is a common phenomenon found in previous explorations of dissipative SIDM (e.g.Essig et al. 2019;Shen et al. 2021;Roy et al. 2023).For example, Shen et al. ( 2021) simulated a DM model that dissipates a constant fraction of a particle's kinetic energy during an interaction.They ran a suite of zoomin simulations of dwarf to Milky Way-mass galaxies with different cross sections in addition to a fiducial CDM model, including the baryonic physics described by the FIRE-2 model (Hopkins et al. 2018).They found that the halo central density increases with cross section and effectively the energy dissipation rate of DM.Similar to the up-scattering in our work, the dissipation causes a particle to lose velocity and fall inwards, thus increasing the central density.
However, increasing χ 2 init leads to a lower core density and higher outer density.The half mass radius r 1/2 increases significantly with χ 2 init , while the virial radius r 200,mean decreases slightly.As shown in the top panel of Figure 5, simulations with higher initial excited state fraction χ 2 init have lower inner densities, varying overall by two orders of magnitude.Exothermic reactions from down-scatterings dominate in these cases and decrease the central density, as previously explored in N-body zoom-in simulations (Vogelsberger et al. 2019).An initially larger fraction of excited particles eventually provides heating of DM particles through downscattering, pushing more particles away from the halo center.This also leads to higher densities at R 200 for simulations with higher χ 2 init , although the difference is less dramatic.The velocity boost also causes virial mass (M 200 ) to decrease with χ 2 init as more particles are driven above escape velocity as shown in Figure 6.
The gradients of the velocity dispersions in the SIDM models vanish at lower radii while the CDM peaks at about 10 kpc.Because heat can transfer between DM particles in a self-interacting model, the velocity dispersion smooths out, The equilibrium velocity dispersion decreases with χ 2 init , and it is higher than the CDM case for low χ 2 init .This implies that there is heat transferred out of the halo when particles start in the excited state.Initial downscattering pushes particles outwards and out of the halo, leaving the halo at a slightly lower mass and with less energy.In the low χ 2 init simulations, on the other hand, up-scattering earlier in the halo's evolution creates a deep gravitational well that binds particles within the halo.The mass of the χ 2 init = 0% is 2.16 × 10 12 M ⊙ and χ 2 init = 100% is 1.70 M ⊙ , and the mass decreases monotonically between the two.
The outwards transfer of energy from down-scattering in the high χ 2 init also increases the core size while the earlier up-scattering in low χ 2 init increases the core density value.When a particle downscatters, it increases its kinetic energy by the energy difference between the ground and excited state δ = 10 keV, so the radius in the gravitational well increases by which comes to ∆r ≈ 0.4 kpc.
In Figure 6 we show M 200,mean and V max for the main halo in

All particles
Figure 3.The projected DM density for the main haloes at z = 0 in the simulations with χ 2 init = 0, 20, 40, 60, 80 or 100% of the particles starting in the excited state.As more particles begin in the excited state, more exothermic reactions occur and the substructure is destroyed.When more particles begin in the ground state, the endothermic reaction plays an important role in the appearance of the main halo at low redshift, and substructure remains.We also note that the halos with low excited state initial fractions are more spherical while the higher initial state fraction simulations have a more oblong shape.each simulation.V max is the maximum circular velocity of the halo.As a result of dominant upscatterings, SIDM models with low χ 2 init values produce more concentrated DM distributions and thus give larger V max values.On the contrary, models with large χ 2 init values feature substantial downscatterings of DM particles that create cored density profiles.The cores not only decrease the V max values in these models but also lead to lower halo virial radius and masses (given the same spherical overdensity threshold).

Excited state fraction
Figure 7 shows the fraction of particles in the excited state at redshift 0 χ 2 (z = 0) as a function of radius r.At small radii, χ 2 (z = 0) shows the opposite trend with respect to χ 2 init , i.e. models with higher initial excited state fraction end up with lower excited state fraction in halo centers at z = 0.This is likely due to the different central velocity dispersions of DM particles reached in different models, which results in different equilibrium states.To quantitatively understand this phenomenon, we can approximate the local equilibrium excited state fraction by equating the upscattering and downscattering rates.The expected time for one ground-state particle to upscatter once is where ⟨...⟩ denotes the thermal average, ρ gr (ρ ex ) denotes the density of ground (excited) state DM particles, and σ up is the upscattering transfer cross section.For simplicity, we assume that the velocities of DM particles locally obey the Maxwell-Boltzmann distribution and we obtain where σ 1d is the local one-dimensional velocity dispersion of DM.Therefore, the expected up-scattering rate in a unit volume can be expressed as And the same rate for downscattering is A local equilibrium occurs when Γ up = Γ down and we obtain the equilibrium excited state fraction as where σ 1d,gr and σ 1d,ex are one-dimensional velocity dispersions of ground and excited state particles, respectively.For simplicity, we assume that both of them are the same as the total DM velocity dispersion σ 1d .On the right y-axis of Figure 5, we show the predicted equilibrium χ 2 as a function of σ 1d , assuming that excited and ground state particles have the same velocity dispersion.At asymptotically high velocities where the probability for up-and down-scattering are consistently similar, we should get 50% of the particles in the excited state.For the velocity dispersion of our simulated haloes, which is σ 3D ∼ 200 − 300 km s −1 so the one-dimensional velocity dispersion is σ 1d ∼ 100 − 150 km s −1 , corresponding to an excited state mass ratio of 0.1 − 0.3, which is consistent with the excited state fractions at lower radii in Figure 7.
On the other hand, the excited state fraction at the outskirt (r ≳ 100 kpc) of the halo in Figure 7 follows a qualitatively dif-ferent trend.χ 2 (z = 0) increases monotonically with χ 2 init and a crossing point of the χ 2 (z = 0) exists around 0.1 − 0.2 R 200,mean .Since the density and velocity dispersion of DM at large radii are low enough such that the up-scattering channel is completely suppressed (see Figure 1), the equilibrium solution cannot be realized.The excited state fraction will simply move monotonically downwards with time due to downscattering, resulting in the trend seen at large radii.The cosmological evolution of excited particle number density when the upscattering channel is turned off would be with n ex (z i ) = χ 2 n dm (z i ) where n dm is the number density of all DM particles.We take z i to be the starting redshift z i = 127 of the simulations.As shown in Figure 1, in the lowvelocity end, the product of particle relative velocities and selfinteraction cross-sections tend to have a stabilized value around 6 × 10 3 cm 2 g −1 km s −1 .Integrating Equation 7 and comparing to the n dm (z = 0), we obtain χ 2 ∼ 0.1, 0.3, 0.4 of excited particles at z = 0 in models with χ 2 init = 0.1, 0.5, 0.9, respectively.This agrees with the excited fraction we found at the outskirts of the simulated .Density profiles and velocity dispersions for the main halo at z=0.Simulations with fewer initially excited particles tend to have denser centers.When all of the particles begin in the ground state, every downscatter must be preceded by an upscatter, so there is more upscattering than downscattering.The right-hand y axis shows the equilibrium excited state fraction, χ 2 equil , at relevant velocity dispersions, using σ 3d = √ 3σ 1d .
haloes.Such a picture also predicts that the excited state fraction in the low-density region of the Universe decreases much faster at earlier times due to higher particle number densities and stabilizes to ≲ 2% level from z = 6 to z = 0, which is consistent with our findings in Section 3.6.

Halo shapes
Here we discuss how the morphology of the halo as a function of χ 2 init .As seen in Figure 3, there is a noticeable difference in the shapes of the haloes as the initial excited state fraction varies.The haloes with a lower χ 2 init are more concentrated and spherical than those with higher χ 2 init .Figure 8 demonstrates this more quantitatively.For each halo, we fit an ellipse at different radii from 0.03 R 200,mean to R 200,mean and divide the middle and minor axes by the major axis to measure how spherical it is.We follow the method used in Chua et al. (2021), which iteratively fits an ellipsoidal shell to the halo while keeping the semi-major axis constant.First, we fix the semi-major axis, R. At each step, the bin is set to be within 0.1 dex of the current ellipse.  .The fraction of main halo particles in the excited state as a function of radius at z = 0.In the outer region (r ≳ 100 kpc), χ 2 increases with χ 2 init .This is due to the lower density at these radii, which decreases the impact of scattering.At smaller radii, χ 2 is not monotonic over χ 2 init , but simulations with higher χ 2 init counterintuitively tend to have lower χ 2 (z = 0).This is likely related to the lower inner velocity dispersion in these simulations, which causes down-scattering to dominate. .We fit an ellipse to each halo at several radii and find the ratio of the middle axis (q, top panel) and minor axis (s, bottom panel) to the major axis.A ratio of s = 1 is spherical, and a ratio of s = 0 is flat.The grey points show the CDM run.Each column of points is the axis ratio at that radius for each simulation halo.We see that low χ 2 init haloes tend to be spherical due to early up-scattering that compresses the halo.High χ 2 init haloes maintain a more oblong shape since down-scattering keeps the halo enlarged.
We begin with a spherical bin, defined by A i j = δ i j .We calculate the shape tensor of the particles inside the bin, defined as Next, we obtain the eigenvectors v 1 , v 2 , v 3 with eigenvalues λ 1 , λ 2 , λ 3 respectively (λ 1 ≥ λ 2 ≥ λ 3 ).These determine the intermediate axis ratios q ≡ λ 2 /λ 1 and s ≡ λ 3 /λ 1 .We also rotate the particles to a coordinate system with axes aligning with the eigenvectors.Finally, we generate a new bin with the new eigenvalues determining q and s.The iteration terminates when q and s change by less than 0.1% between steps.For more details, see Zemp et al. (2011).A ratio of s = 1 is perfectly spherical, and smaller ratios are increasingly oblong.At each radius, we plot q and s for each χ 2 init in different colours.Across all radii, haloes with a lower χ 2 init are more spherical than haloes with a higher χ 2 init .This is consistent with previous results in O'Neil et al. (2023), where the endothermic model resulted in a more spherical halo than the exothermic model.Here, higher χ 2 init haloes are dominated by exothermic reactions since the downscattering reaction is the only available reaction when all particles are in the excited state.In elastic scattering, transfer cross sections of approximately σ T /m χ = 0.1 cm 2 g −1 drive spherical centres in the cores of galaxy clusters (e.g.Miralda-Escudé 2002;Peter et al. 2013;Brinckmann et al. 2018;Shen et al. 2022).Here, we show that significant endothermic and elastic scattering produces similar results, but the exothermic scattering mitigates this, allowing larger elastic transfer cross sections while maintaining a certain amount of asymmetry.To fully understand the constraints in galaxy clusters, a dedicated set of simulations is needed to accurately predict how this model would behave in such environments.
Haloes do not accrete symmetrically, so it is expected that there is a certain amount of asymmetry as seen in the high χ 2 init simulations.This asymmetry is wiped away by the endothermic reactions.Endothermic reactions cause particles to fall deeper into the potential well by reducing their kinetic energy during the reaction.This compresses the inner regions of the halo and causes it to become more spherical, especially when a large amount of particles are up-scattered.This is consistent with findings that low-concentration halos are generally less symmetrical than highconcentration halos (Hashimoto et al. 2007) Even when these particles are later down-scattered and kicked further out in the halo, they cannot regain the asymmetrical orbits they previously may have had and maintain the spherical shape of the halo.Chua et al. (2019) studied halo shapes in Illustris and Illustris-Dark using methods similar to ours and focused on the difference between results with and without radiation.In CDM haloes in Nbody simulations, they found, similar to our results, that the halo is least spherical in the centre, with a median axis ratio of 0.6.They also found that radiation increases how spherical the halo is.We find that SIDM, both up-and down-scattering dominant, increase the how spherical a halo is.However, unlike radiation, we find that SIDM alters the shape of the axis ratio as a function of radius.
In particular, our SIDM haloes become much more spherical in the centre compared to the CDM.This is likely due to the up-scattering reaction in our SIDM model.Even in the high χ 2 init haloes, particles that have up-scattered are more likely to be found closer to the centre of the halo.Up-scattering is also the reaction that increases the axis ratio, which leads to our trend of increasing axis ratio towards the centre of the haloes.

Kinematic properties
The bottom panel of Figure 5 shows the velocity dispersion for all main halo particles as a function of radius.The velocity dispersion indicates the relative velocities between nearby particles, which determines the frequency of each type of scattering.Simulations with higher χ 2 init tend to have lower velocity dispersions in the inner regions, which corresponds to a higher frequency of downscattering and a lower frequency of upscattering.
In Figure 9, we show the anisotropy parameter of the particle velocities in each simulation's main halo.To calcu- late the anisotropy, we first calculate the radial and two components of angular velocity.We then take the variance in the radial component to calculate σ r and then sum the variance of the two angular components to get the tangential variance σ t .
The CDM curve follows the anistropy as a function of radius that is typically found (e.g.Navarro et al. 2010).Once selfinteractions are introduced, the orbits become much more radial.As inelastic reactions occur, the total velocity changes and the particles are scattered isotropically.However, because of the change in kinetic energy, the particles must also move either inwards or outwards in the gravitational potential.This causes the orbits to become more radial compared to any tangential change the particles experience.Hansen & Moore (2006) noted that there is a strong correlation between β and the slope of the density profile, and (Navarro et al. 2010) found that this holds primarily for the inner regions of the halo and not the outer regions.In particular, they find that β and d log ρ d log r are negatively correlated, with higher slopes, i.e. more cored, corresponding to lower β .In Figure 9, we find that the SIDM models all have significantly lower β compared to the CDM model, corresponding to a higher slope that is closer to zero.Although the densities of the SIDM models vary, they all flatten to a core in the centre and their anisotropy values are correspondingly similar at low radii.At larger radii, the SIDM anistropy peaks as the CDM does, though the peak is much wider.This is likely due to the energy transfer that is allowed between SIDM particles, which causes the velocities to be more uniform.The anisotropy is also generally lower at large radii for lower χ 2 init .

Effects on substructures
Figure 10 shows the cumulative number of subhaloes with mass M sub > 10 7 M ⊙ within a radius r to the main halo centre.In the bottom panel, we show the same quantity divided by that of the The number of subhaloes with mass > 10 7 M ⊙ within a radius r of the main halo at z = 0.As χ 2 init increases, the number of haloes decreases, although the shape of the curve remains the same.This is due to down-scattering that gives a velocity kick to particles and makes it more difficult to keep them within a gravitational well.Even at χ 2 init = 0, we see a decrease in substructure due to tidal forces.Bottom: The fractional difference in the subhalo number profile compared to CDM.The χ 2 init = 0 halo has approximately half as many subhaloes at each radius as the CDM halo, and this fraction decreases as χ 2 init increases.
CDM simulation.In models with increasing χ 2 init , there are significantly less substructures.The down-scatterings increase the kinematic energy of the subhalo and make it more difficult for subhalo with small gravitational wells to stay bound.On the other hand, up-scattering decreases the velocity but the velocity dispersion in satellite haloes is generally not high enough to significantly alter the substructure in up-scattering dominant models.The simulations with higher χ 2 init therefore have a much more altered substructure compared to CDM compared to simulations with a higher χ 1 init .However, there is still a decrease by a factor of 2 in even the χ 2 init = 0 case since the main halo increases in central density, which increases the tidal forces and strips the subhaloes.The tidal radius r t of the subhalo decreases as the enclosed mass of the main halo increases where m sub is the mass of the subhalo, M main is the mass of the main halo and r is the distance of the subhalo from the centre of the main halo (King 1962;Pace et al. 2022).Thus, an increased central density of 10 will result in approximately a factor of 2 decrease in a satellite's tidal radius.
Another important effect is the SIDM ram pressure stripping due to self-interactions between satellite and host halo DM parti-cles (e.g.Kummer et al. 2018;Nadler et al. 2020;Jiang et al. 2021), also known as SIDM evaporation (Zeng et al. 2022).Thermalized cores generated by self-interactions can make satellites more prone to this.The relative importance of ram pressure stripping can vary when the cross-section has velocity dependence (e.g.Banerjee et al. 2020;Nadler et al. 2020;Zeng et al. 2022).In our model, the cross-section of the dominant down-scattering process decreases by orders of magnitude with increasing velocities until reaching ∼ 1000 km s −1 , which is larger by a factor of few than the maximum circular velocities of the simulated Milky Way-mass halos.The heat generated from scatterings of satellite particles dominates the heat injected from scattering host halo particles.
Figure 11 shows the cumulative subhalo mass function for the main halo in each solution (top) and the same quantity divided by that of CDM (bottom).We find that simulations with higher χ 2 init tend to have fewer subhaloes at all mass scales.We can do a simple estimate of the typical subhalo mass where the heat generated from down-scattering exceeds the gravitational binding energy.For a halo with an NFW profile, the gravitational binding energy is (Mo et al. 1998) where c h is the concentration of the subhalo, and Meanwhile, the total energy of excited state particles in the subhalo E ex is This represents the maximum possible heating energy to the subhalo through exothermic self-interactions.Therefore, the smallest substructure we expect to survive this heat injection is This is generally above the subhalo masses in simulations and, therefore, the numbers of subhaloes at all mass scales display significant suppression.We note that this estimation assumes all of the energy in excited states is transformed into the kinetic energy of DM particles and should be interpreted as a limiting case.It also assumes that the down-scattering of DM particles within the satellite is the major channel of heat generation, which may not hold for other types of SIDM models or more massive host halos.
In observations of satellite galaxies around the Milky Way and Andromeda (e.g.McConnachie 2012;Lovell & Zavala 2023), galaxies are measured with dynamical masses well below 10 7 M ⊙ .In addition, McConnachie (2012) measures over 60 galaxies with masses larger than 10 7 M ⊙ .Tentatively speaking, these observations exclude our model with χ 2 init ≳ 80%.init due to downscattering kicking particles out of the gravitational wells.At high masses, this effect is less pronounced since the gravitational well is deeper and it is more difficult to escape it.At low χ 2 init , tidal stripping plays a larger role in decreasing the subhalo number than scattering effects.Bottom: The fractional difference in subhalo number compared to the CDM halo.The χ 2 init = 0 halo has approximately half as many subhaloes at each mass as the CDM halo, and this fraction decreases as χ 2 init increases.

Redshift evolution
decreasing in χ 2 and low χ 2 init haloes increasing in χ 2 .The lines do not cross each other, indicating a smooth transition between dominant exothermic and endothermic reactions.The final equilibrium state for the haloes is set by the velocity dispersion and the density of particles as described in Equation 7. Thus, although down-scattering is more likely than up-scattering for the higher χ 2 init haloes with lower velocity dispersions, and lower density, particularly on the outskirts of the haloes, prevents significant scattering of either kind.Therefore, the equilibrium χ 2 remains higher for haloes with higher χ 2 init .Simulations with lower χ 2 init see a jump in up-scattering near z = 2, corresponding to a major halo merger at that time.The merger increases the density and the velocity dispersion of particles in the main halo.The increase in velocity dispersion is particularly important, especially given that we do not see a jump in the higher χ 2 init haloes where down-scattering is more dominant.This is because down-scattering occurs without a velocity threshold so does not need the merger to provide a velocity boost.

CONCLUSIONS
We introduce a suite of zoom-in simulation using a SIDM model that incorporates a mix of elastic and inelastic scattering.The model consists of a two state DM particle, and a particle can move between states during inelastic scatters by exchanging mass energy with kinetic energy.In particular, the model includes an endothermic reaction that is just as likely as exothermic and elastic reactions at relative velocities above 400 km/s.
Past work showed that this new model resulted in cores, with the core density set by the onset time of endothermic reactions (O'Neil et al. 2023).Here, we investigated the effects of the initial conditions, which heavily influence this onset time.We ran a suite of 11 SIDM simulations of the same zoom-in halo varying the initial fractions of ground and excited state particles, with the initial excited state fraction χ 2 init between 0-100% in increments of 10% along with a fiducial CDM model.We summarise our results as follows: • We show in Figure 5 the density profile and velocity dispersion as a function of radius for each of our models.The SIDM model consistently produces a core in the density profile.The size of this core increases with χ 2 init , corresponding to more downscattering reactions.Similarly, the density of the inner regions decreases with χ 2 init and the velocity dispersion also flattens out at low radii.The central density changes by a factor of 40 between the χ 2 init = 0 and the χ 2 init = 1 simulations.• Figure 6 shows the maximum circular velocity as a function of M 200,mean and χ 2 init .V max and M 200,mean steadily decrease as a function of χ 2 init .Lower χ 2 init haloes are more concentrated due to increased up-scattering.
• For each simulation, we fit an ellipse at various radii and compare the major, middle, and minor axes to estimate the shapes of the haloes.Figure 8 shows the middle-to-major and minor-to-major axis ratios (q and s respectively)at several radii for each halo, so a ratio of s = 1 is circular and q = 0 is flat.Higher χ 2 init consistently produces more oblong haloes across all radii with axis ratios as low as q = 0.6 and s = 0.55, whereas low χ 2 init haloes are close to spherical.This is due to early up-scattering compressing the halo such that it loses any oblong shape, and it then remains spherical throughout its evolution.With early down-scattering, the halo remains puffy and there is not a strong enough influence from scattering to compress the halo into a sphere.Figure 9 shows the anisotropy parameter β as a function of radius for each simulation, with higher β corresponding to more radial orbits.In the presence of self-interactions, β decreases across all radii.
• The substructure is significantly decreased as χ 2 init increases.For the χ 2 init = 0, the substructure is reduced by a factor of 2 compared to the CDM model, and the χ 2 init = 1 model decreases the number of satellites by over an order of magnitude.Figures 10  and 11 show the number of subhaloes enclosed within a given radius and the subhalo mass function respectively.Both of these plots show the same trend of the shape of the functions remaining consistent but moving to lower numbers.With higher χ 2 init , down-scattering occurs frequently early in the halo's formation.This gives particles a velocity kick, which makes it more difficult to hold them in gravitational wells.This makes small subhaloes without much gravitational pull difficult to form, and even larger subhaloes lose some of the particles early in their evolution.This effect occurs across radius and mass ranges, indicating that there is not a large environmental dependence on the effect but that the main effect is the initial states and how it alters the early evolution of the halo.The exception to this is in large subhaloes ∼ 10 10 M ⊙ , where the gravity is stronger and the effect is slightly less pronounced.
• Figure 12 shows the fraction of excited state particles as a function of redshift for each simulation.These lines do not cross, so if a halo starts with a higher or lower χ 2 than another, this remains the case throughout the halo's evolution.However, χ 2 for each halo converge to more similar values at low redshift, so higher χ 2 init decrease χ 2 with time and lower χ 2 init increase χ 2 with time.The haloes converge to a value that is set by the velocity dispersion of the halo as described in Equation 6.
There are several other factors that influence halo properties that we did not consider here.The overall environment and mergers may cause certain reactions to dominate or provide an influx of new particles and substructures.A dedicated set of simulations using the same early universe conditions but simulating a variety of halos would shed light on the diversity problem in particular.In addition, baryons could alter several of the properties, for example by deepening the central gravitational potential.This could cause a difference in, for example, whether subhaloes go through core collapse or not.Dissipative dark matter has been shown to accelerate gravothermal collapse (Essig et al. 2019), which mimics the behaviour of the endothermic reactions.However, the two-state dark matter particle includes exothermic reactions, which causes haloes to become more diffuse.Thus, some haloes may experience accelerated core collapse while others may not, depending on the histories of the subhaloes.Currently, most constraints exist for elastic SIDM (Adhikari et al. 2022), while other types of SIDM are not well constrained.Because the two-state model introduces several types of scattering, certain reactions will dominate in different situations.Studying how the halo properties like shape and density profile change depending on mass scale and merger history could help constrain models like this one.Larger haloes with large velocity dispersions may be more dominated by endothermic reactions since particles can more easily cross the velocity threshold.For this, it would be interesting to run a set of cluster simulations using this model.
The initial conditions of the halo result in significantly altered halo properties.In particular, high χ 2 init results in a halo with properties consistent with those from exothermic SIDM, and low χ 2 init results in properties consistent with endothermic SIDM.To accurately predict the expected halo, it is therefore important to have a solid foundation in the initial conditions of the model.The sensitivity of the resulting halo to the initial conditions indicate that this model can produce a range of shapes and density profiles from compact to diffuse, potentially alleviating the diversity of shapes problem.The importance of this work lies in highlighting the sensitivity of the halo to early universe conditions while controlling for these other variables.

Figure 2 .
Figure 2. The projected DM density for the main halo of the CDM simulation at z = 0.The main halo is surrounded by many subhaloes of various sizes and is slightly oblong.

Figure 4 .
Figure 4.The projected DM density for particles in the ground (left) or excited (right) states in the simulations with 0, 20, 40, 60, 80 or 100% of the particles starting in the excited state at z = 0. Most partices are in the ground state, and the halo does not change much in appearance between the ground state projections and the total particle projections.The excited state particles occur mostly in the centre of the halo, where the velocities are high enough for up-scattering to occur.The excited state particles also follow the trend of decreased density as the initial excited state fraction increases.
Figure5.Density profiles and velocity dispersions for the main halo at z=0.Simulations with fewer initially excited particles tend to have denser centers.When all of the particles begin in the ground state, every downscatter must be preceded by an upscatter, so there is more upscattering than downscattering.The right-hand y axis shows the equilibrium excited state fraction, χ 2 equil , at relevant velocity dispersions, using σ 3d =

Figure 6 .
Figure6.M 200,mean and V max for the main halo in each simulation.M 200,mean decreases with χ 2 init due to the kinetic energy boost from downscattering, which can push particles out of the halo.V max also decreases with χ 2 init , reflecting trends in both halo mass and concentration.Haloes with low χ 2 init experience more up-scattering, driving particles inwards.
Figure7.The fraction of main halo particles in the excited state as a function of radius at z = 0.In the outer region (r ≳ 100 kpc), χ 2 increases with χ 2 init .This is due to the lower density at these radii, which decreases the impact of scattering.At smaller radii, χ 2 is not monotonic over χ 2 init , but simulations with higher χ 2 init counterintuitively tend to have lower χ 2 (z = 0).This is likely related to the lower inner velocity dispersion in these simulations, which causes down-scattering to dominate.
Figure8.We fit an ellipse to each halo at several radii and find the ratio of the middle axis (q, top panel) and minor axis (s, bottom panel) to the major axis.A ratio of s = 1 is spherical, and a ratio of s = 0 is flat.The grey points show the CDM run.Each column of points is the axis ratio at that radius for each simulation halo.We see that low χ 2 init haloes tend to be spherical due to early up-scattering that compresses the halo.High χ 2 init haloes maintain a more oblong shape since down-scattering keeps the halo enlarged.

Figure 9 .
Figure9.The anisotropy parameter as a function of radius from the halo centre.A larger β indicates a more radial orbit.The SIDM haloes have a lower β than the CDM halo due to the energy transfer between particles, which causes the orbits to be more isotropic.In particular, the SIDM anisotropy curves reach similar values at low radii, where they have similar cored slopes.

Figure 10 .
Figure10.Top: The number of subhaloes with mass > 10 7 M ⊙ within a radius r of the main halo at z = 0.As χ 2 init increases, the number of haloes decreases, although the shape of the curve remains the same.This is due to down-scattering that gives a velocity kick to particles and makes it more difficult to keep them within a gravitational well.Even at χ 2 init = 0, we see a decrease in substructure due to tidal forces.Bottom: The fractional difference in the subhalo number profile compared to CDM.The χ 2 init = 0 halo has approximately half as many subhaloes at each radius as the CDM halo, and this fraction decreases as χ 2 init increases.

Figure 12 Figure 11 .
Figure12shows the fraction of main halo particles in the excited state by redshift, with color denoting χ 2 init .The simulations converge to more similar values as time goes on, with high χ 2 init haloes

Figure 12 .
Figure12.The fraction of main halo particles that are excited by redshift, in the low-resolution simulations.By z = 0, each line is approximately flat, indicating that the densest regions of the halo have reached equilibrium.At z ≈ 2, simulations with low χ 2 init have a sharp increase in χ 2 , corresponding to a halo merger (grey dashed line) that increases the velocity dispersion and the rate of upscattering.