A galaxy's accretion history unveiled from its integrated spectrum

TL;DR

This paper introduces a novel spectral fitting method to determine a galaxy's accretion history, including the total ex-situ mass fraction and the distribution of accreted satellite masses, using integrated spectra.

Contribution

The paper presents a new technique that accurately derives a galaxy's accretion history from its integrated spectrum, linking stellar populations to merger events and progenitor masses.

Findings

Method recovers total accreted fraction within ~12%.

Accreted satellite stellar mass can be estimated within ~26%.

Slope of satellite mass function is recoverable within ~16%.

Abstract

We present a new method of quantifying a galaxy's accretion history from its integrated spectrum alone. Using full spectral fitting and calibrated regularization techniques we show how we can accurately derive a galaxy's mass distribution in age-metallicity space and further separate this into stellar populations from different chemical enrichment histories. By exploiting the fact that accreted lower mass galaxies will exhibit an offset to lower metallicities at fixed age compared to the in-situ stellar population, we quantify the fraction of light that comes from past merger events, that are long since mixed in phase-space and otherwise indistinguishable. Empirical age-metallicity relations (AMRs) parameterized for different galaxy masses are used to identify the accreted stellar populations and link them back to the progenitor galaxy's stellar mass. This allows us to not only measure the host galaxy's total ex-situ mass fraction ($f_{acc}$), but also quantify the relative amount of accreted material deposited by satellite galaxies of different masses, i.e. the accreted satellite mass function in analogy to the subhalo mass function. Using mock spectra of simulated, present-day galaxies from the EAGLE suite we demonstrate that our method can recover the total accreted fraction to within $\approx 12 \%$, the stellar mass of the most massive accreted subhalo to within $\approx 26 \%$ and the slope of the accreted satellite mass function to within $\approx 16 \%$ of the true values from the EAGLE merger trees. Future application of this method to observations could potentially provide us accretion histories of hundreds of individual galaxies, for which deep integrated light spectroscopy is available.

Figures

Tables

Table 1. Table 1: Overview of the nine EAGLE simulated galaxies. From left to right column : Nr., GalaxyID , total stellar mass inside 100 kpc as well as the mass-weighted mean age and metallicity.

Nr.	GalaxyID	$M_{\star }$	$⟨\mathrm{Age}⟩$	$⟨\left[M/H\right]⟩$
		$\left[M_{\odot }\right]$	$\left[\mathrm{Gyr}\right]$	$\left[\mathrm{dex}\right]$
1	21109760	$6.77\cdot {10}^{11}$	10.43	-0.068
2	20163968	$4.69\cdot {10}^{11}$	9.83	-0.052
3	16971377	$6.72\cdot {10}^{10}$	9.62	0.135
4	9944856	$1.46\cdot {10}^{10}$	8.24	0.050
5	5368409	$1.11\cdot {10}^{10}$	9.26	0.008
6	12092377	$9.09\cdot {10}^9$	9.54	-0.136
7	11487828	$6.82\cdot {10}^9$	7.20	-0.014
8	6031050	$3.68\cdot {10}^9$	8.88	-0.096
9	5443067	$3.29\cdot {10}^9$	8.62	-0.187

Equations

$$Z(t)=-p({\mathrm{M}}_{\star })\ln \mu (t)$$

$${\log }_{10}p({\mathrm{M}}_{\star })=p_0+{\log }_{10}\left(1-\exp {\left(-\frac{{\mathrm{M}}_{\star }}{{\mathrm{M}}_0}\right)}^{{\alpha }_p}\right),$$

$$[\mathrm{M}/\mathrm{H}]=[\mathrm{Fe}/\mathrm{H}]+{\log }_{10}\left(0.694\cdot 10^{[\alpha /\mathrm{Fe}]}+0.306\right).$$

$$\mu (t)=\frac{t-\left(13.5-t_{\mathrm{form}}\right)}{t_{\mathrm{form}}},$$

$$t_{\mathrm{form}}=\min \left[{\left(14-{\log }_{10}{\mathrm{M}}_{\star }\right)}^{{\alpha }_t},13.5\right]$$

Full text

††pubyear: 2018††pagerange: A galaxy’s accretion history unveiled from its integrated spectrum–A

1 Introduction

A galaxy’s total stellar mass budget can be split into two categories: in-situ star formation that originates from cooling of the galaxy’s native gas content and acquisition of ex-situ stars through merging with satellite galaxies. The importance of the ex-situ or accreted material relative to in-situ star formation is expected to be very stochastic, but also strongly mass dependent. For example, cosmological simulations predict that higher mass galaxies with $\mathrm{M}_{\star}\gtrsim 10^{11}\ \mathrm{M}_{\odot}$ have higher ex-situ fractions between 50 and 90 %, while Milky Way (MW) mass haloes typically only exhibit ex-situ mass fractions of 10-20 % (Cooper et al., 2013; Cooper et al., 2015; Rodriguez-Gomez et al., 2016; Qu et al., 2017; Pillepich et al., 2018).

