A model for core formation in dark matter haloes and ultra diffuse galaxies by outflow episodes

TL;DR

This paper introduces a simple analytical model describing how feedback-driven outflows can transform dark matter haloes and ultra-diffuse galaxies into core-like profiles, validated against cosmological simulations.

Contribution

It generalizes previous models to include continuous density and potential profiles, providing a practical tool for understanding core formation due to outflows.

Findings

Successfully predicts dark matter profile evolution in 75% of simulated cases

Models core formation in dark matter haloes and UDGs driven by supernova feedback

Fails mainly during merger events in simulations

Abstract

We present a simple model for the response of a dissipationless spherical system to an instantaneous mass change at its center, describing the formation of flat cores in dark matter haloes and ultra-diffuse galaxies (UDGs) from feedback-driven outflow episodes in a specific mass range. This model generalizes an earlier simplified analysis of an isolated shell into a system with continuous density, velocity and potential profiles. The response is divided into an instantaneous change of potential at constant velocities due to a given mass loss or gain, followed by energy-conserving relaxation to a new Jeans equilibrium. The halo profile is modeled by a two-parameter function with a variable inner slope and an analytic potential profile (Dekel et al. 2017), which enables determining the associated kinetic energy at equilibrium. The model is tested against NIHAO cosmological zoom-in simulations, where it successfully predicts the evolution of the inner dark-matter profile between successive snapshots in about 75% of the cases, failing mainly in merger situations. This model provides a simple understanding of the formation of dark-matter halo cores and UDGs by supernova-driven outflows, and a useful analytic tool for studying such processes.

Tables

Table 1. Table 1: Fiducial parameters used when implementing the model and assessing its success.

Parameter	Description	Value
$t_{\min }$	earliest time considered	3 Gyr
$b$	fixed profile parameter	2
$g$	fixed profile parameter	3
$r_0/R_{\mathrm{vir}}$	radius for evaluating $s_0$ and ${\overline{\rho }}_0$	0.01
$\Delta \log {(r/R_{\mathrm{vir}})}_{\mathrm{fit}}$	radius range for the profile fits	$[-2,0]$
$\Delta \log {(r/R_{\mathrm{vir}})}_{\mathrm{model}}$	radius range for the model	$[-1.75,0]$
$\Delta \log {(r/R_{\mathrm{vir}})}_{\mathrm{success}}$	radius range to define $\delta$	[-2,-1.5]
$\log {(\overline{\rho }/{\overline{\rho }}_0)}_{\min }$	floor to define $\delta$	-1
$r_{\mathrm{merger}}/R_{\mathrm{vir}}$	radius for $f_{\mathrm{merger}}$	$0.15$
$f_{\mathrm{merger},\min }$	threshold for $f_{\mathrm{merger}}$	$10\%$
$f_{\min }$	threshold for ${\left|f\right|}_{\mathrm{RMS}}$	7%
${\delta }_{\min }$	minimum distinguishable ${\delta }_{\mathrm{sim}}$	3%
${\delta }_{\max }$	maximum $\delta$ for a success	10%
$\Delta s_{\max }$	maximum $\left|\Delta s\right|$ for a success	0.10

Table 2. Table 2: Success rate of the instant mass change model for the 33 NIHAO simulated galaxies whose stellar mass lies between 5 × 10 7 5 superscript 10 7 5\times 10^{7} and 5 × 10 9 ​ M ⊙ 5 superscript 10 9 subscript 𝑀 direct-product 5\times 10^{9}~{}M_{\odot} at z = 0 𝑧 0 z=0 , excluding the earliest phase of galaxy evolution before 3 Gyr: p all subscript 𝑝 all p_{\rm all} for all such outputs; p nomergers subscript 𝑝 nomergers p_{\rm nomergers} excluding mergers as indicated in section  4.2 ; p > f min subscript 𝑝 absent subscript 𝑓 min p_{>f_{\rm min}} and p < f min subscript 𝑝 absent subscript 𝑓 min p_{<f_{\rm min}} further focusing either on cases where we expect a significant profile change ( | f | RMS > f min subscript 𝑓 RMS subscript 𝑓 min |f|_{\rm RMS}>f_{\rm min} ) or on cases where we expect negligible changes ( | f | RMS < f min subscript 𝑓 RMS subscript 𝑓 min |f|_{\rm RMS}<f_{\rm min} ). The simulation ID follows the nomenclature used by Wang et al. ( 2015 ) but the indicated virial and stellar masses M vir subscript 𝑀 vir M_{\rm vir} and M ⋆ subscript 𝑀 ⋆ M_{\star} as well as the virial radius R vir subscript 𝑅 vir R_{\rm vir} at z = 0 𝑧 0 z=0 slightly differ from their values since we determined the virial radius R vir subscript 𝑅 vir R_{\rm vir} following Bryan & Norman ( 1998 ) instead of using R 200 subscript 𝑅 200 R_{200} and calculated M ⋆ subscript 𝑀 ⋆ M_{\star} within 0.15 × R vir 0.15 subscript 𝑅 vir 0.15\times R_{\rm vir} instead of 0.2 × R vir 0.2 subscript 𝑅 vir 0.2\times R_{\rm vir} . The stellar half-light radius r e subscript 𝑟 𝑒 r_{e} and the V-band central surface brightness μ 0 , V subscript 𝜇 0 𝑉 \mu_{0,V} at z = 0 𝑧 0 z=0 are also indicated. Daggers ( † † \dagger ) highlight UDGs. Galaxies are ordered by halo mass; fiducial galaxies g6.12e10 and g1.08e11 are highlighted in bold.

ID	$M_{\mathrm{vir}}$	$R_{\mathrm{vir}}$	$M_{\star }$	$r_e$	${\mu }_{0,V}$	$p_{\mathrm{all}}$	$p_{\mathrm{no}\mathrm{mergers}}$	$p_{>f_{\min }}$
	$\left[M_{\odot }\right]$	$\left[\mathrm{kpc}\right]$	$\left[M_{\odot }\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{mag}/{\mathrm{as}}^2\right]$
g2.94e10†	$3.80\times {10}^{10}$	$88.8$	$5.86\times {10}^7$	$1.96$	$24.26$	22/27 = 81 %	22/27 = 81 %	06/09 = 67 %	16/18 = 89 %
g5.05e10†	$4.74\times {10}^{10}$	$95.6$	$9.47\times {10}^7$	$2.09$	$24.28$	27/38 = 71 %	26/33 = 79 %	09/12 = 75 %	17/21 = 81 %
g6.12e10†	$5.50\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}$	$100.5$	$9.13\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟕}}$	$2.67$	$24.03$	27/33 = 82 %	27/30 = 90 %	10/11 = 91 %	17/19 = 89 %
g4.27e10†	$5.54\times {10}^{10}$	$100.7$	$6.15\times {10}^7$	$2.83$	$24.21$	15/20 = 75 %	15/19 = 79 %	09/11 = 82 %	06/08 = 75 %
g4.99e10†	$5.71\times {10}^{10}$	$101.8$	$1.22\times {10}^8$	$3.01$	$24.16$	31/43 = 72 %	31/42 = 74 %	13/23 = 57 %	18/19 = 95 %
g4.86e10†	$5.79\times {10}^{10}$	$102.2$	$1.22\times {10}^8$	$2.36$	$24.37$	45/50 = 90 %	45/50 = 90 %	00/02 = 0 %	45/48 = 94 %
g4.94e10	$5.93\times {10}^{10}$	$103.0$	$1.11\times {10}^8$	$2.42$	$23.89$	25/38 = 66 %	25/33 = 76 %	04/10 = 40 %	21/23 = 91 %
g3.44e10†	$6.65\times {10}^{10}$	$107.0$	$6.32\times {10}^7$	$2.54$	$24.43$	21/33 = 64 %	19/26 = 73 %	13/15 = 87 %	06/11 = 55 %
g4.48e10	$6.77\times {10}^{10}$	$107.7$	$1.37\times {10}^8$	$3.62$	$23.92$	16/31 = 52 %	16/24 = 67 %	06/11 = 55 %	10/13 = 77 %
g6.91e10	$7.82\times {10}^{10}$	$113.0$	$2.50\times {10}^8$	$2.54$	$23.49$	39/50 = 78 %	39/49 = 80 %	07/12 = 58 %	32/37 = 86 %
g6.96e10	$1.08\times {10}^{11}$	$125.9$	$3.64\times {10}^8$	$3.49$	$23.34$	34/47 = 72 %	34/42 = 81 %	14/17 = 82 %	20/25 = 80 %
g6.77e10	$1.09\times {10}^{11}$	$126.1$	$4.83\times {10}^8$	$4.37$	$23.22$	24/42 = 57 %	21/35 = 60 %	19/27 = 70 %	02/08 = 25 %
g8.89e10	$1.10\times {10}^{11}$	$126.6$	$4.02\times {10}^8$	$3.12$	$23.10$	42/50 = 84 %	42/48 = 88 %	07/10 = 70 %	35/38 = 92 %
g9.59e10†	$1.13\times {10}^{11}$	$127.9$	$2.75\times {10}^8$	$4.91$	$24.48$	15/50 = 30 %	13/39 = 33 %	13/33 = 39 %	00/06 = 0 %
g3.23e11†	$1.13\times {10}^{11}$	$127.7$	$3.60\times {10}^8$	$4.85$	$24.05$	25/49 = 51 %	24/43 = 56 %	22/36 = 61 %	02/07 = 29 %
g1.05e11	$1.32\times {10}^{11}$	$134.6$	$5.66\times {10}^8$	$4.99$	$23.39$	25/44 = 57 %	24/38 = 63 %	15/26 = 58 %	09/12 = 75 %
g6.37e10†	$1.35\times {10}^{11}$	$135.7$	$2.11\times {10}^8$	$4.19$	$24.10$	17/41 = 41 %	16/30 = 53 %	10/16 = 62 %	06/14 = 43 %
g1.08e11†	$1.36\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟏}}$	$135.9$	$8.47\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟖}}$	$4.41$	$24.10$	41/50 = 82 %	41/47 = 87 %	10/12 = 83 %	31/35 = 89 %
g2.19e11	$1.52\times {10}^{11}$	$140.9$	$9.27\times {10}^8$	$4.52$	$23.01$	27/47 = 57 %	27/41 = 66 %	20/31 = 65 %	07/10 = 70 %
g1.37e11	$1.69\times {10}^{11}$	$146.2$	$2.02\times {10}^9$	$3.46$	$22.51$	44/50 = 88 %	44/50 = 88 %	03/05 = 60 %	41/45 = 91 %
g1.57e11	$1.78\times {10}^{11}$	$148.6$	$1.15\times {10}^9$	$5.42$	$23.62$	29/50 = 58 %	24/36 = 67 %	18/23 = 78 %	06/13 = 46 %
g1.52e11†	$1.81\times {10}^{11}$	$149.5$	$7.90\times {10}^8$	$5.81$	$24.29$	32/50 = 64 %	32/50 = 64 %	02/07 = 29 %	30/43 = 70 %
g1.59e11†	$2.00\times {10}^{11}$	$154.6$	$6.69\times {10}^8$	$6.21$	$24.97$	39/50 = 78 %	39/48 = 81 %	16/19 = 84 %	23/29 = 79 %
g2.04e11	$2.43\times {10}^{11}$	$164.9$	$4.70\times {10}^9$	$3.07$	$20.18$	43/50 = 86 %	43/46 = 93 %	18/21 = 86 %	25/25 = 100 %
g1.64e11	$2.45\times {10}^{11}$	$165.4$	$9.12\times {10}^8$	$5.72$	$22.41$	10/33 = 30 %	08/21 = 38 %	07/16 = 44 %	01/05 = 20 %
g2.41e11	$2.85\times {10}^{11}$	$173.8$	$4.10\times {10}^9$	$3.97$	$21.48$	36/49 = 73 %	34/43 = 79 %	12/15 = 80 %	22/28 = 79 %
g2.54e11	$3.00\times {10}^{11}$	$177.0$	$3.50\times {10}^9$	$3.60$	$19.84$	28/46 = 61 %	27/38 = 71 %	17/25 = 68 %	10/13 = 77 %
g4.90e11	$3.85\times {10}^{11}$	$192.3$	$3.43\times {10}^9$	$5.77$	$22.95$	37/50 = 74 %	34/43 = 79 %	07/08 = 88 %	27/35 = 77 %
g5.46e11	$3.87\times {10}^{11}$	$192.5$	$3.77\times {10}^9$	$5.29$	$22.47$	37/48 = 77 %	36/43 = 84 %	07/09 = 78 %	29/34 = 85 %
g3.21e11	$3.95\times {10}^{11}$	$193.9$	$3.67\times {10}^9$	$4.88$	$21.29$	25/48 = 52 %	25/39 = 64 %	18/29 = 62 %	07/10 = 70 %
g3.59e11	$4.44\times {10}^{11}$	$201.6$	$4.36\times {10}^9$	$5.10$	$21.76$	26/46 = 57 %	24/33 = 73 %	14/19 = 74 %	10/14 = 71 %
g3.49e11	$5.04\times {10}^{11}$	$210.4$	$3.96\times {10}^9$	$5.23$	$22.02$	39/50 = 78 %	39/49 = 80 %	07/13 = 54 %	32/36 = 89 %
g3.55e11	$5.04\times {10}^{11}$	$210.3$	$3.85\times {10}^9$	$6.48$	$22.44$	19/43 = 44 %	16/33 = 48 %	11/24 = 46 %	05/09 = 56 %
All						962/1446 = 67 %	932/1268 = 74 %	364/557 = 65 %	568/711 = 80 %

