Baryons Shaping Dark Matter Haloes

We analyse galaxies selected from the 100 Mpc sized box reference run of the eagle Project, a suite of hydrodynamical simulation that follows the structure formation in cosmological representative volume. All of them are consistent with the current favoured $\Lambda$-CDM cosmology (Crain et al., 2015; Schaye et al., 2015). These simulations include: radiative heating and cooling (Wiersma et al., 2009), stochastic star formation (Schaye & Dalla Vecchia, 2008), stochastic stellar feedback (Dalla Vecchia & Schaye, 2012) and AGN feedback (Rosas-Guevara et al., 2015). The AGN feedback is particularly important for the evolution of SF activity in massive early-type galaxies (ETGs). An Initial Mass Function (IMF) of Chabrier (2003) is used. A more detailed description of the code and the simulation can be seen in Crain et al. (2015) and Schaye et al. (2015).

The adopted cosmological parameters are $\Omega_{m}=0.307$, $\Omega_{\Lambda}=0.693$, $\Omega_{b}=0.04825$, $\mathrm{H_{0}=100\>\textit{h}\,km\,s^{-1}\,Mpc^{-1}}$, with $h=0.6777$ (Planck Collaboration et al., 2014). The 100 Mpc sized box reference simulation, so called L100N1504, is represented by $1504^{3}$ dark matter particles and the same initial number of gas particles, with an initial mass of $\mathrm{9.70\times 10^{6}M_{\odot}}$ and $\mathrm{1.81\times 10^{6}M_{\odot}}$, respectively. The maximum gravitational softening of $0.7\ \mathrm{kpc}$ is adopted.

The energy injected into the gas corresponds to $10^{51}$ erg per supernovae event times a dimensionless factor, $f_{E}$, that depends on the local gas metallicity and density. For haloes with more than 1500 DM particles (approximately $\mathrm{M}>10^{10}h^{-1}\rm M_{\odot}$), black hole sink particles are placed at the center of the halos. Feedback due to AGN activity is implemented in a similar way to the feedback from star formation. For haloes, at $\rm z=0$, with masses higher than $\mathrm{>10^{13}M_{\odot}}$ the stellar feedback becomes less effective but AGN feedback can still expel baryons (Schaller et al., 2015a).

The halo catalogue ¹¹1We use the publicly available data by McAlpine et al. (2016) http://icc.dur.ac.uk/Eagle/database.php was constructed by using a Friends-of-Friends algorithm, while the substructures were identified using the SUBFIND algorithm (Springel et al., 2001). For the eagle simulations the match of the haloes in the Hydro and DMo simulations was done following Schaller et al. (2015b).

We use the cosmological simulation S230D from the Fenix Project suite. Several properties of galaxies in this simulation have been thoroughly studied such as the morphological properties (Pedrosa et al., 2014), the size-mass relation and the angular momentum evolution (Pedrosa & Tissera, 2015), the stellar and gaseous metallicity gradients of the disc components (Tissera et al., 2016a; Tissera et al., 2016b; Tissera et al., 2017), the chemical abundance of the circumgalactic medium (Machado et al., 2018) and the fundamental properties of the elliptical galaxies, such as Faber-Jackson relations and the Fundamental Plane (Rosito et al., 2018).

This simulation is consistent with a $\Lambda$-CDM universe with $\mathrm{\Omega}_{\mathrm{m}}=0.3$, $\mathrm{\Omega}_{\mathrm{\Lambda}}=0.7$ and $\mathrm{\Omega}_{\mathrm{b}}=0.04$, $\mathrm{H}_{0}=100\ h\mathrm{\ km\ s^{-1}\ Mpc^{-1}}$ with $h=0.7$, and a normalization of the power spectrum of $\mathrm{\sigma}_{8}=0.9$. The simulated box is $14\mathrm{\ Mpc}$ a side. The initial conditions have $2\times 230^{3}$ total particles and a mass resolution of $4.3\ \times{}10^{6}\mathrm{M}_{\odot}$ and $6.4\times{}10^{5}\mathrm{M}_{\odot}$ for the DM particle and initial gas particle, respectively. The maximum gravitational softening is $0.35\ \mathrm{kpc}$. The initial condition has been chosen to correspond to a typical region of the Universe with no massive group present (the largest haloes have virial masses smaller than $\sim\mathrm{10^{13}M_{\odot}}$). We acknowledge the fact that the initial conditions (ICs) represent a small volume of the Universe. Nevertheless, De Rossi et al. (2013) showed that the growth of the simulated haloes are well-described in these simulations by confronting them with those from the Millenium Simulation (Fakhouri et al., 2010).

