Observing the stellar halo of Andromeda in cosmological simulations: the Auriga2PAndAS pipeline

The mock stellar catalogues used as inputs of this pipeline are computed from the Auriga simulations (Grand et al., 2017) in a very similar manner to what was previously presented by Grand et al. (2018b) to produce the Aurigaia mock stellar catalogues. To generate these stellar mocks, two methods were presented, the HITS-mocks and the ICC-mocks. These two methods differ in how the “stars,” generated from the stellar particles of the simulations are distributed in phase space, as well as by the choice of stellar evolution models. For the rest of this paper, we use exclusively the ICC-mocks, based on Lowing et al. (2015). The reason is that this method produces smoother distributions of “stars” and avoids discrete clumps at the coordinates of the parent stellar particles, while still preserving the phase-space distribution. Therefore, with this method, the presence of artificial features, such as fake clumps, in the distribution of “stars”²²2In the rest of the paper, the mock stars generated from a simulation’s star particles will be referred as “stars” to differentiate them from the real stars of the PAndAS observations. from the stellar halo in the Auriga simulations is minimized. With this method, each parent stellar particle is split into $N$ “stars,” assuming a Chabrier IMF, and drawing from a model stellar population of a given age and metallicity according to the parsec isochrones (Bressan et al., 2012), which take into account the evolution of the massive end of the IMF for stellar particles with old ages. The stellar parameters and the absolute magnitudes in the CFHT $g$ and $i$ bands³³3Here, we use the pre-2014 CFHT/MegaCam photometric system, before the current set of filters was built. It is worth noting that the fields obtained before 2007 have been observed with a slightly different $i$-band filter (see Figure 3 of McConnachie et al., 2018) but this change is subtle enough that it should not affect our results. are computed for each mock “star” from the parsec isochrones (Bressan et al., 2012), using the age and metallicity of the parent stellar particle.

The foreground MW extinction was added at a letter step, as will be described later in Section 2.2. Indeed, our study is focused on the region of M31 observed by the PAndAS survey, i.e. the stellar halo, where the extinction caused by the interstellar medium (ISM) of M31 is negligible. The large majority of the dimming of stars observed in the stellar halo of M31 is due to the foreground extinction caused by the ISM of the MW, and is included later in the pipeline.

It is important to note here that for the analysis of the simulations, detailed in Section 3, the central projected 20 kpc (or $1.47\arcdeg$ at the chosen distance of M31) are not taken into account. This region is severely incomplete in the PAndAS survey and a study of this region would require dedicated work (e.g. HST PHAT survey Dalcanton et al., 2012). Therefore, from this stage and for each mock, we decided to completely avoid “stars” within the a central sphere of 20 kpc radius. This action considerably reduces the size of the mocks and the computation time to apply the Auriga2PAndAS pipeline since the large majority of the simulation “stars” reside in this region. By applying this cut to the three-dimensional distances instead of the projected distances at this step, we keep the possibility to project the galaxy using a random point of view around the galaxy (see below).

From these raw mocks that include all “stars” down to an absolute magnitude of $M_{i}=2$ ($\sim 6$ magnitude fainter than the tip of the RGB at the distance of M31), we degrade the data so as to be as close as possible to the observed PAndAS data. The different steps we perform are:

1. 1.

   placing the simulation at the distance of M31, orienting them similarly to M31 and projecting it on the sky;

2. 2.

   masking out “stars” that are not in a PAndAS field or behind saturated foreground stars;

3. 3.

   computing the apparent magnitude of the “stars” and their associated uncertainties;

4. 4.

   making the data incomplete following the observed, field-specific completeness functions;

5. 5.

   selecting the RGB “stars” with a color-magnitude cut;

6. 6.

   adding the contamination from foreground MW stars and background unresolved galaxies.

The effects of these steps on the distribution of “stars” for one of the halos is visible in Figure 1. We now discuss each of these stages in turn.

Following M13, we keep stars that have a color and magnitude compatible with those of M31 RGB stars and remove “stars” outside of the selection box whose $(i_{0},(g-i)_{0})$ vertices are $(21,0.7),(21,2.3),(23.5,0.4),(23.5,1.6)$. For halo H23, this step yields the distribution shown in panel (d) of Figure 1.

