# Halo heating from fluctuating gas in a model dwarf

Mahmoud Hashim,1 Amr A. El-Zant,1 Jonathan Freundlich,2 Justin I. Read3 and
Françoise Combes4,5
1Centre for Theoretical Physics, The British University in Egypt, Sherouk
City, 11837 Cairo, Egypt
2Observatoire Astronomique, Universite de Strasbourg, CNRS, 11 rue de
l’Universite, 67000 Strasbourg, France
3Department of Physics, University of Surrey, Guildford, GU2 7XH, United
Kingdom
4LERMA, Observatoire de Paris, CNRS, PSL Univ., Sorbonne Univ., F-75014 Paris,
France
5Collège de France, 11 place Marcelin Berthelot, F-75005 Paris, France E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The cold dark matter (CDM) structure formation scenario faces challenges on
(sub)galactic scales, central among them being the ‘cusp-core’ problem. A
known remedy, driving CDM out of galactic centres, invokes interactions with
baryons, through fluctuations in the gravitational potential arising from
feedback or orbiting clumps of gas or stars. Here we interpret core formation
in a hydrodynamic simulation in terms of a theoretical formulation, which may
be considered a generalisation of Chandrasekhar’s theory of two body
relaxation to the case when the density fluctuations do not arise from white
noise; it presents a simple characterisation of the effects of complex
hydrodynamics and ‘subgrid physics’. The power spectrum of gaseous
fluctuations is found to follow a power law over a range of scales,
appropriate for a fully turbulent compressible medium. The potential
fluctuations leading to core formation are nearly normally distributed, which
allows for the energy transfer leading to core formation to be described as a
standard diffusion process, initially increasing the velocity dispersion of
test particles as in Chandrasekhar’s theory. We calculate the energy transfer
from the fluctuating gas to the halo and find it consistent with theoretical
expectations. We also examine how the initial kinetic energy input to halo
particles is redistributed to form a core. The temporal mass decrease inside
the forming core may be fit by an exponential form; a simple prescription
based on our model associates the characteristic timescale with an energy
relaxation time. We compare the resulting theoretical density distribution
with that in the simulation.

###### keywords:

dark matter – galaxies: haloes – galaxies: evolution

††pubyear: 2022††pagerange: Halo heating from fluctuating gas in a model
dwarf–C.4

##  1 Introduction

In the standard cosmological scenario, CDM particles decouple while non-
relativistic from the primordial thermal bath. Their consequently small speeds
allow for a clustering that agrees well with an impressive array of
observations on large scales (Frenk & White, 2012). But this same ‘coldness’
leads to excessive central densities in the nonlinear regime. As soon as they
could resolve the central regions, simulations showed that CDM haloes were
endowed with central density cusps, with their density diverging towards the
centre (Dubinski & Carlberg, 1991; Warren et al., 1992). It subsequently
transpired that this reflected nearly ‘universal’, self-similar, behaviour;
with the density profiles of haloes at different masses, identified at various
epochs, well fit with by simple empirical formulae (Navarro et al., 1996b).
But the central density cusps seemed too dense to match those inferred for
dark matter dominated galaxies, like low-surface-brightness and dwarf galaxies
(Moore, 1994; Flores & Primack, 1994).

This ‘cusp-core’ problem is central to the low-redshift small-scale issues
associated with CDM structure formation, which have attracted much attention
during the past decades, as its resolution may simultaneously alleviate other
puzzles associated with small scale structure in CDM simulations, such as the
so-called ‘too-big-to-fail’ problem (Read et al., 2006; Boylan-Kolchin et al.,
2011; Ogiya & Burkert, 2015; for general reviews see, e.g., Del Popolo & Le
Delliou, 2017; Bullock & Boylan-Kolchin, 2017; Salucci, 2019). CDM-based
simulations also predict a general over-abundance of small haloes compared to
small galaxies. However, down to the resolution level of current cosmological
simulations and expected minimal galaxy mass, this issue has been claimed to
be less severe as small Milky Way satellite galaxies are being discovered and
counted in (e.g., Nadler et al., 2020).

Assuming that the data has not been misinterpreted (e.g. Oman et al., 2019),
solutions to the aforementioned problems can be principally divided into two
sets; those modifying the underlying dark matter particle physics model and
those invoking gravitationally mediated CDM interaction with baryons. The most
popular of the first category have been warm dark matter and self interacting
dark matter, and the more recently topical fuzzy dark matter. The simplest
models of the latter type have problems explaining core scaling relations in
galaxies (Deng et al., 2018; Burkert, 2020; Bar et al., 2021), and warm dark
matter does not generally succeed in producing central density cores at all,
as cold collapse occurs unless the dark matter is ‘warm’ enough to prevent the
halo at hand from forming in the first place (Macciò et al., 2012). Warm dark
matter may make baryonic core formation easier, due to lower halo
concentrations, but this is debated since warm dark matter also acts to
suppress early star formation (Governato et al., 2015). In general particle
physics based modifications to the standard structure formation scenario
invoking warm dark matter or fuzzy dark matter, have been more successful in
alleviating the apparent over-abundance of small scale structure in CDM-based
structure formation. But this, instead of being a strong point, has now joined
Lyman-$\alpha$ constraints (e.g., Iršič et al., 2017a, b) and strong lensing
(e.g., Gilman et al., 2020) in setting limits on the efficacy of such models
(Simon & Geha, 2007; Koposov et al., 2008; Kim et al., 2018; Read & Erkal,
2019; Newton et al., 2021; Nadler et al., 2020; Nadler et al., 2021). In
addition, models where small haloes are invariably expected to come with
significant cores, such as in the fuzzy dark matter scenario, may be
constrained by the survival of Milky Way satellites with eccentric
trajectories that take them to the host halo’s central region (Errani et al.,
2023). Interacting dark matter-based models on the other hand generically
suffer from eventual gravothermal contraction of cores, which are not easily
avoided for all relevant spatial and time scales (Burkert, 2000; Kochanek &
White, 2000; Turner et al., 2021), ( although, fine tuning aside, this may
actually be considered a phenomenological strong point of a subset of such
models; Yang et al., 2022).

A non-negligible scattering crossection between standard model particles and
dark matter has also been proposed (Famaey et al., 2020; Salucci et al.,
2020), although this option has not been as fully investigated so far,

On the other hand, gravitational coupling between dark matter and baryons has
long been known to be effective in transforming CDM cusps into cores. Its main
shortcoming however comes from the complex physics and the associated
uncertainties in the parameters involved. Indeed, cusp-core transformation
through such coupling come in three different forms: one time mass blowout due
to a single burst of energy feedback (Navarro et al. 1996a, and more recently
Freundlich et al. 2020a and Li et al. 2022); the pumping of energy from
(gaseous or stellar) baryonic clumps to CDM via dynamical friction (El-Zant et
al., 2001; El-Zant et al., 2004); and density and potential fluctuations in
feedback-driven gas during galaxy formation (Read & Gilmore, 2005; Pontzen &
Governato, 2014). Despite the apparent similarity of the first and third
mechanisms (as they both invoke feedback), it is the last two that are in fact
more closely related at a deeper level — as they both involve a long-lived
fluctuations progressively heating the CDM on timescales much longer than that
of a single feedback starburst, as discussed in El-Zant et al. (2016).

Dynamical friction heating is nevertheless expected to be observationally
distinct from supernova heating since the latter correlates with star
formation (e.g. Read et al., 2019) while the former does not, and both may act
in tandem during galaxy formation (Orkney et al., 2021; Dekel et al., 2021;
Ogiya & Nagai, 2022). In the particular case of dwarf galaxies, there is
mounting observational evidence for dark matter heating from impulsive,
‘bursty’, star formation, ubiquitous in dwarfs below a stellar mass of
$M_{*}\sim 10^{8}$ M⊙ (e.g. Collins & Read, 2022). Indeed, there is an
observed anticorrelation between the inner dark matter density of dwarfs and
their stellar-to-halo mass ratio $M_{*}/M_{200}$ (Read et al., 2019; Bouché et
al., 2022), which is a proxy for the amount of star formation – and therefore
dark matter heating – that has taken place (e.g. Peñarrubia et al., 2012; Di
Cintio et al., 2014; Freundlich et al., 2020b). This may favour a scenario of
feedback-driven core formation in such galaxies.

As conjectured in El-Zant et al. (2016), it may be possible to understand the
process of feedback-driven core formation with repeated bursts from first
principles, bypassing the uncertain complexities of ‘gastrophysics’ and its
various subgrid implementations. Indeed the stochastic dynamical model
presented there suggested that the process of core formation depended
primarily on just two parameters, namely the gas mass fraction and the
strength of the fluctuations characterised by the normalisation of a power-law
power spectrum. It is our purpose here to study the properties of feedback-
driven fluctuations in a full hydrodynamic simulation where a cups-core
transformation occurs, with the aforementioned model forming an interpretive
framework and guide. Thus testing, in the process, its assumptions and
predictions.

##  2 Modelling halo heating from gas fluctuations

###  2.1 Physical setting

The general picture envisioned, in El-Zant et al. (2016), and here, is that of
a gas settling into a CDM halo. As it contracts, a critical density is
reached, leading to star formation and consequent starburst. The star
formation process assumed is akin to that described in Teyssier et al. (2013)
and Read et al. (2016). Namely, the threshold is considered high enough so
that most star formation does not occur in a few bursts, but instead
repetitively over a long timespan. This leads to sustained density and mass
fluctuations in the gas that appear amenable to a description in terms of a
stationary stochastic process over the timescale of the simulation.

If, to a first approximation, the gas density is assumed to be isotropic and
homogeneous when averaged over large spatial or time scales, with average
density $\rho_{0}$, then the mass fluctuations within a spatial scale $R$ can
(as in characterising cosmic structure) be characterized by a dispersion

$\sigma_{R}^{2}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}W^{2}(k,R)\mathcal{P}(k)k^{2}{\rm
d}k,$ (1)

where $\mathcal{P}(k)$ is the equal time power spectrum of the density
fluctuations $\delta\rho({\bf r},t)/\rho_{0}-1$, and $W$ is a Fourier filter
function. If the fluctuations furthermore constitute a stationary Gaussian
random process, this is all what one needs to know to completely characterise
them. The stochastic dynamics can likewise be described completely in terms of
the first and second moment statistics of the force field, which are easily
obtainable from the Poisson equations: in particular, the $k$-modes of the
potential fluctuations are related to those of the density by
$\phi_{\boldsymbol{k}}=-4\pi G\rho_{0}\delta_{\boldsymbol{k}}k^{-2}$, as
discussed in El-Zant et al. (2016) and, in another context, in El-Zant et al.
(2020). As shown in the latter work, this formulation can also describe
standard two-body relaxation. As such, it can be used to calculate the effect
of fluctuations arising from the presence of a collection of massive particles
on a system of lighter ones. Relaxation in this context ‘heats’ the light
particles, by increasing their velocity dispersion, while the heavier
particles lose energy via dynamical friction. The general theoretical
framework used here thus also applies to dynamical friction heating of the
halo in the presence of massive clumps. The relevant power spectrum of density
fluctuations, in that case, is flat (white noise) over scales larger than the
maximum size of the (monolithic) clumps. Thus the description in terms of
dynamical friction heating is more relevant when the gas fluctuations can be
described in terms of a system of long lived distinct clumps with sizes
significantly smaller than the region where the core forms.

When feedback is strong, and the gas fully turbulent on scales relevant to
core formation, the spectrum of density fluctuations may be assumed to be
approximated by a power law over such scales. As in the white noise case, the
fluctuations lead to energy transfer from the gas to the halo particles, which
drives the transformation of the latter’s problematic central cusp into a
core. The energy transfer may still be described in particularly simple terms,
reminiscent of Chandrasekhar’s theory of two body relaxation; for a given halo
and power spectrum $\mathcal{P}\propto k^{-n}$, the associated relaxation time
$t_{\rm relax}$ is predicted to principally depend just on the normalization
of that spectrum at the minimal wavenumber $\mathcal{P}(k_{m})$ (at which the
power law breaks), and the average gas density $\rho_{0}$. More concretely,
for a halo particle moving with an unperturbed speed $v_{p}$ in a field of gas
fluctuations carried by bulk flows moving with characteristic speed $v_{r}$
relative to it, one finds

$t_{\rm relax}=\frac{nv_{r}v_{p}^{2}}{8\pi(G\rho_{0})^{2}\mathcal{P}(k_{m})}.$
(2)

This is the timescale for the random velocity gained by a CDM particle, as a
result of its motion in the stochastic force field born of gas fluctuations,
to reach the unperturbed characteristic speed $v_{p}$. It is analogous to the
standard relaxation time. In fact equation (2) reduces to Chandrasekhar’s
formula for a white noise spectrum ($n=0$), with Coulomb logarithm determined
by the ratio of maximal to minimal wavenumumber cutoffs (El-Zant et al., 2020,
Section 2). Also, as equation (2) is derived using Fourier transformation,
which is transitionally invariant and therefore necessarily leads to loss of
spatial information, it is (like Chandrasekhar’s formula) strictly valid only
for an infinite homogeneous system with a constant level of fluctuations. In
applying it to a real system, we will thus assume an approach that is akin to
that widely invoked with the use of Chandrasekhar’s theory in an inhomogeneous
system, where it is customary to assume that the variation of the relaxation
time with position is accounted for by the density variation. In our
formulation of Chandrasekhar’s theory, referred to above, this would mean
keeping $n=0$ and the ratio of maximal to minimal wavenumber fixed, while
adjusting the density to account for radial variations in the relaxation time.
Analogously, we will keep the power spectrum parameters in equation (2) fixed
and assume that the spatial variation of the energy input to the dark matter
from the gas fluctuations can be accounted for by variation in the density.