Equations

$$\begin{array}[]{@{}l}\displaystyle E_{i}=U_{i}+K_{i}\\ \displaystyle\text{with }~{}U_{i}=-\frac{GM}{r_{i}}~{}\text{ and }~{}K_{i}=\frac{1}{2}\frac{GM}{r_{i}},\end{array}$$

$$\begin{array}[]{@{}l}\displaystyle E_{t}=U_{t}+K_{t}\\ \displaystyle\text{with }~{}U_{t}=U_{i}-\frac{Gm}{r_{i}}~{}\text{ and }~{}K_{t}=K_{i}\end{array}$$

$$\begin{array}[]{@{}l}\displaystyle E_{f}=U_{f}+K_{f}\\ \displaystyle\text{with }~{}U_{f}=-\frac{G(M+m)}{r_{f}}~{}\text{ and }~{}K_{f}=\frac{G(M+m)}{2r_{f}}.\end{array}$$

$$U(r)={\int }_r^{+\infty }F(y)dy=-{\int }_r^{R_t}\frac{GM(y)}{y^2}dy-\frac{G{M}_t}{R_t},$$

$$U(r)=-\frac{GM(r)}{r}-G{\int }_r^{R_t}\frac{1}{y}\frac{dM}{dy}dy=-\frac{GM(r)}{r}-{\int }_r^{R_t}4\pi yG\rho (y)dy;$$

$$K_{\mathrm{isolated}}(r)=\frac{1}{2}\frac{GM}{r}$$

$$K(r)=\frac{dK(<r)}{dM}=\frac{dK(<r)}{dr}\frac{dR}{dr}$$

$$E_i(r_i)=U_i(r_i)+K_i(r_i),$$

$$\begin{array}[]{@{}l}\displaystyle E_{t}(r_{i})=U_{t}(r_{i})+K_{i}(r_{i})\\ \text{with }~{}\displaystyle U_{t}(r_{i})=U_{i}(r_{i})-\frac{Gm}{r_{i}},\end{array}$$

$$\begin{array}[]{@{}l}\displaystyle E_{f}(r_{f})=U_{f}(r_{f})+K_{f}(r_{f},m)\\ \text{with }~{}\displaystyle U_{f}(r_{f})=U_{f}(r_{f})-\frac{Gm}{r_{f}}.\end{array}$$

$$U_i(r_i)-\frac{Gm}{r_i}+K_i(r_i)=U_f(r_f)-\frac{Gm}{r_f}+K_f(r_f,m).$$

$$\overline{\rho }(r)=\frac{\overline{{\rho }_c}}{x^a(1+x^{1/b}{)}^{b(g-a)}}$$

$$\rho (r)=\frac{3-a}{3}\frac{c^3\mu \overline{{\rho }_{\mathrm{vir}}}}{x^a(1+x^{1/2}{)}^{2(3.5-a)}}$$

$$\displaystyle U_{\rm DM}(r)=-V_{\rm vir}^{2}~{}\Bigg{(}1+2c\mu\left[\frac{\chi_{\rm vir}^{2(2-a)}-\chi^{2(2-a)}}{2(2-a)}-\frac{\chi_{\rm vir}^{2(2-a)+1}-\chi^{2(2-a)+1}}{2(2-a)+1}\right]\Bigg{)}$$

$$s(r)=-\frac{d\log \rho }{d\log r}=\frac{a+3.5x^{1/2}}{1+x^{1/2}}$$

$$s_0=s(0.01R_{\mathrm{vir}})$$

$$c_{\max }=\frac{c}{(2-a{)}^2}$$

$$K(r)=\frac{3-2\beta }{2}{\sigma }_r^2,$$

$${\sigma }_r^2(r)=\frac{G}{\rho (r)}{\int }_r^{\infty }\rho (r^{\prime })M(r^{\prime })r^{-2}\mathrm{d}{r}^{\prime }$$

$$K_{\rm DM}(r)=c\mu\frac{GM_{{\rm vir}}}{R_{{\rm vir}}}\frac{(3-a)\overline{\rho_{c}}}{\rho(r)}\Big{[}\mathcal{B}(4-4a,9,\zeta)\Big{]}_{\chi}^{1}$$

$$\frac{M_{\mathrm{tot}}}{M}=X_M{\left(\frac{r}{R_{\mathrm{vir}}}\right)}^{-n}$$

$$K_{\rm multi}(r)=\mu X_{M}c^{n+1}\frac{GM_{{\rm vir}}}{R_{{\rm vir}}}\frac{(3-a)\overline{\rho_{c}}}{\rho(r)}\Big{[}\mathcal{B}(4-4a-2n,9+2n,\zeta)\Big{]}_{\chi}^{1}$$

$${\sigma }_m^2(r)=\frac{G}{\rho (r)}{\int }_r^{\infty }\rho (r^{\prime })m(r^{\prime }){r^{\prime }}^{-2}\mathrm{d}{r}^{\prime }.$$

$${\int }_r^{\infty }\rho (r^{\prime })m(r^{\prime }){r^{\prime }}^{-2}\mathrm{d}{r}^{\prime }=m(r){\int }_r^{\infty }\rho (r^{\prime }){r^{\prime }}^{-2}\mathrm{d}{r}^{\prime },$$

$$K_{m}(r)=1.5\sigma_{m}^{2}(r)=\frac{Gmc}{R_{\rm vir}}\frac{(3-a)\overline{\rho_{c}}}{{\rho(r)}}\Big{[}\mathcal{B}(-2-2a,9,\zeta)\Big{]}_{\chi}^{1}$$

$$f_{\rm merger}=\frac{m}{M_{\rm tot}}\Bigg{|}_{~{}r~{}=~{}0.15R_{\rm vir}},$$

$$\delta =2\left|\frac{{\mathcal{A}}_{\mathrm{model}}-{\mathcal{A}}_{k+1}}{{\mathcal{A}}_{\mathrm{model}}+{\mathcal{A}}_{k+1}}\right|$$

$$(\delta \le {\delta }_{\max })\text{AND}(\delta \le {\delta }_0\text{OR}{\delta }_{\mathrm{sim}}\le {\delta }_{\min }),$$

$$\Delta s=s_{0,\mathrm{model}}-s_{0,\mathrm{f}}$$

$$\Delta c=\log (c_{\max ,\mathrm{model}}/c_{\max ,\mathrm{f}}).$$

Full text

††pubyear: 2019††pagerange: A model for core formation in dark matter haloes and ultra diffuse galaxies by outflow episodes–C

1 Introduction

1.1 The cusp-core discrepancy

The Lambda cold dark matter (LCDM) model of structure formation is extremely successful at describing the evolution of the Universe, from the nearly homogeneous state shown by the cosmic microwave background to the present day clustered distribution of matter with galaxies, clusters of galaxies and voids between them. But while it agrees with observations on the large-scale structure of the Universe, it has faced different challenges at galactic scales. In particular, while LCDM numerical simulations predict steep, ‘cuspy’ density profiles for dark matter halos, observations of dark matter dominated dwarf, low-surface-brightness and dwarf satellite galaxies as well as clusters favor shallower ‘cores’ (e.g., Flores & Primack, 1994; Moore, 1994; McGaugh & de Blok, 1998; van den Bosch & Swaters, 2001; de Blok et al., 2008; de Blok, 2010; Kuzio de Naray & Spekkens, 2011; Oh et al., 2011, 2015; Newman et al., 2013a; Newman et al., 2013b; Adams et al., 2014). Proposed solutions to this ‘cusp-core’ discrepancy and to potentially related challenges of LCDM cosmology such as the ‘too big to fail’ problem111LCDM simulations predict a population of subhaloes around Milky Way-sized galaxies that have higher central surface densities than inferred from observed satellites of the Milky Way and its neighbours (e.g., Boylan-Kolchin et al., 2011) can be broadly divided into solutions considering fundamental changes in the physics of the model and solutions focusing on the baryonic processes at stake during galaxy formation and evolution. The first category comprises alternatives to LCDM such as warm dark matter, self-interacting dark matter and models that fundamentally change the gravitational law (e.g., Milgrom, 1983; Burkert, 2000; Colín et al., 2000; Goodman, 2000; Hu et al., 2000; Spergel & Steinhardt, 2000; Bode et al., 2001; Gentile et al., 2011; Famaey & McGaugh, 2012; Macciò et al., 2012a; Destri et al., 2013; Peter et al., 2013; Marsh & Silk, 2014). Solutions invoking baryonic processes within the LCDM framework are motivated by the fact that the discrepancies between model and observations precisely occur at the scale at which baryons start to play an important role, notably through powerful stellar feedback processes that generate significant movements of the gas. Moreover, hydrodynamical simulations with different feedback implementations are able to reproduce dark matter cores (e.g., Governato et al., 2010, 2012; Macciò et al., 2012b; Martizzi et al., 2013; Teyssier et al., 2013; Di Cintio et al., 2014a; Chan et al., 2015; Peirani et al., 2017).

Baryons mostly affect dark matter halos through their gravitational potential. When baryons cool and contract, they accumulate at the center of the halo, which steepens the potential and causes the dark matter to contract as well through adiabatic contraction (Blumenthal et al., 1986). When a massive object such as a satellite galaxy or a clump of gas moves with respect to the dark matter background, it can transfer part of its kinetic energy to the background through dynamical friction, as the concentration of particles increases in its wake and generates a drag force (Chandrasekhar, 1943). Dynamical friction can thus ‘heat’ the dark matter halo and contribute to remove the central cusp (El-Zant et al., 2001; El-Zant et al., 2004; Tonini et al., 2006; Romano-Díaz et al., 2008; Goerdt et al., 2010; Cole et al., 2011; Nipoti & Binney, 2015). Alternatively, repeated gravitational potential fluctuations induced by stellar winds, supernova explosions and active galactic nuclei (AGNs) can also dynamically heat the dark matter and lead to the formation of a core (Dekel & Silk, 1986; Read & Gilmore, 2005; Mashchenko et al., 2006, 2008; Pontzen & Governato, 2012, 2014; Governato et al., 2012; Zolotov et al., 2012; Martizzi et al., 2013; Teyssier et al., 2013; Madau et al., 2014; Dutton et al., 2016b; El-Zant et al., 2016; Peirani et al., 2017). In this case, variations in the baryonic mass distribution induce violent potential fluctuations which progressively disperse dark matter particles away from the center of the halo.

1.2 Ultra Diffuse Galaxies

Deep imaging observations reveal the existence of a population of low central surface brightness ($\rm\mu_{g,0}>24~{}mag.arcsec^{-2}$) galaxies with stellar masses of dwarfs and effective radii comparable to that of the Milky Way ($\rm r_{eff}>1.5\rm~{}kpc$). These Ultra Diffuse Galaxies (UDGs) appear to be ubiquitous in groups and clusters (e.g., van Dokkum et al., 2015a; van Dokkum et al., 2015b; Janowiecki et al., 2015; Koda et al., 2015; Muñoz et al., 2015; Mihos et al., 2015, 2017; Merritt et al., 2016; van der Burg et al., 2016; Yagi et al., 2016; Janssens et al., 2017; Shi et al., 2017) but are also observed in the field (e.g., Román & Trujillo, 2017; Leisman et al., 2017). Possible formation scenarii include them being (i) failed Milky Way-like galaxies that lost their gas after forming their first stars (van Dokkum et al., 2015a; van Dokkum et al., 2015b, 2016; Yozin & Bekki, 2015; Peng & Lim, 2016; Lim et al., 2018), possibly due to the denser environment in groups and clusters, (ii) the high-spin tail of the dwarf galaxy population (Amorisco & Loeb, 2016; Rong et al., 2017; Shi et al., 2017), (iii) tidal debris from mergers or tidally disrupted dwarfs (Beasley & Trujillo, 2016; Greco et al., 2018; Jiang et al., 2019; Carleton et al., 2019) or (iv) galaxies whose spatial extend is due to episodes of inflows and outflows from stellar feedback (Di Cintio et al., 2017; Chan et al., 2018; Jiang et al., 2019).