Despite this selection, a large fraction of PAndAS stars present in this region of the CMD are actually contaminant objects. The source of this contamination is mostly due to the foreground Milky Way dwarf stars and to unresolved background galaxies that appear as point sources. Therefore, in our attempt to produce realistic “observed” stellar halo mocks, it is important to add this contamination component, especially since the density of contaminants is not constant over the survey because of the increasingly dense MW disk towards the North. The addition of this component is extremely important, in order to treat the observations and the "realistic" mocks with the same methods, since most analyses will remove, or at least take into account, the contamination to enhance the signal produced by the stars from the stellar halo of M31.

The model of the contamination used by the pipeline is based on the 4-dimensional (spatial and color-magnitude) model of M13, for which the density of contaminant objects at a position ($\xi,\eta$) and at a given color and magnitude ($g-i,i$) follow an exponential such that

Practically, and following M13, the number of contaminant objects per spatial pixel of 15 $\times$ 15 arcmin², $N_{j}$, is computed using Eq.(5) for the value of $\alpha_{(g-i,i)}$, $\beta_{(g-i,i)}$ and $\gamma_{(g-i,i)}$ and integrating the model counts in the CMD region delineated by the M31 RGB selection box mentioned above. This number is multiplied by $\mathcal{F}=0.93$, since I14 show that $\simeq 7\%$ of the objects of the external region over which the contamination model was constructed, and which are in a color-magnitude region similar to the M31 RGB box, are actually stars from the halo of M31. Moreover, to incorporate statistical fluctuations to the contamination model, the actual number of contaminants for each pixel is drawn randomly, assuming a Poissonian distribution centered on $N_{j}$. The spatial locations of the $N_{j}$ contaminant particles of a given pixel are then randomly distributed following a uniform distribution over the pixel.

The $g$ and $i$ magnitudes of the contaminant particles are then randomly distributed following the probability distribution function (PDF) of the contaminant CMD at the location of the spatial pixel of the “stars” (see M13). It is important to note here that the magnitudes given by the model are already corrected for the extinction. The photometric uncertainties of each contaminant “star” are then computed using Eq.(2). Finally, only the contaminant “stars” that are contained in the M31 RGB box are kept. This removes the few particles from pixels cut by the M31 RGB box that are outside this box. Depending on the science goal with the mocks, the M31 RGB box criterion can be removed from the overall pipeline (i.e. for particles from the stellar mock and for contaminants). However, in that case the number of contaminants will be slightly overestimated, as mentioned earlier, but by less than $7\%$, since this number was obtained in the M31 RGB box, which contains most of the M31 halo stars.

The final distribution of the “stars” for mock H23, including “stars” from the stellar mock and from the contamination model, is shown in panel (e) of Figure 1.

The density of stars for a full realization of this halo is presented in Figure 6. The lower left panel shows the total density of “stars” that are in the color-magnitude space delimited by the M31 RGB box. This includes the “stars” from the stellar mock, whose distribution is shown in the upper left panel, and contaminants, shown in the upper right panel. While the details of the substructures are of course different because of the difference in the accretion history of the simulated galaxy and Andromeda, we note qualitative consistency between the mock and the PAndAS data.

Figure 7 shows a mapping of the structures present in the stellar halo of M31 and in the 6 mocks. Each of these maps is a red-green-blue image, made from the combination of matched filter (MF) maps whose filters are CMD models of old RGB stars at the distance of M31, convolved by the photometric uncertainties. The background model of the MF technique is the contamination model of M13. The blue image corresponds to a signal from RGB stars with a metallicity of $[\mbox{Fe/H}]=-2.3$, the green to $[\mbox{Fe/H}]=-1.4$, and the red to $[\mbox{Fe/H}]=-0.7$. The maps are made of spatial pixels of $3^{\prime}\times 3^{\prime}$ and have been smoothed by a Gaussian kernel of one pixel width. This is similar to the MF technique used by M13 to produce the map of their Figure 2. The holes visible at the very center of the maps for the mocks are caused by the cut made to mask the stars within a three-dimensional radius of 20 kpc, as described in Section 2.1.