To account for the radial change in the average gas density in a realistic
system, we define the average density, at a given time, within a sphere with
radial coordinate $r$ by

$\rho_{0}(r,t)=\langle\rho_{g}({\bf r},~{}t)\rangle_{|{\bf r}|<r},$ (3)

where $\rho_{g}({\bf r},t)$ is the local gas density and the average is
evaluated within the volume enclosed by $r=|{\bf r}|$. A further average over
time results in $\rho_{0}(r)$. If the averaging is not evaluated over a sphere
centered at the origin, but over spheres with centres at radial coordinate $r$
and radius $r_{\rm av}$, an analogous average may also be defined (as in
equation 9). In Section 4.2, we examine three different approximations for an
effective $\rho_{0}$, in order to account for changes of relaxation time with
radius while using equation (2): one associated with what we will call global
energy transfer, i.e. with $\rho_{0}$ a constant corresponding to an average
inside an energy transfer saturation radius $r_{\rm sat}\gtrsim r_{s}$, which
turns out to be plausible, given the slow variation of average gas density
inside $\sim r_{s}$ and subsequent sharp drop (Fig. 3); a local approximation,
which does not work; and one involving averages over radius $r_{\rm
av}=r_{s}\simeq{\rm 1kpc}\approx 1/k_{m}$, which we will find works best (Fig.
8). This is perhaps not surprising given the aforementioned gas density
variation, and given that equation (2) implies that modes near $k_{m}$
contribute most to the relaxation process. (It in fact arises from a low-$k$
‘effective theory’, as discussed further in the penultimate paragraph of
Appendix A).

The characteristic velocity $v_{r}$ appearing in equation (2) is estimated
through the ‘sweeping’ approximation, long used in turbulence theory (Taylor,
1938; Kraichnan, 1964; Tennekes, 1975) and invoked by El-Zant et al. (2016).
It assumes that the equal-time spatial statistics of the fluctuation field are
swept (’frozen in’) into the time domain through large scale fluid flows. In
this picture $v_{r}$ consists of random and regular velocity components $U$
and $V$, such that $v_{r}=\sqrt{U^{2}+V^{2}}$ that are characteristic of those
flow. In our case, since the halo particles are also moving with respect to
the gaseous flows, the characteristic $v_{r}$ may be considered to include
such motions. 111In principle, each fluctuating mode may have its own sweeping
speed and the velocity distribution of the particles can be taken into
account, as in El-Zant et al. (2020). We do not consider this more complex
case here. We test this assumption in Section 3.3.

In the context of the formulation outlined in this section, the timescale of
significant central halo transformation, from cusp to core, should scale as in
equation 2, with parameters as defined above. Controlled simulations (using
the Hernquist-Ostriker code (Hernquist & Ostriker 1992), whereby Gaussian
random noise was applied to particles of live CDM haloes, produced cores on
the expected timescale, when the haloes were kept strictly spherical (El-Zant
et al., 2016). 222When this condition was relaxed, the timescales for core
formation were found to be about an order of magnitude shorter. But the
scaling with the normalisation of the power spectrum and average density, as
reflected in equation (2), remained. As noted in Section 4.3, we find no
evidence of such accelerated (relative to the relaxation time) core formation
rate here.

###  2.2 Numerical implementation

Figure 1: Gas projected density map inside a box of length 5 kpc (2.5 kpc on a
side from the centre of the halo). The gas contracts from the initial NFW
profile, then much of it is expelled from the inner region during the first
few hundred Myr. It gradually settles into a disk-like configuration, but its
distribution fluctuates over time.

We wish to examine feedback-driven gas fluctuations in a full hydrodynamic
simulation, to find out whether their characteristics and their effects may
fit the picture outlined above. To facilitate the isolation of the principal
features, we focus on the case of an isolated galaxy. Its parameters (Table 1)
are chosen to correspond to those of dwarf galaxies, where dark matter
dominates even in the central regions and the discrepancy between the inferred
density and a centrally concentrated halo profile seems particularly clear.

The simulation of the isolated dwarf is set up and run with the RAMSES
adaptive mesh refinement code (Teyssier, 2002), and already presented in Read
et al. (2016). We use ‘Run M9c224e6’ from that paper, which has 281 time
outputs regularly spaced over a timespan of 13.7 Gyr. The initial condition
for this simulation assumes a NFW dark matter halo (Navarro et al., 1996b)
with total mass inside the virial radius $M_{200}=10^{9}$ M⊙ and concentration
parameter $c_{200}=22.23$. A fraction $f_{b}=0.15$ of that mass ($M_{g}(r_{\rm
vir})=0.15\times 10^{9}\,M_{\odot}$), is in gas that is initially in
hydrostatic equilibrium, with a metallicity of $10^{-3}$ Z⊙ and some angular
momentum set to match median expectations in a $\Lambda$CDM cosmology. While
the simulation assumes that the dwarf galaxy contains the Universal baryon
fraction out to its virial radius, the ratio between the gas mass initially
within the region of interest for core formation ($\sim 1~{}{\rm kpc}$) and
the mass within the virial radius is much lower, and thus consistent with
observational constraints (e.g. Read & Trentham, 2005). At the start of the
simulation, the gas rapidly cools and collapses, then re-expands under the
influence of feedback. In the process, the central gas mass fraction further
decreases (cf. Section 3.1).

The sub-grid model for star formation and feedback is described in detail in
Read et al. (2016). Briefly, star formation follows a Schmidt relation with a
star formation efficiency per free fall time of $\epsilon_{\rm ff}=0.1$ and a
density threshold for star formation of $\rho_{*}=300$ atoms cm-3. The stellar
feedback model is as in Agertz et al. (2013) and includes a model for Type II
and Tyle Ia supernovae, stellar winds and radiation pressure. The simulation
resolution was chosen to capture momentum driving from individual supernovae
events, with a gas spatial resolution (at highest mesh refinement) of $\Delta
x\sim 4$ pc and a mass resolution in gas, stars and dark matter of $m_{\rm
g}=60\,{\rm M}_{\odot}$ and $m_{*}=M_{\rm DM}=250\,{\rm M}_{\odot}$,
respectively. This allows individual star formation events to be resolved,
injecting stellar feedback to the interstellar medium in the correct locations
at the correct times.

Table 1: Initial profile parameters. NFW scale length | $r_{s}$ | $0.88\;{\rm kpc}$
---|---|---
NFW characteristic density | $\rho_{s}$ | $5.34\times 10^{7}~{}{\rm M_{\odot}~{}kpc^{-3}}$
Concentration parameter | c | $22.23$
Mass within the virial radius | $M_{\rm vir}$ | $10^{9}\;{\rm M_{\odot}}$
Baryonic mass fraction | $f_{b}$ | $0.15$

##  3 Characterizing the gas fluctuations

Figure 2: Time variation of the gas mass enclosed within different radii $r$
(in kpc), suggesting a quasi-stationary stochastic process.

Figure 3: Left: Average gas density inside radius $r$, $\rho_{0}(r,t)$, for
all snapshots, obtained by using equation (3). The dashed line is a time
average over the indicated interval, which corresponds to the simulation time
beyond an initial non-equilibrium phase through which much of the gas is blown
out of the region inside $r_{s}\simeq 1~{}{\rm kpc}$. The time averaged
profile is only mildly varying inside $r_{s}$. Right: Ratio of the gas to
combined gas and dark matter masses $f_{g}=M_{g}/(M_{g}+M_{d})$, with
$f_{g}(r_{s})$ denoting the fraction inside $r_{s}$; thus here
$M_{g}=M_{g}(<r_{s},t)$ and $M_{d}=M(<r_{s},t)$, while $f_{g0}$ denotes the
ratio of the gas mass to the initial total mass enclosed:
$f_{g0}(r_{s})=M_{g}(<r_{s},t)/[M_{g}(r_{s},0)+M_{d}(<r_{s},0)]$. Similarly
$f_{g}(r_{1/2})$ and $f_{g0}(r_{1/2})$ refer to those ratios inside the half
mass radius of the stellar component that forms from the repeated starbursts,
$r_{1/2}\simeq 0.5{\rm kpc}$, and which physically correlates with the scale
of the relatively steady gas density associated with strong fluctuations.

In order to apply the physical model described above, leading to equation (2),
we need to know the average gas density $\rho_{0}$. We also need to evaluate
the power spectrum of the gas fluctuations and find out whether it may indeed
be adequately described by a power-law, in which case we need to estimate
values for its normalization $\mathcal{P}(k_{m})$ and index $n$, and verify
that the spatial spectrum is relevant in describing the temporal fluctuations.
This is the objective of this section.

###  3.1 Gas density and mass over time

Fig. 1 shows gas density contrast maps, within a 5 kpc box, 2.5 kpc on a side
from the halo centre333Unless otherwise stated, all centering here refers to
the halo centre located using the shrinking sphere method.. The gas initially
contracts from an NFW profile with scale length $r_{s}$ until the critical
threshold for star formation is reached. The fluid is then driven by the
nascent feedback, much of it driven completely out of the region delineated by
$r_{s}$; the gas mass inside $\sim$$r_{s}$ decreases to less than half its
initial value (but the dark matter mass is largely unaffected, decreasing only
by about $20\%$ within $r_{s}$, and thus the ratio decreases significantly, as
shown in the right hand panel of Fig. 3). After that, the gas mass inside
$\sim$$r_{s}$ is relatively well conserved. Variations in the gas mass
enclosed within a given radius then principally arise from a pattern of
fluctuations (gaseous blobs) materialising on a large range of scales. The gas
progressively settles into a disk-like configuration, but the mass
fluctuations remain sustained and steady.

The assumption that the fluctuations form a quasi-steady stochastic process is
suggested by Fig. 2, showing the mass variation within different radii over
time. The mass fluctuations persist in time but diminish with radius, as
higher mass scales are reached. Significant variations on larger time-scales
persist however. They correspond to large scale flows. In the context of the
model of El-Zant et al. (2016) (recapped in Section 2.1), such large scale
motions ‘sweep’ the smaller scale fluctuations with the characteristic speed
$v_{r}$ (equation 2). For the description of the effect of the gas
fluctuations on the halo particles in terms of a standard diffusion process
(leading to relaxation timescales in the form of equation 2) to be complete,
the statistics of the stationary stochastic process should also be entirely
described by averages and dispersions, as in a Gaussian random process. This
is examined in Appendix A.

The gas density inside spheres of radius $r$ is shown in Fig. 3, as a
functions of time. After the initial outflow during the first few hundred Myr,
the averaged (over time) gas profile is only mildly varying with radius inside
the initial NFW scale length ($r_{s}=0.88~{}{\rm kpc}$), and then starts
decreasing rapidly beyond that. Any core that forms is expected to be of the
order of the scale length of the cusp, which renders the assumptions in El-
Zant et al. (2016) (discussed in Section 2.1.2 therein), plausible. One may
for instance attempt to estimate the relaxation time using equation (2) with
mean density $\rho_{0}$ evaluated at radius $\gtrsim r_{s}$, as the sharp
decrease in density beyond that scale suggests that fluctuations in the gas at
larger radii would contribute little to the overall variations in the gas
potential and associated energy input to the central halo.

In Section 4, we estimate that this may indeed be a good approximation. Here
we note that there are two theoretically independent reasons why, given the
parameters of the simulation, a core is expected to form on a scale of order
$1~{}{\rm kpc}$. For, even if the gas density does not drop beyond $r\gtrsim
r_{s}$, and strong fluctuations are present much beyond that radius, the
resulting stochastic perturbations are expected to have a major dynamical
effect only up to scale of order $r_{s}$. This was in fact found to be the
case even when the amplitude of the fluctuations (as fixed by the
normalization of the power spectrum) were increased well beyond the minimum
level required to significantly affect the inner halo profile (cf. Fig. 6 of
El-Zant et al. (2016)). The phenomenon is likely due to the transition to an
isothermal profile in an NFW halo at $r\simeq r_{s}$, which is relatively
stable against fluctuations, as the differential energy distribution tends to
an exponential form (cf. Appendix C.4). This phenomenon relates to the effect
of the fluctuations.The other reason why a core that forms out of a cusp here
should have radial scale of order $1{\rm kpc}$ has to do with the existence of
fluctuations in a realistic system: the extent of the spatial scale where the
gas density and density contrasts are large is characteristic of the half mass
radius $r_{1/2}$ of the stellar component that forms from the repeated
starbursts. In our case about $0.5{\rm kpc}$. Because of this, we will be able
to define (in relation to Fig. 7) an energy input saturation radius, of order
$\simeq 2r_{1/2}$, which results from the decreasing gas density and level of
fluctuations. Beyond it, there is negligible energy transfer from gas to halo.

The right panel of Fig. 3 shows the fraction of gas inside both $r_{s}$ and
$r_{1/2}$. We note that due to the small variation of the gas density averaged
over radius, the lines closely correspond to one another. In addition, after
the gas loss associated with the initial blowout phase, the gas fraction
decreases to about half its initial value. Thus, any effect due to core
formation will result from fluctuations from a rather small gas fraction in
the inner regions. The efficacy of this perhaps surprisingly small central gas
fraction in modifying the dark matter dynamics is in fact a generic prediction
of our theoretical framework, provided the strength of the gas fluctuations is
at the level of the fiducial model in El-Zant et al. (2016) (cf. Section 2.1.2
and 4.4.2 there). Below, we will find that the fluctuation levels in the
present simulation, as measured by the normalisation of the power spectrum,
are in indeed consistent with those assumed in that fiducial model.

###  3.2 Power spectra

Figure 4: Power spectrum of gas density fluctuations, averaged over 1 Gyr
intervals. The shaded region corresponds to variation of a power law $\sim
k^{-n}$, with $2\leq n\leq 3$; the horizontal line corresponds to
$\mathcal{P}(k_{m})=4~{}{\rm kpc^{3}}$. The spectrum follows a power-law form
over a range of scales, with a best-fit exponential index $n_{\rm
bestfit}=2.31\pm 0.02$ (black dashed line).

We now proceed to determine the power spectrum of the density fluctuations
qualitatively probed above. We evaluate it in the following way. For a given
time snapshot, we consider gas cells within a cubical box of extent 2.5 kpc on
each side from the origin.444Taken here to be the centre of the halo located
through the shrinking sphere method. We have verified that the results are
robust to displacements of that centre $\lesssim{\rm 1~{}kpc}$, and similar
results are obtained when the gas centre of mass or shrinking sphere centres
are used. We assign each cell, depending on its radial location, to a
spherical shell. The average gas density inside the shell $\rho_{\rm sh}$ is
then subtracted from that at a given cell $\rho_{\rm cl}$. The density
contrast in each cell $\delta_{\rm cl}=\rho_{\rm cl}/\rho_{\rm sh}-1$ is
evaluated, and the Fourier transform is taken over the cells and squared in
the usual manner, to obtain the power spectrum (using the Fourier convention
employed in El-Zant et al., 2016). This is done for each time snapshot, and
the results are averaged over a chosen set of snapshots. The outcome is shown
in Fig. 4. The spectrum follows a power-law form over a range of scales, with
an index consistent with that assumed in the fiducial model of El-Zant et al.
(2016), where it was assumed to correspond to 2.4. If we take $k_{m}\simeq
2{\rm kpc^{-1}}$, then $\mathcal{P}(k_{m})\simeq 4~{}{\rm kpc}^{3}$ (again
close to what was assumed in the fiducial model of El-Zant et al. 2016, namely
$\mathcal{P}(k_{m})=4.6~{}{\rm kpc}^{3}$).

Figure 5: Kinetic energy power spectrum of gas motion with least squares best
fit (black dashed line) and Kolmogorov spectrum ($\propto k^{-5/3})$, blue
dashed line). The spectrum is calculated through the square of the Fourier
transform of $\rho^{1/2}_{g_{i}}{\bf v}_{i}$, taken over gas cells $i$ inside
a 5 kpc box around the halo (shrinking sphere) centre. The agreement with the
Kolmogorov form supports feedback-driven fluctuations having characteristics
of fully developed turbulence.