From deep images and stacks, van Dokkum et al. (2015b) and Mowla et al. (2017) find no evidence for the majority of Coma clusters UDGs to be tidally disrupted. Although the abundance of globular clusters in UDGs compared to dwarf galaxies of similar masses displays a significant scatter, its often large value may indicate intense star formation episodes in the past (Beasley et al., 2016; Peng & Lim, 2016; van Dokkum et al., 2016, 2017; Lim et al., 2018; Amorisco et al., 2018). Some UDGs also seem to harbour signs of recent or ongoing star formation (Martínez-Delgado et al., 2016; Trujillo et al., 2017; Shi et al., 2017; Pandya et al., 2018) while both Janowiecki et al. (2015) and Shi et al. (2017) report group UDGs rich in atomic HI gas. Martín-Navarro et al. (2019) further interpret the high [Mg/Fe] abundance ratio of a field UDG with an extended star formation history as a possible result of supernova-driven gas outflows. These observations would fit with UDGs being formed by outflows resulting from stellar feedback. Indeed, the feedback-induced gravitational potential fluctuations leading to dark matter core formation can also cause the expansion of the stellar component (Teyssier et al., 2013; Dutton et al., 2016b; El-Badry et al., 2016) and this process is precisely expected to be most efficient in the mass range of UDGs (Di Cintio et al., 2014a, b; Tollet et al., 2016; Dutton et al., 2016b). In the NIHAO simulations (Wang et al., 2015) studied by Di Cintio et al. (2017), prolonged and persistent bursty star formation histories lead to the simultaneous expansion of both the dark matter and stellar distributions of UDGs.

1.3 Motivation for this work

Although hydrodynamical simulations are able to reproduce the observed cores of low-mass dark matter haloes and a UDG population by adding prescriptions for baryonic processes such as stellar winds, stellar radiation, and supernova feedback, they do not describe nor quantify the exact mechanisms through which baryons affect the dark matter and stellar collisionless spatial distributions. Dekel & Silk (1986) show that a significant supernova-driven gas loss can explain the low surface brightness of dwarf ellipticals, Pontzen & Governato (2012, 2014) that a fluctuating gravitational potential can irreversibly transfer energy to collisionless particles. Attempting to isolate further the physical mechanism at stake during dark matter core formation, El-Zant et al. (2016) propose a theoretical model comparable to two-body relaxation in stellar systems. In this model, core formation occurs from dark matter particles experiencing successive ‘kicks’ resulting from feedback-induced stochastic potential fluctuations, which cumulatively lead the particles to deviate from their trajectories as in a diffusion process. In the context of a fuzzy dark matter halo, Bar-Or et al. (2018) similarly describe the effect on stars and black holes of stochastic density fluctuations due to the large de Broglie wavelength of the dark matter particles they assume while Marsh & Niemeyer (2018) and El-Zant et al. (2019) use the El-Zant et al. (2016) model to place constraints on the fuzzy dark matter particle mass. In yet another context, Penoyre & Haiman (2018) study how particles on Keplerian orbits respond to a sudden drop in the central mass, notably following how the density of the system varies as they move to higher-energy orbits.

In this article, we present and test a simple theoretical model in which core formation results from cycles of inflows and outflows. In particular, we test how an instantaneous inflow or outflow episode followed by relaxation to a new equilibrium can affect the density profile of a spherical dissipationless mass distribution. The model we develop builds upon that presented by Dutton et al. (2016b) in their section 4.3 and recalled in Section 2, which focuses on an isolated spherical shell of collisionless matter, neglects the full halo mass distribution and approximates the shell kinetic energy by its expression derived from the equation of virial equilibrium. As shown in Section 2, these assumptions are not sufficient to describe the whole structure of a collisionless halo and its evolution. We provide here a consistent analytical framework to understand and predict the evolution of a collisionless spherical matter distribution when mass is added or removed at different radii and we test this model with the NIHAO simulations. We show that it accurately predicts the evolution of the dark matter density profile between two successive simulation snapshots in a large majority of cases, failing mostly in violent perturbed situations of major mergers.

This paper unfolds progressively as follows. In section 2, we recall the isolated shell model proposed by Dutton et al. (2016b) and highlight its main limitations. In section 3, we present our new spherical collisionless halo response model for an instant mass change, which enables retrieving the evolution of the halo mass distribution using accurate analytic parametrizations of the potential and kinetic energy stemming from the Dekel et al. (2017) density profile. In section 4, we successfully test this model on NIHAO (Wang et al., 2015) cosmological zoom-in simulations. In section 5, we discuss the reasons for the model’s failures and its ability to predict the evolution of the halo mass distribution over multiple episodes of inflows and outflows. In section 6, we summarize our results and conclude.

2 Isolated shell model and its limits

2.1 The isolated shell model of Dutton et al. (2016b)

In their Section 4.3, Dutton et al. (2016b) consider the evolution of a spherical halo shell enclosing a dark matter mass $M$ initially at radius $r_{i}$ when a baryonic mass $m$ is instantaneously added (or removed, in which case $m<0$) at the shell center. They assume that the halo is spherical, that orbits are circular, that shells do not cross and that energy is conserved between shells enclosing a given collisionless mass before and after the mass inflow (outflow). However, they ignore variations in the mass outside the shell and assume that the expression of the kinetic energy stemming from the virial theorem holds for an isolated shell.

Since the mass variation occurs on a sufficiently short time scale compared to the dynamical time of the system, Dutton et al. (2016b) assume an intermediate transitional state in which the gravitational potential adapts instaneously while the velocities remain frozen at their initial values, following the impulse approximation. They write the initial energy associated with the shell enclosing an isolated dark matter mass $M$ at radius $r_{i}$

[TABLE]

which becomes

[TABLE]

right after the mass variation in this intermediate transitional state. The potential energy per unit mass $U$ corresponds to an isolated spherical system of mass $M$, while the expression for the kinetic energy $K$ is assumed to be that stemming from the virial theorem, even if it is applied to a single shell and not to the whole halo. The system described by Eq. (2) right after the mass variation is out of equilibrium but the shell eventually settles in a new equilibrium state at radius $r_{f}$, where the kinetic energy is again given by the virial theorem. The final energy is then written

[TABLE]

Assuming energy conservation for the shell after the mass variation, i.e., $E_{f}=E_{t}$, they derive an expression for the ratio $r_{f}/r_{i}$ as a function of the ratio $f=m/M$ (their Eq. (16)).

Dutton et al. (2016b) further consider the case in which a slow inflow of mass $m$ is followed by an instantaneous outflow of mass $\eta m$ as well as a succession of multiple such episodes. Expressing the conservation of gas mass by the bathtub model (Bouché et al., 2010; Lilly et al., 2013; Dekel & Mandelker, 2014) and assuming that all baryons accreted at the virial radius penetrate through the halo into the central galaxy, they relate $\eta$ to the integrated star formation efficiency $\epsilon_{\rm SF}=M_{\star}/f_{\rm b}M_{\rm vir}$ with $\eta=1-\epsilon_{\rm SF}$, where $M_{\star}$ is the stellar mass, $f_{\rm b}$ the average baryon fraction in the Universe and $M_{\rm vir}$ the virial mass. Adjusting the number of cycles $N$ to match the simulations given $\epsilon_{\rm SF}$, they obtain simulated halo responses comparable to those expected from the model (cf. their Fig. 14). This agreement supports the general idea according to which the cusp-core transformation occurs through cycles of inflows and outflows and further justifies trying to develop a more precise and more predictive model.

2.2 Gravitational potential of a non-isolated shell

Dutton et al. (2016b) only take into account the enclosed mass in the potential energy although the mass distribution in the outer shells matters because its response to the central mass change affects the potential at the shell in question. Indeed, the potential energy per unit mass $U(r)$ is minus the energy needed to extract a unit mass from $r$ to infinity, or equivalently the work done by the gravitational force between $r$ and infinity:

[TABLE]

assuming that the halo is truncated at a certain radius $R_{t}$ enclosing a mass $M_{t}$. Integrating by parts leads to

[TABLE]

hence $U(r)$ is the sum of a contribution corresponding to the shells inside $r$, $U_{\rm isolated}=-GM(r)/r$, and of a contribution of the shells outside $r$ given by the remaining integral.

Except when applied to the outer layer at $R_{t}$, the gravitational potential is not equal to $U_{\rm isolated}$ as was assumed by Dutton et al. (2016b), and the contribution of the outer shells is always negative: adding the outer shells leads to a deeper potential as it becomes more difficult to extract a particle from any radius owing to the additional mass. As an example, Fig. 1 shows how taking into account the outer layers affects the potential energy per unit mass of test particles at different radii in the cases of a singular isothermal sphere (SIS), of an NFW profile and in a simulated galaxy from the NIHAO sample (g1.08e11; Wang et al., 2015). The figure focuses on simulated outputs after 3 Gyr in order to put aside the more erratic first outputs corresponding to the building up of the halo. As can be seen in the figure, neglecting the outer layers induces a shallower gravitational potential. More importantly, changes in the outer mass distribution change the potential at the radius in question. In fact, the bottom panel of Fig. 1 also shows that most of the changes in the simulated shell potential $U(r)$ come from changes in the outer mass distribution. Indeed, $U_{\rm isolated}(r)$ is more than an order of magnitude smaller than $U(r)$ in absolute value below $0.1R_{\rm vir}$, where changes in the potential are the most significant. In Section 3.2, we indicate the expression of the potential derived from Eq. (4) for a density profile parametrized following Dekel et al. (2017).

2.3 Shell kinetic energy and virial theorem

The scalar virial theorem (e.g., Binney & Tremaine, 2008, section 4.8.3) states that the total kinetic energy of a collisionless system in a steady equilibrium state without surface pressure is half its total potential energy. Dutton et al. (2016b) assume that this is still valid for a shell at radius $r$, which leads them to write the kinetic energy per unit mass of that shell as

[TABLE]

where $M$ is the total mass enclosed within $r$ and where they neglect the contribution of the outer shells when expressing the potential. Fig. 2 compares the local kinetic energy per unit mass in the simulated NIHAO galaxy g1.08e11 with $K_{\rm isolated}$ and $K_{\rm virial}=U/2$, where the potential energy per unit mass $U$ stems from a numerical integration of Eq. (5) as in Fig. 1. The local kinetic energy per unit mass associated with a shell at radius $r$ was computed as

[TABLE]

where $K(<r)$ is the enclosed kinetic energy and $M(r)$ the enclosed mass222 More precisely, $K$ is determined using Eq. (7) with ${K(<r)}=\sum m_{i}v_{i}^{2}/2$ and $M(r)=\sum m_{i}$ summing over all particles within radius $r$; the derivatives $dK(<r)/dr$ and $dM/dr$ are numerically obtained with a Savitzky-Golay smoothing filter (Savitzky & Golay, 1964, implemented as scipy.signal.savgol_filter in scipy) with a polynomial order $n=3$ and a window size $w=11$ (optimized to smooth out the radial oscillations of the profile and to capture its average shape). . As can be seen in the figure, neither $K_{\rm isolated}$ nor $K_{\rm virial}$ follow the simulated kinetic energy $K$. The discrepancy between $K$ and $K_{\rm virial}$ shows that the expression of the kinetic energy derived from the equation of virial equilibrium is not valid for shells embedded in a larger halo, while the discrepancy between $K$ and $K_{\rm isolated}$ shows the limits of this latter expression, in particular towards the center of the halo at $r<0.05R_{\rm vir}$, where we explore the cusp-core evolution. The discrepancy between $K$ and $K_{\rm virial}$ can not be explained by non-isotropic velocities, since the anisotropy parameter is close to zero at most radii and in particular towards the centre (cf. Appendix B.1, Fig. 16). In section 3.3, we carefully adress the parametrization of the kinetic energy, notably by obtaining an analytical expression deriving from the Jeans equation for a density profile parametrized following Dekel et al. (2017).

3 Halo response model

3.1 Instant mass change

We aim at describing the response of a spherical collisionless halo with an initial mass profile $M_{i}(r)$ when mass is added (or removed) instantaneously about the center with a profile $m(r)$, and we wish to determine the final collisionless halo mass profile $M_{f}(r)$. Positive values of $m$ correspond to inflows, negative values to outflows. As in Dutton et al. (2016b), we consider the evolution of a spherical halo shell enclosing a given collisionless mass $M$, initially at radius $r_{i}$ and eventually at radius $r_{f}$ so that $M=M_{i}(r_{i})=M_{f}(r_{f})$.

We assume that the shell is initially at equilibrium (stage 1), its total energy per unit mass being

[TABLE]

where $U_{i}$ and $K_{i}$ are respectively the gravitational potential and kinetic energy per unit mass. Contrarily to Dutton et al. (2016b), we neither neglect the outer shells in the expression of the potential nor assume that the kinetic energy is that derived from the equation of halo virial equilibrium. We keep functional forms for $U_{i}$ and $K_{i}$ for the moment, which will become explicit in the following sections 3.2 and 3.3. The potential $U_{i}(r_{i})$ can be derived from the mass profile $M_{i}(r)$ while the kinetic energy $K_{i}(r)$ is set by Jeans equilibrium.

When a given mass $m$ is instantaneously added (or removed, if $m<0$) inside $r_{i}$, we assume an intermediate transitional state (stage 2) in which the potential immediately adapts to the mass variation while the velocities are frozen to their initial values and the shell still at radius $r_{i}$. This stage notably corresponds to what happens to particles on Keplerian orbits subjected to a sudden change in the central mass, as in Penoyre & Haiman (2018). The total energy of the shell right after the mass variation is

[TABLE]