The measurement of the ex-situ contribution to galaxy evolution is crucial to understand how mergers impact the dynamical, structural and, as we will show here, element abundances of a galaxy population. While the mass distribution of surviving and merged subhaloes is a robust prediction of $\Lambda$CDM in the absence of baryons, there exists an intrinsic degeneracy between dark and baryonic accretion histories due the effect of feedback processes on the baryons. Quantifying past accretion events in observed galaxies is therefore important, but intrinsically difficult for several reasons. For one, the remnant signatures of these events are located in observationally challenging regimes, as ex-situ dominated regions are mainly in the outskirts (> 30 kpc) of galaxies, where the surface brightness becomes very low ($\mu_{g}>25\ \mathrm{mag}/\mathrm{arcsec}^{2}$). Secondly, the dissolution of the merger in the host potential can smoothly distribute stellar material over a large fraction of the virial radius, which after several Gyrs, becomes difficult to identify as projected overdensities or kinematically coherent structures in external galaxies.

For these reasons, the most accurate constraints on accretion histories come from deep photometric surveys of the stellar haloes of Local Group galaxies (e.g. MW: Bell et al., 2008; Carollo et al., 2010; Iorio et al., 2018 and M31 McConnachie et al., 2009; Courteau et al., 2011) and beyond (e.g. Radburn-Smith et al., 2011; Monachesi et al., 2016; Crnojević et al., 2016), which resolve the stellar haloes into individual stars. The resolved photometry can provide information on discrete stellar populations of ex-situ stars in the diffuse outer haloes of galaxies. For example, the observed relation between the metallicity and stellar mass of the halo (Harmsen et al., 2017) can be used to understand the mass assembly history by comparing to correlations found from cosmological simulations (Font et al., 2011; D’Souza & Bell, 2018; Monachesi et al., 2018).

Other studies such as Tonini (2013); Leaman et al. (2013); Kruijssen et al. (2018); Beasley et al. (2018); Mackey et al. (2018); Hughes et al. (2018) use globular clusters (GCs) as bright tracers of past accretion events. By identifying groups of GCs as coherent structures in phase-space or with similar stellar population properties, they can be quantitatively linked back to the stellar mass of the accreted host or the accretion time.

Beyond ($\gtrsim 10$ Mpc) there exist also many deep photometric integrated light studies of stellar haloes, which uncover very low surface brightness substructures and streams (e.g. Martínez-Delgado et al., 2010; Spavone et al., 2017; Huang et al., 2017; Merritt et al., 2016; Hood et al., 2018). In order to estimate an accreted fraction, the measured surface brightness profiles are usually structurally decomposed into multiple components (i.e. an in-situ dominated central component and an accretion dominated "halo" component, see e.g. Huang et al., 2016; Iodice et al., 2016; Spavone et al., 2017) .

These observational methods are restricted in the sense that they provide a single value for the ex-situ contribution or estimate bulk properties of the stellar halo. They also make structural assumptions about where the in-situ and ex-situ contributions should dominate the surface brightness spatially. However, cosmological simulations show that massive galaxies ($\mathrm{M}_{\star}\gtrsim 10^{11}\ \mathrm{M}_{\odot}$) can still have up to 60% of accreted material in the innermost 10 kpc from the galaxy’s centre (Pillepich et al., 2018). Furthermore, the age-metallicity degeneracy, which is especially prominent in colors from integrated photometry, means that the chemical properties are difficult to uniquely interpret.

Given these observational challenges, and the expected variety of merger histories in the galaxy populations, there is a clear need for an observational method that measures the accretion history of galaxies in more detail - including in their central regions as well as in a fashion that can be compared directly with simulations. If such a method can probe the luminous distribution of long disrupted substructures, which are otherwise indistinguishable, it would represent a novel test of galaxy formation models. In addition, if we want to understand the impact of the merger history on other galaxy properties like their morphology or dynamics, we need to go beyond the local universe, as only there we will acquire the needed statistics.

In this paper we present a novel technique to quantify a galaxy’s accretion history from its integrated spectrum, and verify the method with the help of the hydrodynamical cosmological EAGLE simulations (Schaye et al., 2015; Crain et al., 2015). In section 2 we introduce the process of constructing the merger histories and mock spectra from the simulations. In section 3 we describe our method. We first model the age and metallicity distribution of stellar populations in a galaxy with the help of full spectral fitting (pPXF: Cappellari & Emsellem, 2004; Cappellari, 2017) using updated and carefully calibrated regularization methods. The ex-situ contributions to the total stellar mass in age-metallicity space are then linked to satellite galaxy masses by empirically derived mass-dependent age-metallicity relations (AMRs). In section 4 we show that this method provides an observational estimate of the abundance of accreted satellite galaxies as a function of their stellar mass prior to accretion - an analogue to the unevolved111Unevolved denotes that the mass is measured before infall meaning that no mass loss due to tidal interactions has occurred yet. subhalo mass function (Giocoli et al., 2008; Jiang & van den Bosch, 2016). In section 5 we end with a conclusion and an outlook of our method to future applications.

2 Simulated Spectra and Merger Histories from the EAGLE Simulation