The picture of fully turbulent feedback-driven fluctuations is supported by
the near Kolmogorov form of the specific kinetic energy spectrum, shown in
Fig. 5, which is followed for a large range in wave numbers. This spectrum
persists even as the gaseous system evolves from a fully three dimensional
configuration to disk-like form. Indeed Grisdale et al. (2017) found that it
is present even in relatively quiescent phases, as long as feedback is also
present (noting that the spectrum flattened when feedback is eliminated). In
principle, turbulence may also be driven by self gravitating instabilities
rather than feedback (e.g., Yu et al., 2021; Nusser & Silk, 2022). A simple
calculation of the gas Toomre parameter in the present case nevertheless
suggests an absence of strong gravitational instabilities (see however Inoue
et al., 2016).

In the context of our theoretical framework, a further (likely related)
physical distinction may also be reflected in the shape of the density
fluctuation power spectrum; for, as mentioned in Section 2.1, a flat (white
noise) power spectrum corresponds to the limit when heating through dynamical
friction coupling between monolithic clumps and the dark matter is the
dominant process. In this case the clumps are long lived and distinct, and
their maximum size is significantly smaller than the scale associated with the
region where the cusp-core transformation takes place (and on which the
spectrum takes the flat form). The steep power-law dependence that continues
up to large scales, as found here, signals a preeminence of feedback driven
fluctuations. Only at the largest scales ($<k_{m}$) is a flattening of the
spectrum recovered. (The flattening at large $k$ denotes the start of the
transition from the turbulent inertial range to the dissipation range, set
here by the numerical resolution, given the small kinetic viscosity, and may
reflect a ’bottleneck phenomenon’, resulting from a pile-up of inefficiently
dissipated energy as the transition is approached; Falkovich, 1994; Schmidt et
al., 2004; Kritsuk et al., 2007.)

Finally, we note that the break in the kinetic energy power spectrum at
$k_{m}\simeq 2~{}{\rm kpc}^{-1}$ is consistent with that inferred from Fig. 4,
suggesting a demarcation between the turbulent scaling regime and the large
scale flows that may carry the fluctuations into the time domain. We examine
this issue of the transposition of the fluctuation characteristics from the
spatial to temporal domains further in the next subsection.

###  3.3 RMS mass fluctuations

Figure 6: Gas mass dispersion at scale $R$,
$\sigma^{2}_{R}=\sum_{i}(M_{i}(R)-\langle M(R)\rangle_{i})^{2}/\langle
M(R)\rangle_{i}^{2}$, where $i$ denotes different time snapshots, and $M(R)$
is the mass inside a sphere of radius $R$. The centre of the sphere is
randomized; chosen from a homogeneous distribution inside radius 1 kpc
($\simeq r_{s}$). The process is repeated for 80 realizations and the average
$\sigma_{R}$ over the realizations (solid line) and dispersion (error bars)
are evaluated. The dotted line refers to a (least squares) best-fit using
equation (1) with a power-law power spectrum and a top-hat filter, while the
shaded region shows the variation with $n$ in the range of the corresponding
shaded region on Fig. 4, with moderately smaller normalization at
$k_{m}=2~{}{\rm kpc}^{-1}$. All fits assume
$\mathcal{P}(<k_{m})=\mathcal{P}(k_{m})$. The power spectrum parameters,
inferred here from simulation data in the time domain, are generally
consistent with those of the equal-time power spectra of Fig. 4.

Our formulation of the dynamical effect of the gas fluctuations makes use of
the sweeping and random sweeping approximations of turbulence theory,
mentioned at the end of Section 2.1. At some level, this is analogous to an
ergodic assumption, in the sense that statistical properties of the random
field in the spatial domain are transferred into the time domain. If such an
approximation is valid in our case, then we should expect the following: if
the variance $\sigma_{R}^{2}$ on the left hand side of equation (1) is
calculated in the time domain over the simulation snapshots, it should
correspond to the result obtained by plugging the equal-time power spectrum
$\mathcal{P}(k)$, with parameters consistent with what is inferred above (Fig.
4), into the right hand side of equation (1).

In a stochastic process that is homogeneous in space and time, with the
sweeping assumptions holding, the $\sigma_{R}$ measured in the time domain
(over sufficiently long time) will be invariant with respect to the centre it
is measured from. This cannot be strictly the case here, however, since the
average gas density is not strictly homogeneous, although the power spectrum
of relative density fluctuations was found to be rather robust against changes
in the centre of the box inside which it is calculated (Footnote 4).
Furthermore, when attempting to match temporal fluctuations to those captured
by spatial power spectra at given snapshots one needs to eliminate effects not
captured by the frozen-in spectra. For example, the effect of fluctuations
carried by rotational flows will be reduced when calculating $\sigma_{R}$ in
spheres anchored close to the centre of rotation, with little mass transiting
through the shells (and those carried by large scale radial flows are
enhanced). To account for such effects, we randomise the centre of the sphere
within which the mass is calculated. The centres are, in practice, sampled
from a homogeneous distribution within the radius of interest for core
formation (1 kpc from the origin). We then conduct eighty different
realizations of such randomisations and evaluate the average of the dispersion
over the realisations.

Fig. 6 shows the result. The best-fit using equation (1) involves a power-law
power spectrum with index $n$ consistent with the lower limit in the shaded
region of Fig. 4; and also with similar normalization $\mathcal{P}(k_{m})$,
but at smaller $k_{m}$. The change in $k_{m}$, by itself, is not expected to
affect the dynamical effect of the fluctuations in the diffusion limit at the
basis of our model (with consequence that the relaxation time in equation 2
does not directly depend on $k_{m}$; see also Fig. 8 of El-Zant et al. 2016).
On the other hand, the shaded area of Fig. 6 suggests that the mass dispersion
is also consistent with larger values of $k_{m}$, but with mildly smaller
normalization (which does moderately affect the relaxation time).

The general consistency of the parameters inferred here, in comparison with
those obtained from the spatial power spectrum (Fig. 4), supports the
contention that the spatial fluctuations are transported to the time domain in
accordance with the sweeping assumptions (although the separation between the
large scale ‘carrier flows’ and the transported fluctuations may not be sharp,
and this may contribute to the smaller best fit $k_{m}$ here). As we will see
below, these parameters are also consistent with the actual energy input, from
the fluctuating gas to the halo component, in the simulation studied here; in
the sense that the inferred energy input is consistent with theoretical
estimates using such values.

##  4 Stochastic energy transfer and core formation

###  4.1 The relaxation time and its parameters

Now that we have estimates of the parameters entering equation (2), we can
apply it in order to check whether it predicts significant effect for the
model galaxy at hand. As a first estimate we simply assume $v_{p}\approx
v_{r}\approx v_{c}$, where $v_{c}=v_{c}(r_{s})$ is the circular speed at
$r_{s}$ (close to the maximal rotation speed). If we furthermore use a value
for $\rho_{0}$ characteristic of the region when the density starts to decrese
rapidly ($r\gtrsim r_{s}$ from Fig. 3), we find

$t_{\rm relax}=13.2~{}{\rm
Gyr}~{}\frac{n}{2.5}\left(\frac{v}{20{\rm~{}km/s}}\right)^{3}\left(\frac{\mathcal{P}(k_{m})}{3{\rm~{}kpc}^{3}}\right)^{-1}\left(\frac{\rho_{0}}{10^{6}~{}{\rm
M}_{\odot}/{\rm kpc}^{3}}\right)^{-2}\\!\\!.$ (4)

This suggests that gas fluctuations may have a significant effect on dark
matter orbits within $r_{s}$ on a timescale of the order of a Hubble time. If
$v_{p}$ and $v_{r}$ are associated with the characteristic velocity dispersion
(rather than $v_{c}$), which is slightly larger (by a factor of about $10\%$
at $\sim$$r_{s}$), the above timescale increases by a factor of about
$(1.1)^{3/2}$.

The timescale also increases if the characteristic speed of the halo particles
relative to the gas $v_{r}$ (cf. Section2.1) is associated with gas velocities
at $\sim$$r_{s}$, which are around ${\rm 30~{}km/s}$ (cf. Fig. 13). Because
flows may take the gas well beyond $r_{s}$, gas velocities are generally
larger than $v_{p}$. Once a rotating disk starts to materialise (Fig. 1), the
regular component would also include gas rotation, in addition to the orbital
speeds of the halo particles. Below, we will generally set $v_{r}\simeq
30~{}{\rm km/s}$, but note that $v_{r}$ should increase with radius inside
$r_{s}$, and that it may generally be larger than the value adopted here when
all the aforementioned motions and flows are taken into account.

The variation of the relaxation time with radius also depends on whether the
energy input from gas fluctuations is assumed to be primarily global or local:
the relaxation mechanism may in principle be global, in the sense of depending
only on some effective $\rho_{0}$, corresponding to a space and time average
over an appropriately chosen region (e.g., the region within which the
feedback driven fluctuations are significant), or it may be local; depending
on the local gas density $\langle\rho_{g}\rangle(r)$ (averaged only over
time). This is an issue we discuss further below, in connection to Fig. 8.

For a global energy input, the relaxation time depends only on variation in
$v_{p}$, if $v_{r}$ is kept fixed. with radius. If $v_{p}$ is associated with
$v_{c}(r)$, then within the initial NFW cusp $v_{p}^{2}\sim r/r_{s}$ increases
by a factor of 6.25 as $r/r_{s}$ goes from 0.1 to 1. If $v_{p}$ is taken to
correspond to the velocity dispersion then $v_{p}^{2}\sim-r/r_{s}\ln r/r_{s}$
deep inside the cusp, and increases by a factor of about 2.25 between
$r/r_{s}=0.1$ and 1.

Thus, particles nearer to the centre are expected to be affected first by the
fluctuations; if only simply because their initial velocities, and therefore
the relaxation times required to affect them significantly, tend to be
smaller. Core formation would thus proceed inside out, with mass distribution
at larger radii affected at later times.

