\section{Introduction}\label{intro}
\input{docs/intro}
\section{Methodology}\label{method}
\input{docs/Methodology}
\section{Results}\label{results}
\input{docs/Results}
\section{Discussion}\label{discussion}
\input{docs/Discussion}
\section{Conclusions}\label{conclusion}
\input{docs/Conclusion}
\section{Acknowledgements}
This research has made use of the services of the ESO Science Archive Facility. Based on observations collected at the European Southern Observatory under ESO programmes: 094.B-0298 (P.I. Walcher, C.J), 296.B-5054 (P.I. Sarzi, M), 0103.A-0447 \& 097.A-0366 (P.I. Hamer, S), 60.A-9303 (SV), 094.B-0321 P.I. (Marconi, A), 094.B-0711 (P.I. Arnaboldi, Magda). Some contextual panel views were based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA).  This work was completed with support from the ESO Studentship Programme, and the Moses Holden Studentship. 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). AB acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -- Project-ID 138713538 -- SFB 881 (``The Milky Way System'', subproject B08). This research made use of Astropy, \url{http://www.astropy.org} a community-developed core Python package for Astronomy \citep{astropy2013, astropy2018}.

\section{Data Availability}
The data underlying this article are available in the ESO science archive facility (\url{http://archive.eso.org/cms.html}).


\bibliographystyle{mnras}


\subsection{Observations and Analysis}
For all binning methods, a `first sweep' examined the signal-to-noise (S/N) of all spaxels in the datacube. The S/N of each spaxel was defined as the mean variance in a wavelength range of 5450\AA\, < $\lambda$ < 5550\AA\, which allowed an accurate S/N estimate free from the strongest absorption or emission lines. A mask was applied such that any single spaxel with mean S/N < 1.5 within the wavelength range was excluded from further binning. This ensured there would be no pollution of useful regions of a cube by particularly noisy pixels. Pixels were infrequently masked with the mean quantity of initially masked pixels per target being 1.4\%, with the most being masked in NGC\,4696 at 6.1\% initially masked, and some targets (e.g. NGC\,1316) requiring only 5 or 6 initially masked pixels.

Following this, the spaxel mask was built upon to include objects not associated with the galaxy, such as foreground stars or bright background objects. A Moffat function \citep{moffat1969theoretical} was applied to the galaxy at V-band which was subtracted from the original V-band slice. The resulting image was scanned using `photoutils' routine included in MPDAF \citep{baconmpdaf2016} with a sigma value defining the threshold for a source to be considered non-continuous with respect to the galaxy light. This threshold varied case by case, depending on aspects such as galaxy brightness, large scale structure, and disky areas of the galaxy. Any segmented regions above the threshold were added to the general mask of pixels to be excluded whilst binning.

\subsubsection{Voronoi Binning}
The MUSE cubes gathered from the sample presented extremely large datasets, with large variations in S/N. The on-source time of each cube varied between objects and surveys but was at a minimum of 590s per spaxel (central region of NGC\,4594) and a maximum of 11400s per spaxel (combined cubes of outer regions of NGC\,4696). To overcome issues associated with low S/N spectra we employed Voronoi binning to ensure the correct treatment of lower S/N areas, whilst maintaining the usefulness of high S/N spectra (usually located in the galaxy core). 

Voronoi binning of S/N (here specifically using the `Vorbin' package of \citealt{Cappellari2003}) gathers data on a 2D spatial plane such that bins reach a minimum summed S/N threshold. Spaxels are optimally arranged into a bin such that the spatial irregularity of a bin is minimised, and simultaneously the sum of the S/N within a bin minimally above the desired S/N threshold (thus high S/N areas are not wasted).

We use a minimum S/N threshold of 100 per \AA\, uniformly for all galaxies in the sample. Due to differences in exposure time, galaxy brightness and other data quality aspects, this results in varying numbers of Voronoi bins per unit area of a given galaxy. The brightest and deepest cubes necessitated many more Voronoi bins per unit area, as many pixels were often left unbinned. Dimmer and shallower targets consist of fewer bins, containing a greater number of binned pixels.

Spectra within a given Voronoi bin were summed, providing a single spectrum representative of the given bin. Final summed spectra were checked to ensure a minimum of S/N=100 was achieved in all bins. This does not however account for spatial co-variance between pixels. Considering the large spatial coverage of the average Voronoi bin (relative to pixel size) spatial correlations are unlikely to impact on derived values in any bin. Furthermore, all Voronoi diagrams are understood to be a general overview of galaxy properties, rather than an exhaustive pixel-perfect map.




\subsubsection{Elliptical Binning}\label{ellip_bin}
When investigating features and trends associated with galactocentric radius, it can be more instructive to bin elliptically. We use elliptical binning to explore features such as ex-situ fraction as a function of radius. This is simply done, with an origin set at the centre of a galaxy. Elliptical rings derived by projected ellipticity are overlayed onto the cube spatial plane with the major axis increasing by defined bin size. Ellipticity and radial increments are derived from 2MASS LGA data, with axes ratios derived from a J+H+K (super) image at the 3-sigma isophote \citep{jarrett20032mass}. This allows all pixels to be binned by galactocentric radius. We performed this using fractions of effective radius (as a standard using 0.1r$_e$ steps) as the increment of major axes, in order to maintain a comparable increment across all galaxies in the sample. The ratio of pixel size to elliptical bin width varies with galaxy size and apparent radius but has a minimum bin width of 11 pixels, occurring in NGC\,1332 due to the elliptical projection as well as the smaller size against the pixel scale. The mean bin pixel width is greatest for NGC\,1316 which has a minimum bin pixel width of 45 pixels.

\subsubsection{Full Spectral Fitting}\label{fsf}
We fit stellar models to binned spectra using the Penalised Pixel-Fitting (pPXF) method \citep{Cappellari2003, cappellari2016improving}. pPXF uses a maximum penalised likelihood approach to extract the stellar kinematics and stellar populations from the spectra of galaxies. Binned spectra are limited to within the wavelength range (4750\,\AA\, < $\lambda$ < 6800\,\AA). This optimally avoids wavelengths that are sensitive to the IMF (this is discussed further in Section \ref{discussion}) as well as avoiding regions particularly affected by sky lines. This wavelength range has been shown to be successful in recovering stellar populations parameters in \cite{guerou2016exploring}.

To ensure the output of pPXF population estimate is physical, we use the regularisation method built into pPXF. Regularisation is intended to smooth over the output weights to provide population estimates that are optimally physical. Weighted templates that have been combined to produce a target spectrum will often be localised to a single nonphysical solution, with many other valid solutions being overlooked, despite their physicality. To produce more representative distributions, regularisation seeks to smooth the solutions to a physical state. The challenge is to smooth the template weights to a solution that most accurately represents observed conditions, whilst not overlooking genuine fluctuations and details present in the model-fit. 

The value of regularisation is adjusted up to a maximum value such that the  reduced $\chi ^2$ of this smoothest solution is no different than 1-$\sigma$ of the unregularised solution \citep{cappellari2016improving}. This has been shown in literature to be an accurate and useful method of galaxy population extraction \cite[see e.g.][]{comeron2015, guerou2016exploring, faifer2017, Ge2019}. 

We use an iterative method to find the most reasonable value of regularisation for each bin in a galaxy cube.  First the noise from the unregularised solution is rescaled so that $\chi^{2}/N = 1$, where N is the number of voxels in the target spectrum.  Next a series of fits using different amounts of regularisation are produced. We record the output $\Delta \chi^2$ for each fit using a particular regularisation, and fit a function to the output  $\Delta \chi^2$ values with respect to the input regularisation guesses. From this function we can find the maximum regularisation parameter which corresponds to $\chi^2 = \chi^2_0+\sqrt{2N}$. We perform this procedure for all bins, obtaining optimal solutions in all cases. The value of the regularisation parameter has little variation for like spaxels of the same galaxy, but can differ between distinct galaxies. Tests were run without regularisation and similar, albeit far noisier, results were obtained.

A total of 552 SSP models, constructed from the MILES spectral library, were used to fit to galaxy spectra \citep{vazdekis2012miuscat}. These models were of Kroupa \citep{kroupa2001variation} revised initial mass function (log slope of 1.3) using BaSTI isochrones, with a metallicity range of -2.27 to +0.4 [M/H] in 12 non-linear steps, and an age range of 0.1 to 14.0\,Gyr in 46 non-linear steps \cite[][]{cassisi2005basti,pietrinferni2006large,falcon2011updated}.

We perform the pPXF routine twice for each bin, once for the robust extraction of stellar kinematics, and a second time for the extraction of stellar populations assuming the fixed kinematics derived during the initial fit. For robust kinematic fitting we use 4th degree additive polynomials. The population analysis uses no additive polynomials, instead using 16th degree multiplicative polynomials and 2nd order regularisation. This has been shown to be effective in kinematic and population analysis \citep[see e.g.][]{guerou2016exploring}.

\subsection{Ex-Situ Fraction Estimation}\label{exfrac}
The mass fractions returned by pPXF in age-metallicity space represent the contribution to the galaxy’s integrated spectrum from different stellar populations. A host galaxy’s chemical evolution will typically progress from older and more metal poor stellar populations, to younger more metal rich stellar populations. The rate of this chemical enrichment is a combination of many complex processes including outflow and recycling of gas, as well as production and removal of heavy elements through stellar evolutionary processes. Many of these processes are dependent on the star formation rate and stellar mass of the host galaxy, and this typically results in lower mass galaxies having lower metallicities (particularly at fixed age). These mass dependent age-metallicity trends have been seen in the Milky Way’s globular cluster system, and those of its satellite galaxies \citep{forbes2010accreted, leaman2013bifurcated, massari2019origin}. 

In \cite{boecker2020galaxy} the authors leveraged these mass dependent chemical evolution trends to identify which regions of the age-metallicity parameter space were likely associated with accreted satellite galaxies of different masses. In that work, the mass fractions in age-metallicity space returned by pPXF’s fit to an integrated spectra were divided into contributions to the light from galaxies of different masses using a simple mass-dependent chemical evolution model. These flexible chemical evolution tracks provided dividing lines in age-metallicity space between galaxies of different masses.

By comparing to mock spectra of EAGLE galaxies where the accretion history was known, \cite{boecker2020galaxy} showed that this approach could largely identify how much mass from low mass satellites the host galaxy had accreted. There are two key limitations of the simple \cite{boecker2020galaxy} models which we improve upon here in order to more robustly compute ex-situ fractions in real galaxies. 

(1) The mass-dependent chemical evolution templates in \cite{boecker2020galaxy} assumed a constant star formation history, simple leaky-box self enrichment, and self-similar internal spreads in metals at fixed age

(2) The association of mass fractions returned by pPXF in age-metallicity space will be non-uniquely attributed to in-situ or ex-situ components of a galaxy when the mass-ratio of a merger is close to unity (and subsequently the chemical evolution pathways are very similar).

The models for estimating ex-situ fractions from full-spectral fitting outputs used here address both of these aspects, and are presented in full in a forthcoming paper (Leaman et al. in prep.). Below we present the relevant ingredients in brief.

\subsubsection{Mass dependent SFH and chemical evolution}
The end goal of these models is to estimate the probabilistic contribution of accreted (‘ex-situ’) material to every age-metallicity bin, for a host galaxy of a given mass. While low mass accretion events ($\mu \leq 1/10$) may separate cleanly in age-metallicity space from the in-situ component of the massive host, more massive mergers will significantly overlap with the ages and metallicities of the progenitor host galaxy. Here we work to quantify this degeneracy and compute the ex-situ contribution for these galaxies when they are dominated in the limiting case by the most massive accretion event. We compute models where the most massive merger has a stellar mass ratio with respect to the present day galaxy of $\mu_{max} = 1/10, 1/5 \text{ or } 1/3$ (e.g., the mass ratios at time of accretion were $\mu_{0} = 1/9, 1/4, 1/2$).

The sole input required for computing the model is an estimate of the present day total stellar mass of the galaxy $M_{*,gal}$. The most massive accretion event is specified in terms of stellar mass as $M_{sat, max}$ = $M_{*,gal}\,\mu_{max}$. The dark matter mass of the host and most massive accreted satellite are assigned stochastically with an empirical stellar to halo mass relation (SHMR) set forth in \cite{leauthaud2012integrated}. This provides an estimate of the likely total virial mass of an observed galaxy or satellite of a given stellar mass. The $1-\sigma$ range of the quoted SHMR fit parameters in that study are used to introduce stochasticity on the adopted DM virial mass. This DM mass for the galaxy is used to specify the likely redshift of infall for the satellite given the DM mass ratio, with functional form and scatter taken from the statistics of subhalo infall in \cite{boylan2011too}. Specifically those authors show there is a mild decrease in the redshift of accretion as the mass ratio increases, which we include to set the accretion redshift when the satellite first enters the virial radius of the host. A coarse prescription for the time for the orbit to decay due to dynamical friction is adopted from Equation 5 of \cite{boylan2011too} with the redshift and mass distribution of orbital circularity taken from \cite{wetzel2010determines}, and is primarily driven by the mass ratio of the host and satellite.

The SFH of the host and most massive satellite are stochastically evolved using empirical relations for galaxy SFRs at different redshifts from \cite{genzel2015combined}. This study characterises the typical specific star formation rate (sSFR) for galaxies of a given mass and redshift. We use this empirical relation recursively starting from redshift zero with the present day mass of the host used to specify the SFR. We proceed back in time in steps of 1 Myr, with the assumption of a constant SFR within those intervals, and subtract off the amount of mass formed in that timestep (e.g., $M_{i-1} = M_{i} - (t_{i}-t_{i-1})SFR|(t_{i})$). A given SFH constructed in this way is obviously approximate, but reproduces the mass dependent evolutionary pathways of galaxies through the star forming main sequence. Each SFH is given correlated noise with a power law slope of $\beta = 2$ to mimic the bursty star formation histories seen in simulated galaxies \citep{iyer2020diversity}. The stochastic assignment of SFHs using both of these aspects typically results in SFRs that vary by $\sim 50\%$ at any given redshift for a single galaxy’s model.

The nominal mass dependent SFHs of the host and most massive satellite are further augmented by a variable quenching time, after which the SF is terminated. In the case of the host, for any given model iteration this quenching time is drawn uniformly between $t_{young} \leq t_{quench} \leq min\left \{ t_{young} + 3.5Gyr, t_{Hubble}\right \}$, where $t_{young}$ is the age of the youngest stellar population detected in the full spectral fit of the observed spectrum. For the satellite galaxy, the quenching time is specified by the time of infall to the host virial radius, plus the orbital decay time as described previously. The 3.5Gyr value acts as a buffer value encompassing uncertainty in the quenching times of galaxies, as well as accounting for uncertainties in t$_{young}$ at high redshift. These models neglect any triggered star formation in the merger, a complex process which may contribute newly formed stars with some metallicity representative of a mixture of the satellite and host gas metallicities at that epoch.


The chemical enrichment in any timestep corresponds to a leaky or accreting box prescription with equal probability. The effective yield is given a mass dependent parameterisation as in \cite{boecker2020galaxy}, with the parameters of the mass scalings chosen stochastically in each model from a range that ensures the final galaxies reproduce the global mass-metallicity relation, and its scatter, observed in SDSS and the Local Volume \citep{gallazzi2005ages, kirby2013universal}. The gas fraction time evolution, which drives any analytic chemical evolution model, is here coupled to the empirical SFHs by re-expressing this quantity in terms of a sSFR ($t_{gas} = (1 + (t_{dep}sSFR)^{-1})^{-1}$, where we adopt the redshift scaling of depletion time from \cite{tacconi2013phibss}. This allows us to self-consistently use the same relations to drive the star formation histories and chemical evolution.

\subsubsection{Model mass fractions and uncertainties}\label{mmr_uncert}
For a single merger history (characterised by the satellite mass ratio and time of accretion) we run 200 iterations of the SFH/chemical model for the host and merging satellite. Due to the stochastic sampling of those model ingredients, the resulting evolutionary tracks in age-metallicity-SFR space span a probable distribution of pathways that a galaxy of that mass might take through this parameter space. These age-metallicity tracks are weighted by the SFR, and the composite density of mass-weights in age-metallicity space are smoothed and binned to the SSP model grid used by pPXF in a metallicity range of -2.27 to +0.4 [M/H] in 12 non-linear steps, and an age range of 0.1 to 14.0 Gyr in 46 non-linear steps.

As this is done for the host and the most massive satellite galaxy which merges, for a given cell in age-metallicity space we can ask what is the plausible contribution from overlapping ex-situ/in-situ material. This grid of ex-situ fractions can be multiplied by the mass-fractions returned by pPXF from a full-spectrum fit to an observed galaxy spectrum, to produce an estimate of the contribution of ex-situ material to that galaxy’s spectrum. To quantify the uncertainty in these ex-situ fraction estimates, we run 150 merger histories (stochastically producing different quenching and accretion times) for a model corresponding to a given mass-ratio merger (an example of these stochastic histories is shown in Figure \ref{stochastic}). The variance of these 200x150 models serves as one estimate of the systematic uncertainty in this exercise. The primary uncertainties in the chemical and star formation histories are encompassed within the variations of a single merger history. For these limiting estimates we assume that any cells in the age-metallicity grid which fall below the lower envelope of the stochastic AMR tracks for the host galaxy are from material which must have come from even lower mass satellites - and hence these cells are set to values with $f_{ex-situ} = 1$.  Similarly we assume that any cells with metallicity above the upper envelope of the most massive satellite galaxy stochastic AMRs are contributed to by pure in-situ material, thus accounting for the difficulty in obtaining precise ages in the oldest age bins with any spectroscopic measurement approach. The models are importantly allowing for probabilistic recovery of ex-situ fractions in the the most important regions between these two limits where the AMR distributions
of the host and satellite overlap, which is extremely important for major mergers. This is shown for one example in Figure \ref{1407_poster}.

Final ex-situ fraction and uncertainty can be calculated from a given population as:
\begin{equation}\label{var_eq}
    f_{ex,tot\pm1\sigma}=\frac{\sum_{t,[M/H]}G_{re}\times(0\leq avg(f_{ex})\pm \sqrt{var(f_{ex})}\leq 1)}{\sum_{t,[M/H]}G_{re}}
\end{equation}
Where $G_{re}$ is the product age-metallicity stellar population map derived from full spectral fitting at a given radial annuli.

We stress that while simple, these analytic models are stochastically exploring plausible evolutionary pathways using empirical and physically motivated mass-dependent galaxy scaling relations and theoretical prescriptions. That these models reproduce additional redshift zero scaling relations they were not informed by (scatter and mean of the mass-metallicity relation, SFH-mass relations), is one simple check that they are reasonable ways to characterise the age-metallicity distributions of galaxies. The forthcoming paper (Leaman et al. in prep) will present more details on the model ingredients and calibration checks.

To ensure that assumptions regarding the metallicity gradients of galaxies were valid, extensive investigation into the possible effects of intrinsic metallicity gradients was performed. Appendix A describes this, and highlights the consistency of our method by exploring metallicity gradients with enhanced models and the EAGLE simulations.

\begin{figure*} 
	\includegraphics[width=0.9\linewidth]{images/posterplot2.pdf}
    \caption{An example process for NGC\,1407. Here we show the analysis and reduction from the collapsed cube, to the final ex-situ estimation. The leftmost panel in the top row shows a full view of the collapsed cube. Colour is inverted, and rings of 1r$_e$ and 2r$_e$ are shown projected onto the cube. The middle panel of the top row shows the stellar velocity as derived by pPXF and applied to a Voronoi binned map of NGC\,1407. The rightmost panel in the top row shows the velocity dispersion for the same map. Other kinematic and chemical properties for NGC\,1407 can be found in Figure \ref{1407_full}. The spectra at 3 different radii are shown in the second row. The r$_e$ value provided is the lower bound of an annulus with an outer boundary of radius (n+0.1)r$_e$. All spectra within these bounds are summed and shown in the panels. The black line represents the original summed spectrum, and the red line shows the final reproduced spectrum derived from full-spectral fitting. A grey shaded region represents a mask used to obscure a strong skyline ([OI]5577\AA). In the third row we show the weighted stellar models in the age-metallicity plane of NGC\,1407 for 3 different annuli in units of effective radius, as derived from the full spectral fitting fit. The left panel of the 4th row shows the ex-situ distribution obtained from the model realisation of many $\mu$=1:3 mergers at all radii, whilst the middle panel shows the resultant variance for all age-metallicity bins in the map. Variance is applied to the total age-metallicity map as per equation \ref{var_eq}. The final Voronoi map of ex-situ stars is shown in the right panel of row 4. This is shown for all galaxies in Figure \ref{ex-sit_vis}. Finally in the 5th row we show the profile of ex-situ fraction in units of r$_e$.}
    \label{1407_poster}
\end{figure*}

\begin{figure*} 
	\includegraphics[width=\linewidth]{images/Model_thread_b.pdf}
    \caption{Example empirical SFHs for one merger realisation of NGC\,1380 (merger mass ratio = 1:3). Each thread shows the path in SFR-z space of a single model, with many hundreds of stochastically varied models producing a variance map of expected z=0 population ages and metallicities. Colour shows the mean metallicity increase within the galaxy over time. The time of accretion and quenching time is stochastically generated, informed by the constraints of the youngest stellar population detected in the spectrum (3.5Gyr for NGC\,1380).}
    \label{stochastic}
\end{figure*}










\subsection*{Appendix B: Galaxy Sample - Interesting Cases and Points of Note}

The galaxy sample was sub-selected from the wider available MUSE data as described in Section \ref{method}. Figure \ref{subsel} shows this selection from available MUSE data.
\begin{figure} 
	\includegraphics[width=\linewidth]{images/sel_dist.pdf}
    \caption{Distribution of available MUSE targeted galaxies (grey dots), shown against final sample galaxies (red star) after sub-selection described in section \ref{method}. Distance for non-selected galaxies was calculated with a simplistic recessional velocity calculation assuming Hubble flow and a H0 value of 67km/s/Mpc.}
    \label{subsel}
\end{figure}
In this section we will lay out some points of interest and note regarding particular sample galaxies.
\subsubsection*{IC1459}\label{1459}
In Figure \ref{1459_full} we see a clear counter-rotating core at the centre of IC\,1459. This is well known and well studied \citep{franx1988counterrotating, cappellari2002counterrotating, prichard2019unravelling}. Counter-rotation has long been associated with galaxy mergers \citep[see e.g.][]{bertola1988counter}. The galaxy appears to show nested rings of distinct populations which are particularly visible in the M/H panel of Figure \ref{1459_full}. This feature was questioned considering the lack of kinematic features, however the ring can clearly be seen as a region of gas experiencing starformation, for instance in H-Alpha or in H-Beta, shown in Figure \ref{1459_hbeta} of the Appendix. 
\begin{figure} 
	\includegraphics[width=\linewidth]{images/HB.pdf}
    \caption{H-Beta in IC\,1459, showing a clear ring and matching to that seen distinctly in metallicity. Colour shows relative magnitude of H-Beta.}
    \label{1459_hbeta}
\end{figure}
\subsubsection*{NGC\,1316}
NGC\,1316 is a giant elliptical galaxy in the Fornax cluster, and one of the brightest radio sources present in the sky \citep{geldzahler1984radio}. The general galaxy morphology was extensively studied in \cite{schweizer1980optical, schweizer1981optical}. \cite{mackie1998evolution} use X-ray data to provide strong evidence from \cite{schweizer1981optical} that the galaxy underwent an intermediate merger, further evidenced by the presence of shells and strong tidal features, and evidence of accretion \citep{horellou2001atomic, serra2019neutral}. 
\subsubsection*{NGC\,1404}
The elliptical galaxy NGC\,1404 is located in Fornax and is one of the brighter members. NGC\,1404 is known to have historical and ongoing interaction with the Fornax central galaxy NGC\,1399, with a recent fly-by at around 1.1-1.3Gyr ago \citep{sheardown2018recent}. The interactions between NGC\,1404 and NGC\,1399 appear to have resulted in the stripping of globular clusters from NGC\,1404 which are now found to be associated with NGC\,1399 \citep{Bekki2003}. MUSE observations of the galaxy resulted in a chance nova outburst observation seen in data taken in November 2017 \citep{smith2020discovering}. In Figure \ref{1404_full} we find a double sigma spike in the centre of NGC\,1404. This matches to a small area of counter-rotation also in the centre, which is the most likely cause of this double peak in dispersion.
\subsubsection*{NGC\,2992}
NGC\,2992 is a constituent of Arp\,245, along with NGC\,2993 and the dwarf galaxy A245N. NGC\,2992 appears to be in the early stages of interaction \citep{guolo2021exploring}. IR imaging shows a very weak stellar stream connecting NGC\,2992 with NGC\,2993. We see the affects of this stream in Figures \ref{2992_full} and \ref{ex-sit_vis}. In Figure \ref{2992_full} we see younger than average populations in the lower left of the 3-6Gyr stellar age panel in the same area as the connecting stellar stream. These younger stars dilute older populations as seen in the 12-14Gyr panel. Furthermore in Figure \ref{ex-sit_vis} we see higher than average ex-situ fractions in the lower left of the panel, in the direction of the stellar stream.
\subsubsection*{NGC\,4594}
NGC\,4594, otherwise known as `The Sombrero Galaxy' is an elliptical galaxy with an inclination of 84 degrees. Its prominent disk has made it a target of both scientific and artistic value. A large tidal stream of stellar material has recently been identified in \cite{martinezdelgado2021feather}, further pointing to evidence of historical accretion. The dust lane presents particular challenges in full spectral fitting. Though pPXF has advanced tools that are well equipped to deal with dust features in galaxies, the thick dust lanes in NGC\,4594 appears to be too thick for pPXF to accurately handle. Residuals from the spectral fitting were on average 30\% higher in the thick dusty regions, compared to residuals in the bulge. This pushed many bins above the threshold residual limit and they were excluded. This also causes issues (to a far lesser degree) in Figures of NGC\,1316 and NGC\,2992.
\subsubsection*{NGC\,4696}
NGC\,4696 is the central galaxy of the Centaurus cluster of galaxies. A clear asymmetric dust feature reaches out, visually disrupting the otherwise seemingly placid galaxy. There are also less visually distinct filaments of gas and dust seen in line-emission studies \citep{fabian1982optical, laine2003hubble, sparks1989imaging}. \cite{crawford2005extended} suggest the filaments are the result of buoyant gas bubbles lifting the gas from the galaxy centre, and find areas of low X-ray pressure supporting this hypothesis. Others suggest the gas features are evidence of the accretion of a gas rich galaxy \citep{farage2010optical}.
\subsubsection*{NGC\,5846}
NGC\,5846 is a massive elliptical galaxy with a broadly spherical shape, with an axis ratio of $\sim$0.95 \citep{jarrett20032mass}. The galaxy is well known to be a near-by example of an AGN galaxy, with extensive studies in X-ray \citep{trinchieri2002peculiar, Machacek_2005}. NGC\,5846 is well studied largely due to both its proximity as an AGN galaxy, and its extensive globular cluster count \citep{forbes1996hubble}. There is strong evidence that NGC\,5846 experienced a recent flyby by its visual companion, NGC\,5846A. Indications hint at tidal stripping of gas and stars from a number of close companions in the associated galaxy group \citep{mahdavi2005ngc, zhu2016discrete}










\subsection*{Appendix A: Verification of Metallicity Effects}
\subsubsection*{A.1 Eagle Simulation Verification}
The methodology used to ascertain the ex-situ fractions of sample galaxies relied on an expectation of metallicity evolution with respect to ex-situ fractions, and largely ignored possible radial metallicity gradients attributed to internal mechanisms independent of the ex-situ fractions. The true interplay between in- and ex-situ populations that result in the final galaxy metallicity gradient are hugely complex and difficult to disentangle. We therefore present robust tests regarding this possible metallicity complication to ensure that our final results are valid irrespective of small underlying internal processes.

We use 718 galaxies from the EAGLE simulations \citep{crain2015eagle,schaye2014eagle} to assess the validity of the metallicity assumptions. Galaxies were selected to be of `red type' (more spheroidal and passively evolving), by selecting only galaxies with co-rotational values of $\kappa _{co} < 0.4$ \citep{correa2017relation}. Galaxies were further chosen to ensure a minimum of 500 stellar particles were present within 3re, and that star formation in the prior 1Gyr remained low. Each galaxy was subsequently divided into five radial bins between 0 and 5 effective radii.

For each radial bin of all galaxies, mock spectra were produced by summing weighted spectra from an AMR grid derived from stellar particles present in a galaxy radial bin. Gaussian noise was applied and spectra were convolved with Gaussian functions to replicate a spectral line width consistent with 100km/s velocity dispersion. These mock spectra were then analysed with pPXF to extract an estimated AMR. This allowed us to also ensure that degeneracy was not affecting our final results.

As it was not feasible to run simulations for each individual galaxy, ex-situ fraction was estimated using the method described in \cite{boecker2020galaxy}. As described in Section \ref{exfrac}, mass dependent chemical evolution trends are used to identify which regions of the age-metallicity parameter space are likely associated with ex-situ populations. A total of 500 Monte-Carlo simulations were produced for each bin of every galaxy, providing a robust estimate to the recovered ex-situ fraction.

This recovered estimate was then measured against the true value from the EAGLE galaxies. The mean absolute difference in ex-situ estimate varied for mass bins. In the lowest and highest mass bins, ex-situ fraction was either very low or very high respectively. This allowed easy recovery of ex-situ fraction. For intermediate masses, the AMR was more complex and showed greater variation between galaxies. At the very most, the mean difference between the recovered and the true ex-situ fraction was 9.5 percentage points. Even if all galaxies in our main MUSE sample were incorrect by 9.5 percentage points, no results would be significantly affected, and a division between the ex-situ fractions of the lowest and most massive galaxies would remain.


\begin{figure} 
	\includegraphics[width=\linewidth]{images/mean_stdv.pdf}
    \caption{Mean absolute difference of ex-situ recovery from real, shown for five stellar mass bins. A total of 718 EAGLE simulation galaxies were used, with mock spectra generated for each in 5 radial bins between 0 and 5 effective radii. The y-axis shows the value of difference between the EAGLE galaxies and the real galaxies studied in this paper. Green circles indicate the mean difference in ex-situ fraction across all radial bins for all galaxies, for 5 different stellar mass bins. Error bars indicate the standard deviation in difference within the mass bin. Points are shown at the centre of their mass bin, and cover $\pm$0.5 dex in stellar mass.}
    \label{subsel}
\end{figure}

\subsubsection*{A.2 Adaptive Metallicity Gradient Modelling}
While galaxies show a wide diversity of radial metallicity gradients, all observations to date quantify average metallicity from the \emph{total} stellar population light at any location (e.g., an unknown blend of in-situ and ex-situ stars). The interplay between the radial variation in ex-situ fraction that we recover, and any intrinsic radial metallicity gradients in the in-situ (and ex-situ) stars are therefore complex. We explore changes in recovery of ex-situ fractions when assuming intrinsic radial metallicity gradients, by adapting the models used as per Section \ref{mmr_uncert}. To examine possible effects in detail we have constructed a version of our analytic ex-situ fraction model grids which include spatial variations in the chemical evolution. We describe below how we use recent simulation and observational results to guide our choices in exploring any potential impacts of radial metallicity gradients of the in-situ and ex-situ stellar components separately.

We proceed by computing our model chemical evolution tracks for the in-situ and ex-situ stellar populations as described in Section \ref{mmr_uncert}. For a given stochastic model of many merger histories and chemical evolution histories, these produce density distributions for the ex-situ and in-situ components within the AMR grid, and these are used to decompose the returned pPXF mass fractions into ex-situ fractions. Here we simply modify this exercise to produce model density distributions of the ex-situ and in-situ components at each measured radius (0.1 R/Re steps) within the galaxy, and modify the mean metallicity of the AMR density distributions according to a particular radial metallicity gradient.

For the in-situ component we adopt a logarithmic metallicity gradient of d[M/H]/dRe = -0.2 dex/Re, which \cite{zhuang2019dynamical} showed is typical for galaxies of similar mass to NGC\,1380. For the ex-situ component we adopt three possible cases:  i) no radial metallicity gradient, ii) the same radial metallicity gradient as the in-situ component and iii) a metallicity gradient found for only the ex-situ particles in galaxies of that mass using TNG-50 simulations (\cite{nelson2019first}, Zhuang, Leaman, Pillepich, in prepaparation). 

In each case, the final model ex-situ fraction grid was used to assign the mass density in the AMR plane (derived from pPXF) an ex-situ fraction. However this was done at every radius with a locally specified ex-situ grid tailored to that radius in each of the three cases for the radial metallicity gradients described above. The attached figure shows how these new model grids (which incorporate reasonable amounts of radial dependence on metallicity informed from observations and simulations), modify the ex-situ fraction profiles for NGC\,1380.  It is evident that the changes introduced by this are important, but within the systematic uncertainties due to degeneracies, regularisation, and resolution metallicity and ages of populations. Figure \ref{comp_mods} shows the results of the adapted models, with coloured lines indicating the newly derived ex-situ profiles for NGC\,1380. All are found to be largely within the intrinsic uncertainty, as demonstrated by the grey shaded regions.

\begin{figure} 
	\includegraphics[width=\linewidth]{images/compare_models3_V.pdf}
    \caption{Ex-situ fraction profile estimations for NGC\,1380, using three different possible intrinsic metallicity gradients. The chemical evolution models were adapted according to three intrinsic metallicity cases, and their predicted ex-situ fractions were computed as described in Section \ref{mmr_uncert}. The three cases were: no radial metallicity gradient (red line); equal in-situ and ex-situ gradient magnitudes (orange line); ex-situ metallicity gradient as predicted by the TNG-50 simulations (blue line). These are overplotted on the predicted ex-situ fractions from the base models with no extra consideration of intrinsic profiles (black dotted line), with uncertainties shown with a grey hilighted region.}
    \label{comp_mods}
\end{figure}




\subsection{Observations and Analysis}
For all binning methods, a `first sweep' examined the signal-to-noise (S/N) of all spaxels in the datacube. The S/N of each spaxel was defined as the mean variance in a wavelength range of 5450\AA\, < $\lambda$ < 5550\AA\, which allowed an accurate S/N estimate free from the strongest absorption or emission lines. A mask was applied such that any single spaxel with mean S/N < 1.5 within the wavelength range was excluded from further binning. This ensured there would be no pollution of useful regions of a cube by particularly noisy pixels. Pixels were infrequently masked with the mean quantity of initially masked pixels per target being 1.4\%, with the most being masked in NGC\,4696 at 6.1\% initially masked, and some targets (e.g. NGC\,1316) requiring only 5 or 6 initially masked pixels.

Following this, the spaxel mask was built upon to include objects not associated with the galaxy, such as foreground stars or bright background objects. A Moffat function \citep{moffat1969theoretical} was applied to the galaxy at V-band which was subtracted from the original V-band slice. The resulting image was scanned using `photoutils' routine included in MPDAF \citep{baconmpdaf2016} with a sigma value defining the threshold for a source to be considered non-continuous with respect to the galaxy light. This threshold varied case by case, depending on aspects such as galaxy brightness, large scale structure, and disky areas of the galaxy. Any segmented regions above the threshold were added to the general mask of pixels to be excluded whilst binning.

\subsubsection{Voronoi Binning}
The MUSE cubes gathered from the sample presented extremely large datasets, with large variations in S/N. The on-source time of each cube varied between objects and surveys but was at a minimum of 590s per spaxel (central region of NGC\,4594) and a maximum of 11400s per spaxel (combined cubes of outer regions of NGC\,4696). To overcome issues associated with low S/N spectra we employed Voronoi binning to ensure the correct treatment of lower S/N areas, whilst maintaining the usefulness of high S/N spectra (usually located in the galaxy core). 

Voronoi binning of S/N (here specifically using the `Vorbin' package of \citealt{Cappellari2003}) gathers data on a 2D spatial plane such that bins reach a minimum summed S/N threshold. Spaxels are optimally arranged into a bin such that the spatial irregularity of a bin is minimised, and simultaneously the sum of the S/N within a bin minimally above the desired S/N threshold (thus high S/N areas are not wasted).

We use a minimum S/N threshold of 100 per \AA\, uniformly for all galaxies in the sample. Due to differences in exposure time, galaxy brightness and other data quality aspects, this results in varying numbers of Voronoi bins per unit area of a given galaxy. The brightest and deepest cubes necessitated many more Voronoi bins per unit area, as many pixels were often left unbinned. Dimmer and shallower targets consist of fewer bins, containing a greater number of binned pixels.

Spectra within a given Voronoi bin were summed, providing a single spectrum representative of the given bin. Final summed spectra were checked to ensure a minimum of S/N=100 was achieved in all bins. This does not however account for spatial co-variance between pixels. Considering the large spatial coverage of the average Voronoi bin (relative to pixel size) spatial correlations are unlikely to impact on derived values in any bin. Furthermore, all Voronoi diagrams are understood to be a general overview of galaxy properties, rather than an exhaustive pixel-perfect map.




\subsubsection{Elliptical Binning}\label{ellip_bin}
When investigating features and trends associated with galactocentric radius, it can be more instructive to bin elliptically. We use elliptical binning to explore features such as ex-situ fraction as a function of radius. This is simply done, with an origin set at the centre of a galaxy. Elliptical rings derived by projected ellipticity are overlayed onto the cube spatial plane with the major axis increasing by defined bin size. Ellipticity and radial increments are derived from 2MASS LGA data, with axes ratios derived from a J+H+K (super) image at the 3-sigma isophote \citep{jarrett20032mass}. This allows all pixels to be binned by galactocentric radius. We performed this using fractions of effective radius (as a standard using 0.1r$_e$ steps) as the increment of major axes, in order to maintain a comparable increment across all galaxies in the sample. The ratio of pixel size to elliptical bin width varies with galaxy size and apparent radius but has a minimum bin width of 11 pixels, occurring in NGC\,1332 due to the elliptical projection as well as the smaller size against the pixel scale. The mean bin pixel width is greatest for NGC\,1316 which has a minimum bin pixel width of 45 pixels.

\subsubsection{Full Spectral Fitting}\label{fsf}
We fit stellar models to binned spectra using the Penalised Pixel-Fitting (pPXF) method \citep{Cappellari2003, cappellari2016improving}. pPXF uses a maximum penalised likelihood approach to extract the stellar kinematics and stellar populations from the spectra of galaxies. Binned spectra are limited to within the wavelength range (4750\,\AA\, < $\lambda$ < 6800\,\AA). This optimally avoids wavelengths that are sensitive to the IMF (this is discussed further in Section \ref{discussion}) as well as avoiding regions particularly affected by sky lines. This wavelength range has been shown to be successful in recovering stellar populations parameters in \cite{guerou2016exploring}.

To ensure the output of pPXF population estimate is physical, we use the regularisation method built into pPXF. Regularisation is intended to smooth over the output weights to provide population estimates that are optimally physical. Weighted templates that have been combined to produce a target spectrum will often be localised to a single nonphysical solution, with many other valid solutions being overlooked, despite their physicality. To produce more representative distributions, regularisation seeks to smooth the solutions to a physical state. The challenge is to smooth the template weights to a solution that most accurately represents observed conditions, whilst not overlooking genuine fluctuations and details present in the model-fit. 

The value of regularisation is adjusted up to a maximum value such that the  reduced $\chi ^2$ of this smoothest solution is no different than 1-$\sigma$ of the unregularised solution \citep{cappellari2016improving}. This has been shown in literature to be an accurate and useful method of galaxy population extraction \cite[see e.g.][]{comeron2015, guerou2016exploring, faifer2017, Ge2019}. 

We use an iterative method to find the most reasonable value of regularisation for each bin in a galaxy cube.  First the noise from the unregularised solution is rescaled so that $\chi^{2}/N = 1$, where N is the number of voxels in the target spectrum.  Next a series of fits using different amounts of regularisation are produced. We record the output $\Delta \chi^2$ for each fit using a particular regularisation, and fit a function to the output  $\Delta \chi^2$ values with respect to the input regularisation guesses. From this function we can find the maximum regularisation parameter which corresponds to $\chi^2 = \chi^2_0+\sqrt{2N}$. We perform this procedure for all bins, obtaining optimal solutions in all cases. The value of the regularisation parameter has little variation for like spaxels of the same galaxy, but can differ between distinct galaxies. Tests were run without regularisation and similar, albeit far noisier, results were obtained.

A total of 552 SSP models, constructed from the MILES spectral library, were used to fit to galaxy spectra \citep{vazdekis2012miuscat}. These models were of Kroupa \citep{kroupa2001variation} revised initial mass function (log slope of 1.3) using BaSTI isochrones, with a metallicity range of -2.27 to +0.4 [M/H] in 12 non-linear steps, and an age range of 0.1 to 14.0\,Gyr in 46 non-linear steps \cite[][]{cassisi2005basti,pietrinferni2006large,falcon2011updated}.

We perform the pPXF routine twice for each bin, once for the robust extraction of stellar kinematics, and a second time for the extraction of stellar populations assuming the fixed kinematics derived during the initial fit. For robust kinematic fitting we use 4th degree additive polynomials. The population analysis uses no additive polynomials, instead using 16th degree multiplicative polynomials and 2nd order regularisation. This has been shown to be effective in kinematic and population analysis \citep[see e.g.][]{guerou2016exploring}.

\subsection{Ex-Situ Fraction Estimation}\label{exfrac}
The mass fractions returned by pPXF in age-metallicity space represent the contribution to the galaxy’s integrated spectrum from different stellar populations. A host galaxy’s chemical evolution will typically progress from older and more metal poor stellar populations, to younger more metal rich stellar populations. The rate of this chemical enrichment is a combination of many complex processes including outflow and recycling of gas, as well as production and removal of heavy elements through stellar evolutionary processes. Many of these processes are dependent on the star formation rate and stellar mass of the host galaxy, and this typically results in lower mass galaxies having lower metallicities (particularly at fixed age). These mass dependent age-metallicity trends have been seen in the Milky Way’s globular cluster system, and those of its satellite galaxies \citep{forbes2010accreted, leaman2013bifurcated, massari2019origin}. 

In \cite{boecker2020galaxy} the authors leveraged these mass dependent chemical evolution trends to identify which regions of the age-metallicity parameter space were likely associated with accreted satellite galaxies of different masses. In that work, the mass fractions in age-metallicity space returned by pPXF’s fit to an integrated spectra were divided into contributions to the light from galaxies of different masses using a simple mass-dependent chemical evolution model. These flexible chemical evolution tracks provided dividing lines in age-metallicity space between galaxies of different masses.

By comparing to mock spectra of EAGLE galaxies where the accretion history was known, \cite{boecker2020galaxy} showed that this approach could largely identify how much mass from low mass satellites the host galaxy had accreted. There are two key limitations of the simple \cite{boecker2020galaxy} models which we improve upon here in order to more robustly compute ex-situ fractions in real galaxies. 

(1) The mass-dependent chemical evolution templates in \cite{boecker2020galaxy} assumed a constant star formation history, simple leaky-box self enrichment, and self-similar internal spreads in metals at fixed age

(2) The association of mass fractions returned by pPXF in age-metallicity space will be non-uniquely attributed to in-situ or ex-situ components of a galaxy when the mass-ratio of a merger is close to unity (and subsequently the chemical evolution pathways are very similar).

The models for estimating ex-situ fractions from full-spectral fitting outputs used here address both of these aspects, and are presented in full in a forthcoming paper (Leaman et al. in prep.). Below we present the relevant ingredients in brief.

\subsubsection{Mass dependent SFH and chemical evolution}
The end goal of these models is to estimate the probabilistic contribution of accreted (‘ex-situ’) material to every age-metallicity bin, for a host galaxy of a given mass. While low mass accretion events ($\mu \leq 1/10$) may separate cleanly in age-metallicity space from the in-situ component of the massive host, more massive mergers will significantly overlap with the ages and metallicities of the progenitor host galaxy. Here we work to quantify this degeneracy and compute the ex-situ contribution for these galaxies when they are dominated in the limiting case by the most massive accretion event. We compute models where the most massive merger has a stellar mass ratio with respect to the present day galaxy of $\mu_{max} = 1/10, 1/5 \text{ or } 1/3$ (e.g., the mass ratios at time of accretion were $\mu_{0} = 1/9, 1/4, 1/2$).

The sole input required for computing the model is an estimate of the present day total stellar mass of the galaxy $M_{*,gal}$. The most massive accretion event is specified in terms of stellar mass as $M_{sat, max}$ = $M_{*,gal}\,\mu_{max}$. The dark matter mass of the host and most massive accreted satellite are assigned stochastically with an empirical stellar to halo mass relation (SHMR) set forth in \cite{leauthaud2012integrated}. This provides an estimate of the likely total virial mass of an observed galaxy or satellite of a given stellar mass. The $1-\sigma$ range of the quoted SHMR fit parameters in that study are used to introduce stochasticity on the adopted DM virial mass. This DM mass for the galaxy is used to specify the likely redshift of infall for the satellite given the DM mass ratio, with functional form and scatter taken from the statistics of subhalo infall in \cite{boylan2011too}. Specifically those authors show there is a mild decrease in the redshift of accretion as the mass ratio increases, which we include to set the accretion redshift when the satellite first enters the virial radius of the host. A coarse prescription for the time for the orbit to decay due to dynamical friction is adopted from Equation 5 of \cite{boylan2011too} with the redshift and mass distribution of orbital circularity taken from \cite{wetzel2010determines}, and is primarily driven by the mass ratio of the host and satellite.

The SFH of the host and most massive satellite are stochastically evolved using empirical relations for galaxy SFRs at different redshifts from \cite{genzel2015combined}. This study characterises the typical specific star formation rate (sSFR) for galaxies of a given mass and redshift. We use this empirical relation recursively starting from redshift zero with the present day mass of the host used to specify the SFR. We proceed back in time in steps of 1 Myr, with the assumption of a constant SFR within those intervals, and subtract off the amount of mass formed in that timestep (e.g., $M_{i-1} = M_{i} - (t_{i}-t_{i-1})SFR|(t_{i})$). A given SFH constructed in this way is obviously approximate, but reproduces the mass dependent evolutionary pathways of galaxies through the star forming main sequence. Each SFH is given correlated noise with a power law slope of $\beta = 2$ to mimic the bursty star formation histories seen in simulated galaxies \citep{iyer2020diversity}. The stochastic assignment of SFHs using both of these aspects typically results in SFRs that vary by $\sim 50\%$ at any given redshift for a single galaxy’s model.

The nominal mass dependent SFHs of the host and most massive satellite are further augmented by a variable quenching time, after which the SF is terminated. In the case of the host, for any given model iteration this quenching time is drawn uniformly between $t_{young} \leq t_{quench} \leq min\left \{ t_{young} + 3.5Gyr, t_{Hubble}\right \}$, where $t_{young}$ is the age of the youngest stellar population detected in the full spectral fit of the observed spectrum. For the satellite galaxy, the quenching time is specified by the time of infall to the host virial radius, plus the orbital decay time as described previously. The 3.5Gyr value acts as a buffer value encompassing uncertainty in the quenching times of galaxies, as well as accounting for uncertainties in t$_{young}$ at high redshift. These models neglect any triggered star formation in the merger, a complex process which may contribute newly formed stars with some metallicity representative of a mixture of the satellite and host gas metallicities at that epoch.


The chemical enrichment in any timestep corresponds to a leaky or accreting box prescription with equal probability. The effective yield is given a mass dependent parameterisation as in \cite{boecker2020galaxy}, with the parameters of the mass scalings chosen stochastically in each model from a range that ensures the final galaxies reproduce the global mass-metallicity relation, and its scatter, observed in SDSS and the Local Volume \citep{gallazzi2005ages, kirby2013universal}. The gas fraction time evolution, which drives any analytic chemical evolution model, is here coupled to the empirical SFHs by re-expressing this quantity in terms of a sSFR ($t_{gas} = (1 + (t_{dep}sSFR)^{-1})^{-1}$, where we adopt the redshift scaling of depletion time from \cite{tacconi2013phibss}. This allows us to self-consistently use the same relations to drive the star formation histories and chemical evolution.

\subsubsection{Model mass fractions and uncertainties}\label{mmr_uncert}
For a single merger history (characterised by the satellite mass ratio and time of accretion) we run 200 iterations of the SFH/chemical model for the host and merging satellite. Due to the stochastic sampling of those model ingredients, the resulting evolutionary tracks in age-metallicity-SFR space span a probable distribution of pathways that a galaxy of that mass might take through this parameter space. These age-metallicity tracks are weighted by the SFR, and the composite density of mass-weights in age-metallicity space are smoothed and binned to the SSP model grid used by pPXF in a metallicity range of -2.27 to +0.4 [M/H] in 12 non-linear steps, and an age range of 0.1 to 14.0 Gyr in 46 non-linear steps.

As this is done for the host and the most massive satellite galaxy which merges, for a given cell in age-metallicity space we can ask what is the plausible contribution from overlapping ex-situ/in-situ material. This grid of ex-situ fractions can be multiplied by the mass-fractions returned by pPXF from a full-spectrum fit to an observed galaxy spectrum, to produce an estimate of the contribution of ex-situ material to that galaxy’s spectrum. To quantify the uncertainty in these ex-situ fraction estimates, we run 150 merger histories (stochastically producing different quenching and accretion times) for a model corresponding to a given mass-ratio merger (an example of these stochastic histories is shown in Figure \ref{stochastic}). The variance of these 200x150 models serves as one estimate of the systematic uncertainty in this exercise. The primary uncertainties in the chemical and star formation histories are encompassed within the variations of a single merger history. For these limiting estimates we assume that any cells in the age-metallicity grid which fall below the lower envelope of the stochastic AMR tracks for the host galaxy are from material which must have come from even lower mass satellites - and hence these cells are set to values with $f_{ex-situ} = 1$.  Similarly we assume that any cells with metallicity above the upper envelope of the most massive satellite galaxy stochastic AMRs are contributed to by pure in-situ material, thus accounting for the difficulty in obtaining precise ages in the oldest age bins with any spectroscopic measurement approach. The models are importantly allowing for probabilistic recovery of ex-situ fractions in the the most important regions between these two limits where the AMR distributions
of the host and satellite overlap, which is extremely important for major mergers. This is shown for one example in Figure \ref{1407_poster}.

Final ex-situ fraction and uncertainty can be calculated from a given population as:
\begin{equation}\label{var_eq}
    f_{ex,tot\pm1\sigma}=\frac{\sum_{t,[M/H]}G_{re}\times(0\leq avg(f_{ex})\pm \sqrt{var(f_{ex})}\leq 1)}{\sum_{t,[M/H]}G_{re}}
\end{equation}
Where $G_{re}$ is the product age-metallicity stellar population map derived from full spectral fitting at a given radial annuli.

We stress that while simple, these analytic models are stochastically exploring plausible evolutionary pathways using empirical and physically motivated mass-dependent galaxy scaling relations and theoretical prescriptions. That these models reproduce additional redshift zero scaling relations they were not informed by (scatter and mean of the mass-metallicity relation, SFH-mass relations), is one simple check that they are reasonable ways to characterise the age-metallicity distributions of galaxies. The forthcoming paper (Leaman et al. in prep) will present more details on the model ingredients and calibration checks.

To ensure that assumptions regarding the metallicity gradients of galaxies were valid, extensive investigation into the possible effects of intrinsic metallicity gradients was performed. Appendix A describes this, and highlights the consistency of our method by exploring metallicity gradients with enhanced models and the EAGLE simulations.

\begin{figure*} 
	\includegraphics[width=0.9\linewidth]{images/posterplot2.pdf}
    \caption{An example process for NGC\,1407. Here we show the analysis and reduction from the collapsed cube, to the final ex-situ estimation. The leftmost panel in the top row shows a full view of the collapsed cube. Colour is inverted, and rings of 1r$_e$ and 2r$_e$ are shown projected onto the cube. The middle panel of the top row shows the stellar velocity as derived by pPXF and applied to a Voronoi binned map of NGC\,1407. The rightmost panel in the top row shows the velocity dispersion for the same map. Other kinematic and chemical properties for NGC\,1407 can be found in Figure \ref{1407_full}. The spectra at 3 different radii are shown in the second row. The r$_e$ value provided is the lower bound of an annulus with an outer boundary of radius (n+0.1)r$_e$. All spectra within these bounds are summed and shown in the panels. The black line represents the original summed spectrum, and the red line shows the final reproduced spectrum derived from full-spectral fitting. A grey shaded region represents a mask used to obscure a strong skyline ([OI]5577\AA). In the third row we show the weighted stellar models in the age-metallicity plane of NGC\,1407 for 3 different annuli in units of effective radius, as derived from the full spectral fitting fit. The left panel of the 4th row shows the ex-situ distribution obtained from the model realisation of many $\mu$=1:3 mergers at all radii, whilst the middle panel shows the resultant variance for all age-metallicity bins in the map. Variance is applied to the total age-metallicity map as per equation \ref{var_eq}. The final Voronoi map of ex-situ stars is shown in the right panel of row 4. This is shown for all galaxies in Figure \ref{ex-sit_vis}. Finally in the 5th row we show the profile of ex-situ fraction in units of r$_e$.}
    \label{1407_poster}
\end{figure*}

\begin{figure*} 
	\includegraphics[width=\linewidth]{images/Model_thread_b.pdf}
    \caption{Example empirical SFHs for one merger realisation of NGC\,1380 (merger mass ratio = 1:3). Each thread shows the path in SFR-z space of a single model, with many hundreds of stochastically varied models producing a variance map of expected z=0 population ages and metallicities. Colour shows the mean metallicity increase within the galaxy over time. The time of accretion and quenching time is stochastically generated, informed by the constraints of the youngest stellar population detected in the spectrum (3.5Gyr for NGC\,1380).}
    \label{stochastic}
\end{figure*}










\subsection*{Appendix B: Galaxy Sample - Interesting Cases and Points of Note}

The galaxy sample was sub-selected from the wider available MUSE data as described in Section \ref{method}. Figure \ref{subsel} shows this selection from available MUSE data.
\begin{figure} 
	\includegraphics[width=\linewidth]{images/sel_dist.pdf}
    \caption{Distribution of available MUSE targeted galaxies (grey dots), shown against final sample galaxies (red star) after sub-selection described in section \ref{method}. Distance for non-selected galaxies was calculated with a simplistic recessional velocity calculation assuming Hubble flow and a H0 value of 67km/s/Mpc.}
    \label{subsel}
\end{figure}
In this section we will lay out some points of interest and note regarding particular sample galaxies.
\subsubsection*{IC1459}\label{1459}
In Figure \ref{1459_full} we see a clear counter-rotating core at the centre of IC\,1459. This is well known and well studied \citep{franx1988counterrotating, cappellari2002counterrotating, prichard2019unravelling}. Counter-rotation has long been associated with galaxy mergers \citep[see e.g.][]{bertola1988counter}. The galaxy appears to show nested rings of distinct populations which are particularly visible in the M/H panel of Figure \ref{1459_full}. This feature was questioned considering the lack of kinematic features, however the ring can clearly be seen as a region of gas experiencing starformation, for instance in H-Alpha or in H-Beta, shown in Figure \ref{1459_hbeta} of the Appendix. 
\begin{figure} 
	\includegraphics[width=\linewidth]{images/HB.pdf}
    \caption{H-Beta in IC\,1459, showing a clear ring and matching to that seen distinctly in metallicity. Colour shows relative magnitude of H-Beta.}
    \label{1459_hbeta}
\end{figure}
\subsubsection*{NGC\,1316}
NGC\,1316 is a giant elliptical galaxy in the Fornax cluster, and one of the brightest radio sources present in the sky \citep{geldzahler1984radio}. The general galaxy morphology was extensively studied in \cite{schweizer1980optical, schweizer1981optical}. \cite{mackie1998evolution} use X-ray data to provide strong evidence from \cite{schweizer1981optical} that the galaxy underwent an intermediate merger, further evidenced by the presence of shells and strong tidal features, and evidence of accretion \citep{horellou2001atomic, serra2019neutral}. 
\subsubsection*{NGC\,1404}
The elliptical galaxy NGC\,1404 is located in Fornax and is one of the brighter members. NGC\,1404 is known to have historical and ongoing interaction with the Fornax central galaxy NGC\,1399, with a recent fly-by at around 1.1-1.3Gyr ago \citep{sheardown2018recent}. The interactions between NGC\,1404 and NGC\,1399 appear to have resulted in the stripping of globular clusters from NGC\,1404 which are now found to be associated with NGC\,1399 \citep{Bekki2003}. MUSE observations of the galaxy resulted in a chance nova outburst observation seen in data taken in November 2017 \citep{smith2020discovering}. In Figure \ref{1404_full} we find a double sigma spike in the centre of NGC\,1404. This matches to a small area of counter-rotation also in the centre, which is the most likely cause of this double peak in dispersion.
\subsubsection*{NGC\,2992}
NGC\,2992 is a constituent of Arp\,245, along with NGC\,2993 and the dwarf galaxy A245N. NGC\,2992 appears to be in the early stages of interaction \citep{guolo2021exploring}. IR imaging shows a very weak stellar stream connecting NGC\,2992 with NGC\,2993. We see the affects of this stream in Figures \ref{2992_full} and \ref{ex-sit_vis}. In Figure \ref{2992_full} we see younger than average populations in the lower left of the 3-6Gyr stellar age panel in the same area as the connecting stellar stream. These younger stars dilute older populations as seen in the 12-14Gyr panel. Furthermore in Figure \ref{ex-sit_vis} we see higher than average ex-situ fractions in the lower left of the panel, in the direction of the stellar stream.
\subsubsection*{NGC\,4594}
NGC\,4594, otherwise known as `The Sombrero Galaxy' is an elliptical galaxy with an inclination of 84 degrees. Its prominent disk has made it a target of both scientific and artistic value. A large tidal stream of stellar material has recently been identified in \cite{martinezdelgado2021feather}, further pointing to evidence of historical accretion. The dust lane presents particular challenges in full spectral fitting. Though pPXF has advanced tools that are well equipped to deal with dust features in galaxies, the thick dust lanes in NGC\,4594 appears to be too thick for pPXF to accurately handle. Residuals from the spectral fitting were on average 30\% higher in the thick dusty regions, compared to residuals in the bulge. This pushed many bins above the threshold residual limit and they were excluded. This also causes issues (to a far lesser degree) in Figures of NGC\,1316 and NGC\,2992.
\subsubsection*{NGC\,4696}
NGC\,4696 is the central galaxy of the Centaurus cluster of galaxies. A clear asymmetric dust feature reaches out, visually disrupting the otherwise seemingly placid galaxy. There are also less visually distinct filaments of gas and dust seen in line-emission studies \citep{fabian1982optical, laine2003hubble, sparks1989imaging}. \cite{crawford2005extended} suggest the filaments are the result of buoyant gas bubbles lifting the gas from the galaxy centre, and find areas of low X-ray pressure supporting this hypothesis. Others suggest the gas features are evidence of the accretion of a gas rich galaxy \citep{farage2010optical}.
\subsubsection*{NGC\,5846}
NGC\,5846 is a massive elliptical galaxy with a broadly spherical shape, with an axis ratio of $\sim$0.95 \citep{jarrett20032mass}. The galaxy is well known to be a near-by example of an AGN galaxy, with extensive studies in X-ray \citep{trinchieri2002peculiar, Machacek_2005}. NGC\,5846 is well studied largely due to both its proximity as an AGN galaxy, and its extensive globular cluster count \citep{forbes1996hubble}. There is strong evidence that NGC\,5846 experienced a recent flyby by its visual companion, NGC\,5846A. Indications hint at tidal stripping of gas and stars from a number of close companions in the associated galaxy group \citep{mahdavi2005ngc, zhu2016discrete}










\subsection*{Appendix A: Verification of Metallicity Effects}
\subsubsection*{A.1 Eagle Simulation Verification}
The methodology used to ascertain the ex-situ fractions of sample galaxies relied on an expectation of metallicity evolution with respect to ex-situ fractions, and largely ignored possible radial metallicity gradients attributed to internal mechanisms independent of the ex-situ fractions. The true interplay between in- and ex-situ populations that result in the final galaxy metallicity gradient are hugely complex and difficult to disentangle. We therefore present robust tests regarding this possible metallicity complication to ensure that our final results are valid irrespective of small underlying internal processes.

We use 718 galaxies from the EAGLE simulations \citep{crain2015eagle,schaye2014eagle} to assess the validity of the metallicity assumptions. Galaxies were selected to be of `red type' (more spheroidal and passively evolving), by selecting only galaxies with co-rotational values of $\kappa _{co} < 0.4$ \citep{correa2017relation}. Galaxies were further chosen to ensure a minimum of 500 stellar particles were present within 3re, and that star formation in the prior 1Gyr remained low. Each galaxy was subsequently divided into five radial bins between 0 and 5 effective radii.

For each radial bin of all galaxies, mock spectra were produced by summing weighted spectra from an AMR grid derived from stellar particles present in a galaxy radial bin. Gaussian noise was applied and spectra were convolved with Gaussian functions to replicate a spectral line width consistent with 100km/s velocity dispersion. These mock spectra were then analysed with pPXF to extract an estimated AMR. This allowed us to also ensure that degeneracy was not affecting our final results.

As it was not feasible to run simulations for each individual galaxy, ex-situ fraction was estimated using the method described in \cite{boecker2020galaxy}. As described in Section \ref{exfrac}, mass dependent chemical evolution trends are used to identify which regions of the age-metallicity parameter space are likely associated with ex-situ populations. A total of 500 Monte-Carlo simulations were produced for each bin of every galaxy, providing a robust estimate to the recovered ex-situ fraction.

This recovered estimate was then measured against the true value from the EAGLE galaxies. The mean absolute difference in ex-situ estimate varied for mass bins. In the lowest and highest mass bins, ex-situ fraction was either very low or very high respectively. This allowed easy recovery of ex-situ fraction. For intermediate masses, the AMR was more complex and showed greater variation between galaxies. At the very most, the mean difference between the recovered and the true ex-situ fraction was 9.5 percentage points. Even if all galaxies in our main MUSE sample were incorrect by 9.5 percentage points, no results would be significantly affected, and a division between the ex-situ fractions of the lowest and most massive galaxies would remain.


\begin{figure} 
	\includegraphics[width=\linewidth]{images/mean_stdv.pdf}
    \caption{Mean absolute difference of ex-situ recovery from real, shown for five stellar mass bins. A total of 718 EAGLE simulation galaxies were used, with mock spectra generated for each in 5 radial bins between 0 and 5 effective radii. The y-axis shows the value of difference between the EAGLE galaxies and the real galaxies studied in this paper. Green circles indicate the mean difference in ex-situ fraction across all radial bins for all galaxies, for 5 different stellar mass bins. Error bars indicate the standard deviation in difference within the mass bin. Points are shown at the centre of their mass bin, and cover $\pm$0.5 dex in stellar mass.}
    \label{subsel}
\end{figure}

\subsubsection*{A.2 Adaptive Metallicity Gradient Modelling}
While galaxies show a wide diversity of radial metallicity gradients, all observations to date quantify average metallicity from the \emph{total} stellar population light at any location (e.g., an unknown blend of in-situ and ex-situ stars). The interplay between the radial variation in ex-situ fraction that we recover, and any intrinsic radial metallicity gradients in the in-situ (and ex-situ) stars are therefore complex. We explore changes in recovery of ex-situ fractions when assuming intrinsic radial metallicity gradients, by adapting the models used as per Section \ref{mmr_uncert}. To examine possible effects in detail we have constructed a version of our analytic ex-situ fraction model grids which include spatial variations in the chemical evolution. We describe below how we use recent simulation and observational results to guide our choices in exploring any potential impacts of radial metallicity gradients of the in-situ and ex-situ stellar components separately.

We proceed by computing our model chemical evolution tracks for the in-situ and ex-situ stellar populations as described in Section \ref{mmr_uncert}. For a given stochastic model of many merger histories and chemical evolution histories, these produce density distributions for the ex-situ and in-situ components within the AMR grid, and these are used to decompose the returned pPXF mass fractions into ex-situ fractions. Here we simply modify this exercise to produce model density distributions of the ex-situ and in-situ components at each measured radius (0.1 R/Re steps) within the galaxy, and modify the mean metallicity of the AMR density distributions according to a particular radial metallicity gradient.

For the in-situ component we adopt a logarithmic metallicity gradient of d[M/H]/dRe = -0.2 dex/Re, which \cite{zhuang2019dynamical} showed is typical for galaxies of similar mass to NGC\,1380. For the ex-situ component we adopt three possible cases:  i) no radial metallicity gradient, ii) the same radial metallicity gradient as the in-situ component and iii) a metallicity gradient found for only the ex-situ particles in galaxies of that mass using TNG-50 simulations (\cite{nelson2019first}, Zhuang, Leaman, Pillepich, in prepaparation). 

In each case, the final model ex-situ fraction grid was used to assign the mass density in the AMR plane (derived from pPXF) an ex-situ fraction. However this was done at every radius with a locally specified ex-situ grid tailored to that radius in each of the three cases for the radial metallicity gradients described above. The attached figure shows how these new model grids (which incorporate reasonable amounts of radial dependence on metallicity informed from observations and simulations), modify the ex-situ fraction profiles for NGC\,1380.  It is evident that the changes introduced by this are important, but within the systematic uncertainties due to degeneracies, regularisation, and resolution metallicity and ages of populations. Figure \ref{comp_mods} shows the results of the adapted models, with coloured lines indicating the newly derived ex-situ profiles for NGC\,1380. All are found to be largely within the intrinsic uncertainty, as demonstrated by the grey shaded regions.

\begin{figure} 
	\includegraphics[width=\linewidth]{images/compare_models3_V.pdf}
    \caption{Ex-situ fraction profile estimations for NGC\,1380, using three different possible intrinsic metallicity gradients. The chemical evolution models were adapted according to three intrinsic metallicity cases, and their predicted ex-situ fractions were computed as described in Section \ref{mmr_uncert}. The three cases were: no radial metallicity gradient (red line); equal in-situ and ex-situ gradient magnitudes (orange line); ex-situ metallicity gradient as predicted by the TNG-50 simulations (blue line). These are overplotted on the predicted ex-situ fractions from the base models with no extra consideration of intrinsic profiles (black dotted line), with uncertainties shown with a grey hilighted region.}
    \label{comp_mods}
\end{figure}



\section{Introduction}\label{intro}
\input{docs/intro}
\section{Methodology}\label{method}
\input{docs/Methodology}
\section{Results}\label{results}
\input{docs/Results}
\section{Discussion}\label{discussion}
\input{docs/Discussion}
\section{Conclusions}\label{conclusion}
\input{docs/Conclusion}
\section{Acknowledgements}
This research has made use of the services of the ESO Science Archive Facility. Based on observations collected at the European Southern Observatory under ESO programmes: 094.B-0298 (P.I. Walcher, C.J), 296.B-5054 (P.I. Sarzi, M), 0103.A-0447 \& 097.A-0366 (P.I. Hamer, S), 60.A-9303 (SV), 094.B-0321 P.I. (Marconi, A), 094.B-0711 (P.I. Arnaboldi, Magda). Some contextual panel views were based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA).  This work was completed with support from the ESO Studentship Programme, and the Moses Holden Studentship. 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). AB acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -- Project-ID 138713538 -- SFB 881 (``The Milky Way System'', subproject B08). This research made use of Astropy, \url{http://www.astropy.org} a community-developed core Python package for Astronomy \citep{astropy2013, astropy2018}.

\section{Data Availability}
The data underlying this article are available in the ESO science archive facility (\url{http://archive.eso.org/cms.html}).


\bibliographystyle{mnras}