In order to test whether our new method recovers a known accretion history, we make use of galaxies from the EAGLE simulation suite (Schaye et al., 2015; Crain et al., 2015), where we have an exact knowledge of a galaxy’s true mass distribution in age-metallicity space, and its accretion history. The EAGLE simulation is a hydrodynamical cosmological simulation and therefore treats the evolution of baryonic and dark matter self-consistently and adopts a series of subgrid models in order to account for baryonic physics below the resolution scale. In particular, the parameters of the subgrid models were calibrated to ensure reproduction of the $\mathrm{z}=0$ galaxy stellar mass function. The data from the EAGLE simulation have been publicly released in McAlpine et al. (2016); The EAGLE team (2017).

We utilize the particle information of nine early-type galaxies at $\mathrm{z}=0$ spanning a stellar mass range from $10^{9}$ to $10^{12}\ \mathrm{M}_{\odot}$ contained in the largest volume simulation ("REFL0100N1504"), which has a baryonic particle mass resolution of $1.81\cdot 10^{6}\ \mathrm{M}_{\odot}$. This sample is not representative and was not picked for any other reason than to explore the method’s success with a range of stellar masses and star formation histories. Therefore, the demonstration of our developed method throughout section 3 and 4 will focus on a single, central galaxy with a total stellar mass of $\mathrm{M}_{\star}\approx 6.77\cdot 10^{11}\ \mathrm{M}_{\odot}$, which we will refer to as Galaxy Nr. 1 according to Table 1. Additional information about basic properties of our full galaxy sample can be found in the appendix A.

For every stellar particle we have the basic information like its mass, age, total metallicity [M/H] and position, but also the SnapNum (or redshift), the GroupNumber and SubGroupNumber222The GroupNumber and SubGroupNumber are integer identifiers of the Friends-of-Friends halo and subhalo respectively, which a given galaxy currently resides in. Those values are not unique across snapshots. at which the particle first appeared as well as a flag, whether the particle was born in-situ or ex-situ. This additional information is obtained by stepping through all the snapshots and identifying particles that later end up in our selected galaxy at $\mathrm{z}=0$. Ex-situ particles are identified as those not born on the main branch of the simulated galaxy (Qu et al., 2017). For this study, we will focus on particles within a 3D aperture of 100 kpc with respect to the galaxy’s centre of mass.

As we want to go one step further and associate each ex-situ particle to a host galaxy of mass $\mathrm{M}_{\mathrm{sat}}$ prior to merging with the primary halo, we step through the merger trees and track those particles. The birth SnapNum, GroupNumber and SubGroupNumber of every ex-situ particle in our simulated galaxy uniquely identifies the GalaxyID333This is a unique integer identifier of a galaxy in the simulation. of the progenitor. By traversing up the merger tree, we can identify when the progenitor (i.e. the galaxy containing the ex-situ stars) merged onto the main branch of the host galaxy. We follow Rodriguez-Gomez et al. (2015, 2016) and quantify the stellar mass of the subhalo at the time $\mathrm{t}_{\mathrm{max}}$, where it reaches its maximum mass prior to merging - as this mass should provide the most representative indication of its chemical enrichment efficiency.

In the left panels of Figure 1 we show the stellar surface mass density in the xy-plane separated into the in-situ and ex-situ components as well as the mass distribution in age-metallicity space for the entire stellar population of a EAGLE simulated galaxy with $\mathrm{M}_{\star}\approx 4.69\cdot 10^{11}\ \mathrm{M}_{\odot}$ (Galaxy Nr. 2, see Table 1 for basic properties). While the ex-situ component is more spatially extended than the in-situ component, there is significant overlap spatially, and the bulk of the ex-situ component is centrally concentrated and exhibits a smooth distribution. In the right panels of Figure 1 we plot the ex-situ age-metallicity distribution further split up into the corresponding stellar mass ranges of the satellite galaxies they came from. This high mass galaxy had a 1:1 merger, but also several lower mass ratio accretion events, which brought in stellar material that is on average more metal poor at fixed age.

In order to test whether our machinery can successfully extract a galaxy’s accretion history from its integrated spectrum alone, we construct a mock spectrum from the simulated galaxy. We represent each stellar particle within those 100 kpc as a single stellar population (SSP). We use the SSP model library MILES (Vazdekis et al., 2010; Vazdekis et al., 2015) based on the BaSTI isochrones (Pietrinferni et al., 2004, 2006). These cover an irregular grid with 53 ages spanning from 30 Myr to 14 Gyr and 12 metallicities from -2.27 dex to 0.4 dex. We choose a bimodal initial mass function (IMF) with a slope of 1.3444The EAGLE simulation uses a Chabrier IMF (Chabrier, 2003), however as we will use the same SSP models for fitting the mock spectrum in section 3.1, this does not influence our results in any way. (Vazdekis et al., 1996) and added no additional dust. In this way every particle is assigned a representative spectrum suitable for its age and metallicity. The integrated spectrum of the simulated galaxy is then obtained by summing up all the SSP spectra weighted by their corresponding particle stellar mass, because the luminosity of the SSP models is expressed per unit solar mass. Naturally, the ages and metallicities of the particles in the simulation have more values than the SSP library. However, we deemed interpolating among the models was unnecessary, as the mean absolute difference (MAD) between the integrated spectrum using interpolated SSPs and non-interpolated ones was less than 0.08%. Finally, we added random Gaussian noise to the integrated spectrum to achieve a SNR of 100 Å*-1* for the test presented here.555In Boecker (2018) we investigated the SNR variations, which had an impact on the recovered distribution in age-metallicity space from regularized full spectral fitting, but was less significant with regard to the accreted satellite mass function.