where we further assume that the mass inflow (or outflow) is spherically symmetric and limited to within $r_{i}$ so that its contribution to the potential is simply $-Gm/r_{i}$. This expression for the potential is rigorously valid when mass is added (or removed) directly at the center of the halo, but it can also describe an inflow/outflow mass profile $m(r_{i})$ provided that the shell of interest is outside the main body of the mass variations, i.e., that the contribution of the mass added (removed) at $r>r_{i}$ to the potential is small compared to $Gm/r_{i}$. In Section 3.3.3, we will use the same assumption to express the contribution of the added (removed) mass to the local kinetic energy.

The system at this stage is out of equilibrium, and we assume that it subsequently relaxes to a new equilibrium state (stage 3) where the total energy per unit mass of the shell, now at $r_{f}$, is

[TABLE]

In this expression, the contribution to the potential of the additional mass has been singled out and the local kinetic energy $K_{f}$ is set by the new equilibrium state.

We now postulate that the total energy is the same for corresponding shells before and after the relaxation, $E_{f}(r_{f})=E_{t}(r_{i})$. If there is no shell crossing following the sudden mass change (e.g., as expected in the case of contraction due to added mass), we can have $r_{f}$ and $r_{i}$ refer to the same Lagrangian shell, and then indeed expect conservation of energy based on the energy conservation in the whole halo and the approximate self-similarity of the halo density profile over a range of radii. When there is shell crossing, we make the ansatz that the equal energy equation is still valid but for shells encompassing the same mass before and after the relaxation even when they do not represent the same Lagrangian shell. This approximation is not formally justified a priori – its validity will be verified by the success or failure of the model in reproducing the response in simulations. We note that in Penoyre & Haiman (2018), shell crossing manifests itself by the presence of caustics in the perturbed radial density profile where shells overlap, which we do not take into account in our parametrisation of the density profile. Combining Eqs. (9) and (10), the energy equality yields

[TABLE]

Given functional forms $U(r;p)$ and $K(r;p,m)$ with parameters $p$, where $K$ stems from Jeans equilibrium, it is possible to solve numerically this equation in order to obtain the final radius $r_{f}$ and the final parameters $p_{f}$ entering in the expressions of $U(r;p)$ and $K(r;p,m)$ from given initial conditions ($r_{i}$, $p_{i}$) and the mass change $m(r_{i})$. Note that the expression of the kinetic energy needs to account for the mass that is added to (or removed from) the system, as we would otherwise be left with the trivial solution $r_{f}=r_{i}$ and $p_{f}=p_{i}$. One can apply Eq. (11) to a set of $N$ values of $r_{i}$ in a certain range of interest (for example towards the center of the halo to describe the cusp-core transition) and thus obtain $N$ constraints that should be approximately fullfilled. By finding the best-fit values for the final parameters $p_{f}$ that minimize the residuals between the two sides of Eq. (11) across all the $N$ radii, the toy model presented here makes a prediction for the evolution of the collisionless density profile when mass is added (or removed) to the halo. The functional forms for the gravitational potential and the kinetic energy can be derived under certain assumptions from the density profile or its parametrisation and from the Jeans equation. In the following section 3.2, we describe one such parametrisation and the resulting expression for the gravitational potential. In section 3.3, we derive the corresponding kinetic energy from the Jeans equation describing the steady-state equilibrium of a collisionless system, assuming complete spherical symmetry and an isotropic velocity dispersion.

To summarize, the key ingredients and main assumptions of the model are: (1) the system is spheri-symmetric with a given profile that is described by a parametric function, and the velocities are isotropic; (2) the initial halo is in Jeans equilibrium; (3) the mass change is instantaneous, with a mass profile $m(r)$ about the halo center; (4) the change puts the halo system in an intermediate out-of-equilibrium state in which the gravitational potential is changed instantaneously while the velocities are frozed at their initial values; (5) the halo relaxes to a new Jeans equilibrium with a final profile that is described by the same function as initially but with new parameters; (6) the energy is the same before and after the relaxation process for shells that enclose a given halo mass.

Fig. 3 illustrates the assumed steps of this model and the evolution of the density profile for a fiducial choice of the initial collisionless mass distribution and a mass $m$ equal to 0.2% of the virial mass directly removed at the center. The initial density profile corresponds to a fit to output 32 ($t=7.14~{}\rm Gyr$) of the NIHAO galaxy g1.08e11 according to Eq. (13) below. The gravitational potential is given by Eq. (14) and the kinetic energy by Eq. (20). As can be seen in the figure, the density profile becomes shallower when the mass is removed: the inner slope at $0.01R_{\rm vir}$ drops from 0.37 to 0.25. A given collisionless mass $M$ initially enclosed within radius $r_{i}$ is enclosed within $r_{f}>r_{i}$ after the relaxation, while the density decreases. Following the model, the equilibrium state of the initial shell enclosing $M$ at $r_{i}$ is described by Eq. (8) (stage 1); the counterpart of this shell during the transitional out-of-equilibrium state verifies Eq. (9) (stage 2), where the potential energy has increased because of the mass removal but the kinetic energy and radius remain those of the initial shell; and the equibrium state of final shell enclosing $M$ at $r_{f}$ is characterized by Eq. (10) (stage 3). The potential energy associated to the final stage is higher than the initial one, but the kinetic energy lower. Minimizing Eq. (11) for shells at different initial radii enables to retrieve both the parameters of the final density profile and the final characteristics of each shell. Fig. 4 proposes a schematic representation of the different stages assumed by the model in the case of a sudden gas outflow towards the center of a dark matter halo. The cartoon is inspired by Fig. 3 but includes a stellar component for illustration.

3.2 Parametrisation of the density profile

Dekel et al. (2017) propose a parametrisation of the density profile of dark matter haloes with a free inner slope and an analytic potential, developed to describe the transition from cusps to cores and alterations of the profile due to environmental effects. This parametrisation relies on expressing the mean density profile within a sphere of radius $r$ as:

[TABLE]

where $\overline{\rho_{c}}$ is a characteristic density, $x=r/r_{c}$ with $r_{c}$ an intermediate characteristic radius, $a$ and $g$ the inner and outer asymptotic slopes and $b$ a middle shape parameter. As the virial radius $R_{{\rm vir}}$ is set by cosmology for a given halo mass, this functional form effectively depends on four shape parameters: $a$, $b$, $g$ and $c=R_{{\rm vir}}/r_{c}$, to be reduced to two free parameters in our analysis. Indeed, the normalisation factor $\overline{\rho_{c}}$ can be expressed as $\overline{\rho_{c}}=c^{3}\mu\overline{\rho_{\rm vir}}$, with $\mu=c^{a-3}(1+c^{1/b})^{b(g-a)}$, and $\overline{\rho_{\rm vir}}$ the mean mass density within $R_{{\rm vir}}$. When $g=3$ and $b=2$, the local density profile is

[TABLE]

with two remaining free parameters ($a$ and $c$). The associated potential energy stemming from Eq. (5) is

[TABLE]

when $a\neq 2$ or $5/2$, with $\chi=x^{1/2}(1+x^{1/2})^{-1}$, $\displaystyle\chi_{\rm vir}=\chi(c)$ and $V_{\rm vir}=GM_{\rm vir}/R_{\rm vir}$, as detailed in Appendix A.

We caution that parameter $a$ is not the slope at the resolution limit ($0.01R_{\rm vir}$ in the case of the NIHAO simulations) and that $c$ does not necessarily reflect the actual concentration of the halo as for an NFW (Navarro et al., 1996, 1997) profile. In fact, these two parameters can differ significantly from the concentration and the inner slope. The logarithmic slope of the density profile when $g=3$ and $b=2$ is indeed

[TABLE]

so negative values of parameter $a$ can be compatible with a positive logarithmic slope at the resolution limit, in particular for large values of $c$ since $x=cr/R_{\rm vir}$. For the NIHAO simulations, we define

[TABLE]

the inner logarithmic slope at the resolution limit. We tried enforcing a positive value for $s_{0}$, but it did not improve the fits so we subsequently prefered to leave both parameters $a$ and $c$ free. Given the expressions for the mass and velocity profiles that can be derived from Eq. (12) (cf. Appendix A), the velocity peaks at $x_{\rm max}=(2-a)^{2}$ when $g=3$ and $b=2$ such that

[TABLE]

defines a more physical concentration parameter than $c$, coinciding with the latter only when $a=1$. General expressions for $s$ and $c_{\rm max}$ are given in Dekel et al. (2017) and recalled in Appendix A. At fixed $b$, $g$, $R_{\rm vir}$ and $M_{\rm vir}$, there is a bijection between the couple $(a,c)$ and $(s_{\rm 0},c_{\rm max})$.

The functional form proposed by Dekel et al. (2017) enables to describe different inner and outer slopes for the density profile, contrarily to the usual NFW and Einasto (1965) profiles, and allows analytical expressions for the density profile, the integrated mass profile as well as for the gravitational potential and the velocity dispersion when $g=3$ and $b$ a natural number. Dekel et al. (2017) show that this functional form yields excellent fits for haloes in simulations with and without baryons, ranging from steep cusps to flat cores. In particular, they show that the parametrisation with $g=3$ and $b=2$ matches simulated profiles at least as well as the usual NFW and Einasto profiles but captures cores better in addition to providing fully analytical expressions for the density, the mass and the gravitational potential. This motivates us to adopt the Dekel profile with $g=3$ and $b=2$ to follow the transition from cusps to cores. As an example, Fig. 5 compares this density profile with the NFW profile for one output of NIHAO galaxy g1.08e11, illustrating how it enables to better account for the inner part of the density profile. In Fig. 6 we show that the expression of the gravitational potential derived from the Dekel et al. (2017) parametrisation of the dark matter density profile (Eq. (14)) reproduces the simulated quantity. There is a small systematic trend for it to underestimate the simulated gravitational potential towards the center due to the presence of stars and gas, but the relative error remains below 10% and its mean RMS below 5%.

3.3 Parametrisation of the kinetic energy

The model presented in section 3.1 relies on analytic functional forms both for the gravitational potential and for the local kinetic energy $K(r)$ under Jeans equilibrium. Assuming complete spherical symmetry and isotropy without global rotation, we derive the local kinetic energy $K_{\rm DM}$ for a single-component halo described by a Dekel et al. (2017) density profile with $b=2$ and $g=3$ in steady-state Jeans equilibrium. Using a simple parametrisation of the total mass profile, we further account for the multiple components (dark matter, stars and gas) of simulated galaxies and derive the corresponding local kinetic energy $K_{\rm multi}$, which is found to be in agreement with the simulations. Additional terms can potentially be added to account for the presence of other components, for example a central mass or uniformly spread baryons.

3.3.1 Spherical symmetry and isotropy

The local kinetic energy per unit mass of a spherical shell at radius $r$ within a spherical collisionless system, as defined in Eq. (7), can be written as

[TABLE]

where $\sigma_{r}$ is the radial component of the velocity dispersion and $\beta$ the anisotropy parameter (e.g., Binney & Tremaine, 2008, Section 4.8.3). This expression assumes complete spherical symmetry, i.e., that not only the potential and density profiles are spherically symmetric but also the distribution function describing the collisionless system such that there is no global rotation of the halo. It yields $\displaystyle K(r)=1.5\sigma_{r}^{2}$ when $\beta=0$, which is what we assume in the following. In Appendix B.1, we indeed show that assuming complete spherical symmetry and isotropy are valid assumptions for dark matter haloes in the NIHAO simulations.

3.3.2 Jeans equilibrium

The equilibrium of a spherical collisionless system can be described by the spherical Jeans equation stemming from the Boltzmann equation (Binney & Tremaine, 2008, Eq. 4.215), which yields

[TABLE]

when $\beta=0$ and the boundary condition is $\lim_{r\rightarrow+\infty}\sigma_{r}^{2}=0$. We show in Appendix B.2 that the Jeans equation describing the steady-state equilibrium of a spherical collisionless system is valid in the NIHAO simulations within a mean RMS difference of 13% for $r$ between $0.02R_{\rm vir}$ and $R_{\rm vir}$. For a single-component spherical halo whose density profile is given by Eq. (13), Eq. (19) yields a local kinetic energy per unit mass

[TABLE]

where $\mathcal{B}(a,b,x)=\int_{0}^{x}t^{a-1}(1-t)^{b-1}dt$ is the incomplete beta function and the outer brackets refer to the difference of the enclosed function between 1 and $\chi$, i.e. $\displaystyle\left[f(\zeta)\right]_{\chi}^{1}\equiv f(1)-f(\chi)$ for any function $f$. Note that although the incomplete beta function is formally only defined when its parameters are positive, we extend the notation to negative parameters since Eq. (20) arises from a finite integral that ensures the difference to be finite too. In fact, Eq. (20) can also be expressed in terms of different finite series, whose expressions can be found in Appendix B.3.1.

In the case of a multi-component halo as in the simulations, the enclosed mass entering the gravitational term in the right-hand side of Eq. (19) should include all components (dark matter, stars and gas). In Fig. 7, we show how the total enclosed mass $M_{\rm tot}$ differs from that of the dark matter component $M$ in the different outputs of NIHAO galaxy g1.08e11 after 3 Gyr. To account for the difference between $M_{\rm tot}$ and $M$, we model their ratio as a power law

[TABLE]

