Anisotropic satellite accretion onto the Local Group with HESTIA

The HESTIA simulation suite (Libeskind et al., 2020) is a set of low-resolution and high-resolution cosmological magnetohydrodynamical simulations of the Local Group (LG).

The Cosmicflows-2 (CF2) peculiar velocities (Tully et al., 2013) are used to construct the initial conditions (within the framework of constrained simulations). The CF2 catalog is grouped (Sorce & Tempel, 2018) and bias-minimized (Sorce, 2015). The simulations are built using the AREPO code to solve the ideal magnetohydrodynamics (MHD) equations (Springel, 2010; Pakmor et al., 2016; Weinberger et al., 2020), and the AURIGA galaxy formation model (Grand et al., 2017), which is based on the Illustris model (Vogelsberger et al., 2013) and implements various physical processes that are the most relevant for the formation and evolution of galaxies. Vogelsberger et al. (2014a, b) uses the same galaxy formation model. The reader can also refer to Vogelsberger et al. (2020) for a review on cosmological simulations.

The entire suite of the HESTIA simulations assumes a $\Lambda$CDM cosmology with the following Planck Collaboration et al. (2014) values: $\sigma_{8}=0.83$, $H_{0}=100h$ km s^-1 Mpc^-1 where $h=0.677$, $\Omega_{\Lambda}=0.682$, $\Omega_{m}=0.270$ and $\Omega_{b}=0.048$.

The publicly available AHF halo finder (Knollmann & Knebe, 2009) identifies haloes, subhaloes, and their properties, at each redshift by detecting gravitationally bound particles. Structures containing less than 20 bound particles are not considered. Halo histories are derived by the MERGER TREE algorithm, a package included in the AHF software.

In this paper, we consider three Local Groups simulated at high resolution, namely 09_18, 17_11, 37_11 (the numbers are based on the random seed used). These three simulations each consists of a region of two 3.7 Mpc spherical volumes centered on the two primary haloes at $z=0$, filled with $8192^{3}$ effective particles, achieving a mass and spatial resolution of $m_{\mathrm{dm}}=1.5\times 10^{5}$ M_⊙, $m_{\mathrm{gas}}=2.2\times 10^{4}$ M_⊙ and $\epsilon=220$ pc. A larger number of lower resolutions simulations are considered in the appendix.

The reader may refer to Libeskind et al. (2020) for more details about the HESTIA suite and the algorithms mentioned above.

To evaluate the cosmic web, we consider the velocity shear tensor, derived from the distribution of the particles in the simulation box. First, a clouds-in-cell (CIC) technique is applied to the simulation to compute the density and velocity fields in a $256^{3}$ grid of side length $100$ Mpc/$h$, hence resulting in a spatial resolution of $0.39$ Mpc/$h$. The grid is then smoothed with a Gaussian kernel, and various smoothing lengths $r_{s}=1,2,5$ Mpc/$h$ are considered.

Finally, the shear in the velocity field is derived from the velocity field as (Hoffman et al., 2012):

where $H_{0}$ is the Hubble constant, and the partial derivatives of the velocity $V$ are derived along the directions $\alpha$ and $\beta$, representing the orthogonal Supergalactic Cartesian axes $x$, $y$ and $z$. In each grid cell, we obtain the eigenvectors $\vec{e_{1}}$, $\vec{e_{2}}$ and $\vec{e_{3}}$, and the corresponding eigenvalues ordered such that $\lambda_{1}>\lambda_{2}>\lambda_{3}$.

We also consider in this work the tidal shear tensor (Zel’Dovich, 1970; Hahn et al., 2007a), defined as the Hessian of the gravitational potential $\phi$:

where the gravitational potential is rescaled by a factor $4\pi G\bar{\rho}$ (where G is the gravitational constant and $\bar{\rho}$ is the mean density of the Universe), and obeys the Poisson equation $\nabla^{2}\phi=\delta$, where $\delta$ is the matter overdensity field.