3 Methodology

3.1 Recovering an extended mass distribution in age-metallicity space from an integrated spectrum

The first step in our method is to use full spectral fitting, where the mock integrated galaxy spectrum is fitted with a linear combination of SSP model templates in order to infer stellar population parameters like age and metallicity. We do this by using pPXF (Cappellari & Emsellem, 2004; Cappellari, 2017), which optimizes the weights for the linear combination of SSP templates and hence provides us with a distribution in age-metallicity space instead of average quantities. Because SSP models are normalized to one solar mass, the recovered weights are mass fractions. For simplicity and to better understand the sources of uncertainty in our final quantity of interest, we fit the same SSP models that were also used when constructing the mock spectrum.

This inverse problem can in principle be solved by a simple linear least squares minimization. However, if we compare the solution from this with the true distribution of the EAGLE particle data for our particular galaxy binned to the SSP model grid in Figure 2, we see that the recovered mass fractions are distributed very sparsely and not necessarily in a physically meaningful way across the age-metallicity plane. This behaviour is due to the nature of the SSP models, as their shapes are determined by stellar physics and evolution models and are hence degenerate. Therefore, the optimal representation of the mock spectrum can be found with an ambiguous and almost arbitrary combination of SSP templates (ill-posed problem). Furthermore, due to the flexibility of optimal weight combinations, any change in the input galaxy spectrum due to noise can change the solution drastically (ill-conditioned problem).

pPXF offers a way to circumvent those problems via regularization. Mathematically speaking this dampens weight solutions that are driven by noise in the data. In an astrophysical context, it ensures a certain smoothness between neighbouring SSP bins in age-metallicity space, which is appropriate as chemical enrichment proceeds smoothly in galaxies as seen in Figure 1. The amount of smoothness in pPXF is controlled by the regularization parameter $\lambda$ and the way the weights are smoothed out is defined by the regularization matrix $\bm{B}$.

Crucially, $\lambda$ needs to be calibrated for every spectrum that is fit, as too low or too high regularization parameters will recover an astrophysically different star formation history and chemical enrichment. We follow the calibration procedure that is described in the source code of pPXF, which was adopted from Press (2007, section 19.4.1) and has been used in astronomical literature (see e.g. Norris et al., 2015; McDermid et al., 2015). The "optimal" regularization parameter represents the maximal amount of smoothness for a given data fidelity, such that the solution is still consistent with the unregularized solution. This is found in the following way: first the noise vector is re-scaled such that the reduced $\chi^{2}$ of the unregularized fit becomes unity. Then a series of fits adopting different regularization parameters is performed until the $\chi^{2}$ has increased by one standard deviation, which corresponds to a $\Delta\chi^{2}$ of $\sqrt{2N}$, where $N$ is the number of pixels included in the fit.

The default regularization matrix that is adopted in pPXF is the second order finite difference operator, $\bm{B}=\operatorname{diag}(1,-2,1)$, while the user also has the option to choose the first order one, $\bm{B}=\operatorname{diag}(1,-1)$. However, we have found that neither of those two options recover the very extended distribution in age-metallicity for our given simulated galaxies well enough. We therefore manually introduced a third order finite difference operator, $\bm{B}=\operatorname{diag}(1,-3,3,-1)$, in the source code of pPXF, which seems to recover the overall shape of the true distribution much better. The optimal regularized solution following the calibration, and using the default as well as our implemented matrix is shown in Figure 2 for comparison. The median absolute deviation (MAD) of the absolute residuals between the true mass distribution in age-metallicity space and the recovered one with the updated regularization matrix is 14 % lower than for the default matrix.

We also show in Figure 2 the best-fit spectrum for the regularized solution using the third order difference operator, however we note that in practice this is indistinguishable from the unregularized and default regularized case. The standard deviation of the residuals are 1 %, which shows that the spectrum is fitted down to the injected noise level with SNR of 100 even for the regularized case.

3.2 Flexible mass-dependent chemical enrichment templates

With a robust recovery of extended mass distributions in age-metallicity space from regularized pPXF fitting, we now need to associate the mass fractions to potential accreted galaxies. We do this by constructing flexible, mass-dependent templates in age-metallicity space, which describe how galaxies of a given mass should (on average) chemically evolve. While the detailed chemical evolution of distant high mass galaxies is not observationally constrained, here we attempt to construct a physically motivated, flexible, mass-dependent chemical framework.