In what follows we evaluate the energy input in the simulation at hand and
compare it with the theoretical picture sketched here; testing the ‘inside
out’ scheme of core formation and the assumption of dependence of the energy
transfer principally on the average density inside $\sim$$r_{s}$ rather than
the local density.

###  4.2 Energy input

Figure 7: Energy input to halo particles from gas fluctuations inside
indicated radii $r$. $E_{r}$ is evaluated via equation (8), and measured from
the temporal zero point $t=t_{\rm diff}$, beyond which the diffusion limit may
be assumed to hold. In line with the discussion following equation (5), we
take $t_{\rm diff}=1~{}{\rm Gyr}$. In order to reduce fluctuations related to
the precise choice of zero point, we calculate in practice $\Delta
E_{r}=E_{r}(t>t_{\rm diff})-\langle E_{r}\rangle_{t_{\rm diff}}$, where the
average is taken in the interval $0.5~{}{\rm Gyr}\leq t\leq 1.5~{}{\rm Gyr}$.
The dashed line indicates energy transfer through a stationary stochastic
process, with parameter values appearing in equation (7). It assumes that
energy input saturates around $r_{\rm sat}=1.5~{}r_{s}\simeq 1.3~{}{\rm kpc}$
(corresponding to the converging lines). For smaller radii, the lines flatten
after a few Gyr due to mass and energy transfer resulting from particle
migration towards larger radii.

Figure 8: Local vs global energy transfer. Left: Radial energy transfer as a
function of time. Right: Energy input rate as a function of radius. The left
hand panel shows that radial energy transfer is slow compared to the energy
input; the straight lines represent steady diffusive energy input, with no
transfer outside the radii enclosed. While this holds, one may estimate the
energy input rate at various radii by simply evaluating the average slope of
the generally linear energy increase. The right hand panel shows those slopes
over the first 1, 2, 3 and 4 Gyr (solid lines with indicated values of $T$).
The dashed line corresponds to estimate assuming global energy transfer; using
equation (6) with $n=2.5$, $v_{r}=30.0~{}{\rm km/s}$,
$\mathcal{P}(k_{m})=3.6~{}{\rm kpc^{3}}$, and average gas density within
$r_{s}$, $\rho_{0}(r_{s})=\langle\rho_{0}(r_{s},t)\rangle_{T}$, with
$\rho(r,t)$ measured directly from the simulation (using equation 3) and
averaged over $T=[0:4]~{}{\rm Gyr}$ (giving $\rho_{0}=1.95\times
10^{6}~{}M_{\odot}/{\rm kpc^{3}}$). The dotted line corresponds to the
assumption of purely local energy input; using the indicated equation, where
the local gas density $\rho_{g}({\bf r},t)$ is time-averaged (in the range
$T=[0:4]~{}{\rm Gyr}$) over spherical shells at $r$. This vastly overestimates
the energy input at all radii. The dashed dotted line represents an
intermediate alternative, whereby $\rho_{g}$ is averaged over spheres of
radius $r_{s}$ and centered at $r$ (as in equation 9).

Equation (2) implies that a particle increases its velocity variance with time
$T$ as

$\langle(\Delta
v)^{2}\rangle=\frac{8\pi(G\rho_{0})^{2}\mathcal{P}({k_{m}})}{nv_{r}}T.$ (5)

The linear time dependence here is characteristic of a diffusion process. It
is valid in the diffusion limit, which means that $T\gg t_{\rm
diff}\equiv(v_{r}k_{m})^{-1}$ (El-Zant et al., 2016). To account for this in
practice, we will generally measure $T$ and $\langle(\Delta v)^{2}\rangle$
from zero points at times $t\gtrsim t_{\rm diff}$. This is necessitated not
just by the applicability of the diffusion limit, but also by the assumption
of a quasi-steady stochastic process, which is clearly invalid during the
first few hundred Myr, which are characterised by a highly evolving phase,
involving the initial gas contraction and the triggering of a starburst that
expels most of it out of $r_{s}$.

Assuming $v_{r}\simeq 30~{}{\rm kpc/Gyr}$, gives $(v_{r}k_{m})^{-1}\lesssim
1/30{~{}\rm Gyr}$, for $1/k_{m}\lesssim 1~{}{\rm kpc}$. Thus the steady state
diffusion limit may be safely assumed to be established beyond 1 Gyr. We will
take it to be our zero point for $T=t-t_{\rm diff}$.

As in standard two body relaxation, one may assume that fluctuations initially
change just the velocities of the halo particles 555This is strictly speaking
valid when the orbital period of halo particles $T_{p}$ is much longer than
the fluctuation timescale $T_{f}$. Here, they are generally of the same order.
However the impulsive approximation is known to remain reasonable even when
$T_{p}$ is moderately smaller than $T_{f}$ (e.g. Binney & Tremaine, 2008;
Aguilar & White, 1985; Peñarrubia, 2019a). For a circular orbit and gas mode
of wavenumber $k$, $T_{p}/T_{f}=krv_{r}/v_{c}(r)$. This is minimized in the
initial cuspy potential, and when $k=k_{m}$, and $v_{r}\simeq 30{\rm km/s}$
includes only the gas speed. For $k_{m}=2{\rm kpc}^{-1}$ (consistent with Fig.
4), this minimal ratio is $T_{p}/T_{f}\simeq 3$ for $r\simeq 1{\rm kpc}$. It
is significantly smaller than unity ($\lesssim 1/3$) only at $r<30{\rm pc}$.
Adiabatic corrections may be included (Peñarrubia, 2019a), though the
corrections are generally modest for our range of $T_{p}/T_{f}$ (Cf. Fig. 2 of
Peñarrubia, 2019a). In our context the impulsive regime also in fact places a
lower limit on $k_{m}$, independent of the flattening of the power spectrum in
Fig. 4, as the effect of fluctuations rapidly decreases beyond it.

The average energy change per unit mass resulting from this is $\langle\Delta
E\rangle=\langle\Delta~{}\frac{1}{2}(\Delta v)^{2}\rangle=v\langle\Delta
v\rangle+\frac{1}{2}\langle(\Delta v)^{2}\rangle$. In principle, the particles
may gain energy through the first term or lose it to the fluctuating field
through dynamical friction, as in any general diffusive dynamical process. But
in practice the mass of the particles is far too small and so the average
energy gain per unit mass is simply $\langle\Delta
E\rangle=\frac{1}{2}\langle(\Delta v)^{2}\rangle$.

By fixing $v_{r}$, and assuming the energy transfer mechanism to be global —
again, in the sense of depending on an effective characteristic density
$\rho_{0}$ rather than on the local $\langle\rho_{g}\rangle_{T}(r)$ — the
expected energy input to halo particles within radius $r$ may be estimated as

$E_{\rm in}(<r)=~{}\langle M(<r)\rangle~{}\langle\Delta E\rangle~{}=~{}\langle
M(<r)\rangle~{}\frac{4\pi(G\rho_{0})^{2}\mathcal{P}_{k_{m}}}{nv_{r}}~{}T,$ (6)

where $\langle M(<r)\rangle$ is the (time) average of the halo mass enclosed
within $r$.

We discuss the assumption of global energy input below. First, we assume it
holds at least approximately, and use it to estimate the total energy transfer
from the gas to the halo. Given the rapidly decreasing gas density beyond
$r_{s}$ (Fig. 3), one may suppose that the energy input saturates at $r_{\rm
sat}\gtrsim r_{s}$. If the energy transfer depends on global properties, then
$\rho_{0}$ may be considered to correspond to the space and time average of
the gas density inside $r_{\rm sat}$, i.e. $\rho_{0}({\rm
sat})=\langle\rho_{0}(r_{\rm sat},t)\rangle_{T}$, where $\rho_{0}(r_{\rm
sat},t)$ is given by (3), with $r=r_{\rm sat}$ (unless otherwise stated,
particularly in relation to Fig. 8, time averages are evaluated from $T=0$,
i.e. $t=t_{\rm diff}=1~{}\rm Gyr$, to the end of the simulation at $t=13.7\rm
Gyr$). For $r_{\rm sat}\gtrsim r_{s}$, $\langle M(<r)\rangle$ is close to the
initial mass at $T=0$, since the dark matter mass inside it is essentially
conserved. We thus write, in terms of the initial mass $M_{0}(<r)$,

$E_{\rm in}(<r_{\rm sat})=10^{9}~{}M_{\odot}~{}{\rm
kpc^{2}~{}Gyr^{-2}}~{}\frac{M_{0}(<r_{\rm
sat})}{10^{8}M_{\odot}}~{}\frac{T}{\rm Gyr}~{}\times\\\
\left(\frac{n}{2.5}\right)^{-1}\left(\frac{v_{r}}{30~{}{\rm
km/s}}\right)^{-1}\left(\frac{\mathcal{P}(k_{m})}{3{\rm~{}kpc}^{3}}\right)\left(\frac{\rho_{0}~{}(r_{\rm
sat})}{10^{6}~{}M_{\odot}~{}{\rm kpc}^{-3}}\right)^{2},$ (7)

where we have chosen $v_{r}$ in accordance with the discussion of the previous
subsection, and inserted values characteristic of the initial halo mass and
average gas density (as defined above) around $r_{s}$. More precisely, the
numbers correspond to values at $r=r_{\rm sat}=1.5r_{s}\simeq 1.3~{}{\rm
kpc}$, but we note that the estimate of the total energy input is not
sensitive to the choice of $r_{\rm sat}$; as the product
$M_{0}(<r)\rho_{0}^{2}(r)$ varies slowly with radius around $r_{s}$ (indeed,
more generally, $\langle M(<r)\rangle\rho_{0}^{2}(r)$ varies by at most a
factor of 2 in the range $0.1~{}r_{s}\leq r\leq 2r_{s}$). The power spectrum
parameters are estimated from the results in figures 4 and 6 (a larger
$\mathcal{P}(k_{m})=4{\rm kpc}^{3}$ gives the same result if $v_{r}=40{\rm
km/s}$, which would take into account enhancement in the speeds of halo halo
particles, relative to the gas, due to their own motion).

We now wish to compare the prediction in equation (7) with the actual energy
input inferred from the simulation. For this purpose, we define the quantity

$E_{r}=\frac{m}{2}\left(\sum_{i}v_{i}^{2}+\Phi_{i}\right),$ (8)

where $m$ is the halo particle mass in the simulation, $v_{i}$ is the speed of
particle $i$ and $\Phi_{i}$ is the Newtonian potential at its location. The
summation is evaluated over all particles within radius $r$. Increase in this
quantity measures the energy input inside radius $r$, provided there is no
dark matter mass or energy outflow from that radius and the mass distribution
beyond it remains constant (so that changes in the potential are solely due to
modification of the mass distribution within $r$) 666The potential $\Phi_{i}$
at particle $i$ is due to other halo particles; we ignore changes in the
potential arising directly from the gas flows, as these constitute a
fluctuating contribution with small average, particularly after the initial
expulsion episode during the first few hundred Myr. We may expect that these
conditions hold, at all radii, for sufficiently small times, as the kinetic
energy acquired by the halo particles is still being converted into changes in
the mass distribution and potential, a process slower than that characterising
the initial energy input. The conditions should hold for longer times as $r$
increases; and we expect that beyond $r_{\rm sat}\gtrsim r_{s}$, the energy
input saturates to a specific value regardless of radius, and little energy is
transferred beyond that radius, so the change in $E_{r}$ from its initial
value (at $T=0$) corresponds to the total energy input from the fluctuating
gas. In particular, if that input may be estimated theoretically using
equation (7), then one expects $\Delta E_{r_{\rm sat}}\simeq E_{\rm
in}(<r_{\rm sat})$.

These expectations are confirmed in Fig. 7. At all radii, we initially find a
general linear increase in $E_{r}$, as expected from a steady diffusive
process, with stochastic variations on this general trend. At smaller radii,
the lines clearly flatten after a few Gyr, as a result of energy and mass
outflow onto larger $r$. A saturation radius $r_{\rm sat}\gtrsim r_{s}$,
beyond which the lines converge, can be defined with value consistent with
$r_{\rm sat}=1.5~{}r_{s}\simeq 1.3~{}{\rm kpc}$ as assumed in equation (7).
The slope derived from that relation is also consistent with that derived
directly from the simulation, as indicated by the dashed line.

The above suggests that the total energy input within $r_{\rm sat}$ may be
adequately described by assuming a global energy transfer mechanism,
parameterised by the average gas density and halo mass at $r_{\rm sat}$
inferred from the simulation; and with $v_{r}$ and density contrast power
spectrum parameters also consistent with those found in the simulated
hydrodynamics.

We now wish to ask to what extent this global approximation holds in general;
for, as the time averaged gas density inside $r_{s}$ is not strictly constant,
energy transfer may, in principle, depend instead on the local density. This
is a possibility akin to invoking the purely local approximation when
evaluating the effect of two body relaxation in stellar dynamics, by using
Chandrasekhar’s formula and plugging in the local stellar density. This is in
principle plausible, but more difficult to interpret in the present case;
whereas in two body relaxation logarithmic intervals in impact parameters (or
spatial scales) equally contribute to the fluctuations leading to relaxation,
here fluctuations at the largest scales (characterized by
$\mathcal{P}(k_{m})$) are more important.