These simulations were run using gadget-3, an update version of gadget-2 (Springel & Hernquist, 2003; Springel, 2005), optimised for massive parallel simulations of highly inhomogeneous systems. It includes treatments for metal-dependent radiative cooling, stochastic star formation, chemical and energetic supernovae (SN) feedback (Scannapieco et al., 2005, 2006). This feedback model is able to reproduce galactic mass-loaded winds without introducing any mass-dependent scaling parameter. It also includes a multiphase model for the ISM that allows the coexistence of the hot, diffuse phase and the cold, dense gas phase (Scannapieco et al., 2006, 2008), where star formation takes place. Part of the stars ends their lives as Type II and Type Ia Supernovae, injecting energy and chemical elements into the ISM. Each SN event releases $7\ \times\ 10^{50}$ erg, which are distributed equally between the cold and hot phases surrounding the stellar progenitor. The adopted code uses the chemical evolution model developed by Mosconi et al. (2001) and adapted to gadget-3 by Scannapieco et al. (2005). This model considers the enrichment by SNII and SNIa adopting the yield prescriptions of Woosley & Weaver (1995) and Iwamoto et al. (1999), respectively and the Initial Mass Function (IMF) of Salpeter (1955). The synthesised chemical elements are distributed between the cold and hot phase ($\mathrm{80\%}$ and $\mathrm{20\%}$, respectively). The lifetimes for SNIa are randomly selected within the range [0.1, 1] Gyr. Albeit simple, this model reproduces well the mean chemical patterns obtained by adopting the single degenerated model (Jiménez et al., 2015). This combination of star formation and feedback parameters produced systems which can reproduce the angular momentum and the mass-size relation of the discs and bulges of galaxies in agreement with the observational trends (Pedrosa & Tissera, 2015).

The halo catalog was constructed by using a Friends-of-Friends algorithm, while the substructures were identified using the SUBFIND algorithm (Springel et al., 2001). As a result 317 galaxies were identified, resolved with more than $2000$ DM particles. The virial masses of the haloes ($\rm M_{200}$) were defined as the mass within a sphere of radius ($\rm R_{200}$) containing$\sim 200$ times the cosmic critical matter density at the corresponding redshift. The hydrodynamic simulation has a DMo run counterpart. We refer to them as "Hydro" and "DMo", respectively. The Hydro and DMo simulations start from identical ICs. The shrinking sphere method proposed by Power et al. (2003) is applied to find the coordinates of the centre of mass of the haloes in each simulation. Haloes are matched between the Hydro and DMo runs by one-by-one as performed for the EAGLE set.

The modification of the potential well due to the gathering of baryons in the inner regions to form a galaxy, induces changes in the mass distribution of the DM haloes. Several studies have addressed this issue adopting different radius to evaluate the effects (Butsky et al., 2016; Zavala et al., 2016; Thob et al., 2019; Chua et al., 2019). In order to estimate the radius that maximizes the signal, we inspect three possible selections defined by the radii that enclose 5, 10 and 20 percent of M₂₀₀ (hereafter measuring radius).

This redistribution can be appreciated in Fig. 1 where we show the ratio between the DM masses measured at the three defined measuring radii in the Hydro and DMo runs, $\rm f_{{}_{200}}^{{}^{i\%}}$= M^Hydro/M${}^{\rm DMo}(\mathrm{i\%R_{200}})$, as a function of the stellar-to-halo mass ratio $\rm log(M_{\rm star}/M_{\rm 200})$. The label $i$ denotes the measuring radii defined to enclose 5, 10, and 20 percent of M₂₀₀.

The results in Fig. 1 show that, for the inner regions, haloes in eagle can contract or expand (i.e. higher or lower DM masses in the Hydro run compared to the DM only). This type of behaviour has also been detected by Dutton et al. (2016) using 100 hydrodynamical cosmological zoom-in simulations of the NIHAO Project. While the contraction is the response to the increase of the potential well due to the accumulation of baryons, the expansion is claimed to be the result of stellar feedback.

An important fraction of the eagle haloes expand when comparing Hydro and DMo runs, while Fenix haloes contract in almost all cases and radii. The expansion in eagle haloes in the central regions has been thoroughly analysed by Schaller et al. (2015a) and is the result of a more effective stellar feedback and a greater loss of baryons (Sawala et al., 2013). Also it can be noticed in Fig. 1 that for higher masses, the contraction is more important as a result of a less effective feedback action. As can be seen, in the case of Fenix simulation, the subsample of haloes studied always contracts when baryons are present within the mass range obtained with this experiment.