For this we use results from Leaman et al. (2013), who derived empirical age-metallicity relations spectroscopically for local group dwarf galaxies from resolved stellar populations. To first order, leaky box chemical evolution models describe the metallicity distribution function (MDF) and age-metallicity relation (AMR) of those galaxies with only a galaxy mass-dependent variation in the effective yield $p(\mathrm{M}_{\star})$. For a given galaxy, the chemical evolution is described as:

[TABLE]

where $Z$ is the metallicity, $t$ (in Gyr) is time since the Big Bang, $\mu$ the galaxy’s gas fraction and $p(\mathrm{M}_{\star})$ the mass-dependent effective yield.

The mass-dependent effective yield is empirically measured for the Local Group galaxies below $\mathrm{M}_{\star}\leq 10^{9}\ \mathrm{M}_{\odot}$ in Leaman et al. (2013) (see also Lee et al., 2006) and yields the observed relation of $p(\mathrm{M}_{\star})\propto\mathrm{M}_{\star}^{\alpha_{p}}$ with $\alpha_{p}\simeq 0.4$. Above this mass, the observed mass-metallicity relation (MMR) is seen to flatten (e.g. Gallazzi et al., 2005), and we take this into account by modifying the functional form of the $p(\mathrm{M}_{\star})$-relation such that a galaxy’s average stellar mass will reproduce the turnover in the MMR (at $\sim 10^{10}\ \mathrm{M}_{\odot}$), while still matching the MDF of low mass galaxies:

[TABLE]

here $p_{0}$ describes the value the relation asymptotes towards for high galaxies masses, $\mathrm{M}_{0}$ is the turn-over mass and $\alpha_{p}$ is the low-mass slope.

In order to convert the iron abundances [Fe/H], as they were derived in Leaman et al. (2013), to total metallicity [M/H], which is used in the SSP model grid, we adopt the following relation from Salaris & Cassisi (2005)

[TABLE]

We further allow each galaxy to have a mass-dependent evolution in [$\alpha$/Fe], by utilizing the empirical [$\alpha$/Fe]-[Fe/H]-relation derived from individual stars on Local Group galaxies in de Boer et al. (2014). They found that the "knee" in the [$\alpha$/Fe]-[Fe/H] diagram ([Fe/H]knee) occurs at higher metallicity for high mass galaxies (see also Kirby et al., 2011). We adopt this mass dependence of the [Fe/H]knee position as well as of the slope of the low alpha sequence. We calibrate the [$\alpha$/Fe]-plateau value to have a mass dependence as well, such that the mean $\alpha$-abundance versus galaxy mass trends seen in observations (Thomas et al., 2010) and EAGLE (Segers et al., 2016) are reproduced.

Lastly, we account for a mass-dependent gas fraction evolution $\mu(t)$, which can be due to gas being consumed in star formation and/or being removed through feedback processes or quenching. Higher mass galaxies ($\gtrsim 10^{10}\ \mathrm{M}_{\star}$) typically exhaust their in-situ gas much quicker than lower mass galaxies (e.g. McDermid et al., 2015; Pacifici et al., 2016a; Pacifici et al., 2016b) and we therefore parameterize the gas fraction as:

[TABLE]

where $t_{\mathrm{form}}$ is an epoch by which the galaxy has formed its in-situ stars. We allow a galaxy mass dependence to enter through this formation time as:

[TABLE]

with $\alpha_{t}$ influencing how long the star formation duration will be for a given galaxy mass. This will result in AMR curves that do not evolve until $\mathrm{z}=0$, but reach their maximum metallicity at $t_{\mathrm{form}}$ for higher galaxy masses ($\gtrsim 10^{10}\ \mathrm{M}_{\star}$).

A set of these AMR templates for different galaxy masses and with parameters $p_{0}=0.1$, $M_{0}=10^{10.5}\ \mathrm{M}_{\odot}$, $\alpha_{p}=0.4$ and $\alpha_{t}=2.0$ is plotted in Figure 3a) as an example. Importantly however, these parameterizations enable us to flexibly vary the shape of the mass dependent chemical evolution tracks allowing for a stochastic assessment of the uncertainties in our final quantities of interest (see section 3.3.1).

3.3 Associating the spectroscopic mass fractions to accreted satellite galaxies

With the galaxy mass-dependent chemical templates described above, we are able to link the mass fractions in age-metallicity space recovered from the regularized pPXF solutions to astrophysical quantities of interest - such as the galaxy’s total fraction of accreted mass as well as the distribution of merged satellite galaxies. The amount of mass in accreted satellites of a given mass can be straightforwardly computed by overlaying our mass-dependent age-metallicity relation templates onto the spectroscopically recovered mass distribution in age-metallicity space, as seen in Figure 3a). Every mass weight $m_{i}$ recovered from pPXF lying below a certain AMR curve is treated as potentially coming from accreted satellite galaxies with chemical evolution representative of that mass or lower, i.e. $\mathrm{f}_{\mathrm{acc}}(\leq\mathrm{M}_{\mathrm{sat}})=\sum_{i}m_{i}(t,[\mathrm{M}/\mathrm{H}]\leq[\mathrm{M}/\mathrm{H}]_{\mathrm{template}}(t|\mathrm{M}_{\ \mathrm{sat}}))$.