A visual inspection of these maps shows that the simulated stellar halos present numerous, well-defined substructures on a similar scale to those observed around M31 (central panel). In the central 50 kpc, the simulations and observations are very similar, with an inner stellar halo dominated by metal rich-stars ([Fe/H] $\sim-0.7$) and a small number of identifiable structures. Most of the mocks, with the exception of H6 and H24, have a large and metal-rich stream or shell, sign of a recent, or on-going, massive accretion. These structures are similar to the Giant stream⁷⁷7We follow Lewis et al. (2013) and McConnachie et al. (2018) for the nomenclature of the halo structures of M31. (Ibata et al., 2001) visible in PAndAS that is the consequence of the accretion of a galaxy with a mass similar to the Large Magellanic Cloud $2-3$ Gyrs ago (Fardal et al., 2013).

On average, the mock halos are redder, and so more metal-rich, than the observations. The halo of M31 presents more structures, such as stellar streams or clouds, at intermediate/low metallicity than the simulations. It is important to note here that the absence of a galaxy similar to M33 in some of the simulations is not surprising. Indeed, M33 would not be visible in the PAndAS footprint if it were located at the same projected distance but at a different angle around M31, as it is likely to be the case in the mocks. Furthermore, not all simulated halos presently have a satellite that is as massive as M33.

For each of the halos, we compute the surface density of RGB stars in spatial pixels of $0.1\arcdeg\times 0.1\arcdeg$ for different ranges of metallicities, in the same way as for the observations. In each of these pixels, the number of contaminant objects have been statistically removed assuming the model of M13. To compare directly the mocks to the observations, the metallicities have been derived photometrically for each star by comparing their colors and magnitudes to the Dartmouth isochrones (Dotter et al., 2008) assuming that the stars are at the distance of M31, have an age of 10 Gyrs and an alpha enhancement of [$\alpha$/Fe]$=+0.2$, typical of the populations in the halo of M31/MW (i.e. Helmi, 2008; Kilic et al., 2019). The surface density maps in four ranges of metallicities (all, metal-poor, intermediate metallicity and metal-rich stars) for the halo of M31 and for mock H$23$ are shown on Figure 8. The information this figure contains is similar to that visible in Figure 7, but it is easier to see the contribution of each structure in the different ranges of metallicity. For instance, it clearly shows that H23 has a similar number of dwarf galaxies at intermediate metallicities compared to M31, but has only about half its number of metal-poor satellite galaxies similar to AndXXI or AndXXIII ($M_{V}\sim-9.5$). This characteristic is present in all mocks, especially for H6 and H21, for which the metal-poor galaxies are far less numerous but also more centrally distributed than observed around M31.

We use these surface-density maps to derive the median surface density radial profile of RGB stars in the four ranges of metallicities for all the halos. The median profile is favored over the mean one since it limits the contribution from compact substructures, such as satellite galaxies. The resulting profiles are shown in Figure 9 and fitted with a single power-law profile in the range of projected distances of $23-120$ kpc ($9\arcdeg$ at the distance of M31). The resulting power-law slope values are listed on Table 2. We purposefully avoid the region beyond 120 kpc since it was used to construct the contamination model of M13, leading to profiles that could be systematically biased at these radii, despite the correction mentioned in Section 2.3.

When considering the surface density profiles for the full metallicity range, all mocks apart from H27 have profile slopes relatively similar to M31. However, most of the mocks are more centrally populated than M31, with higher surface densities, up to a factor $\sim 4$, until $\sim 75$ kpc. Beyond $\sim 100$ kpc, most of the mocks are less densely populated than the observed M31 halo. It is also interesting to note here that the profile of H21 shows a clear drop at $\sim 75$ kpc, also visible in the metal-rich range. This is caused by the 2 major accretions that are ongoing and largely dominate the rest of the halo population up to $\sim 75$ kpc.