Throughout this paper and for both tensors, the eigenvectors are computed following the methodology described in Wang et al. (2020a), which differs from the usual Nearest Grid Point method, as in this case the Hessian matrix is computed and solved for each halo in the simulation, instead of doing it once for the entire grid. In practice the tensor fields are first evaluated on a gird. Instead of then begin evaluated on the grid, the method of Wang et al 2020 computes a single tensor associated to each halo before diagonalisation. The tensor field associated with each halo is constructed based on the 7 adjoining cells in a CIC inspired way, namely weighted by distance to the halo in consideration. Each halo thus has a unique velocity or tidal tensor associated with it which is then diagonalized and whose eigenvectors are then used.

AHF uses the parameter $R_{200}$¹¹1The parameter $R_{200}$ is defined as the radius at which the mean enclosed matter density is $\rho(<R_{200})=200\times\rho_{\mathrm{crit}}$, where $\rho_{\mathrm{crit}}$ corresponds to the critical density for closure. to determine masses and radii of objects. In this case, the actual virial radius $R_{\mathrm{vir}}$ of halos sits between $R_{200}$ and $2\times R_{200}$. Hence, except when specifically mentioned, we consider the moment of infall at $2\times R_{200}$, as $R_{\mathrm{vir}}$ is larger than $R_{200}$ and we do not expect much signal at that threshold (for reasons that will be explained and examined below).

Satellite galaxies are identified as being within $R_{200}$ of their main host at $z=0$. These are then tracked backwards in time by following the main branch of the merger tree, which tracks the location of the most massive progenitor (sub)halo. A moment of accretion is identified as when the satellite first crossed $2R_{200}$. This is termed $z_{inf}$. For clarity we repeat that satellites are identified as those haloes that reside within $R_{200}$ at $z=0$ but their infall time is when they first crossed $2R_{200}$. This is because we are interested in examining the origin of the classical satellite population (not all LG dwarfs) and if their accretion happened at earlier times.

This is done by examining two adjacent snapshots. A subhalo is deemed to have been accreted if in the later snapshot it is located within a sphere of radius $2\times R_{200}$ centered on a main progenitor of one of the two LG haloes (MW or M31), where $R_{200}$ is defined by the corresponding main progenitor. The progenitor of the subhalo has to be located outside of the so defined sphere at the earlier snapshot. By going backwards from $z=0$ through all available snapshots, two at a time, we can follow the trajectories of satellites and identify the place and time of accretion. As mentioned earlier, only the subhaloes that survive up to $z=0$ (i.e the ones that are present within $R_{200}$ in the $z=0$ snapshot), are considered. Moreover, we only look at the first time the subhaloes cross the $2\times R_{200}$ threshold.

We refer to the time of infall as the last snapshot, corresponding to the redshift $z_{\mathrm{inf}}$, before a given subhalo crosses for the first $2\times R_{200}$. We can then define the direction of infall $\vec{r}_{\mathrm{inf}}$ as the coordinates of the subhaloes at infall, re-centered on the respective main host (MW or M31). We also note the total mass of a subhalo at infall $M_{200}^{\mathrm{inf}}$.

The distribution of the redshifts $z_{\mathrm{inf}}$, i.e when the accreted subhaloes cross $2\times R_{200}$ of the MW (red lines) or M31 (blue lines), is shown in Figure 1. The three LGs considered are represented by thin solid, dashed or dotted lines (09_18, 17_11 and 37_11 respectively). The bold lines show the total redshift distribution of all satellites taken all together, for each main halo. Two peaks can be observed: late infall and early infall. The majority of satellites are accreted recently, after $z_{\mathrm{inf}}=0.7$. Other subhaloes get accreted much earlier $z_{\mathrm{inf}}>0.7$ and survive up to $z=0$. Also, we observe that the three individual LGs show a similar trend. This demonstrates the quality and power of the constrained simulations, as such a behavior would not be expected from randomly selected LGs. In fact it has been previously established (Sorce et al., 2016; Libeskind et al., 2010; Carlesi et al., 2020) that one of the main traits of constrained simulations is that they limit the assembly histories of specific objects. We also note that both peaks of the distributions (total and for each LG) of the two main hosts overlap. This indicates that it is a feature of the entire LG, and that the haloes of the MW and M31 have likely had similar accretion histories.