This method of recovering the contributions from accreted satellite galaxies of a given mass is only considered valid in a cumulative sense, as there is an astrophysical degeneracy inherent to the assignment of mass fractions of low-mass satellites. We do not know a priori, whether these recovered stellar populations with lowest [M/H] at fixed age are coming from low-mass satellites directly accreted to the host or if they were first accreted to an intermediate mass satellite, which then merged with the host.

The chemical evolution templates extend up to an AMR associated with some most massive accreted satellite galaxy, $\mathrm{M}_{\mathrm{sat,\ max}}$. At this mass, the method has provided an estimate of the total accreted fraction for the host galaxy, i.e. $\mathrm{f}_{\mathrm{acc},\ \mathrm{tot}}=\mathrm{f}_{\mathrm{acc}}(\leq\mathrm{M}_{\mathrm{sat,\ max}})$. Every mass fraction lying above the AMR template associated with $\mathrm{M}_{\mathrm{sat,\ max}}$ is considered in-situ according to our method, but naturally there will be some overlap with the ex-situ contributions, as the age and metallicity properties are very similar in that mass regime.

In Figure 3 we show a schematic of our method to associate the spectroscopically recovered mass fractions with accreted satellite galaxies of different masses. AMR templates corresponding to accreted galaxies of masses between $10^{6}$ and $10^{11}\ \mathrm{M}_{\odot}$ are overlaid on the spectroscopically recovered mass fractions in Figure 3a). The mass fractions below an AMR curve represents the contribution to the galaxy’s merger history from satellite galaxies of this mass as seen in Figure 3b). The resultant cumulative accretion fractions versus the associated satellite galaxy masses are then produced by summing up the mass fractions lying in the shaded regions respectively and are plotted in Figure 3c). In Figure 3d) we show an analogue to the (unevolved) subhalo mass function, which can be calculated by dividing the recovered cumulative accreted mass by its associated satellite galaxy mass (i.e. $\mathrm{N}_{\mathrm{sat}}(\geq\mathrm{M}_{\mathrm{sat}})=\mathrm{M}_{\mathrm{acc}}(\leq\mathrm{M}_{\mathrm{sat}})/\mathrm{M}_{\mathrm{sat}}$, where $\mathrm{M}_{\mathrm{acc}}(\leq\mathrm{M}_{\mathrm{sat}})=\mathrm{f}_{\mathrm{acc}}(\leq\mathrm{M}_{\mathrm{sat}})\cdot\mathrm{M}_{\mathrm{host}}$).

Figure 3 shows the recovered accreted mass fractions for one realization of the chemical evolution templates ($p_{0}=0.1$, $M_{0}=10^{10.5}\ \mathrm{M}_{\odot}$, $\alpha_{p}=0.4$ and $\alpha_{t}=2.0$), however in the final results (section 4) we show them as the median of many realizations (see section 3.3.1 for more details).

3.3.1 Stochastically assessing the systematic uncertainties

As the shape of the chemical evolution templates is not well constrained by observations, especially for galaxy masses higher than found in the Local Group, we introduce flexibility and perform a Monte-Carlo (MC) simulation. We randomly draw for 1000 trials the parameters $p_{0}$, $\mathrm{M}_{0}$, $\alpha_{p}$ (equation 2) and $\alpha_{t}$ (equation 4) from a uniform distribution in a range of [-0.2,0.2], $10^{[9.5,11]}$, [0.4,0.6] and [1.7,2.5] respectively. We also add a scatter to the sampled AMR curve, which is drawn from a Gaussian distribution with a standard deviation of $\pm 0.2\ \mathrm{dex}$ in order to account for the intrinsic scatter of a galaxy’s metallicity at fixed stellar age (see e.g. Leaman et al., 2013). The choices of the parameter variation ranges were made such that the calculated MMR from the median of the mass-dependent AMR curves lie in the scatter of the observed MMR (Gallazzi et al., 2005; Kirby et al., 2013). Uncertainties in the derived [$\alpha$/Fe]-[Fe/H]-mass relation are not accounted for, as their impact is negligible in comparison to the other parameter variations. Hence, the recovered cumulative accretion fraction as a function of satellite galaxy mass is the median of the 1000 MC trials and the uncertainty is expressed by the 16th and 84th percentiles of the trials.

While this takes into account systematic uncertainties, random errors due to signal-to-noise variations in the integrated spectrum or the exact nature of the regularization calibration procedure are not considered. We however tested the variation in the shape of the recovered satellite mass function for SNR of 50, 100, 200 and 500 as well as the first, second and third order difference operator as the regularization matrix. The scatter was found to be much smaller than the uncertainties of the chemical evolution templates.

3.3.2 Defining the most massive accretion event

In principle, the cumulative accretion fraction corresponding to the highest satellite galaxy mass marks the total accretion fraction, $\mathrm{f}_{\mathrm{acc},\ \mathrm{tot}}$, of the host galaxy. However, as we do not have any prior knowledge about the most massive satellite ever accreted by the host, $\mathrm{M}_{\mathrm{sat,\ max}}$, our derived curve extends arbitrarily higher than the true total accretion fraction.