Global energy transfer implies the validity of equation (6). Local energy
transfer, on the other hand, requires the evaluation of the time average
$\langle\rho_{g}(r,t)\rangle^{2}(r)$ over spherical shells centered at $r$,
and integrating it over the dark matter mass in those shells, $4\pi
r^{2}\rho_{d}(r)$ (assuming a fixed $v_{r}$). The right panel of Fig. 8
(dotted line) shows that, when assuming power spectrum parameters and velocity
$v_{r}$ consistent with the simulation, this vastly overestimates the energy
input. The non-local assumption — with energy input rate within radius $r$
scaling simply as $\rho_{0}^{2}M(<r)$ as implied by equation (6) — results in
much better agreement, with a moderate discrepancy at small radii that may at
least partly be accounted for by an increase of $v_{r}$ with radius (as
discussed above and expected from Fig. 13), which we do not take into account
here for simplicity.

The global approximation, with total energy input within radius $r$ simply
proportional to $M(<r)$, clearly cannot remain valid as $r\rightarrow r_{\rm
sat}$, when the gas density and fluctuations rapidly decrease. Indeed, the
difference between $\rho_{0}(r_{s})$ that fits the energy transfer rate in
Fig. 8 and $\rho_{0}(r_{\rm sat})$ used in Fig. 7, principally reflects the
gradual saturation process; the fit in Fig. 7 effectively assumes sudden
saturation, while in reality the process is gradual. To take this into account
one may invoke an intermediate regime, between the purely local and purely
global energy transfer limits. We do this by defining

$\rho_{0}(r,r_{\rm av})=\langle~{}\rho_{g}({\bf r},~{}{\bf r}-{\bf
r}_{g})~{}\rangle_{r,|{\bf r}-{\bf r}_{g})|<r_{\rm av},T},$ (9)

where the average over the gas density is evaluated over time at points ${\bf
r}_{g}$ inside spheres with centres at radial coordinate $r$ and radii $r_{\rm
av}$. Thus in this context, $\rho_{0}(r_{\rm sat})=\rho_{0}(0,r_{\rm sat})$
and $\rho_{0}(r_{s})=\rho_{0}(0,r_{s})$. The results for $\rho_{0}(r,r_{s})$
are shown by the dashed dotted line in the right hand panel of Fig. 8. They
suggest that the energy transfer process is best considered as non-local, with
a range $\sim$$r_{s}$.

Finally, we note that although the energy is assumed to be transferred to halo
particles initially as kinetic energy, due to modification in their velocities
(as given by equation 5), the changes eventually affect the self consistent
potential. The resulting average gain in total energy per unit mass turns out
to be amenable to estimation from a low energy cutoff that appears in the
phase space distribution function. In Appendix B, we show how this can be
related to the energy input calculated here, and connected to the change in
the potential that accompanies core formation.

Figure 9: Mass migration and core formation. The left panel shows the ratio
between the halo mass within radii $r$ (in kpc) and the corresponding initial
mass at $T=0$ ($t=t_{\rm diff}=1{\rm Gyr}$, as discussed in relation to Fig.
7). The dashed lines are exponential fits, $\propto\exp(\alpha T)$, with
numbers corresponding to $\alpha(r)$ in ${\rm Gyr}^{-1}$. The exponential
decay may be derived from a simple theoretical model for the mass transfer; in
its context $\alpha$ is predicted to scale with the inverse of an energy
relaxation time (obtained using Eqs. 14. 10, and 16). The right panel shows
the corresponding evolution in density. Theoretical expectations are shown by
the dashed lines. They are obtained by differentiating equation (15), starting
from an NFW fit to the dark matter density in the simulation at $T=0$. To
reduce uncertainties arising from fluctuations in the density profile around
$T=0$, we average simulation outputs over the range $0.5~{}{\rm Gyr}\leq t\leq
1.5~{}{\rm Gyr}$ (as in Fig. 7). Model predictions are shown at $T=0$ (fitting
the averaged profile), and at $2,4,6,8,10$ and 12 Gyr, corresponding to times
$1,3,5,7,9,11$ and 13 Gyr.

###  4.3 Mass migration and core formation

The energy input from the fluctuating gas leads to the migration of halo
particles from the inner radii, which decreases the enclosed mass, as shown in
Fig. 9, left panel. Straight lines on the log linear scale suggest a general
exponential decrease in mass with time. A full examination of its origin in
the context of a diffusion model would require a full Fokker Planck
formulation, using the full (first and second order) diffusion coefficients
and explicitly including changes in the potential due to the evolving mass
distribution. Here we proceed heuristically, leaving the examination of the
reasons why the drastic simplifications invoked seem to work to Appendix C.

We suppose that the mass flux across energy surface $E$ changes the mass
within it through the first order energy diffusion coefficient

$D[\Delta E]=\frac{1}{M(<r_{\rm sat})}\frac{E_{\rm in}(r_{\rm sat})}{T},$ (10)

describing the average rate of change of halo particle energy per unit mass
due the gas fluctuations. At the energy input rate inferred from Fig. 7, this
is of the order of $10~{}{\rm kpc^{2}~{}Gyr^{-3}}$ in the simulation at hand
(which is consistent with theoretical expectations, given the power spectrum
and density inferred from simulation, cf. equation 7, but is actually
empirically independent of the model).

The change in mass of particles with energy less than $E$ is

$\frac{\partial M(<E)}{\partial t}~{}\approx-\frac{\partial M(<E)}{\partial
E}~{}D~{}[\Delta E],$ (11)

where $M(<E)$ is the mass in halo particles with energy less than $E$ and
$\partial M(<E)/\partial E$ is the mass-weighed differential energy
distribution (Binney & Tremaine, 2008). Furthermore, we approximate this
latter quantity by its average within $E$, such that

$\frac{\partial M(<E)}{\partial
E}\approx~{}a_{p}~{}\frac{M(<E)-M(0)}{E-E(0)},$ (12)