Figure 2 shows the distribution of the cosine of the angle between the infall direction $\vec{r}_{\mathrm{inf}}$, vector at infall of the accreted satellite galaxy and the velocity vector at infall $\vec{v}_{\mathrm{inf}}$. Both of these vectors are centered on the parent halo’s frame. Thin solid, dashed and dotted lines correspond to a single realization of the Local Group. The bold solid lines gives the total distribution, combining all three simulations together. We can observe that the distribution peaks around $\cos{(\vec{r}_{\mathrm{inf}}.\vec{v}_{\mathrm{inf}})}=-1$, i.e, the two vectors are aligned with opposite direction, meaning that the satellites are accreted radially, along their velocity.

We can approximate the pre-infall distance travelled by satellites by considering the co-moving coordinate distance between their position at the first snapshot recorded (i.e, at birth) and their position at $z_{\mathrm{inf}}$, both re-centered on the respective host halo. This approximation gives the minumum distance travelled by satellites before accretion, as we assume that their path from birth to infall is a straight line, and do not consider their actual trajectories. The distribution of the distance travelled within the host rest frame, combining all 3 simulations, is shown in Figure 3. We observe that subhaloes have travelled a distance up to 4 Mpc before getting accreted. More importantly, we notice a peak at 1 Mpc, showing that on average most $z=0$ subhaloes traveled a distance of roughly 1 Mpc before crossing the $2\times R_{200}$ threshold of the MW or M31. That said, the median of the distribution is greater than this, implying that the LG accretes from, roughly, a sphere of up to 4 Mpc around the two main hosts.

In Figure 4, the minimum co-moving distance travelled described above is plotted as a function of the redshift at time of infall $z_{\mathrm{inf}}$. Dark matter only haloes are represented by black dots, while haloes containing baryons are highlighted by red triangles. Firstly, one can easily notice the two peaks in the redshift distribution as seen in Figure 1, before and after $z=0.7$. Moreover, the decreasing trend shows that the satellites accreted recently (recent $z_{\mathrm{inf}}$) are the one that travelled the most, i.e with the largest minimum travelled distance. This is because the haloes that were recently accreted have had the most amount of time to travel across the Universe. Those accreted at early times had less time at their disposal and thus were not able to cross large distances.

In this section, we analyze the alignment of the infall direction $\vec{r}_{\mathrm{inf}}$ with the three eigenvectors of the velocity shear and tidal tensors, by considering the distribution of the cosine of the angle between $\vec{r}_{\mathrm{inf}}$ and $\vec{e_{1}}$, $\vec{e_{2}}$ and $\vec{e_{3}}$. As eigenvectors are non-directional, $0\leq\cos{(\vec{r}_{\mathrm{inf}}.\vec{e}_{1,2,3})}\leq 1$. A cosine of 1 means that the two directions are parallel while 0 implies they are perpendicular. To quantify the strength of the distribution of $\cos{(\vec{r}_{\mathrm{inf}}.\vec{e}_{1,2,3})}$, we calculate its statistical significance as the difference between the actual distribution and the median of 10,000 uniform distributions of the same size, averaged over all bins and calculated in units of the Poisson error. The uniform distributions are derived from 10,000 random samples of points uniformly distributed on a sphere. High values of the statistical significance indicate a strong deviation from a uniform distribution, hence a strong signal, while low values show a statistically weak signal, close to or consistent with a uniform distribution.