We approximate the "stopping point" by calculating the stellar mass of the accreted satellite galaxy, where the number of accreted satellite galaxies, $\mathrm{N}_{\mathrm{sat}}$, becomes unity on our spectroscopically recovered accreted satellite mass function. In return, we can truncate the cumulative accretion fractions at that point and derive the total accretion fraction of the host galaxy. Similarly, errors are computed by the intersection of the 16th and 84th percentile uncertainty of the accreted satellite mass function with unity.

This procedure works quite well even though our method can only compute a lower limit of the accreted satellite mass function, as it does not take into account any higher orders of subhalo-subhalo mergers. We found that this gives better results than using theoretical predictions of cosmological simulations, which can provide us with a statistical relation of the stellar mass of the most massive accreted satellite for a given host galaxy (see e.g. D’Souza & Bell, 2018). In addition, it allows us to implicitly characterize $\mathrm{M}_{\mathrm{sat,\ max}}$ as an intrinsic measure of our method.

4 Results

Having optimized the recovery of extended age-metallicity distributions, and formulated a way to link these to accreted galaxies of different masses, we can proceed with an example on a massive EAGLE galaxy (Galaxy Nr. 1) - one where we can independently verify the spectroscopically recovered accretion history by comparing to the known merger tree.

In the left panel of Figure 4 we show the derived cumulative accretion fraction as a function of accreted satellite galaxy mass as well as the uncertainties as calculated from the section above. We computed the cumulative accretion fraction with AMRs corresponding to galaxy masses between $10^{6.2}$ and $10^{12.4}\ \mathrm{M}_{\odot}$ in 0.2 dex steps.

The agreement with the true accretion fractions obtained from the EAGLE merger trees is remarkable given that we obtained the "observed" quantity purely from a simulated integrated spectrum. This result suggests that the assumptions behind the mass-dependent chemical evolution templates are reasonable and provide a novel way to recover signatures of otherwise unobservable ancient merger events. For this example galaxy, we recover a total accreted fraction of $\mathrm{f}_{\mathrm{acc},\ \mathrm{tot}}=0.461^{+0.104}_{-0.113}$, while the "true" total accreted fraction from the ex-situ particle classification is 0.421. For the most massive progenitor we find a stellar mass of $\log_{10}\mathrm{M}_{\mathrm{sat,\ max}}=11.50^{+0.11}_{-0.21}$ dex, whereas the actual value is 11.06 dex. The errors have been calculated by the intersections of the scatter of the MC trials of the recovered satellite mass function and where $\mathrm{N}_{\mathrm{sat}}=1$.

In the right panel of Figure 4 we compare the true satellite mass function of the EAGLE galaxy with our estimate, which also shows excellent agreement. We measured the slope ($\alpha_{\mathrm{sat}}$) of the recovered and true satellite mass function by fitting a power law between $10^{7}\ \mathrm{M}_{\odot}$ and $\mathrm{M}_{\mathrm{sat,\ max}}$ in log space. We found the values for $\alpha_{\mathrm{sat}}$ to be $-0.65^{+0.08}_{-0.17}$ and $-0.60$ respectively. The errors where estimated by calculating the slope of the extreme points of the recovered satellite mass function by taking the scatter of the MC trials as well as the uncertainty of $\mathrm{M}_{\mathrm{sat,\ max}}$ into account. The agreement of our result with the true satellite mass function is somewhat surprising as our method is formally providing a lower limit to the satellite mass function, as it cannot differentiate mergers that happened prior to a galaxy merging to the primary halo. For example, if a late time merger of high mass had its own accretion history, this would be degenerate with our solutions resulting in a flattening of the slope of the $\mathrm{f}_{\mathrm{acc}}$-$\mathrm{M}_{\mathrm{sat}}$ relation and steepening of the $\mathrm{N}_{\mathrm{sat}}$-$\mathrm{M}_{\mathrm{sat}}$ relation. However, we expect that this effect will be very small and well inside our uncertainties (see e.g. Jiang & van den Bosch, 2016, Figure 9). Nevertheless, the validation of this method using the EAGLE simulations, suggests it is a powerful new way to recover the accreted satellite mass function and distribution of merger mass ratios in galaxies.

We have also performed the same analysis on eight more EAGLE simulated galaxies, which all lie in a mass range of $10^{9}$ to $10^{12}\ \mathrm{M}_{\odot}$. This is by no means a statistical or comprehensive sample, but we wanted to illustrate that this method works on more than one handpicked galaxy. In Figure 5 and 6 we show a comparison between the true parameters of $\mathrm{f}_{\mathrm{acc,tot}}$, $\mathrm{M}_{\mathrm{sat,\ max}}$ and $\alpha_{\mathrm{sat}}$ and those retrieved from an integrated spectrum with our method. They overall lie on the 1:1 relation within the uncertainties suggesting that this method performs well on different galaxies. The mean accuracy averaged across all nine analyzed EAGLE galaxies is 12 %, 26 % and 16 % for $\mathrm{f}_{\mathrm{acc,tot}}$, $\mathrm{M}_{\mathrm{sat,\ max}}$ and $\alpha_{\mathrm{sat}}$ respectively.