where $E(0)$ is the energy of a particle at rest at the centre of the
potential (so $E(0)=\Phi(0)$ and $M(0)=0$), and $a_{p}$ is a constant
numerical factor of order 1. This holds exactly inside pure power-law cusps:
e.g. as derived in Appendix C, one finds in such cases that
$a_{p}=(3+\gamma)/(2+\gamma)$ for $\rho\propto r^{\gamma}$ (thus $a_{p}=2$ for
a $1/r$ cusp, tends towards $1.5$ for flatter profiles, and is larger for
steeper ones (finally diverging in the case of the singular isothermal sphere,
where a potential zero point at the centre cannot be fixed as above due to
logarithmic divergence).

With this in mind we will assume that

$\frac{\partial M(<E)}{\partial
t}~{}=-a~{}\frac{M(<E)-M(0)}{E-E(0)}D~{}[\Delta E],$ (13)

where the parameter $a$ absorbs all of the major approximations involving the
use of the first order diffusion coefficient (equation 11), the association of
the differential energy distribution with its average value (equation 12), and
the neglect of the effect of the evolution of the self consistent potential.

In Appendix C, we discuss how such a formulation may be justified if the
differential energy distribution may be approximated by a power law; and how,
when one includes both diffusion coefficients, a zero Fokker Planck flux may
be attained if the distribution takes an exponential form.

We now define an energy relaxation time as

$t_{\rm relax}(E)=D[\Delta E]^{-1}~{}[E-E(0)].$ (14)

For timescales small compared to the initial relaxation time (with $E$ and
$E_{0}$ defined at the initial time), one may take this timescale to be
constant, evaluated by inserting the initial values of $E$ and $E_{0}$ in the
above equation. We will assume that this is the case here (as we will see
below, the energy changes are small enough, during timescales of interest, to
render this approximation reasonable). Then, assuming that equation (11) is
applicable in the diffusion limit, solving it (starting at time $T=0$) results
in

$M(<E)=M_{0}(<E)~{}\exp\left[-a~{}T/t_{\rm relax}(E)\right],$ (15)

where $M_{0}$ refers to the mass at time $T=0$.

Within this picture, the numbers on the lines in the left panel of Fig. 9,
denoted by $\alpha$, should correspond to $a/t_{\rm relax}(E)$, if $E$ is
associated with particle energies inside the chosen radii. To make this
correspondence, we fix $E=\langle E\rangle(r)$ to be the average specific
energy at radius $r$ in the initial profile. As a first approximation, we
simply associate this with the initial NFW profile of the configuration. The
specific potential energy $\Phi(r)-\Phi(0)$ in this case is $4\pi
G\rho_{s}r_{s}^{2}[1-\ln(1+x)/x]$, with $x=r/r_{s}$, and $r_{s}$ and
$\rho_{s}$ as in Table 1. To obtain a simple form for the kinetic energy we
make use of established empirical relations according to which the pseudo
phase space density varies with radius approximately as $\rho/\langle
v^{2}\rangle^{3/2}\sim r^{-1.875}$ (Taylor & Navarro, 2001), and normalize the
velocity variance to its value around $r=r_{s}$, such that $\langle
v^{2}\rangle\simeq 470~{}{\rm km^{2}/s^{2}}$. One may then write

$E-E(0)\simeq 4\pi
G\rho_{s}r_{s}^{2}\left[1+\frac{2^{4/3}}{10}\frac{\left(x^{0.875}\right)^{2/3}}{(1+x)^{4/3}}-\frac{\ln(1+x)}{x}\right].$
(16)

Using this formula in conjunction with equation (14), and using $a=3.4$, we
find $a/t_{\rm relax}=0.034,0.047,0.071,0.11,0.17~{}{\rm Gyr}^{-1}$ at
$r=1,0.5,0.25,0.125,0.0625~{}{\rm kpc}$ respectively. These numbers agree to
better than $20\%$ with the values of $\alpha$ indicated on the fitting lines
in Fig. 9 (left panel). The corresponding energy relaxation times are about
$100,72.3,47.9,30.9$ and $20~{}{\rm Gyr}$.

A more careful analysis, taking into account that the quasi-steady state
diffusive process is well established only after $t\gtrsim t_{\rm diff}=1\rm
Gyr$, suggests $a\simeq 2.5$ (cf. below). The corresponding energy relaxation
times are therefore lower by a factor of about 1.4 than those quoted above.
Even then, these relaxation times are still significantly larger than what is
obtained from the velocity variance (equation 4). This is because $E-E(0)$ is
generally larger than the average kinetic energy (from equation 16, by a
factor of about four at $r_{s}$). Indeed, changes in kinetic energy lead to
relatively little change in total energy. In fact the relative change in
$E_{r}$ at 1 kpc (as calculated from equation 8, starting at $T=0$), is
$\Delta E_{r}/E_{r}\simeq 0.2$, over 12.7 Gyr. Initially the rate of relative
change in energy is much larger at smaller radii (as illustrated by the
differences in $E_{\rm in}/T$ inferred from Figs. 7 and 8). But it
subsequently saturates, as the energy is redistributed in the system, with
mass and energy flowing towards outer radii (Fig. 8, left panel). When thus
redistributed, the modest changes in total energy lead to a modification of
the self consistent potential, including the minimal possible energy in it, as
discussed in Appendix B.

This leads to core formation, as observed in the right panel of Fig. 9, There,
we note that we find no evidence of any accelerated core formation, relative
to the velocity relaxation time (equation 4), as found in El-Zant et al.
(2016), when non-spherical modes were used in conjunction with the Hernquist-
Ostriker code.

Having obtained the evolution of the enclosed mass within a given radius, it
is also possible to define a theoretical density $\rho=r^{-2}dM/dr$ and derive
a closed form formula for it using Eqs. (15), (14), and (16). The result may
then be compared with the dark matter density evolution in the simulation.
More care is however required here in defining the initial energies in
equation (16), as we must start our comparison at $T=t-t_{\rm diff}=0$, i.e.,
when the diffusion limit at the basis of our model is valid (cf. the
discussion in relation to Fig. 7). For this purpose, we average simulation
outputs in the time range $0.5~{}{\rm Gyr}\leq t\leq 1.5~{}{\rm Gyr}$, and
then fit the resulting dark matter density with an NFW profile. The
corresponding parameters are $r_{s}=1.17~{}{\rm kpc}$ and $\rho_{s}=2.7\times
10^{7}~{}M_{\odot}/{\rm kpc^{3}}$. When the latter is adjusted by adding a gas
fraction of about $0.075$ (assumed for simplicity to be also NFW with same
$r_{s}$, but taking into account the initial central gas expulsion), the
multiplicative factor $\rho_{s}r_{s}^{2}$ in front of the bracket in equation
(16) remains approximately the same as in the case of the $t=0$ profile. After
some trials, we found that values of $a$ in the range 2.5 to 2.8 provide
reasonable approximations to the density evolution inferred from the
simulation (perhaps surprisingly so, given the various simplifying assumptions
of our mass transfer model). This is illustrated in the right panel of Fig. 9,
where we compare the density evolution expected from our model (fixing
$a=2.5)$ with the results from the simulation.

Note that, other than the value of $a$, which is theoretically fairly tightly
constrained by the arguments of Appendix C (cf. equation 39), there are no
free parameters in these fits. The energy distribution of halo particles at
$T=0$, and the energy transfer rates from gas to halo, are of course
characteristic of a given system. They will thus vary from simulation to
another, but they are entirely fixed for any one system (within numerical
uncertainties in their inference from the simulation).

##  5 Conclusion

Potential fluctuations from feedback driven gas can ‘heat’ halo cusps, turning
them into cores. This work aimed at quantifying this, by measuring the gas
fluctuations and tracking the way they transfer energy to the central halo,
forcing the outward migration of the dark matter. The interpretive framework
we use is a model, first outlined in El-Zant et al. (2016), which predicts
that these processes principally depend on the amplitude of the fluctuations,
as measured by the normalisation of their power spectrum, and the average gas
density. The result being a standard diffusion process, characterised by a
linear temporal increase in velocity variance and energy of halo particles as
a result of their interaction with the fluctuating gas field. As shown
previously (El-Zant et al., 2020), in this picture the effect of the
fluctuations reduces to standard two body relaxation in the case of a white
noise power spectrum. It may thus, in this limit, also describe halo heating
via dynamical friction from a system of compact, monolithic massive clumps
moving among much lighter dark matter particles.

To test this interpretative framework, we measure the density fluctuations in
feedback-driven gas from a full hydrodynamic simulation of a model dwarf
galaxy (Figs. 1-3), and obtain their density contrast power spectrum (Fig. 4).
We find that the spectrum follows a power-law in wave number ($\propto
k^{-n}$), with an exponent $2\lesssim n\lesssim 3$, for the whole period of
the simulation (spanning a Hubble time). Although the time-averaged density of
the driven gas varies with radius, the variation is modest inside the initial
NFW scale length, relative to a sharp drop outside. This suggests a
characteristic density that may be used as input for the model. We also
examine the velocity distribution of the gas and find that it is approximately
fit by Maxwellians at larger radii (Appendix A). The kinetic energy power
spectrum is close to the Kolmogorov form over large range of scales (Fig. 5).
This reinforces the picture of a fully turbulent medium.

With the input parameters directly measured from the simulation, we use our
model to calculate the total energy transfer rate from fluctuating gas to the
central halo. The result is compared with the actual energy input to the
system of halo particles, as directly inferred from the simulation (Fig. 7).
This is found to generally agree — displaying the general linear increase,
expected of a steady diffusion process — over almost a Hubble time. We also
examine the radial distribution of the energy input rate (Fig. 8). We find
that the energy transfer is indeed much better approximated as a global rather
than local process; in the sense that it depends on the average gas energy
density within the core region rather than the local density at each radius.

The energy is initially transferred from the gas as kinetic energy to
individual halo particles, but it is then redistributed through the self
consistent gravity of the system, a process through which the core replaces
the cusp. The process comes with a low energy cutoff in the halo phase space
distribution function, as particles migrate to higher energy levels. It is
straightforward to link the level of this cutoff to the total energy input
from the gas, and to the resulting change in halo gravitational potential that
comes with core formation (Appendix B).

The energy flow upwards is accompanied by a mass flow outwards. We empirically
find an exponential decrease in halo mass with time, within a given radius, in
the initial cusp. We then devise a simple approximate description of the mass
flow, based on our model, from which the exponential form may be inferred. In
this context, the exponential decay time scales with (and is of the order of)
the local (energy) relaxation time, and the evolution of the corresponding
theoretical density profile mimics that in the simulation (Fig. 9).

Strictly speaking, energy transfer via a standard diffusion process requires
the gas force fluctuations to be normally distributed. We have verified this
to be the case to a good approximation (even if the larger density
fluctuations can be lognormal rather than Gaussian at small radii, cf.
Appendix A). Another assumption of the model is the use of the ‘sweeping’
approximations of turbulence theory, whereby the statistics of the spatial
fluctuations are transferred to the time domain through large scale flows. The
general consistency of the spatial power spectrum parameters with those that
fit the mass dispersion in the time domain (Fig. 6), confirm that this may be
an applicable approximation.

In general, the assumptions and predictions of the theoretical framework seem
vindicated for the model galaxy studied here. This suggests a remarkably
concise description of halo core formation from gaseous fluctuations,
summarizing the effect of much complex ‘gastrophysics’ in terms of the two
principal parameters of gas density and fluctuation levels. Obvious extensions
include considering different galaxy masses, as well as verifying that the
model also predicts a lack of core formation in simulations that do not
produce them. Success there would help delineate particular elements in the
physical ‘subgrid’ input that are most crucial in producing the required
conditions for core formation, and examining how they compare with
observations.

It may also be possible to test the predictions regarding the level of the
fluctuations required for core formation directly from observations. In
particular, it has already been possible to derive surface density power
spectra for larger,and relatively quiescent galaxies. Estimating these for
actively star forming galaxies, and calibrating them with three dimensional
density contrast power spectra entering into our model calculations, may
provide a direct test of our picture of core formation through gaseous
fluctuations.

## Acknowledgements

We would like to thank the referee for a careful reading, thorough report and
constructive criticism that helped improve our manuscript.

## Data Availability

The simulation data underlying this article will be shared on reasonable
request.

## References

  * Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  * Aguilar & White (1985) Aguilar L. A., White S. D. M., 1985, ApJ, 295, 374
  * Bar-Or et al. (2013) Bar-Or B., Kupi G., Alexander T., 2013, ApJ, 764, 52
  * Bar et al. (2021) Bar N., Blum K., Sun C., 2021, arXiv e-prints, p. arXiv:2111.03070
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic dynamics
  * Bouché et al. (2022) Bouché N. F., et al., 2022, A&A, 658, A76
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, Annual Review of Astronomy and Astrophysics, 55, 343
  * Burkert (2000) Burkert A., 2000, ApJl, 534, L143
  * Burkert (2020) Burkert A., 2020, arXiv e-prints, p. arXiv:2006.11111
  * Coles & Jones (1991) Coles P., Jones B., 1991, MNRAS, 248, 1
  * Collins & Read (2022) Collins M. L. M., Read J. I., 2022, Nature Astronomy, 6, 647
  * Dekel et al. (2021) Dekel A., et al., 2021, MNRAS, 508, 999
  * Del Popolo & Le Delliou (2017) Del Popolo A., Le Delliou M., 2017, Galaxies, 5, 17
  * Deng et al. (2018) Deng H., Hertzberg M. P., Namjoo M. H., Masoumi A., 2018, Phys. Rev. D, 98, 023513
  * Di Cintio et al. (2014) Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415
  * Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
  * El-Zant (2008) El-Zant A. A., 2008, ApJ, 681, 1058
  * El-Zant et al. (2001) El-Zant A., Shlosman I., Hoffman Y., 2001, ApJ, 560, 636
  * El-Zant et al. (2004) El-Zant A. A., Hoffman Y., Primack J., Combes F., Shlosman I., 2004, ApJ, 607, L75
  * El-Zant et al. (2016) El-Zant A. A., Freundlich J., Combes F., 2016, MNRAS, 461, 1745
  * El-Zant et al. (2020) El-Zant A. A., Freundlich J., Combes F., Halle A., 2020, MNRAS, 492, 877
  * Errani et al. (2023) Errani R., Navarro J. F., Peñarrubia J., Famaey B., Ibata R., 2023, MNRAS, 519, 384
  * Evans (1994) Evans N. W., 1994, MNRAS, 267, 333
  * Falkovich (1994) Falkovich G., 1994, Physics of Fluids, 6, 1411
  * Famaey et al. (2020) Famaey B., Khoury J., Penco R., Sharma A., 2020, J. Cosmology Astropart. Phys., 2020, 025
  * Flores & Primack (1994) Flores R. A., Primack J. R., 1994, ApJl, 427, L1
  * Frenk & White (2012) Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507
  * Freundlich et al. (2020a) Freundlich J., Dekel A., Jiang F., Ishai G., Cornuault N., Lapiner S., Dutton A. A., Macciò A. V., 2020a, MNRAS, 491, 4523
  * Freundlich et al. (2020b) Freundlich J., et al., 2020b, MNRAS, 499, 2912
  * Gilman et al. (2020) Gilman D., Birrer S., Nierenberg A., Treu T., Du X., Benson A., 2020, MNRAS, 491, 6077
  * Governato et al. (2015) Governato F., et al., 2015, MNRAS, 448, 792
  * Grisdale et al. (2017) Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093
  * Hernquist & Ostriker (1992) Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375
  * Inoue et al. (2016) Inoue S., Dekel A., Mandelker N., Ceverino D., Bournaud F., Primack J., 2016, MNRAS, 456, 2052
  * Iršič et al. (2017a) Iršič V., et al., 2017a, Phys. Rev. D, 96, 023522
  * Iršič et al. (2017b) Iršič V., Viel M., Haehnelt M. G., Bolton J. S., Becker G. D., 2017b, Phys. Rev. Lett., 119, 031302
  * Kim et al. (2018) Kim S. Y., Peter A. H. G., Hargis J. R., 2018, Phys. Rev. Lett., 121, 211302
  * Kochanek & White (2000) Kochanek C. S., White M., 2000, ApJ, 543, 514
  * Koposov et al. (2008) Koposov S., et al., 2008, ApJ, 686, 279
  * Kraichnan (1964) Kraichnan R. H., 1964, Physics of Fluids, 7, 1723
  * Kritsuk et al. (2007) Kritsuk A. G., Norman M. L., Padoan P., Wagner R., 2007, ApJ, 665, 416
  * Li et al. (2022) Li Z., Dekel A., Mandelker N., Freundlich J., François T., 2022, arXiv e-prints, p. arXiv:2206.07069
  * Macciò et al. (2012) Macciò A. V., Paduroiu S., Anderhalden D., Schneider A., Moore B., 2012, MNRAS, 424, 1105
  * Moore (1994) Moore B., 1994, New Astron.t, 370, 629
  * Nadler et al. (2020) Nadler E. O., et al., 2020, arXiv e-prints, p. arXiv:2008.00022
  * Nadler et al. (2021) Nadler E. O., Birrer S., Gilman D., Wechsler R. H., Du X., Benson A., Nierenberg A. M., Treu T., 2021, ApJ, 917, 7
  * Navarro et al. (1996a) Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72
  * Navarro et al. (1996b) Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, 563
  * Newton et al. (2021) Newton O., et al., 2021, J. Cosmology Astropart. Phys., 2021, 062
  * Nusser & Silk (2022) Nusser A., Silk J., 2022, MNRAS, 509, 2979
  * Ogiya & Burkert (2015) Ogiya G., Burkert A., 2015, MNRAS, 446, 2363
  * Ogiya & Nagai (2022) Ogiya G., Nagai D., 2022, MNRAS, 514, 555
  * Oman et al. (2019) Oman K. A., Marasco A., Navarro J. F., Frenk C. S., Schaye J., Benítez-Llambay A., 2019, MNRAS, 482, 821
  * Orkney et al. (2021) Orkney M. D. A., et al., 2021, MNRAS, 504, 3509
  * Peñarrubia (2019a) Peñarrubia J., 2019a, MNRAS, 484, 5409
  * Peñarrubia (2019b) Peñarrubia J., 2019b, MNRAS, 490, 1044
  * Peñarrubia et al. (2012) Peñarrubia J., Pontzen A., Walker M. G., Koposov S. E., 2012, ApJ, 759, L42
  * Pontzen & Governato (2014) Pontzen A., Governato F., 2014, Nature, 506, 171
  * Read & Erkal (2019) Read J. I., Erkal D., 2019, MNRAS, 487, 5799
  * Read & Gilmore (2005) Read J. I., Gilmore G., 2005, MNRAS, 356, 107
  * Read & Trentham (2005) Read J. I., Trentham N., 2005, Philosophical Transactions of the Royal Society of London Series A, 363, 2693
  * Read et al. (2006) Read J. I., Wilkinson M. I., Evans N. W., Gilmore G., Kleyna J. T., 2006, MNRAS, 367, 387
  * Read et al. (2016) Read J. I., Agertz O., Collins M. L. M., 2016, MNRAS, 459, 2573
  * Read et al. (2019) Read J. I., Walker M. G., Steger P., 2019, MNRAS, 484, 1401
  * Salucci (2019) Salucci P., 2019, Astron. Astrophys. Rev., 27, 2
  * Salucci et al. (2020) Salucci P., Turini N., di Paolo C., 2020, Universe, 6, 118
  * Schmidt et al. (2004) Schmidt W., Hillebrandt W., Niemeyer J. C., 2004, arXiv e-prints, pp astro–ph/0407616
  * Simon & Geha (2007) Simon J. D., Geha M., 2007, ApJ, 670, 313
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters
  * Taylor (1938) Taylor G. I., 1938, Proceedings of the Royal Society of London Series A, 164, 476
  * Taylor & Navarro (2001) Taylor J. E., Navarro J. F., 2001, ApJ, 563, 483
  * Tennekes (1975) Tennekes H., 1975, Journal of Fluid Mechanics, 67, 561
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
  * Turner et al. (2021) Turner H. C., Lovell M. R., Zavala J., Vogelsberger M., 2021, MNRAS, 505, 5327
  * Warren et al. (1992) Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405
  * Widrow (2000) Widrow L. M., 2000, ApJS, 131, 39
  * Yang et al. (2022) Yang D., Nadler E. O., Yu H.-b., 2022, arXiv e-prints, p. arXiv:2211.13768
  * Yu et al. (2021) Yu X., Bian F., Krumholz M. R., Shi Y., Li S., Chen J., 2021, MNRAS, 505, 5075

##  Appendix A Statistics of the gas fluctuations

For the description of the effect of the gas fluctuations on the halo
particles in terms of a diffusion process, leading to relaxation timescales of
the form of Eq (2), to be complete, their statistics should be entirely
described by averages and dispersions, as in a Gaussian random process. We
wish to examine whether this is the case with the gas fluctuations considered
in this work.

Figure 10: Probability distribution function (PDF) of the mass fluctuations
within different radii $r$. The fluctuations are well-fit by Gaussians for
$r\geq 0.50.5{\rm kpc}$ and by a lognormal distribution below. The fits are
shown by the red dashed-dotted lines Figure 11: PDF of the fluctuations in the
gas gravitational potential and the magnitude of the fluctuating force inside
$r=1.0\;{\rm kpc}$. The fits (dashed-dotted red lines) are Gaussian for the
potential fluctuations and Maxwellian (i.e., assumed to emanate from a
Gaussian distribution in each degree of freedom, as suggested by Fig. 12) for
the absolute value of the force. Figure 12: PDF of X,Y,Z components of force
fluctuations with Gaussian fits (red dashed-dotted lines). Figure 13: PDF of
gas speeds, with Maxwellian fits at larger radial bins (red dashed-dotted
lines).

Fig. 10 shows the statistics of the mass fluctuations in three different bins
in the region where the core forms. For $r\gtrsim 0.5r_{s}$ the fluctuations
are well fit by Gaussians. Inside $0.5r_{s}$ however they are better described
by a lognormal distribution. This is not surprising, given that the mass
fluctuations at these scales are quite large (as may also be probed from Fig.
6), and the density is bounded (by zero) from below but not necessarily from
above. A similar situation, with analogous description, also appears in the
description of large scale structure at large density contrasts (Coles &
Jones, 1991).

The density and mass fluctuations do not however affect the dynamics directly.
Rather what affects CDM trajectories are potential and force fluctuations. As
the potentials and forces are related (through the Poisson equation) to the
density through integrals, their fluctuations are generally milder in
gravitational systems. We now examine these.

In order to obtain the gravitational field exerted by the gas perturbations on
DM particles, we first sample random points inside a spherical region centred
at the DM halo’s centre of potential. Then we calculate gas mass fluctuations
defined as $\delta M=M-\langle M\rangle$, where $\langle M\rangle$ is the
average gas mass profile from the gas centre of mass. The gravitational
potential of the gas perturbations is then obtained using:

$\phi=\sum\frac{\delta M}{d},$ (17)

where $d$ is the distance between a DM particle and gas points.

The gravitational force is directly obtained using the relations:

$\displaystyle F_{x}$ $\displaystyle=$ $\displaystyle\sum\frac{\delta
M}{d^{2}}\frac{x-x_{c}}{d},$ (18) $\displaystyle F_{y}$ $\displaystyle=$
$\displaystyle\sum\frac{\delta M}{d^{2}}\frac{y-y_{c}}{d},$ (19)
$\displaystyle F_{z}$ $\displaystyle=$ $\displaystyle\sum\frac{\delta
M}{d^{2}}\frac{z-z_{c}}{d},$ (20)

where $x_{c},y_{c},z_{c}$ are the coordinates of the gravitational field
centre. The force amplitude is given as
$F=\sqrt{F_{x}^{2}+F_{y}^{2}+F_{z}^{2}}$.

The results are shown in figures 11 and 12. Apart from modest departures —
some long tail effects, as well as some skewness in the distribution of
potential fluctuations — the fluctuating potential and force field is indeed
Gaussian. This is particularly clear in the distribution of the force
components shown in Fig. 12.

We note nevertheless that the modest but non-negligible non-Gaussianities
imply that higher order moments of velocity perturbations ($\langle(\Delta
v)^{m}\rangle$, $m>2$) may contribute to the effect of gas fluctuations; the
departure of halo particle trajectories from their collisionless counterparts
may involve corrections beyond the increase in velocity variance in the
standard diffusion limit as described by (5). In standard two body relaxation,
the long velocity tails arise from strong, close encounters, and are
alleviated when the force perturbations are caused by extended structures
(Peñarrubia, 2019b). In our case, for power spectra $\mathcal{P}(k)\sim
k^{-n}$ with $n>0$), the relaxation process is expected to be dominated by
large scale modes (analogous to distant encounters the standard case); as
equation (2) in fact emanates from an effective theory in which a higher scale
$k_{x}$ is eliminated when $n>0$ and $k_{x}\gg k_{m}$ (cf. El-Zant et al.,
2016 Section 2). The long tails here are thus likely to be due to halo
particle interactions with rare large scale fluctuations (with large
perturbing forces; cf. equation 33 of El-Zant et al., 2016). Their effect may
be limited however by the flattening of the power spectrum at $k<k_{m}$ (Fig.
4), and slow temporal variation of such modes (cf. Footnote 5). Verifying
this, and comparing with theoretical formulations of anomalous diffusion (e.g.
Bar-Or et al., 2013), may be of interest for further study of the dynamical
effects of stochastic fluctuations with power law spectra considered here.