We also estimate the galaxy compactness factor, $\rm e_{\rm 1/2}\left[\%\right]$ (Dutton et al., 2016), defined as the ratio between the galactic half-mass radius³³3The half-mass radius, $\rm r_{\rm 1/2}$, is defined as the one that enclosed 50 percent of the baryons., $\rm r_{\rm 1/2}$, and the virial radius, i.e. $\mathrm{(100*r_{1/2})/R_{200}}$. For convenience it is expressed as a percentage. The colours in Fig. 1 map $\rm e_{\rm 1/2}\left[\%\right]$. In the most inner regions, within 5 percent of $\mathrm{R_{200}}$, haloes hosting more compact galaxies tend to be more concentrated as can be seen in the top panel of Fig. 1. When these ratios are measured at $\mathrm{0.10R_{200}}$ ($\rm f_{{}_{200}}^{{}^{10\%}}$) and $\mathrm{0.20R_{200}}$ ($\rm f_{{}_{200}}^{{}^{20\%}}$), this trend gradually disappears. This is expected since the larger effects on baryons in the density distribution will be in the very central region where they are located. The concentration of the most inner regions of the haloes depends not on the global baryonic mass gathered in the very centre of haloes but also on their degree of compactness (e.g. Pedrosa et al., 2010; Tissera et al., 2010; Dutton et al., 2016).

The change of the DM concentration in the inner regions can be also quantified through the contraction/expansion parameter $\mathrm{\Delta_{\rm v/2}}$, defined as the ratio of the mean DM density to the critical density $\rho_{crit}$ within the radius at which the DM reaches half of its maximum circular rotational velocity, $\mathrm{R_{\rm V/2}}$.

This parameter was previously implemented in Pedrosa et al. (2010) in the study of the halo contraction evolution with the redshift and originally introduced by Alam et al. (2002). By using the radius "half way up", $\mathrm{R_{\rm V/2}}$, the rising part of the rotation curve, Alam et al. (2002) focus on the region where conflicts between predicted and observed DM halo densities are more severe (and where the observations are still typically robust against resolution uncertainties).

In Fig. 2 we compare the concentration parameters from the Hydro and DMo runs as a function of $\mathrm{R_{\rm V/2}}$ and the stellar-to-virial-mass ratio. Haloes with higher contraction (i.e., $\Delta_{\rm v/2}^{\rm Hydro}/\Delta_{\rm v/2}^{\rm DMo}\gg 1$) are populated with more compact galaxies (i.e., lower $\rm e_{\rm 1/2}\left[\%\right]$). The compactness of galaxies appears as a key parameter to trace the degree of concentration of the dark matter haloes. When haloes expand (i.e., $\Delta_{v/2}^{\rm Hydro}/\Delta_{v/2}^{\rm DMo}\ll 1$), they tend to host more extended galaxies. In the top panel it can be seen that compact galaxies are hosted by haloes with shorter $\mathrm{R_{\rm V/2}}$, indicating that the DM maximum rotational velocity $V_{\rm max}$ moves to the inner part of the halo with respect to their DMo counterparts. The lower panel shows that both Fenix and eagle subsamples present a slight correlation (Spearman rank coefficient $\rho=0.15$) with $\rm log(M_{\rm star}/M_{\rm 200})$, in the sense that more massive galaxies with respect to halo mass tend to gather also larger DM masses within the very inner regions when baryons are present.

These results confirm the discussed above findings, showing that less concentrated galaxies are located in haloes with expanded central regions with respect to their DMo counterparts. This is clear for the eagle halo sample, which in turn spans on a wide range of stellar-to-dark-matter ratio.

Our main interest is to dig into the interdependence of halo shape and galaxy morphology and try to quantitatively correlate them.

We describe the shapes using the semi-axes of the triaxial ellipsoids, $\rm a>b>c$, where $\rm a$, $\rm b$ and $\rm c$ are the major, intermediate and minor axis respectively of the $shape$ $tensor$ $\rm S_{ij}$ (e.g. Bailin & Steinmetz, 2005; Zemp et al., 2011). An iterative method is used, starting with particles selected in a spherical shell (i.e. $\rm q=s=1$ Dubinski & Carlberg, 1991; Curir et al., 1993).