By comparing the density profiles in the different metallicity intervals, it is clear that the global surface density profiles of the stellar halos are dominated by metal-rich and intermediate metallicity stars, both in M31 and in the mocks. However, the difference of density between the observations and the simulations are most prominent for metal-rich stars ($-0.7\leq$[Fe/H]$<0.0$), for which the density of stars in the mocks is higher than the observations at all the radii. In this metallicity range, H16 is the exception, with surface density profiles that are similar to the ones observed for M31 beyond 30 kpc. At intermediate metallicities ($-1.7\leq$[Fe/H]$<-0.7$), the differences are smaller than in the metal-rich regime, with profile of the mocks being, on average, $\simeq$ 1.5 more populated that of M31 at intermediate metallicities, against $\simeq$ 2.2 at metal-rich metallicities. Overall, the mocks are still more densely populated than observed until a projected radius of $\sim 100$ kpc. Beyond that distance, the density in most of the simulated halos drops below the observed density.

For the metal-poor star selection ($-2.5\leq$[Fe/H]$<-1.7$), half of the mocks are $\sim 1.7$ times more populated than M31. The other half have a similar inner density than observed in M31. However, in this metallicity range, the surface density profile of the mocks is systematically steeper than observed, and for most of them, the distant halo ($>60$ kpc) shows a clear deficit of metal-poor stars compared to the observation. H27 is an exception compared to the other simulations since, for all metallicity intervals, the surface density of the mock is always higher than observed in M31 at any radius.

Given these observations, we can conclude that the halos of the simulated galaxies are, in general, more populated than the M31 halo, especially in the central region. In particular, the simulations are more populated by metal-rich stars than observed around M31 at any radius. Moreover, the metal-poor stars in the mocks are more centrally concentrated and are less extended than observed by PAndAS around M31.

In the previous section, the analysis is partly based on photometric metallicities derived using the Dartmouth isochrones, while the Padova isochrones are used to give an absolute magnitude to the “stars” in the simulations. In this section, we show that the choice of the set of isochrones has a negligible impact on our analysis, and that this choice does not alter the conclusions that we draw from the analysis of the mocks.

For the comparison between the mocks and the observations, the Dartmouth isochrones were preferred over the Padova ones for mainly two reasons. The first is that, historically, the previous works done by the PAndAS collaboration were used mainly the Dartmouth isochrones (e.g. McConnachie et al., 2009; Collins et al., 2010; Martin et al., 2013; Ibata et al., 2014; McConnachie et al., 2018). Thus, using the same set of isochrones facilitates a direct comparison with these previous analyses. The second is that using a set of isochrones that is independent from the set used to generate the stellar particles in the simulations removes any suspicion of possible biases in the analysis.

Figure 10 displays the median surface density profile for mock H23 in the four previous ranges of metallicities, for the case where the photometric metallicities are derived using the 10 Gyr old Dartmouth isochones with [$\alpha$/Fe]=+0.2 (purple lines), like in the previous section, and with the Padova isochrones (blue dashed lines) corresponding to a population of 8 Gyr, the median age of the “stars” in this specific halo. This figure also include the median surface profile of the un-contaminated “stars” (i.e. that does not include the “stars” from the foreground model), using the metallicities directly provided by the simulations (red dots). These different profiles are very similar in all metallicities ranges, except for the metal-poor regime. In this range, the two profiles obtain using the photometric metallicities are still very similar to each other but are systematically higher than the“true” profile by a factor $\sim 2$. This is explained by the fact that isochrones of metal-poor stars are very close to each other in a CMD. Thus, with the photometric uncertainties and the intrinsic scatter of a stellar population, metal-poor stars overlap each other in a CMD, leading to a lower accuracy of the photometric metallicities in that regime. Actually, it seems that a photometric metallicity range of $-2.5\leq$[Fe/H]${}_{\mathrm{photo}}<-1.7$ corresponds to a range of $-2.5\leq$[Fe/H]${}_{\mathrm{sim}}<-1.5$ for the metallicities of the mocks, for which the profiles are represented by the green dots on the lower-right panel.

However, it is important to keep in mind that our analysis is based on relative metallicities, with four broad ranges. Therefore, whether a given “star” has [Fe/H]=-1.7 or [Fe/H]=-1.5 does not significantly impact the fact that metal-poor stars are less numerous than the intermediate or metal-rich ones. In addition, we perform the exact same analysis on the PAndAS observations and on the Auriga mocks. The conclusions that we draw are therefore based on a relative comparison of these results and do not rely of the absolute metallicity values.