Finally we show in Fig. 13 the distribution of speeds of the gas fluid
elements. This is calculated from the speeds in all gas cells within the
indicated radii and the result is averaged over all time snapshots (excluding
the first Gyr). At larger radial bins the speeds closely follow a Maxwellian
distribution, (but with long tail at high speeds), with an average value
saturating at $\simeq 30~{}{\rm km/s}$. The speeds are smaller as one moves
towards the centre.

##  Appendix B Energy input and cutoff in the halo distribution function

Figure 14: Dimensionless phase space distribution function $F$ in terms of
dimensionless energy $\mathcal{E}$, measured from the minimal energy in the
initial NFW potential. In the initial NFW cusp the phase space density of halo
particles continuously increases as the zero energy point of the NFW system
(at $\mathcal{E}=-1$) is approached. As energy is transferred from the
fluctuating gas, however, a low energy cutoff in the distribution emerges, as
particles migrate towards progressively higher energies. Figure 15: The ratio
$R=T_{r}/\Phi_{r}$, where $T_{r}=M\langle T_{i}\rangle$,
$\Phi_{r}=\frac{M}{2}\langle\Phi_{i}\rangle$, and $T_{i}$ and $\Phi_{i}$ are
halo particle specific kinetic energy and potential energies, and the averages
are taken inside radii $r$. Figure 16: Comparison of total energy transferred
from gas fluctuations to the halo (as in Fig. 7) with the average energy
gained by dark matter particles evaluated via the cutoff in their distribution
function $\mathcal{E}_{c}$. The latter is inferred from Fig. 14 and
transformed through the factor $\mathcal{R}$ from equation 25. Figure 17:
Comparison of the minimal dimensionless specific particle energy
$\mathcal{E}_{c}$ (corresponding to the energy cutoffs in Fig. 14) with that
of the dimensionless potential $\Psi=\Phi/4\pi G\rho_{s}r_{s}^{2}$ at $r=0$
(as in Fig. 14 the normalisation corresponds to the initial NFW parameters
listed in Table 1).

We obtain the phase space distribution $f$ of halo particles by using the
Eddington formula for a spherical density distribution with isotropic
velocities (Binney & Tremaine, 2008), by inserting numerical values for the
dark matter potential and density inferred from the simulation. In Fig. 14 we
plot the dimensionless distribution function of halo particles,
$F(\mathcal{E})=f/(4\pi G)^{3/2}r_{s}^{3}\rho_{s}^{1/2}$, where
$\mathcal{E}=E/4\pi G\rho_{s}r_{s}^{2}$, against $1+\mathcal{E}$, which
measures excess in energy above the minimum energy in the initial NFW system
(i.e with zero point at $\Phi_{\rm NFW}(r=0)$). Here one probes a cutoff in
the distribution function at progressively larger values of $1+\mathcal{E}$;
as the evolution progresses, halo particles migrate to higher energies as a
result of their interaction with the fluctuating gas. Lower energy levels are
emptied in the process.

We wish to relate this effect to the total energy input examined in Section 4,
which is initially pumped into the particle trajectories in the form of
kinetic energy. For this purpose, we may rewrite equation (8) in terms of the
average kinetic and potential energies of the halo particles,

$E_{r}=M\left(\langle T_{i}\rangle+\frac{1}{2}\langle\Phi_{i}\rangle\right),$
(21)

where the averages are taken over all particles within radius $r$ and
$M=M(<r)$ refers to the mass enclosed therein. Then we define $T_{r}=M\langle
T_{i}\rangle$, $\Phi_{r}=M\langle\Phi_{i}\rangle$ and
$R_{v}(r)=-T_{r}/\Phi_{r}$ such that

$E_{r}=\frac{M}{2}\left(1-R_{v}\right)\langle\Phi_{i}\rangle.$ (22)

Also, using $\langle E_{i}\rangle=\langle T_{i}\rangle+\langle\Phi_{i}\rangle$
and $\langle T_{i}\rangle/\langle\Phi_{i}\rangle=-R_{v}/2$,

$\langle E_{i}\rangle=\left(1-\frac{R_{v}}{2}\right)\langle\Phi_{i}\rangle.$
(23)

And thus,

$E_{r}=\frac{M}{2}~{}\left(\frac{1-R_{v}}{1-\frac{R_{v}}{2}}\right)~{}\langle
E_{i}\rangle.$ (24)

Fig. 15 shows the ‘virial ratio’ $R_{v}$ of the collection of halo particles
within radius $r$. Within $r_{s}$, it is seen to vary from $\sim$$0.1$, in the
inner region, to $\sim$$0.4$ at $r\simeq r_{s}$. The latter is the relevant
value if we are interested in evaluating the average energy input within the
initial NFW cusp.

To compare with the energy input inferred from Fig. 7, we first estimate the
cutoff $\mathcal{E}_{c}$ through the intersection of the distribution function
with a horizontal line drawn around the maximum values in the distribution
function in Fig. 14 (we don’t explicitly draw this; it can be considered to
correspond to the upper horizontal limiting line of the figure panel at
$F=10^{4}$).

As with Fig. 7, we subtract the average over $0.5\leq t\leq 1.5~{}{\rm Gyr}$
from the numbers thus inferred and start our comparison at $t=t_{\rm
diff}=1~{}{\rm Gyr}$ (which is again our $T=0$ point). We denote the result by
$\Delta\mathcal{E}_{r}$, and associate it with the average energy input per
particle $\langle\Delta E_{i}\rangle=(4\pi
G\rho_{s}r_{s}^{2})~{}\Delta\mathcal{E}_{r}$. The full transformation is then
$\Delta E_{r}=\mathcal{R}~{}\Delta\mathcal{E}_{r}$, with

$\mathcal{R}=2\pi
G\rho_{s}r_{s}^{2}~{}M~{}\left(\frac{1-R_{v}}{1-\frac{R_{v}}{2}}\right).$ (25)

Fig. 16 shows the result of matching the total energy input at the saturation
radial level (converging lines in Fig. 7), with
$\mathcal{R}~{}\Delta\mathcal{E}_{r}$, obtained as just described, assuming
$R_{v}=0.4$ and $M=M(<1{\rm kpc})$ in equation (25).

Finally we show that the dimensionless minimal cutoff energy scale
$\mathcal{E}_{c}$ in fact corresponds to the corresponding dimensionless
potential ($\Psi$) at $r=0$ (Fig. 17). This illustrates the process of energy
redistribution; even if the initial energy is transferred from the gas to the
halo particles in the form of kinetic energy, particles may still end up with
zero kinetic energy, but still with increased energy as a result of a
shallower self consistent potential well, corresponding to the formed core.

##  Appendix C Fokker Planck flux for power law and exponential differential
energy distributions

In this appendix we examine the three major simplifications made to arrive at
equation (13). Namely, that the differential energy distribution may be
approximated by its average value (equation 12); that the energy space flux of
dark matter particles, migrating into higher energy levels under the influence
of stochastic gas fluctuations, may be described in terms of the first order
diffusion coefficient (as in equation 11); and that the flux through energy
surface $E$ due to the self consistent evolution of the system may likewise be
absorbed in the parameter $a$.

In particular, we show how the assumption of absorbing the complexity of the
situation in a single parameter $a$ may be reasonable when the differential
energy distribution may be approximated by a power law (which we will find to
be the case in the region where the halo cusp is transformed into a core). We
also show that this parameter is positive, allowing for outward mass flow, for
systems with density decreasing with radius.

We finally discuss the qualitatively different behaviour at intermediate
radii, where the differential energy distribution may be approximated instead
by an exponential. In this latter case, the competing effects of the fluxes
associated with the first and second order diffusion coefficients may cancel
out, and the Fokker Planck flux may in fact vanish.

###  C.1 Approximating the differential energy distribution by its average

Consider a spherical density distribution $\rho_{p}\propto r^{\gamma}$, with
$-2<\gamma<0$. If the distribution is infinite in radial extent, the potential
$\Phi(r)$, defined in terms of the work done to take a particle to infinity,
diverges. But $\phi(r)\propto r^{-\beta}$, associated with the work done in
moving a particle from the centre to radius $r$, is well defined, and from the
Poisson equation $\gamma=-(\beta+2)$. One can also define the specific energy
$E_{p}=\phi(r)+\frac{1}{2}v^{2}$. In the notation used in this work, $E_{p}$
will correspond to $E-E_{0}$ when $\Phi(r)$ is well defined — as in the
central region of a density distribution that, though well approximated by
$\rho_{p}\sim r^{\gamma}$ in the central region, eventually falls off
sufficiently steeply with increasing radius. In this case
$\phi(r)=\Phi(r)-\Phi(0)$.