In order to obtain these ratios $\rm q\equiv b/a$ and $\rm s\equiv c/a$, we diagonalize the reduce inertia tensor to compute the eingervectors and eigenvalues, as described in Tissera & Dominguez-Tenreiro (1998). Traditionally the $\rm s$ shape parameter has been used as a measure of halo sphericity (e.g. Allgood, 2005; Vera-Ciro et al., 2014; Chua et al., 2019).

We adopt the triaxiality parameter, defined as $\rm T\equiv(1-q^{2})/(1-s^{2})$, which quantifies the degree of prolatness or oblatness: $\rm T=1$ describes a completely prolate halo ($\rm a>b\approx c$) while $\rm T=0$ describes a completely oblate halo ($\rm a\approx\rm b>c$). Haloes with $\rm T>0.67$ are considered prolate and haloes with $\rm T<0.33$ oblates, while those with $\rm 0.33<T<0.67$ are considered triaxials (Allgood, 2005; Artale et al., 2019).

In the next sections we inspect the shape parameters of the haloes in relation with the measuring radius of the inner region, the D/T fraction and the stellar mass of the hosted galaxy.

The velocity anisotropy of individual haloes presents a variety of behaviours depending on their particular formation history (see e.g., Tissera et al., 2010). In this section we analyse the change in velocity anisotropy of the DM haloes as a function of radius and the morphology of the hosted galaxy. The velocity structure of the DM particles are closely related with the resulting shapes of the haloes. Additionally, it has important implications for the predicted scattering rates in direct detection experiments (Kuhlen et al., 2010).

We inspect the velocity anisotropy parameter defined as:

where $\sigma_{r}$ and $\sigma_{t}$ are the radial and tangential velocity dispersions averaged over concentric spherical shells. This anisotropy parameter provides some indication of the velocity distribution of the haloes. Isotropic distributions have $\beta\approx 0$ whereas those that more radially biased have $\beta>0$.

In Fig. 9 we show the median distributions of $\beta(r)$ as a function of the radius normalized by the virial one, $\rm r/R_{\rm 200}$, for galaxies grouped according to their morphology: $\rm D/T<0.30$, $0.30\leq\rm D/T<0.50$, $0.50\leq\rm D/T<0.70$, $0.70\leq\rm D/T<0.90$, and $0.90\leq\rm D/T$, for the Fenix (left panel) and eagle (right panel) haloes. Both eagle and Fenix simulations show a noticeable trend for radius smaller than approximately $\mathrm{0.20R_{200}}$: haloes hosting almost bulgeless galaxies present higher values of $\beta$. The velocity structure within the central regions of DM haloes hosting the highest D/T fraction depart from isotropy. For radius larger than $\rm\sim 0.20R_{200}$ this trend disappear.

When the DMo is inspected (lower panels), $\beta^{DMo}(r)$ show a slight trend in the same sense albeit weaker. As we have found for the shapes, there is an intrinsic behaviour indicating that haloes that host important stellar disc structures tend to have a distinct, slightly less isotropic, DM velocity pattern. To better quantify this behaviour, we estimate $(\beta-\bar{\beta})/\bar{\beta}$ as a function of $\rm r/R_{\rm 200}$ for both the Hydro and DMo runs for galaxies with different morphologies.

Fig.  10 show the relative change of $\beta(r)$ with respect to the median relations (i.e $(\beta-\bar{\beta})/\bar{\beta}$) for both Fenix (panels a,b) and eagle (panels c, d) haloes as a function of $\rm r/R_{\rm 200}$. $\beta^{DMo}(r)$, panels b and d, show a slight trend in the same sense albeit weaker. As we have found for the shapes, there is an intrinsic behaviour indicating that haloes that host more important stellar disc structures tend to have a distinct, slightly less isotropic, DM velocity pattern. The error regions were computed by using a bootstrap method with a resampling of 200 for each D/T interval. As can be seen, the larger velocity anisotropies are found in haloes with almost bulgeness galaxies within $\sim\rm 0.20R_{200}$. The correspondent DMo counterpart haloes show similar estimations but using the $\beta^{\rm DMo}$.

A greater velocity anisotropy could be related with modification in the orbit structure of DM velocity. We have also looked for DM discs in the central regions of the selected Fenix haloes without any detection. Schaller et al. (2016) concludes that the presence of a dark disc is unlikely in the eagle haloes. The DM configuration could be related to more tube-type orbits, as found by Zhu et al. (2017). So this change in orbit configuration towards higher central velocity anisotropies together with intrinsic DMo haloes shapes could foster preferentially the formation of stellar discs.

The presence of a similar excess in the same radial range is consistent with the claim of an intrinsic characteristics of the haloes that might contribute to the formation of the disc-dominated galaxies.