Figure 5 shows the distribution of $\cos{(\vec{r}_{\mathrm{inf}}.\vec{e}_{1,2,3})}$ at $R_{200}$ (top row) and $2\times R_{200}$ (bottom row). Note this is the only plot in the paper where we have examined the accretion at $R_{200}$. The three eigenvectors $\vec{e_{1}}$, $\vec{e_{2}}$ and $\vec{e_{3}}$ are represented respectively by blue, orange and green lines. Both tidal and shear tensors have been checked, however, the obtained distributions being similar, only the tidal tensor eigenvectors are shown for clarity. Similarly, as all smoothing lengths applied to the tensors show similar trends, only the distributions obtained with a smoothing length of $r_{s}=1$ Mpc are displayed. From left to right columns depict all satellites (combining all three simulations and both host haloes), MW satellites only, and M31 satellites only, respectively. The statistical significance of each distribution is noted on the top of each panel, with the same color code as their respective eigenvector: blue for $\vec{e_{1}}$, orange for $\vec{e_{2}}$, and green for $\vec{e_{3}}$. The grey shaded regions depict the Poisson error, i.e the standard deviation $\sigma$ of 10,000 uniform distributions at each bin. The dark grey shows the $\pm 1\sigma$ interval, while the lighter grey shows the $\pm 2\sigma$ interval.

The distributions in top row panels of Figure 5 being mostly located within the grey regions, and the values of the statistical significance being mostly less or close to unity, indicate there is no significant alignment of the infall direction with any eigenvector of the tidal tensor. We remind the reader that the top row corresponds to the distribution of accreted subhaloes at $R_{200}$. However, when looking at the distributions at $2\times R_{200}$ in the bottom row panels of Figure 5, we notice that the infall direction $\vec{r}_{\mathrm{inf}}$ is strongly aligned with the axis of slowest collapse $\vec{e_{3}}$. We also notice that $\vec{r}_{\mathrm{inf}}$ is orthogonal to $\vec{e_{1}}$, while no significant alignment is observed with $\vec{e_{2}}$. At $R_{200}$, i.e in the non-linear regime, the signal is extremely weak due to non-linear dynamics. Beyond $R_{200}$, we transition from the virial regime to the quasi-linear (QL) regime, which explains why the signal is much stronger at $2\times R_{200}$.

Besides, when comparing the middle (MW satellites) and right (M31 satellites) columns of Figure 5, corresponding to satellites accreted by each host separately, no significant difference can be observed in the top row, i.e at $R_{200}$, for the same reasons described above. At $2\times R_{200}$, the distributions corresponding to both hosts look similar, however the statistical significance corresponding to the satellites accreted by M31 is $\sim 1.7$ times higher than the one corresponding to the MW satellites. This shows that, although the alignment with $\vec{e_{3}}$ is observed for both hosts, it is stronger for M31, and weaker for the MW satellites. Such a difference may be explained by the fact that M31 is more massive than the MW ²²2The more massive halos are called M31 by definition in HESTIA simulations..

In the next figures, we look at the distribution of $\cos{(\vec{r}_{\mathrm{inf}}.\vec{e}_{1,2,3})}$ at $2\times R_{200}$ split by mass and redshift, in order to check for any dependence.

First, Figure 6 shows the distribution of $\cos{(\vec{r}_{\mathrm{inf}}.\vec{e}_{1,2,3})}$ at $2\times R_{200}$, where $e_{1,2,3}$ are represented by blue, orange and green lines, respectively. For clarity purposes, only the results obtained from the tidal tensor, smoothed by 1 Mpc, are shown, however both shear and tidal tensors show similar trends. Similarly, the results obtained with other smoothing lengths are not displayed as the distributions are very similar. The distributions are split by total mass of satellites at accretion, namely $M_{200}^{\mathrm{inf}}$, using the following bins: $M_{200}^{\mathrm{inf}}<10^{7}$ M_⊙ (first column), $10^{7}\leq M_{200}^{\mathrm{inf}}<10^{8}$ M_⊙ (second column) and $M_{200}^{\mathrm{inf}}>10^{8}$ M_⊙ (third column). The statistical significance and the standard deviation of a uniform distribution are shown on each panel. The reader may refer to the previous figures for a detailed description. All three mass bins show alignment of the accretion however we observe a mild mass dependence wherein the smallest infall masses ($M_{200}^{\mathrm{inf}}>10^{7}M_{\odot}$) exhibit a weaker signal than the more massive haloes. The $10^{7}\leq M_{200}^{\mathrm{inf}}<10^{8}$ mass bin is stronger.