where $X_{M}$ and $n$ are ajustable parameters, and we replace $M$ by $M_{\rm tot}$ in the integral of Eq. (19). In the case of the galaxy shown in Fig. 7, the best-fit parameters to the median ratio are $X_{M}=1.03$ and $n=0.07$. Other approaches would be possible (e.g., solving Eq. (20) numerically or running the model on the three halo components simultaneously), but this parametrisation keeps the model simple, analytic and cost-effective in terms of computational time. With this parametrisation, the local kinetic energy per unit mass becomes

[TABLE]

(cf. Appendix B.3.2). Fig. 8 shows that this parametric expression of the local kinetic energy reproduces the simulated local kinetic energy $K$ with a relative error lower than 10%, the agreement being particularly good towards the center. The small negative deviation in the relative difference around $0.1R_{\rm vir}$ comes from an accumulation of gas at the outskirts of the stellar disk in this galaxy, which results in a relative decrease of the dark matter fraction and hence in a slight overestimation of the total mass and subsequently also of the local kinetic energy.

3.3.3 Presence of an additional mass

The key element of our model is the introduction of an additional mass $m$ (negative if mass is removed) to the system, which results in an additional term to the velocity dispersion $\sigma_{r}^{2}$ of Eq. (19),

[TABLE]

This term depends on how $m$ is spread over the halo. Assuming that it is dominated at each radius $r$ by the contribution of the mass entering (outflowing) at $r$ itself, i.e., that

[TABLE]

the corresponding contribution to the local kinetic energy per unit mass is

[TABLE]

(cf. Appendix B.3.3). The assumption explicited by Eq. (24) is valid outside the main body of $m(r)$ and in particular for the introduction (removal) of a point mass. If mass is added (removed) only up to a certain radius $r_{m}$, Eq. (24) is indeed rigourously valid for $r\geq r_{m}$. If $m(r)$ further does not oscillate between negative and positive values between the center and $r_{m}$, i.e. if the mass variation is a pure inflow (outflow), the right-hand side of Eq. (24) provides an upper limit to the left-hand integral and hence also to the contribution of the additional (removed) mass to the local kinetic energy. The assumption of Eq. (24) is motivated by a desire to obtain a simple analytic formulation. Alternative assumptions that also achieve this goal include assuming that the mass change is uniformly distributed up to a radius $r_{m}$, that it follows a singular isothermal sphere profile, or that if follows a Dekel et al. (2017) profile as the collisionless halo. Expressions for the corresponding contributions to the local kinetic energy $K_{m}$ are given in Appendices B.3. Another possibility would be to integrate Eq. (23) numerically at the expense of the analytical formulation, which we avoid in the present work.

In brief, the local kinetic energy per unit mass for a single-component spherical collisionless halo described by a Dekel et al. (2017) density profile in Jeans equilibrium in the presence of an additional mass $m(r)$ can be parametrized analytically as $K=K_{\rm DM}+K_{m}$; for a multi-component halo where $M_{\rm tot}/M$ is modeled as a power law in $r/R_{\rm vir}$, it can be parametrized as $K=K_{\rm multi}+K_{m}$. The different expressions for $K_{m}$ considered in Appendix B.3 further enable to describe other mass distributions combining central masses, Dekel et al. (2017) haloes, singular isothermal spheres and uniform spheres.

4 Testing the model with simulations

4.1 Description of the simulation suite

We test our theoretical model using simulated galaxies taken from the Numerical Investigation of a Hundred Astrophysical Objects project (NIHAO; Wang et al., 2015), which provides a set of about 90 cosmological zoom-in hydrodynamical simulations ran with the improved Smoothed Particle Hydrodynamics (SPH) code gasoline2 (Wadsley et al., 2004; Keller et al., 2014; Wadsley et al., 2017). For each hydrodynamical simulation, a corresponding dark matter only simulation was run at the same resolution to study how baryonic processes affect the structure of dark matter haloes. The simulations assume a flat $\Lambda$CDM cosmology with Planck Collaboration et al. (2014) parameters, namely $\Omega_{m}=0.3175$, $\Omega_{r}=0.00008$, $\Omega_{\Lambda}=1-\Omega_{m}-\Omega_{r}=0.6824$, $\Omega_{b}=0.0490$, $H_{0}=67.1~{}\rm km.s^{-1}.Mpc^{-1}$, $\sigma_{8}=0.8344$ and $n=0.9624$. They include a subgrid model to account for the turbulent mixing of metals and thermal energy (Wadsley et al., 2008), cooling via hydrogen, helium and other metal lines in a uniform ultraviolet ionizing and heating background (Shen et al., 2010) and star formation according to the Kennicutt-Schmidt relation when the temperature is below $15000~{}\rm K$ and the density above $10.3~{}\rm cm^{-3}$ (Stinson et al., 2013). Stars inject energy back to their surrounding interstellar medium (ISM) both through ionizing feedback from massive stars (Stinson et al., 2013) and through supernova feedback ejecting both energy and metals into the ISM (Stinson et al., 2006). The pre-supernova feedback consists in ejecting 13% of the total stellar luminosity, which is typically $2.10^{50}~{}\rm erg$ per $\rm M_{\odot}$ of the entire stellar population over the 4 Myr preceding the explosion of high-mass stars as supernovae, into the surrounding gas. The relatively high 13% efficiency is set to match the Behroozi et al. (2013) abundance matching results despite the increased mixing of gasoline2. During the second supernova feedback phase, stars whose mass is comprised between 8 and 40 $\rm M_{\odot}$ eject 4 Myr after their formation both an energy $E_{\rm SN}=10^{51}~{}\rm erg$ and metals into their surrounding ISM according to the blast-wave formalism described in Stinson et al. (2006). In this formalism, the total feedback energy is distributed amongst gas particles within the maximum blastwave radius from Chevalier (1977) and McKee & Ostriker (1977). In order to prevent the energy from supernova feedback to be radiated away, cooling is delayed for 30 Myr inside the blast region. Without cooling, the added supernova energy heats the surrounding gas, which both prevents star formation and models the high pressure of the blastwave. The NIHAO sample comprises isolated haloes chosen from dissipationless cosmological simulations (Dutton & Macciò, 2014) with halo masses between $\rm\log(M_{{\rm vir}}/M_{\odot})=9.5-12.3$ but without considering their halo merging histories, concentrations and spin parameters. Most importantly, the particle masses and force softening lengths are chosen to resolve the dark matter mass profile below 1% of the virial radius at all masses in order to resolve the half-light radius of the galaxies. Stellar masses range from $5.10^{4}$ to $2.10^{11}~{}\rm M_{\odot}$, i.e., from dwarfs to Milky Way sized galaxies, with morphologies, colors and sizes that correspond well with observations (e.g., Wang et al., 2015; Stinson et al., 2015; Dutton et al., 2016a).

We test our model on the simulated NIHAO galaxies whose stellar mass at $z=0$ lies between $5\times 10^{7}$ and $5\times 10^{9}~{}\rm M_{\odot}$, which is the range where we expect the most important changes in the dark matter density profile according to Tollet et al. (2016) and Dutton et al. (2016b). Amongst these galaxies, we further use g1.08e11 and g6.12e10 as a fiducial cases given their intermediate stellar and halo masses at $z=0$, $M_{\star}=8.47\times 10^{8}~{}\rm M_{\odot}$ and $M_{\rm vir}=1.20\times 10^{11}~{}\rm M_{\odot}$ for the former and $M_{\star}=9.13\times 10^{7}~{}\rm M_{\odot}$ and $M_{\rm vir}=5.50\times 10^{10}~{}\rm M_{\odot}$ for the latter, their relatively quiet merging history, and the fact that they belong to the UDGs sample studied by Di Cintio et al. (2017) and Jiang et al. (2019). Their stellar half light radius and V-band central surface brightness at $z=0$ are indeed respectively $r_{e}=4.4~{}\rm kpc$ and $\mu_{0,V}=24.1~{}\rm mag.arcsec^{-2}$ for g1.08e11 and $r_{e}=2.6~{}\rm kpc$ and $\mu_{0,V}=24.03~{}\rm mag.arcsec^{-2}$ for g6.12e10. As they are UDGs, we expect both the dark matter and the stellar systems to expand and develop a core.

4.2 Implementation of the model in the simulations

For each simulation output, we first derive global properties of the central galaxy identified with the Amiga Halo Finder (AHF; Gill et al., 2004; Knollmann & Knebe, 2009). The virial radius $R_{\rm vir}$ is defined as the radius within which the average total density is $\Delta_{c}$ times the critical density of the Universe, where $\Delta_{c}$ is defined according to Bryan & Norman (1998). The virial mass $M_{\rm vir}$ is the total mass enclosed within $R_{\rm vir}$. The stellar mass $M_{\star}$ is calculated within $0.15R_{\rm vir}$; the stellar half-light radius $r_{1/2}$ corresponds to the sphere enclosing half of $M_{\star}$. The central V-band magnitude $\mu_{V,0}$ is computed within $0.25R_{1/2}$ using Padova Simple stellar populations (Marigo et al., 2008; Girardi et al., 2010) as implemented in pynbody (Pontzen et al., 2013). Galaxies with $r_{1/2}>1.5~{}\rm kpc$ and $\mu_{V,0}>24~{}\rm mag.arcsec^{-2}$ are considered as UDGs as in Jiang et al. (2019). The timestep between two successive simulation outputs is 216 Myr.

We further derive radial profiles of mass ($M$), density ($\rho$), mean density ($\overline{\rho}$), radial and tangential velocity dispersions ($\sigma_{r}$ and $\sigma_{t}$) and anisotropy parameter ($\beta=1-\sigma_{t}^{2}/2\sigma_{r}^{2}$) using the pynbody.analysis.profile module for the different components constituting the halo, i.e., gas, stars and dark matter. The profile radii $r$ are spaced logarithmically with $\sim$100 radii between $0.01R_{\rm vir}$ and $R_{\rm vir}$. For each radius $r$, we compute the total mass including all components flowing in ($m_{\rm in}$) and out ($m_{\rm out}$) between all successive simulation outputs as well as the corresponding net mass change $m=m_{\rm in}-m_{\rm out}$, which is the main input of our theoretical model. We also define $f=m/M$ the fraction between this mass change and the enclosed mass of each component. We fit the average dark matter density profiles between $0.01R_{\rm vir}$ and $R_{\rm vir}$ according to the Dekel et al. (2017) parametrization (Eq. (12)) through a least-square minimization with python’s lmfit package, imposing $R_{\rm vir}$, $M_{\rm vir}$, $b=2$ and $g=3$ while leaving parameters $a$ and $c$ free. The inner slope at the resolution limit $s_{0}=s(0.01R_{\rm vir})$ and the effective concentration parameter $c_{\rm max}$ can be derived from Eqs. (15) and (17) for each set of parameters.

We implement the instant mass change theoretical model presented in section 3.1 in all successive simulation outputs excluding the first 3 Gyr of evolution, which are relatively perturbed. As mergers and fly-bys are expected to break the assumed spherical symmetry and to affect the halo matter distribution through processes that are not accounted for in the model, we define a merger-free subsample of successive outputs by excluding cases where the ratio between the total mass change including all components (dark matter, stars and gas) and the total enclosed mass at $0.15R_{\rm vir}$,

[TABLE]

is above $f_{\rm merger,min}=10\%$ during the timestep at stake. Amongst the processes potentially induced by mergers and fly-bys, dynamical friction can heat the background particles (e.g. El-Zant et al., 2001; El-Zant et al., 2004; Tonini et al., 2006; Goerdt et al., 2010) and gravitational torques can drive large amounts of gas towards the center in a non-spherical way (Barnes & Hernquist, 1991; Mihos & Hernquist, 1996; Fensch et al., 2017).

As presented in section 3.3, we account for the difference between the total enclosed mass $M_{\rm tot}$ and the collisionless component mass $M$ through a power-law radial dependence according to Eq. (21). We derive the best-fit parameters $X_{M}$ and $n$ at each output and assume that this parametrization holds until the next output time. The principle of our model was outlined in section 3.1: for a given best-fit parametrization of the collisionless density profile, a power-law parametrization of $M_{\rm tot}/M$ and the net mass change $m$ from the simulation, we make a prediction for the new collisionless density profile through a least square minimization of the energy equation $E_{f}(r_{f})=E_{t}(r_{i})$ explicited in Eq. (11). The minimization is carried out simultaneously at all radii between $0.02R_{\rm vir}$ and $R_{\rm vir}$ (or $\log(r/R_{\rm vir})$ between $-1.75$ and [math]), since this is the domain where the Jeans equation is valid, as shown in Appendix B.2. The final density profile is parametrized according to the Dekel et al. (2017) functional form. As for the profile fits, its parameters $a$ and $c$ are allowed to vary while $g=3$ and $b=2$. $R_{\rm vir}$ and $M_{\rm vir}$ are fixed to their initial values.

Fig. 9 shows an example evolution of the dark matter average density profile between two successive outputs of the NIHAO g1.08e11 simulation compared to the instant mass change model prediction. The bottom panel further shows the radial profile of the mass ratio $f=m/M$ for the dark matter component stemming from the simulation that is implemented in the model. As can be seen in the figure, the dark matter average density profile evolves towards a flattening of its central part. This evolution agrees very well with the model prediction (red line in the upper panel) and is consistent with a mass change dominated by outflows (red portions of the curve in the bottom panel).