It has been shown with semi-analytical and hydrodynamical simulations of L$\star$-galaxies that their stellar halo can be decomposed into two populations with very different origins, the in-situ and the accreted components (also referred to as the ex-situ component in the literature, e.g. Bullock & Johnston, 2005; Cooper et al., 2010; Font et al., 2011; Purcell et al., 2011; Pillepich et al., 2014; Monachesi et al., 2019). The accreted component is made of stars initially hosted by dwarf galaxies and globular clusters that have since been disrupted by tidal effects and and which have deposited their stars into the main halo. The in-situ component is made of stars from the galactic disk that have been ejected into the halo due to the dynamical heating produced by the interactions with sub-halos, galaxies or massive molecular clouds, and by stars that were born in the halo of the proto-galaxy.

To differentiate the contribution of these two populations in the mocks, we redo the analysis done in Section 3.1, but select only the accreted “stars” in the mocks. Mock “stars” are classified as accreted or in-situ according to their parent stellar particle, following the definition presented in Section 2.3 of Monachesi et al. (2019). In this definition, stellar particles bound to the host galaxy at birth are considered as in-situ stars.

As visible in Figure 11, considering only the accreted stars, the surface density profiles of the simulated galaxies are generally in better agreement with the profile of M31. In all the different metallicity intervals, the shape of the surface density profiles of the accreted “stars” alone are closer to the profiles observed in M31 than if we consider the full stellar halo populations (in situ and ex situ), especially for simulation H23. The strongest deviations between the surface density profiles of the accreted populations and the overall population (accreted+in-situ) are visible with the most metal-rich stars, followed by the intermediate metallicities range. This is not surprising since the accreted stars are on average more metal-poor than the stars formed in-situ, as seen in the Milky Way (i.e. Carollo et al., 2007, 2010; Belokurov et al., 2019) and in different suites of cosmological simulations (i.e. Purcell et al., 2011; Pillepich et al., 2014; Monachesi et al., 2019).

As visible on Figure 12, the accreted stars make up the majority of the stellar halos (at least beyond 23 kpc) for all metallicities ranges, as noted by Monachesi et al. (2019). However, the radial profile of the fraction of accreted stars is very different in the different metallicity intervals. Indeed, for metal-rich stars, $\simeq 50-65\%$ of them are accreted at 23 kpc and this fraction increases slowly to reach $\simeq 80-90\%$ at 120 kpc, while for the intermediate and metal-poor stars the fraction of accreted stars is almost constant across all radii, accounting for $\sim 70-80\%$ and $\sim 90\%$ of the stars in these respective ranges. This confirms that the in-situ component is mostly composed of centrally concentrated metal-rich stars, similar to what is observed in the Milky Way or in different simulations (i.e. Purcell et al., 2010; Cooper et al., 2015; Pillepich et al., 2015; Monachesi et al., 2019; Sanderson et al., 2018; Belokurov et al., 2019). Note that the ratio shown in this figure corresponds to the ratio of stars from the smooth halo component (i.e. not in sub-structures). Taking into account the overall population, the ratios are similar for all mocks, except for H16 and H24 (referred as the Full populaiton of Figure 12), for which the fraction of in-situ stars is significantly higher, up to $\sim 35$ and $25$, kpc respectively. This is caused by the presence of very extended galactic disk in these simulations, as noticed by Grand et al. (2018a). However, the goal of this section is to compare the fraction of accreted stars in the halos of the different simulated galaxies and so we purposefully do not take into account the stars in the disk or those present in sub-structures.

By comparing the surface density profile observed in M31 in the different ranges of metallicity and the profiles of the accreted and of the accreted+in-situ populations in the mocks, it seems that the mocks show a density of in-situ stars that is about twice as large as observed in M31. Indeed, by reducing roughly by a factor 2 the contribution of the in-situ component, the surface density profiles will be closer to the profile observed in PAndAS, especially for H16, H21 and H23. This is most visible when comparing the accreted and overall metal-rich surface density profile in the inner-halo ($<30$ kpc). This conclusion is in agreement with the observation of Monachesi et al. (2019) using the Auriga simulations, and similar results have recently been found by Merritt et al. (2020) comparing the Illustris TNG100 simulations (Springel et al., 2018; Pillepich et al., 2018b) to the Dragonfly Nearby Galaxies Survey (Merritt et al., 2016).