5 Conclusion and Outlook

In this work we have presented and validated a new method of measuring a galaxy’s accretion history from its integrated spectrum alone. Not only are we able to quantify the total accreted fraction of a galaxy, but also measure relative amounts of accreted material coming from merged galaxies of different masses.

We hence provide, for the first time, an observational method to asses the fraction of accreted material from completely disrupted satellites as a function of their original stellar mass prior to merging. This is possible, because we exploit the fact that the chemical evolution of a galaxy varies for different masses, and hence accreted material does not lie in the same region in age-metallicity space as in-situ created material. Other studies such as Tonini (2013); Leaman et al. (2013); Kruijssen et al. (2018); Beasley et al. (2018) have also successfully exploited similar arguments to quantitatively link globular clusters to the total accreted fraction of their host. Importantly, identifying these signatures from a spectrum is only possible, if the full spectral fitting code can recover a distribution simultaneously in age and metallicity.

The advantage of this method is that it provides a more detailed understanding of the accretion history of a galaxy than existing methods (Spavone et al., 2017; Huang et al., 2016; Crnojević, 2017; Harmsen et al., 2017; D’Souza & Bell, 2018; Monachesi et al., 2018), while simultaneously being applicable to potentially every galaxy at low or high redshift, for which an integrated spectrum exists. In particular, it is able to make predictions about completely disrupted past merger events, which are otherwise indistinguishable, because they are no longer evident in phase-space in the form of streams or shells. Because this method also links chemical signatures recovered in the spectrum back to the original stellar mass of the accreted satellite galaxy prior to merging, it is so far the closest observational equivalent to a simulation’s merger tree.

In future work, we will perform more tests that will especially focus on comparing the recovery of the accretion history when only considering the light contribution of the centre or the outskirts of a galaxy. This will grant us insight as to whether there will be enough line-of-sight integration of ex-situ material in the halo, such that a central high signal-to-noise spectrum is sufficient for our method. If that is the case, this method could also be applied to higher redshift galaxies. Next, we will focus on observationally verifying this technique on targets, where there is an alternative measure of stellar population properties from resolved stars available. This is important in the sense that we can determine, whether this method is robust against difficult-to-model observational complications, as well as being able to fine tune the fitting technique. The necessary high signal-to-noise spectrum can observationally be achieved with today’s modern integral field units such as MUSE (Bacon et al., 2010).

In summary, our presented method provides a novel opportunity to validate aspects of hierarchical structure formation, as well as to make a connection to fundamental physical processes on galactic scales such as quenching or dynamical transformations. Together these will help us understand the different evolutionary pathways galaxies can experience in their mass assembly, eventually leading to the large diversity in galaxies we observe today.

Acknowledgements

RL acknowledges support from the Natural Sciences and Engineering Research Council of Canada PDF award, and by Sonderforschungsbereich SFB 881 "The Milky Way System" (subproject A7and A8) of the German Research Foundation (DFG) and DAAD PPP exchange program Projekt-ID 57316058. GvdV acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 724857 (Consolidator Grant ArcheoDyn). RAC is a Royal Society University Research Fellow. This work made use of high performance computing facilities at Liverpool John Moores University, partly funded by the Royal Society and LJMU’s Faculty of Engineering and Technology.

Appendix A EAGLE Galaxy Sample

In total we have applied this method to nine EAGLE simulated galaxies in a mass range between $10^{9}$ and $10^{12}\ \mathrm{M}_{\odot}$, as summarized in Table 1. We followed the exact same procedure as described throughout this paper. In Figure 5 we compare the parameters of the satellite mass function, $\mathrm{f}_{\mathrm{acc,\ tot}}$, $\mathrm{M}_{\mathrm{sat,\ max}}$ and $\alpha_{\mathrm{sat}}$, derived with our spectroscopic method and the truth from the EAGLE merger trees. Almost all lie on the 1:1 relation.

In addition, we derived the mass distribution in age-metallicity space for a mock spectrum with SNR of 50, 100, 200 and 500 as well as for the first, second and third order difference operator as a regularization matrix, which in total makes 12 derived satellite mass functions for every simulated galaxy. The averages of the derived parameters, $\mathrm{f}_{\mathrm{acc,\ tot}}$, $\mathrm{M}_{\mathrm{sat,\ max}}$ and $\alpha_{\mathrm{sat}}$, of those 12 variations for our galaxy sample are summarized in Figure 6. Consequently, the errorbars in this plot do not only encompass the systematic uncertainties in the chemical evolution, but also random uncertainties, when deriving the age-metallicity distribution from full spectral fitting. The mean accuracy averaged across all analyzed EAGLE galaxies and the 12 realizations is 12 %, 26 % and 16 % for $\mathrm{f}_{\mathrm{acc,tot}}$, $\mathrm{M}_{\mathrm{sat,\ max}}$ and $\alpha_{\mathrm{sat}}$ respectively.

This sample of nine EAGLE galaxies is by no means comprehensive, but it shows that this method performs well on different galaxies with different accretion histories and across a wide range of masses.