4.3 Assessing the success of the model

To quantitatively assess the success of the model in retrieving the inner density profile, we define $\mathcal{A}_{k}$, $\mathcal{A}_{k+1}$ and $\mathcal{A}_{\rm model}$ the areas between the horizontal line $\log(\overline{\rho}/\overline{\rho}_{0})=-1$ and the Dekel et al. (2017) fit to the average density profiles $\log(\overline{\rho}/\overline{\rho}_{0})$ corresponding respectively to the initial output ($t=6.7$ Gyr in the example shown in Fig. 9), the final output ($t=6.9$ Gyr) and to the model, $\log(r/R_{\rm vir})$ being between $-2$ and $-1.5$. The shaded areas of Fig. 9 help visualizing these areas, with $\mathcal{A}_{k}$ from the bottom to the gray dashed line, $\mathcal{A}_{k+1}$ from the bottom to the black dashed line and $\mathcal{A}_{\rm model}$ from the bottom to the red plain line. The absolute ratio

[TABLE]

between $\mathcal{A}_{\rm model}-\mathcal{A}_{k+1}$ and $(\mathcal{A}_{\rm model}+\mathcal{A}_{k+1})/2$ provides a relative measurement of the difference between the model prediction and the final output, which we use to assess the success of the model. Although rather arbitrary, the choice of the horizontal line $\log(\overline{\rho}/\overline{\rho}_{0})=-1$ corresponds well with our visual assessment of the successes.

We similarly define $\delta_{0}$ to measure the difference between $\mathcal{A}_{\rm model}$ and $\mathcal{A}_{k}$ (instead of $\mathcal{A}_{k+1}$), i.e. between the model prediction and the initial output, and $\delta_{\rm sim}$ between $\mathcal{A}_{k}$ and $\mathcal{A}_{k+1}$, i.e. between the initial and the final outputs. If $\delta>\delta_{0}$, the model prediction is closer to the initial density profile than to the final one and hence fails to retrieve the simulated evolution. However, if $\delta_{\rm sim}$ is particularly small, the comparison between $\delta$ and $\delta_{0}$ loses its relevance. We thus consider the model successful at predicting the evolution of the average collisionless density profile when (i) $\delta$ is below a certain threshold $\delta_{\rm max}$ and (ii) $\delta<\delta_{0}$ unless $\delta_{\rm sim}$ is too small to differenciate between $\delta$ and $\delta_{0}$. In terms of logical operators, the success condition we choose can be expressed as

[TABLE]

where $\delta_{\rm min}$ corresponds to the minimum discrepancy between the two outputs we deem physically meaningful. In the following, we take $\delta_{\rm max}=10\%$ as success threshold and $\delta_{\rm min}=3\%$. This value of $\delta_{\rm max}$ enables to consider the case shown in Fig. 9 as a clear success. A different threshold would change our results quantitatively (a lower value of $\delta_{\rm max}$ leading to fewer reported successes and vice-versa) but would not change the qualitative trends reported in the next sections.

The Dekel et al. (2017) parametrization of the average density profile further enables to measure the success of the model in terms of the errors in inner logarithmic slope

[TABLE]

and effective concentration parameter

[TABLE]

Since the success measure $\delta$ aims at capturing changes in the inner part of the average density profile, it correlates strongly with $\Delta s$ (with a Pearson correlation coefficient $r=0.71$ between $\delta$ and $\log|\Delta s|$), less with $\Delta c$ ($r=0.43$ between $\delta$ and $\log|\Delta c|$), while $\Delta s$ and $\Delta c$ are themselves moderately correlated ($r=0.53$ between the logarithms of their absolute values). Figures showing these correlations can be found in Appendix C (Fig. 21).

Table 1 summarizes the fiducial parameters used for fitting the average density profiles, implementing the model and assessing its success.

4.4 Model success results

Fig. 10 shows the evolution of the inner part of the average dark matter average density profile between successive outputs of the fiducial NIHAO simulations g1.08e11 and g6.12e10 after the perturbed earliest gigayears. Only cases where the RMS value of the mass ratio between $0.01R_{\rm vir}$ and $R_{\rm vir}$, $|f|_{\rm RMS}$, is above a certain threshold $f_{\rm min}=7\%$ are shown, since these are the cases for which we expect a significant change in the collisionless density profile. With the success criterion indicated above (Eq. (28)), the success rate of the model amongst the cases with $|f|_{\rm RMS}>7\%$ shown in Fig. 10 is 79%, with 22 successes out of 28. We note that the model systematically fails during the three merger cases. Excluding mergers, the success rate when $|f|_{\rm RMS}>7\%$ becomes 88% (22/25). Amongst the cases with $|f|_{\rm RMS}\leq 7\%$ that are not shown in Fig. 10, the success rate is 89% (48/54), but these cases display little evolution. Including all successive outputs of both g1.08e11 and g6.12e11, the overall success rate is 82% (68/83) and 88% (68/77) without mergers. The mass change $m$ is the main input to the model with the initial parametrization of the density profile and the ratio $M_{\rm tot}/M$ between the total mass and that of the dark matter component. Fig. 19 in the appendices shows the radial profiles of the ratio $f=m/M$ between the total mass change including all components and the enclosed dark matter mass together with its RMS average $|f|_{\rm RMS}$ for the cases shown in Fig. 10. We note that mergers coincide with large values of $|f|_{\rm RMS}$.

For the whole sample of simulated galaxies whose stellar mass lies between $5\times 10^{7}$ and $5\times 19^{9}~{}\rm M_{\odot}$ at $z=0$, the success rate given the chosen success criterion is $p_{\rm all}=67\%$ overall and $p_{\rm no~{}mergers}=74\%$ excluding mergers. In this latter case (excluding mergers), the success rate is $p_{>f_{\rm min}}=65\%$ when $|f|_{\rm RMS}\geq f_{\rm min}$ and $p_{<f_{\rm min}}=80\%$ when $|f|_{\rm RMS}<f_{\rm min}$. Table 2 in the appendices summarizes these trends galaxy per galaxy and for the whole sample. The standard deviation of the success rate is about 20%, with significant variations from one galaxy to another. Successful evolutions according to Eq. (28) systematically correspond to

[TABLE]

which can thus provide a relatively equivalent success criterion (cf. also Fig. 21). In contrast, $\Delta c$ fluctuates more since it captures the overall dark matter density profile and not specifically its inner part. The success criterion defined by Eq. (31) leads to similar success rates as indicated above, with 68% cases with $|\Delta s|\leq 0.10$ overall and 73% excluding mergers.

Since the halo dynamical time after 3 Gyr spans over more than 3 time steps, the merger criterion may not exclude all the outputs affected by mergers. Considering that mergers are galaxies whose $f_{\rm merger}$ is above 10% both during the evolution at stake (outputs $k\rightarrow k+1$) and during the previous one (outputs $k-1\rightarrow k$) leads to $p_{\rm no~{}mergers}=76\scalebox{1.0}{\%}$. Further extending this criterion to another previous evolution (outputs $k-2\rightarrow k-1$) leads to $p_{\rm no~{}mergers}=77\scalebox{1.0}{\%}$. These alternative merger criteria remove some failures of the model that coincide with mergers but their incidence on the success rates is small and we thus prefer to retain our initial, simpler merger criterion. We note that the success rate for the fiducial galaxies g1.08e11 and g6.12e10 are higher than average since they were selected to have a relatively quiet merging histories. In the next section, we indeed show that the success rate of the model primarily depends on the merger indicator $f_{\rm merger}$.

Fig. 11 compares the changes in $s_{\rm 0}$ and $c_{\rm max}$ between two successive simulation outputs with those predicted by the model. The relatively high Pearson correlation coefficients show that the model captures the overall evolution of the dark matter density profile, but the non-negligible scatter hinders the capacity of the model to yield precise, quantitative predictions for the profile parameters. Both points are already illustrated in Fig. 10 for two of the simulated galaxies, since the predicted average density profiles is in most cases comparable to the final ones but not fully overlaping. In this figure, the correlations between $\delta$, $\Delta s$ and $\delta c$ (cf. Appendix C) are apparent, since $|\Delta s|>0.10$ always coincides with failures while all successes but one have $|\Delta c|<0.20$. For merger cases, the variations in $s_{\rm 0}$ and $c_{\rm max}$ are on average bigger than for non-merger cases, but they are still correlated to the actual variations. For non-merger cases, the relation between $s_{\rm 0,model}-s_{\rm 0,i}$ and $s_{\rm 0,f}-s_{\rm 0,i}$ is significantly tighter than for the whole sample (and for merger cases) as indicated both by the correlation coefficients and the scatter. Although the correlation coefficient of the relation between $\log(c_{\rm max,model}/c_{\rm max,i})$ and $\log(c_{\rm max,f}/c_{\rm max,i})$ is higher for merger cases, the scatter is significantly higher as well. We note that the quantities plotted against each other in Fig. 11 are not formally independent (they both involve $s_{\rm in,i}$ in the left panel and $\log(c_{\rm max,i})$ in the right panel). If we consider directly the relation between $s_{\rm 0,model}$ and $s_{\rm 0,f}$ (and between $\log(c_{\rm max,model})$ and $\log(c_{\rm max,f})$), the Pearson correlation coefficient yields 0.85 (0.95) but this is mostly due to the fact that the change between two successive outputs is relatively small compared to the initial slope (concentration). The null hypothesis of no change between two successive outputs would indeed yield a strong correlation between $s_{\rm 0,model}$ and $s_{\rm 0,f}$ (or between $\log(c_{\rm max,model})$ and $\log(c_{\rm max,f})$) since $s_{\rm 0,f}$ and $s_{\rm 0,i}$ are correlated with a correlation coefficient 0.73 (0.88 between $\log(c_{\rm max,f})$ and $\log(c_{\rm max,i})$). In the next section, we characterize successes and failures in terms of $\delta$, $\Delta s$ and $\Delta c$ and try to identify the main causes of failure of the model to characterize its validity domain and suggest possible future improvements.

5 Discussion

5.1 Mergers as the main cause of failure

To identify systematic causes for the failures of the model after the first 3 Gyr of evolution, we follow how the success measure $\delta$ (and the errors in inner slope and concentration, $\Delta s$ and $\Delta c$) depends on the following quantities relevent to disk and halo dynamics: time, masses of the dark matter, stellar and gas components, stellar size, gas fraction, as well as the merger indicator $f_{\rm merger}$ and a measure of the deviation from Jeans equilibrium, $\Delta_{\rm Jeans}$. This latter quantity corresponds to the RMS average of the relative difference between the two sides of the Jeans equation (Eq. 43) for the dark matter component between 2% and 60% of the virial radius, namely $\log(r/R_{\rm vir})$ between $-1.75$ and $-0.25$. This radius range covers most of the halo but avoids regions where the anisotropy parameter $\beta$ and the logarithmic slopes $\alpha$ and $\gamma$ entering the Jeans equation fluctuate more significantly (cf. Fig.16 and the right panel of Fig. 18). Given the success condition and threshold, we further derive the probability $p$ for $\delta\leq\delta_{\rm max}$ in equally-spaced bins of these quantities to highlight the systematic trends: in each bin, $p$ is the ratio of the number of cases with $\delta\leq\delta_{\rm max}$ to the total number of cases in the bin333The probability $p$ differs from the success rates defined in Section 4.4 since it does not take into account whether the model prediction is closer to the final density profile than to the initial one, i.e., it does not take into account the second part of the success criterion explicited in Eq. (28).. The figure with the different correlations is shown in Appendix C (Fig. 20). We find that $\delta$ correlates primarily with $f_{\rm merger}$ (with a Pearson correlation coefficiant $r=0.56$) and then with $\Delta_{\rm Jeans}$** ($r=0.39$)**, time ($r=0.28$), and the gas-to-virial mass ratio $\mu_{\rm gas/vir}=M_{\rm gas}/M_{{\rm vir}}$ ($r=0.25$), which relates to the gas fraction. In contrast, it correlates weakly with the gas mass ($r=0.11$) and the stellar mass ($r=0.15$), and very weakly with the virial mass ($r=0.03$) and the stellar radius ($r=0.03$). Similarly, the errors on the inner slope and concentration correlate first with $f_{\rm merger}$ ($r=0.38$ for $\log|\Delta s|$ and $0.25$ for $\log|\Delta c|$), and then with $\Delta_{\rm Jeans}$ ($r=0.31$ and $r=0.21$), time ($r=0.21$ and $0.17$), $\mu_{\rm gas/vir}$ ($r=0.20$ and $0.14$) and $M_{\rm star}$ ($r=0.13$ and $0.16$). We note that the scatter in the different success measures increases significantly with $\mu_{\rm gas/vir}$. Given these correlations, the main causes for the failures of the model appear to be mergers, deviation from Jeans equilibrium, gas content, and time.