In Figure 7, the distributions of $\cos{(\vec{r}_{\mathrm{inf}}.\vec{e}_{1,2,3})}$ at $2\times R_{200}$ are split into two bins of redshift at infall $z_{\mathrm{inf}}$: late infall and early infall. Subhaloes accreted recently at $z_{\mathrm{inf}}<0.7$ are shown in the left column while the ones accreted earlier at $z_{\mathrm{inf}}\geq 0.7$ are shown in the right column. The top row shows the distribution with the tidal tensor’s eigenframe while the bottom rows show this distribution for the shear tensor. The reader may refer to the texts of the previous figures for a detailed description of the color code and quantities displayed on the panels. We notice that the alignment with $\vec{e_{3}}$ is dominated by early infall and is more isotropic at later times.

Finally, Figure 8 depicts the location of the MW (left) and M31 (right) satellites entry points as they cross the $2\times R_{200}$ of their respective main halo. The entry points are shown as a “heat map” and plotted on an Aitoff projection of the $2R_{200}$ sphere, oriented in the eigenframe of the tidal tensor, smoothed at 1 Mpc. As the eigenvectors are non-directional, only a single octant of the $2R_{200}$ sphere is shown, instead of an entire full sky projection. The north pole of the projection corresponds to the direction of $\vec{e_{1}}$, the right point on the horizontal axis at $+90\deg$ corresponds to the direction of $\vec{e_{2}}$, and the midpoint corresponds to the direction of $\vec{e_{3}}$. Areas within $15\deg$ of the eigenvectors $\vec{e_{1}}$, $\vec{e_{2}}$ and $\vec{e_{3}}$ are defined respectively by blue, orange and green circles. The yellow pixels correspond to a high density of satellite entry points, while dark blue pixels correspond to a low density of entry points. The location of the simulated Virgo in all three realizations is shown as white scattered points (circle: 09_18, triangle: 17_11, square: 37_11). In both cases, we observe a high density of points around $\vec{e_{3}}$, which is in accordance with the results shown above. Besides, the high density around the axis of slowest collapse looks more spherical in the case of M31 than for the MW, confirming the different trends between the two main hosts observed in the middle and right panels of Figure 5.

We remind the reader that Figure 8 represents a single octant of the $2R_{200}$ sphere as eigenvectors are non-directional. However one can give directions to the eigenvectors in order to get a full sky representation of this sphere by choosing positive direction. We establish an oriented eigenframe with respect to Virgo, by choosing as positive, the direction of each eigenvector $\vec{e_{1}}$, $\vec{e_{2}}$ and $\vec{e_{3}}$ as the direction closest to the orientation of Virgo (effectively putting Virgo into the 1st octant of the coordinate system). This new setting of the eigenframe is shown in Figure 9. Similarly to the previous Figure 8, Figure 9 shows the location of MW (left) and M31 (right) subhaloes entry points as they cross the $2\times R_{200}$ threshold of their respective main host. The colormap and colored circles are the same as in 8. The red point and error bars correspond respectively to the mean and standard deviation of the locations of Virgo in the three simulations. As in Figure 8, the infall pattern indicates that subhaloes appear to be accreted from the general direction of Virgo. Put another way, the $\vec{e}_{3}$ vector points towards Virgo.