For such power law systems, the mass weighed distribution function of a
configuration with isotropic velocities can be written as
$g=g_{0}E_{p}^{2/\beta-1/2}$ and the density of states
$p=p_{0}E_{p}^{1/2-3/\beta}$, where $g_{0}$ and $p_{0}$ are positive constants
(Evans, 1994; El-Zant, 2008). In general, the density of states corresponds to
the area that encloses a phase space volume

$q(<E)=\frac{16\pi^{2}}{3}\int_{0}^{r_{\rm max}}[2(E-\Phi)]^{3/2}r^{2}dE,$
(26)

where $r_{\rm max}(E)$is the maximum radius that a particle with specific
energy $E$ can reach. Thus $p(E)=\partial q(<E)/\partial E$.

The differential energy distribution is thus $\partial M(<E_{p})/\partial
E_{p}=pg=p_{0}g_{0}E_{p}^{-1/\beta}$, which for a power law distribution
embedded in a larger halo tends to

$\frac{\partial M(<E)}{\partial E}=pg=p_{0}g_{0}(E-E_{0})^{-1/\beta},$ (27)

from where it also follows that

$\frac{\partial M(<E)}{\partial E}=(1-1/\beta)\frac{M(<E)}{E-E_{0}}.$ (28)

Therefore, for such systems, equation (12) becomes exact with
$a_{p}=1-1/\beta=(3+\gamma)/(2+\gamma)$, as mentioned in Section 4.3.

###  C.2 Use of first order diffusion coefficient

As the process of core formation proceeds on a significantly larger timescale
than the dynamical time, one may use the orbit averaged Fokker Planck equation
to describe the cumulative effect of small perturbations emanating from
stochastic fluctuations (e.g. Spitzer, 1987; Binney & Tremaine, 2008). In flux
conservation form, the mass change inside energy level $E$ can then be written
as (13) for an isotropic system (i.e., with $g=g(E)$) may then be written as
(El-Zant, 2008)

$\frac{\partial M(<E)}{\partial E}=-F_{E}+g\frac{\partial q}{\partial E},$
(29)

which says that a decrease in total mass of (in our case dark matter)
particles with energy less than $E$ is due to the Fokker Planck flux out of
the volume enclosed by $E$ in addition to changes in that volume resulting
from the evolution in the structure of the system.

A major simplification leading to equation (13) involves the use of the first
order energy diffusion coefficient $D[\Delta E]=\langle\Delta E\rangle/T$, the
value of which is inferred directly from the simulation, instead of $F_{E}$.
Theoretically however, terms associated with the second order coefficient
$D[(\Delta E)^{2}]=\langle(\Delta E)^{2}\rangle/T$, constitute a comparable
contribution, as both coefficients involve the velocity variance
$\langle(\Delta v)^{2}\rangle$. Nevertheless, as we now show, if $D[\Delta E]$
can be considered an energy independent constant, as assumed in equation (10),
the total flux is simply proportional to $D[\Delta E]$ 777In the terminology
of Sections 4.1 and 4.2, equation (10) effectively assumes a global energy
input, with effective gas density $\rho_{0}$ corresponding to the average
density within the saturation radius $r_{\rm sat}$, and parameters
corresponding to those in equation 7. Furthermore, for power law systems, the
inclusion of $D[(\Delta E)^{2}]$ in the evaluation of the total flux leads to
a simple rescaling that can be absorbed in a parameter $a$.

In general, for a spherical system with isotropic velocities, the orbit
averaged Fokker Planck flux can be written as (e.g. El-Zant, 2008) 888Note
that $\langle\Delta E\rangle$ and $\langle(\Delta E)^{2}\rangle$ are defined
there directly as the rates of changes of energy and its variance, and
therefore are themselves the diffusion coefficients referred to here by
$D[\Delta E]$ and $D[(\Delta E)^{2}]$.

$F_{E}=\left(D_{E}-D_{EE}\frac{\partial}{\partial E}\right)g,$ (30)

where

$D_{E}=pD[\Delta E]-\frac{1}{2}\frac{\partial}{\partial E}\left(pD[(\Delta
E)^{2}]\right),$ (31)

and

$D_{EE}=\frac{1}{2}pD[(\Delta E)^{2}].$ (32)

As dark matter particles are too light to exchange energy through dynamical
friction with the fluctuating gas field, $\langle\Delta
E\rangle=\frac{1}{2}\langle(\Delta v)^{2}\rangle$. Also, to leading order in
perturbation $\Delta{\bf v}/{\bf v}$, $\langle(\Delta
E)^{2}\rangle=\langle({\bf v}.\Delta{\bf
v})^{2}\rangle=\frac{1}{3}v^{2}\langle(\Delta v)^{2}\rangle$ (e.g. Peñarrubia,
2019a). If we assume that $D[\Delta E$] is independent of energy so will
$\langle(\Delta v)^{2}\rangle$. In this case, the first order energy diffusion
coefficient is already orbit averaged, and averaging the second order one
simply involves evaluating the mean of $v^{2}$ over the energy shell of
surface area $p$, which we denote by $\langle v^{2}\rangle_{E}$. The Fokker
Planck flux across that shell is then given by

$F_{E}=D[\Delta E]\left[pg-\frac{1}{3}\frac{\partial}{\partial E}\left(\langle
v^{2}\rangle_{E}pg\right)\right].$ (33)

The total flux is thus simply proportional to the first order coefficient
$D[\Delta E]$.

Furthermore, for a power law region embedded in a larger halo, the density of
states may be approximated by $p=p_{0}(E-E_{0})^{1/2-3/\beta}$, and the
corresponding $\langle v^{2}\rangle_{E}=6\beta(3\beta-6)^{-1}(E-E_{0})$, and
thus

$F_{E}=p_{0}g_{0}D[\Delta
E]\left(\frac{6+\beta}{6-3\beta}\right)(E-E_{0})^{-1/\beta},$ (34)

From equations (27), and (28) this can be rewritten as

$F_{E}=\left(\frac{(6+\beta)(\beta-1)}{\beta(6-3\beta)}\right)\frac{M(<E)}{E-E_{0}}D[\Delta
E].$ (35)

Thus, for power law systems, the total flux is proportional to to the flux
associated with $D[\Delta E]$ through a constant scaling factor (in brackets),
that depends on the steepness of the power law. It is also always positive for
systems with density decreasing with radius ($\beta>-2$), allowing for mass
outflow from stochastic fluctuations.

###  C.3 Self consistent evolution and evaluation of parameter $a$

Figure 18: Dimensionless potential (defined as in Fig. 17) as a function of
radius and time, as the central halo region evolves under the influence of the
gas flcutuations.

To estimate the second term on the right hand side of equation (29) we note
that, from equation (26),

$\frac{\partial q}{\partial t}=-16\pi^{2}\int_{0}^{r_{\rm
max}}\frac{\partial\Phi}{\partial t}[2(E-\Phi)]^{1/2}r^{2}dE.$ (36)

Therefore the orbit averaged change of potential due to the self consistent
evolution at energy $E$ is then

$\langle\dot{\Phi}(r)\rangle_{E}=-\frac{1}{p}\frac{\partial q}{\partial t}.$
(37)

Recall that in Fig. 17 we found that the change of potential at the centre of
the halo ($r=0$), resulting from the gas fluctuations, corresponds to the
average energy gained per particle, taking into account the self consistent
evolution of the system. Fig. 16, on the other hand, showed that that this
corresponds to the average energy per unit mass pumped into the central halo
multiplied by a factor $\frac{2-R_{v}}{1-R_{v}}$, with $R_{v}=0.4$; i.e., by a
factor factor of about $2.7$. Thus $\dot{\Phi}_{E=E_{0}}=2.7\times D[\Delta
E]$. But, as can be seen Fig. 18, the change in the potential with time is
significantly smaller at larger radii. Thus $\dot{\Phi}_{E>E_{0}}<2.7\times
D[\Delta E]$; at $r=1{\rm kpc}$ (i.e., about the radius inside which the core
forms) the total change is a bit more than a quarter of the value at $r=0$. As
the radial variation is not very steep we will approximate the change in
energy due to evolution of the self consistent potential with an effective
value, such that $\dot{\Phi}_{\rm eff}=\eta D[\Delta E]$, where from the above
considerations $0.7<\eta<2.7$. As we will see below, the mass particles within
energy shell $E$ increases steeply with energy away from the zero point in the
central power law regions, with $M(<E)\sim(E-E_{0})^{2.3}$. An effective
value, replacing $\langle\dot{\Phi}(r)\rangle_{E}$, would thus correspond to
more weight over larger energies, and therefore larger radii. We may thus
expect $\eta$ to be closer to its lower, rather than upper, limit.

Valid values of $\eta$ are actually tightly constrained by the mass and
density evolution in our model galaxy (Fig. 9), and the differential energy
distribution in the inner halo. To see this we replace $\partial q/\partial t$
by $-p\eta D[\Delta E]$, and use equation (35), to rewrite equation (29) as

$\frac{\partial M(<E)}{\partial
t}=-\left(\frac{(6+\beta)(\beta-1)}{\beta(6-3\beta)}+\eta\right)\frac{M(<E)}{E-E_{0}}D[\Delta
E].$ (38)

This is equivalent to equation (13) if $a$ is equal to the term in brackets.
That is,

$a=\frac{(6+\beta)(\beta-1)}{\beta(6-3\beta)}+\eta.$ (39)

To determine $\beta$ we need to show that the differential energy distribution
of the inner halo may be fit by a power law and determine its index, in which
case $\beta$ is fixed by equation (27). Fig. 19 shows that in the region where
the core replaces the cusp the differential energy distribution of an NFW
profile may indeed be fit by a power law, with index of about $1.35$, which
corresponds to $\beta=-0.74$. Furthermore, as can also be seen from Fig. 19
the differential energy slopes of models with the same asymptotic physical
density distribution for $r\gg r_{s}$, but have milder inner cusps
($\gamma=-1/2$), or even cores ($\gamma=0$), can be approximated by similar
power laws (with $\beta\simeq-0.8$), suggesting that the power law
approximation may hold as the cusp transforms into core.

Furthermore, in Section 4.3 we found that a value of $2.5\lesssim a\lesssim
2.8$ can be used to obtain solutions to equation (13) that describe the mass
and density evolution fairly well (Fig. 9, right hand panel). For the
aforementioned values of $\beta$, this corresponds to $1\lesssim\eta\lesssim
1.4$, which is within the expected range inferred above.

###  C.4 The vanishing of the flux for exponential differential energies (and
constant dispersion)

Figure 19: Differential energy distributions in the inner regions of three
density profiles with $\rho=\rho_{s}r^{\gamma}(1+r/r_{s})^{-3-\gamma}$, where
$\gamma=-1$ for the NFW profile. The dimensionless energy $\mathcal{E}=E/4\pi
G\rho_{s}r_{s}^{2}$ is scaled such that the three profiles converge to the
same densities as $r\gg r_{s}$ (thus $C=1$ for the NFW, $C=3/2$ for
$\gamma=-1/2$ and $C=2$ for $\gamma=0$; cf. Widrow, 2000 for more details
regarding these models). The dashed dotted lines are power law fits. The one
fitting the NFW has slope $1.35$. It approximates the differential energy
distribution over energy scales $0.03\lesssim\mathcal{E}+1\lesssim 0.33$,
which enclose zero velocity spheres with radii of about $50{\rm pc}$ and
$1{\rm kpc}$ in our simulation (with $r_{s}=0.88{\rm kpc}$ as in Table 1),
which correspond, respectively, to about the smallest radius inside which the
halo density may be reliably resolved and the region inside which the gas
fluctuations have significant effect on the halo profile. The dashed dotted
line fitting the shallower profiles ($\gamma=-1,2,0$) has only slightly
flatter slope of $1.25$. Figure 20: The differential energy distributions of
the models of Fig. 19 in the intermediate energy range, where the quasi-linear
variation on the log-linear plot suggests an exponential scaling, particularly
for $0.35\lesssim|\mathcal{E}|\lesssim 0.65$, corresponding to zero velocity
spheres with radii of about 1.4 and 4.7 kpc in our model galaxy.

At larger energies, corresponding to intermediate radii beyond $r_{s}$, the
differential energy distribution is better approximated by an exponential,
rather than a power law, as Fig. 20 illustrates. At such radii, the halo is
nearly isothermal, with density $\rho_{p}\sim 1/r^{2}$, and $\langle
v^{2}\rangle_{E}$ approximately constant. Assuming this, and inserting
$\partial M(<E)/\partial E\propto e^{\beta E}$ into (33), it is easily seen
that the flux is nullified if $\beta=1/\sigma^{2}$, where
$\sigma^{2}=\frac{1}{3}\langle v^{2}\rangle$ is the one dimensional velocity
variance. This may explain the formation of cores up to the region where the
differential energy distribution may be approximated by a power law, where the
flux is positive, but that the response to fluctuations is more robust when
the exponential is a better approximation.