Fig. 12 projects the distribution of $\delta$ onto the planes involving to these four main causes of failure in order to disentangle their effect and highlight their relative importance. Firstly, the figure clearly shows that amongst the four quantities at stake, the main drivers for failures (high values of $\delta$, in red) are mergers and deviations from Jeans equilibrium. Secondly, we note that the four quantities are not uniformly sampled. In particular, the two upper left panels further indicate that mergers are more prevalent at high gas fraction and early time, which contributes to the dependences on $\mu_{\rm gas/vir}$ and $t$. The median value of $\delta$ at $\log\mu_{\rm gas/vir}>-1.75$ is notably significantly more affected by mergers than below. Finally, Fig. 12 shows secondary trends as a function of time and gas content: at fixed $f_{\rm merger}$ or $\Delta_{\rm Jeans}$, $\delta$ slightly decreases with time and increases with $\mu_{\rm gas/vir}$. Since $\mu_{\rm gas/vir}$ is inversely correlated with time (cf. Fig. 20, with $r=-52$), these two secondary trends may have the same physical origin. We further note in Fig. 20 that while the scatter in $\delta$ is fairly uniform with time, it increases significantly with increasing gas content, hinting at a specific role for the gas. Plots comparable to Fig. 12 for $\log|\Delta s|$ and $\log|\Delta c|$ are shown in Appendix C, showing similar trends.

We conclude from this analysis that amongst the considered physical properties, (i) mergers and deviations from dynamical equilibrium are the main markers of failure of the model, (ii) that there are secondary trends in gas content and time, and (iii) that part of these secondary trends stem from the first two markers. As can be seen in Fig. 20, $f_{\rm merger}$ and $\Delta_{\rm Jeans}$ are themselves significantly correlated, with a Pearson correlation coefficient $r=0.47$, and it is thus difficult to fully disentangle their effects. Nevertheless, since mergers and fly-bys not only deviate the system from equilibrium but are also accompanied by significant non-sphericity that breaks another core assumptions of the model, they may constitute a better marker of failure – which is notably reflected in the tighter correlation of $\delta$ with $f_{\rm merger}$ than with $\Delta_{\rm Jeans}$ ($r=0.56$ instead of $0.39$). Mergers are furthermore likely to invalidate energy conservation within shells enclosing a given collisionless mass during the relaxation. In the model, the sudden change of potential indeed imparts mean motions to the collisionless particles that are subsequently dissipated as the system relaxes to its new equilibrium. The change in the mass distribution during the relaxation may thus be envisioned to result from initial oscillations of the shells due to the imparted motions. If these oscillations are large, which would be the case during mergers, shell crossing would prevent energy conservation within shells enclosing the same collisionless mass. Since Pontzen & Governato (2012) and El-Zant et al. (2016) invoke small successive mass changes to explain core formation, the model we propose here may thus also dynamically describe the situations they study and, in principle, provide a simple means to derive the evolution of the corresponding density profiles. The secondary dependences of the success measure on time and gas fraction can be assumed to stem from the gas content, since the dissipationless nature of the gas makes it likely to violate the assumption of energy conservation during the relaxation phase. We also note in Fig. 8 that an accumulation of gas around $0.1R_{\rm vir}$ at the outskirts of the stellar disk induces a deviation in the ratio between the total mass and the dissipationless mass introduced in Section 3.3.2. Such deviations due to the gas distribution affect the parametrization of the kinetic energy and may thus also affect the ability of the model to predict the evolution of the dark matter density profile. Better accounting for the gas component, for example by parametrizing $M_{\rm tot}/M$ directly from the initial mass distribution instead of parametrizing it as in Eq. (21), may thus help improve the model in the future. We also further expect deviations from spherical symmetry to affect the success rate of the model independently of mergers.

5.2 Contribution to the global change in inner density slope

Since our model describing the evolution of the density profile from inflow and outflow episodes fails during mergers, we assess the contribution of time steps devoid of mergers to the overall evolution of the inner slope $s_{\rm 0}=s(0.01R_{\rm vir})$. In this effect, we consider the inner slope variation between two successive time steps $\Delta s_{\rm 0}(t)=s_{\rm 0}(t)-s_{\rm 0}(t-\Delta t)$, where $\Delta t$ is the time step, and its absolute value $|\Delta s_{\rm 0}|$.

We find that merger-free time steps contribute to 71% of $\Sigma\Delta s_{\rm 0}$ the total evolution of the inner slope for all galaxies of the sample (55% for g1.08e11) and 81% after 3 Gyr (83% for g1.08e11). Given these numbers and the high success rate of the model in the non-merger cases, our relatively simple model may thus be relevent to explain about 80% of the inner slope variation after 3 Gyr of galaxies whose stellar mass at $z=0$ lies between $5\times 10^{7}$ and $5\times 10^{9}~{}M_{\odot}$ and which are making cores. If instead of $\Sigma\Delta s_{\rm 0}$ we consider $\Sigma|\Delta s_{\rm 0}|$ the sum of all inner slope variations in absolute value, we find that merger-free time steps contribute to 66% for all galaxies of the sample (71% for g1.08e11) and 70% after 3 Gyr (94% for g1.08e11).

However, in the current stage of the model, successful non-merger cases after 3 Gyr correspond only to 50% of $\Sigma\Delta s_{\rm 0}$ (36% for g1.08e11, where all the failures after 3 Gyr correspond to a decrease of $s_{0}$) and to 59% of $\Sigma|\Delta s_{\rm 0}|$ (78% for g1.08e11). In the next section, we discuss the ability of the model to predict the evolution of the dark matter density profile over multiple time steps. We leave the description of the global time evolution of the inner logarithmic slope and concentration to future work, together with a sytematic study of the Dekel et al. (2017) parametrisation of the dark matter density profile in the NIHAO simulations.

5.3 Multiple episodes

The model developed in Section 3 was tested so far on successive outputs of the NIHAO simulations, separated by a timestep $\Delta t=216~{}\rm Myr$. We showed in Section 4.4 that this model was successful at predicting the evolution of the dark matter density profile in about 74% of non-merger cases given the success criterion in terms of $\delta$ stated in Eq. (28) (or similarly that in terms of $\Delta s$ indicated in Eq. (31)) and in Section 5.1 that mergers were the main cause of failure. In the previous Section 5.2, we further showed that a large part of the evolution of the inner slope of the dark matter density profile (about 70%) occured in the absence of mergers and hence in the validity domain of the model. However, to what extend is the model relevent to predict the evolution of the density profile over times larger than one timestep? Although we leave the full analysis to future work, we explore here the relevance of the model over multiple timesteps.

In this effect, we consider all suites of outputs within the sample, differentiating those that are interupted by mergers and those that are not. For each of these suites, we first apply the model to the initial dark matter distribution with the first mass change over $\Delta t$, then to the resulting prediction with the second mass change, and so on to the successive predictions. If $t_{0}$ is the initial time, we compare the dark matter density profiles parametrised by the successive parameters $p_{n}$ at $t=t_{0}+n\Delta t$ predicted by the model with the actual profiles.

Fig. 13 shows the time evolution of the inner logarithmic slope $s_{0}$ for simulated galaxies g1.08e11 and g6.12e10 together with two suites of successive predictions of the model for each galaxy, one uninterrupted by mergers and one including a merger. The correlations between the different curves and their slopes highlight the ability of the model to predict the evolution of the inner dark matter density profile between two successive outputs, as already shown in Section 4.3, while the deviations between the curves relates to the scatter shown in Fig. 11. The suites of successive predictions uninterupted by mergers follow closely the actual evolution of the inner logarithmic slopes over a large number of outputs. In contrast, the stochasticity introduced by the merger in the other suites leads the model predictions to deviate significantly from the actual values. Nevertheless, we note the consistency of the inner slope changes between the different suites.

Fig. 14 shows how the scatters in $\Delta s$ and $\Delta c$, namely the scatter in Fig. 11 or the difference between the curves in Fig. 13, vary with the number of successive outputs $N$ taken into account. Similarly, Fig. 15 shows how the success measure $\delta$ varies with $N$. As already seen in Fig. 13, suites of outputs interrupted by mergers first display much more scatter and a much larger success measure than those that are not interrupted by mergers. Scatters and success measures all increases with $N$ as a result of the accumulated uncertainties, however the evolution is significantly milder for the merger-free suites. Indeed, while the evolution of the scatter in the merger cases is compatible with an evolution in $\sqrt{N}$, it is flatter in the non-merger cases, in particular for $s_{0}$. Finally, the median $\delta$ remains below $\delta_{\rm max}=10\%$ until $N\sim 15$, so the probability associated to $\delta\leq\delta_{\rm max}$ remains above 50% until this $N$, i.e., until about 3 Gyr. We thus conclude that while our simple theoretical model is unable to predict the global evolution of the dark matter density profile over more than a few Gyr because of the cumulative effect of mergers, of the scatter and of the progressive increase of $\delta$, it is able to predict with a non-negligible scatter the evolution of the dark matter density profile from episodes of inflows and outflows up to a few Gyr in the absence of mergers.

6 Conclusion

In this article we presented a simple model describing the response of a dissipationless spherical system to an instantaneous mass change, intended to describe the formation of flat cores in dark matter haloes and UDGs by outflow episodes induced by stellar feedback. The model applies rigorously to a mass added or removed at the center of the system but can be extended to mass changes at different radii under certain assumptions. In addition to describing the effect of bulk mass changes, this model may further describe core formation from small successive mass variations, as envisioned by Pontzen & Governato (2012) and El-Zant et al. (2016). The response of the system is divided into an instantaneous change of the gravitational potential to adapt to the mass variation while the velocities are frozen to their initial values, followed by relaxation to a new equilibrium. The key assumption of the model is that the local energy is conserved during the relaxation phase for shells enclosing a given dissipationless mass. This assumption is not formally justified in case of shell crossing, so its validity is tested by the success of the model in reproducing the simulation results. Other assumptions include spherical symmetry, isotropic velocities, and Jeans equilibrium for the initial and final systems. In order to follow the evolution of the dissipationless density profile with time, we parametrize it through a two-parameter function with a variable inner slope and an analytic gravitational potential profile following Dekel et al. (2017) and we derive analytic expressions for the resulting local kinetic energy. The contribution of the added or removed mass to the kinetic energy is determined assuming that most of it is located inside the spherical dark matter shells that are considered, which is notably valid when mass is added or removed directly at the center. To account for the different components of multi-component systems, such as in hydrodynamical simulations, we further use a simple parametrisation of the ratio between total and dissipationless enclosed masses.

The model is tested against NIHAO cosmological zoom-in simulations (Wang et al., 2015) in the context of dark matter core formation. Firstly, we build upon Dekel et al. (2017) showing that their parametrisation of the dark matter density profile matches the simulated profiles significantly better than the usual NFW parametrisation. Secondly, we test the assumptions of spherical symmetry, velocity isotropy, Jeans equilibrium and the simple parametrization of the ratio between total and dissipationless masses. We show that they enable to recover the local kinetic energy with a relative error lower than 10%. We note however that the assumption of energy conservation in shells encompassing a given collisionless mass during the relaxation phase can not be tested directly in the simulations since it is not possible to separate the two phases of evolution assumed in the model. Finally, we compare the model predictions for the evolution of the dark matter density profile given the mass change during each 216 Myr timestep with the actual simulation outputs for the 31 galaxies whose stellar mass at $z=0$ is comprised between $5\times 10^{7}$ and $5\times 19^{9}~{}\rm M_{\odot}$. This is indeed the range where we expect dark matter core formation according to Tollet et al. (2016) and Dutton et al. (2016b). Focusing on times after the relatively perturbed first 3 Gyr of evolution and excluding mergers, we obtain a $\sim$74% success rate, whether success is defined by a relative measure of the proximity between the inner density profile predicted by the model and that of the simulation output ($\delta\leq 10\%$) or by the accuracy of the predicted inner slope variation ($|\Delta s|\leq 10\%$).

We show that mergers and deviations from Jeans equilibirum are the main cause of failure for the model, together with a small influence of the gas mass fraction – higher mass fractions leading to more failures. Mergers are indeed expected to break the assumed spherical symmetry and perturb the equilibrium of the system, in particular through dynamical friction and gravitational torques. The presence of a large dissipative gas component makes the assumption of energy conservation during the relaxation phase less valid. Nevertheless, the overall success of the model in reproducing the halo response in the simulations justifies a posteriori the assumption of energy conservation during the relaxation phase. Since about 70% of the evolution of the inner dark matter density slope occurs without mergers, our model is relevent to describe a large part of this evolution given its success rate. We further test the model over times larger than 216 Myr, spanning multiple mass change episodes. We find that although the scatter increases with time, the success rate remains above 50% over time intervals below 3 Gyr in the absence of mergers.

The model we propose provides a simple understanding of the formation of dark-matter cores and UDGs by supernova-driven outflow episodes. We note that the time step over which the model was tested (216 Myr) is above the dynamical time at the center of the halo, which may contribute to the success of the model since the system has enough time to relax between two successive outputs. Instead of complex hydrodynamical simulations, idealized dissipationless simulations may further help testing this model and its assumptions in fully controlled experiments both for single mass changes and for multiple episodes. Possible improvements and extensions of the model were suggested throughout the text and are left to future studies. They include extending the model to cylindrical systems, adding a disk component to the spherical halo, taking into account the anisotropy of the velocity distribution, better accounting for the difference between total and dark matter mass and for the radial distribution of the mass change. Since simulations suggest that bursty star formation histories can lead both to dark matter core formation and to the expansion of the stellar distribution in UDGs, such galaxies may provide future observational tests for the model proposed here.