However, even by reducing the number of the stars formed in-situ by a factor of 2, most of the mocks have halos that are overall more populated than M31’s and their profiles present a great diversity. Therefore, outside the inner halo, where the majority of the in-situ stars are located, this scenario is not the only explanation for the observed difference between M31 and the simulated galaxies. Because galaxies like H$23$ have a very similar profile to the observed one, this overpopulation of most of the galaxies and the diversity of profiles is possibly driven by the stochastic process on the merger histories of each galaxies.

From all the simulated halos analysed here, the halo H$23$ is the most representative of the halo of M31. Not only does this halo have a higher number of dwarf galaxies similar to AndXXI or AndXXIII than the other simulations, but its surface density profile is the closest to the profile constructed from PAndAS in all metallicities ranges, especially for the accreted “stars”. However, as for all the simulated halos, the in-situ component is $\sim 2$ times more populated compared to M31.

As we have seen in the previous sections, the Auriga-PAndAS mocks are, overall, reasonable approximations of the stellar halo of M31 and its diversity of structures. Qualitatively, the inner halo of the simulations are very similar to the inner halo of M31. Moreover, most of the halos are populated by a number of intermediate metallicity dwarf galaxies that are similar to what we see around M31. They cover a broad range of luminosities, from galaxies with a size and a metallicity similar to NGC 147 and NGC 185 (corresponding to a stellar mass of $\simeq 5\times 10^{7}$ M_⊙), to galaxies similar to AndXXI and AndXXIII (corresponding to stellar masses of $\sim 1\times 10^{6}$ M_⊙). However, in certain aspects, the mocks differ from what is observed around M31. Although the goal of the paper is not to analysis the reasons of the (overall) small differences between the simulations and the observations, we mention hereafter a few points that could be explored in future works.

For instance, we can notice that the simulated stellar halos are on average more metal-rich than observed in M31. As explained in the previous section, this is likely a consequence of an over representation of the in-situ stars in the simulations, these stars being more metal-rich than the accreted ones. As mentioned by Monachesi et al. (2019), this over-population of stars formed in-situ could be a consequence of the disks of the host galaxy, which, in the Auriga simulations, are typically larger than observed for M31 (see Grand et al. (2017) and Monachesi et al. (2019) for a discussion on the size of the disk). However, even with considering only the accreted stars, the simulations tends to be more metal rich than observed for all radii, which might be a consequence of the properties of the accreted dwarf galaxies that formed these structures but also, potentially, an under-quenching of the star formation feedback in the last billion years, especially in the most massive satellites, as proposed by Monachesi et al. (2019). Moreover, the surface density profile is lower that the majority of the simulation out to $\sim$ 70-80 kpc, even when considering only the accreted population. This is surprising given that the simulated galaxies have a similar or lower stellar mass than observed in M31. This could also be a consequence of the under-quenching of the star formation feedback, of the way in which the stellar particles are split into “stars” in the simulations, and/or of the IMF choice or the stellar populations models.

One of the other noticeable differences between the simulations and the observations of M31 is the relatively small number of metal-poor dwarf galaxies present around the simulated galaxies, compared to the numbers observed in PAndAS. Galaxies similar to And XIV or And XII are not visible at all in the simulations. Simpson et al. (2018) show that the number of dwarf galaxy Auriga simulations converge down to $2-3\times 10^{4}$ M_⊙ for different resolutions. However, it is possible that the small number of metal-poor dwarf galaxies in the simulation might be due to an incomplete sampling of the total phase-space of these galaxies, despite the fact that, with the ICC method, “stars” conserve the phase space of the initial stellar particle. It will be interesting to quantify these apparent differences by searching and characterizing the dwarf galaxies in the mocks in the same way they were found and characterized in PAndAS.