Acknowledgements

We thank A. El-Zant for his careful reviewing, F. Combes, A. Nusser, F. van den Bosch, O. Ginzburg, T. Nussbaum, Y. Birnboim and M. Kretschmer for stimulating discussions and physical insights. This work was partly supported by the grants France-Israel PICS, I-CORE Program of the PBC/ISF 1829/12, BSF 2014-273, NSF AST-1405962, GIF I-1341-303.7/2016, and DIP STE1869/2-1 GE625/17-1. NIHAO simulations were carried out at the Gauss Centre for Super-computing e.V. (www.gauss-centre.eu) at the GCS Supercomputer SuperMUCat Leibniz Supercomputing Centre (www.lrz.de) and on the High Performance Computing resources at New York University Abu Dhabi.

Appendix A Density and potential for a Dekel et al. (2017) mean density profile

The functional form for the average density proposed by Dekel et al. (2017) and presented in Eq. (12) yields the following analytical expressions for the mass, velocity and force per unit mass profiles:

[TABLE]

[TABLE]

and

[TABLE]

where $x=r/r_{c}$, $c=R_{\rm vir}/r_{c}$, $\mu=c^{a-3}(1+c^{1/b})^{b(g-a)}$, $M_{\rm vir}$ is the virial mass, $V_{\rm vir}^{2}=GM_{\rm vir}/R_{\rm vir}$ and $F_{\rm vir}=-GM_{\rm vir}/R_{\rm vir}^{2}$. We recall that $\overline{\rho_{c}}=c^{3}\mu\overline{\rho_{\rm vir}}$ with $\overline{\rho_{\rm vir}}$ the mean mass density within $R_{\rm vir}$.

The density profile is obtained by derivating the expression of the enclosed mass:

[TABLE]

This expression is in most cases very different from that of the average density profile (Eq. (12)), but does yield a similar expression in the case $g=3$, with an outer asymptotic slope equal to $3+1/b$.

The potential energy at radius $r$ is the work done by the gravitational force on a unit mass from $r$ to infinity. To derive the expression of the potential energy per unit mass from the force per unit mass $F$, we assume that it vanishes at infinity and that the halo density profile is truncated at the virial radius, in which case Eq. (4) yields

[TABLE]

from which Eq. (14) is derived.

The logarithmic slope of the density profile expressed in Eq. (35) is

[TABLE]

leading to Eq. (15) when $g=3$ and $b=2$, while the velocity given by Eq. (33) peaks at

[TABLE]

This defines a concentration parameter

[TABLE]

which coincides with $c$ when $a+g=4$ and leads to Eq. (17) when $g=3$ and $b=2$. The logarithmic slope $s$ at the resolution limit ($0.01R_{\rm vir}$ for the NIHAO simulations) and $c_{\rm max}$ can be used as effective inner slope and concentration when describing the density profile.

Appendix B Parametrizing the local kinetic energy with the Jeans equation

B.1 Spherical symmetry and isotropy

Complete spherical symmetry yields Eq. (18) expressing the local kinetic energy per unit mass as a function of the anisotropy parameter $\beta$ and the radial velocity dispersion $\sigma_{r}$ while we subsequently assume $\beta=0$ to write $K(r)=1.5\sigma_{r}^{2}$. As shown in Fig. 16, $\beta$ is indeed close to zero at most radii, although with a small offset that tend to increase with radius. In the left panel of Fig. 17, we test Eq. (18) for the outputs of one NIHAO galaxy, computing independently the local kinetic energy per unit mass $K=dK(<r)/dM$ on one side and $K_{\sigma}=(1.5-\beta)\sigma_{\rm r}^{2}$ on the other. The relative root mean square (RMS) difference between the two terms is below $5\%$, hence validating the assumption of complete spherical symmetry: we can neglect the global rotation of the halo when parametrizing the kinetic energy. In the right panel of Fig. 17, we show that further assuming $\beta=0$ in Eq. (18) still leads to an accurate estimate of the kinetic energy, with a 10% relative difference.

B.2 Jeans equilibrium

The steady-state equilibrium of a spherical collisionless system can be described by the spherical Jeans equation stemming from the Boltzmann equation (Binney & Tremaine, 2008, Eq. 4.215),

[TABLE]

where $\phi$ is the gravitational potential. Since $\sigma_{r}^{2}$ is related to the pressure support, this relation is analogous to hydrostatic equilibrium in hydrodynamics. When the anisotropy parameter $\beta$ is a constant, the Jeans equation corresponds to a differential equation in $\rho\sigma_{r}^{2}$ whose solution is

[TABLE]

assuming that it satisfies the boundary condition $\lim_{r\rightarrow+\infty}\sigma_{r}^{2}=0$. As the gravitational force is given by ${\rm d}\phi/{\rm d}r=GM(r)/r^{2}$, this yields

[TABLE]

leading to Eq. (19) when $\beta=0$. Introducing the logarithmic slopes

[TABLE]

the spherical Jeans equation can be rewritten

[TABLE]

where $V_{c}^{2}\equiv rd\phi/dr=GM/r$. In Fig. 18, we compare the two terms of Eq. (43) independently of any assumptions on the anisotropy parameter $\beta$. The figure shows that the Jeans equation describing the steady-state equilibrium of a spherical collisionless system is valid within a mean RMS difference of 13% for $\log(r/R_{\rm vir})$ between $-1.75$ ($0.02R_{\rm vir}$) and $-0.25$ ($0.56R_{\rm vir}$). As $\beta\ll 1.5$ (Fig. 16), the sum $\alpha+\gamma-2\beta$ is dominated by the first two terms and hence barely affected by $\beta$. The relative difference between the two terms of Eq. (43) in a given radius range, and in particular in the range $0.02-0.60~{}R_{\rm vir}$, can be used to characterize the deviation from Jeans equilibrium, as mentioned in Section 5.1 and Appendix C.

B.3 Analytical expressions of the local kinetic energy

Eq. (19) enables to obtain analytical expressions for the local kinetic energy per unit mass for a spherical, isotropic halo in Jeans equilibrium which is not globally rotating given the mass and density profiles.

B.3.1 Single-component halo

For a single-component halo whose density profile follows the Dekel et al. (2017) parametrization with $g=3$ and $b=2$ as in Eq. (13), the enclosed mass $M(r)$ is given by Eq. (32), which yields

[TABLE]

where $\mathcal{B}(a,b,x)=\int_{0}^{x}t^{a-1}(1-t)^{b-1}dt$ is the incomplete beta function and the outer brackets denote the difference of the enclosed function between 1 and $\chi$, i.e., $\left[f(\zeta)\right]_{\chi}^{1}\equiv f(1)-f(\chi)$. As mentioned in section 3.3, we extend the definition of the incomplete beta function appearing inside the outer brackets to negative parameters since the integral of Eq. (45) is well-defined as long as $\chi>0$. Following Eqs. 19 and A.9-11 of Zhao (1996), Eq. (20) (Eq. (46)) can also be expressed (see also Dekel et al., 2017, Eq. A.10) as the following sum:

[TABLE]

Alternatively, noticing that $\mathcal{B}(a,b,x)=a^{-1}x^{a}(1-x)^{b-1}+a^{-1}(b-1)\times\mathcal{B}(a+1,b-1,x)$, we can deduce that

[TABLE]

where $\Gamma$ denotes the usual gamma function. This enables to write the local kinetic energy as

[TABLE]

The different expressions for $K_{\rm DM}$ are equivalent.

B.3.2 Power-law multi-component halo

For a multi-component halo in which the ratio between the enclosed total mass $M_{\rm tot}$ and the enclosed collisionless mass $M$ is parametrized as a power-law according to Eq. (21), substituting $M_{\rm tot}$ to $M$ in Eq. (19) yields the local kinetic energy per unit mass

[TABLE]

As shown in Fig. 8 this expression reproduces very well the kinetic energy profiles of the NIHAO simulations. For such a multi-component halo, the gravitational potential from Eq. (36) further becomes

[TABLE]

which reverts to Eq. (14) when $n=0$ and $X_{M}=1$.

B.3.3 Additional central point mass

The contribution $\sigma_{m}$ of an additional mass $m$ to the velocity dispersion is expressed in Eq. (23); its contribution to the local kinetic energy is $K_{m}=1.5\sigma_{m}^{2}$ for an isotropic velocity field. For an additional point mass at the center, this yields

[TABLE]

This expression is approximately valid for an additional mass distribution dominated by the regions inside radius $r$, as notably assumed in Eq. (24).

B.3.4 Additional uniform sphere

For an additional (removed) mass $m$ uniformly spread between the halo center and a given radius $r_{m}$, the additional mass enclosed within a radius $r$ is $m(r/r_{m})^{3}$ for $r<r_{m}$ and $m$ for $r\geq r_{m}$. If $r\geq r_{m}$, its contribution to the kinetic energy corresponds to that of a point mass $m$ as in Eq. (59). However, if $r\leq r_{m}$, Eq. (23) yields

[TABLE]

where the second term corresponds to the contribution of a point mass $m$ evaluated at $r_{m}$. Hence if $r\leq r_{m}$,

[TABLE]

where $x_{m}=r_{m}/r_{c}$ and $\chi_{m}=\chi(x_{m})$.

B.3.5 Additional singular isothermal sphere

The mass distribution of a singular isothermal sphere of total mass $m$ and radius $r_{m}$ is $m(r/r_{m})$ for $r\leq r_{m}$ and $m$ for $r\geq r_{m}$. If $r\geq r_{m}$, its contribution to the kinetic energy corresponds to that of a point mass $m$ as in Eq. (59). If $r\leq r_{m}$, Eq. (23) yields

[TABLE]

B.3.6 Additional Dekel at al. (2017) sphere

If the additional mass is rather spread according to a Dekel et al. (2017) density profile $\rho_{m}(r)$ as in Eq. (13) with a total mass $m$ within $R_{\rm vir}$, the additional enclosed mass within $r$ is $\mu mx^{3}\overline{\rho(r)}/\overline{\rho_{c}}$ and the calculation of $K_{m}$ follow that of $K_{\rm DM}$, replacing $M_{{\rm vir}}$ by $m$ so that

[TABLE]

Appendix C Additional figures and table

To complement Section 4.4, Fig. 19 shows the radial profiles of the mass ratio $f=m/M$ between the successive outputs shown in Fig. 10. Together with the initial density profiles, these profiles are used to predict the evolution of the mass distribution according to the model. Table 2 indicates the success rate of the model galaxy per galaxy and for the whole sample, differenciating the overall success rate ($p_{\rm all}$), that excluding mergers ($p_{\rm no~{}mergers}$), that excluding mergers and with $|f|_{\rm RMS}\geq f_{\rm min}$ ($p_{>f_{\rm min}}$), as well as that excluding mergers but with $|f|_{\rm RMS}<f_{\rm min}$ ($p_{<f_{\rm min}}$).

Fig. 20 shows the correlations between different quantities relevent to disk and halo dynamics and how the success measures $\Delta s$, $\Delta c$, and $\delta$ as well as the probability $p$ associated to $\delta<\delta_{\rm max}$ depend on them. The quantities we consider are the time $t$, the virial mass $M_{\rm vir}$, the stellar mass $M_{\rm star}$, the gas mass $M_{\rm gas}$, the gas-to-virial mass ratio $\mu_{\rm gas/vir}=M_{\rm gas}/M_{\rm vir}$, the stellar half-light radius $R_{1/2}$, the merger indicator $f_{\rm merger}$, and the deviation from Jeans equilibrium $\Delta_{\rm Jeans}$ defined from the relative difference shown in the lower left panel of Fig. 18. We note that a correlation coefficient $r\gtrsim 0.15$ between the success measures and any given quantity is sufficient to induce a noticeable trend of the associated probability $p$. The figure notably shows that the main quantities on which the model success depends are $f_{\rm merger}$, $\Delta_{\rm Jeans}$, $t$, and $\mu_{\rm gas/vir}$. Fig 21 further shows the correlations between the three success measures $\delta$, $\Delta s$ and $\Delta c$. As expected from the fact that both $\delta$ and $\Delta s$ assess the accuracy of the model to describe changes in the inner part of the average density profile while $\Delta c$ relates to the dark matter distribution at larger scales, $\delta$ is more closely related to $\Delta s$ (Pearson correlation coefficient $r=0.71$) than to $\Delta c$ ($r=0.43$). In Fig. 22, we show the counterparts of Fig. 12 for $\Delta s$ and $\Delta c$ the errors in inner logarithmic slope and concentration. The plots display similar trends as in Fig. 12. The secondary trends as a function of the deviation from Jeans equilibrium $\Delta_{\rm Jeans}$, time $t$ and gas-to-virial mass ratio $\mu_{\rm gas/vir}$ are visible.