Differences in the properties of disrupted and surviving satellites of Milky-Way-mass galaxies in relation to their host accretion histories

From the chemo-dynamical properties of tidal debris in the Milky Way, it has been inferred that the dwarf satellites that have been disrupted had different chemical abundances from their present-day counterparts of similar mass that survive today, specifically, they had lower [Fe/H] and higher [Mg/Fe]. Here we use the ARTEMIS simulations to study the relation between the chemical abundances of disrupted progenitors of MW-mass galaxies and their stellar mass, and the evolution of the stellar mass - metallicity relations (MZR) of this population with redshift. We find that these relations have significant scatter, which correlates with the accretion redshifts ($z_{\rm acc}$) of satellites, and with their cold gas fractions. We investigate the MZRs of dwarf populations accreted at different redshifts and find that they have similar slopes, and also similar with the slope of the MZR of the surviving population ($\approx 0.32$). However, the entire population of disrupted dwarfs displays a steeper MZR, with a slope of $\approx 0.48$, which can be explained by the changes in the mass spectrum of accreted dwarf galaxies with redshift. We find strong relations between the (mass-weighted) $\langle z_{\rm acc} \rangle$ of the disrupted populations and their global chemical abundances ($\langle$[Fe/H]$\rangle$ and $\langle$[Mg/Fe]$\rangle$), which suggests that chemical diagnostics of disrupted dwarfs can be used to infer the types of merger histories of their hosts. For the case of the MW, our simulations predict that the bulk of the disrupted population was accreted at $\langle z_{\rm acc} \rangle \approx 2$, in agreement with other findings. We also find that disrupted satellites form and evolve in denser environments, closer to their hosts, than their present-day counterparts.


INTRODUCTION
In the standard ΛCDM cosmological model, large galaxies like the Milky Way (MW) grow hierarchically, by accreting many smaller galaxies over time (White & Frenk 1991).In MW-like galaxies, signatures of accreted and disrupted dwarf galaxies can be recovered from the properties of the tidal debris found at present-day in the stellar halo (Bullock & Johnston 2005;Cooper et al. 2010;Monachesi et al. 2019;Helmi 2020).Many observations support this hierarchical paradigm of galaxy formation, numerous tidal streams having being already discovered in the MW (e.g., Helmi et al. 1999;Newberg et al. 2003;Majewski et al. 2004;Belokurov et al. 2007;Myeong et al. 2019;Li et al. 2022;Naidu et al. 2020;Yuan et al. 2020), in M31 (e.g., Martin et al. 2014), and in more distant spiral galaxies (e.g., Martínez-Delgado et al. 2010).
★ Corresponding author:A.S.Font@ljmu.ac.uk Combining kinematical and chemical abundance data for stars in the tidal debris is deemed to be a powerful tool for pinning down the properties the progenitors, and hence, to reconstruct the merger history of the host galaxy (McWilliam 1997;Freeman & Bland-Hawthorn 2002).Previous numerical simulations have shown the validity of this approach and were able to put some constraints on the merger history of the Milky Way (e.g., Font et al. 2006b;Johnston et al. 2008;Gilbert et al. 2009;Lee et al. 2015;Cunningham et al. 2022;Deason et al. 2023).However, to date, the full reconstruction of MW's merger history is limited by the availability of the observational data.Reconstructing the full merger history will also help us place our Galaxy within the larger cosmological context, and determine whether this history is "typical" among galaxies of similar mass (e.g., De Rossi et al. 2009).
For the Milky Way there are now exquisite kinematical and chemical abundance data for many stellar tidal streams, thanks to precise astrometric measurements from Gaia (Gaia Collaboration et al. 2018) and high-resolution spectroscopy from the APOGEE (Majewski et al. 2017), GALAH (De Silva et al. 2015) and H3 (Conroy et al. 2019) surveys.This wealth of data has enabled the reconstruction of many of the past accretion events.For example, it is now well established that a massive dwarf galaxy, called Gaia Enceladus/Sausage (GES), was accreted by the MW about 8−9 Gyr ago (Belokurov et al. 2018;Helmi et al. 2018).It is estimated that, at the time of accretion, GES had a stellar mass similar to that of the Large Magellanic Cloud (LMC) today.Recently, streams from other, less massive, disrupted satellites have been uncovered by observations.These include the Kraken (Kruĳssen et al. 2019), Sequoia (Myeong et al. 2019), Thamnos (Koppelman et al. 2019), Wukong and I'toi (Naidu et al. 2020).However, although the main branches of the MW's "merger tree" have been identified, many other tidal streams may still remain hidden.Cosmological models predict that the disrupted dwarf galaxies make up the vast majority (∼ 70−80%) of all the satellites ever accreted onto the MW (Fattahi et al. 2020;Santistevan et al. 2020).Given that the current census of surviving satellites is > ∼ 50 (Drlica-Wagner et al. 2020), this suggests that the merger history of the MW is still to be determined.
With a growing number of tidal streams from disrupted satellites, we begin to have a better picture about the properties of these early accreted dwarf galaxies, and how these compare with the properties of present-day satellites.The emerging picture is that disrupted dwarf galaxies (i.e., the early progenitors, or the MW "building blocks") had different evolutionary paths than their present-day counterparts.For example, by using data from Gaia and the H3 survey, Naidu et al. (2020) were able to assign a good portion of the stellar halo to individual tidal streams from disrupted galaxies.With these data, Naidu et al. (2022) found that at fixed stellar mass, disrupted satellites have lower metallicities than the surviving ones, ([Fe/H] lower by ≃ 0.4 dex).They also found that, at fixed [Fe/H], the disrupted satellites are more -enriched than present-day satellites ([Mg/Fe] ratios higher by ∼ 0.2 − 0.3 dex).
These results are consistent with the well known chemical abundance differences between the stars in the halo of the MW and stars in present-day satellites.At fixed [Fe/H], the stellar halo shows higher [/Fe] values than those in surviving dwarf galaxies, typically by 0.2 dex (Fulbright 2002;Shetrone et al. 2003;Venn et al. 2004;Tolstoy et al. 2009).This suggests that the surviving population is unlike the MW's building blocks and that it has followed different star formation and chemical evolution paths than the disrupted population.Simulations run in ΛCDM models are able to explain these differences, showing that the early accreted dwarfs tend to have more bursty star formation histories than the later accreted ones (Robertson et al. 2005;Khoperskov et al. 2023).Consequently, the accreted component of the MW stellar haloes in these simulations have higher [/Fe] than surviving satellites which, having been accreted later, had more time to form stars and increase their Fe content, and therefore have lower [/Fe] abundances (see, for example, Font et al. 2006a).
The observed differences in [Fe/H] between the disrupted and surviving populations prompt the question whether they also display different stellar mass -metallicity relations (MZR).The lower [Fe/H] at fixed stellar mass displayed by disrupted satellites implies that they obey a MZR with a lower normalisation than corresponding relation for surviving satellites.Determining the offset in normalizations between the MZRs of the two populations, and also, if there are any differences in the slopes of these relations, allows us to understand further the intrinsic differences between the two populations, and even to constrain the accretion history of the MW.
The MZR of surviving population is well determined.Present-day dwarf galaxies follow the same global scaling relation for other, more massive, galaxies (Grebel et al. 2003;Kirby et al. 2013Kirby et al. , 2020)), while the relation is also observed to continue down to the mass regime of the ultra-faint dwarfs, with the same slope and normalisation as measured for the classical ones (Simon 2019).Moreover, the MZR has been observed for galaxies outside the Local Group (e.g., Gallazzi et al. 2005) and at also at higher redshifts (e.g., Cullen et al. 2019).State-of-the-art cosmological simulations obtain that the present-day stellar MZR extends over several decades in galaxy mass (e.g., Schaye et al. 2015;De Rossi et al. 2015, 2017).
In contrast, the MZR for the disrupted population of MW satellites was characterized only recently.By reconstructing the properties of disrupted dwarfs from the information encoded in the tidal debris, Naidu et al. (2022) found that their MZR has a similar slope to that of the MZR of surviving population (∼ 0.3).The MZR of the disrupted population is offset lower in [Fe/H], by ≈ 0.4, compared with the MZR at present-day from Kirby et al. (2013).As previously mentioned, Naidu et al. (2022) also find that disrupted satellites are more enriched in [Mg/Fe] than surviving dwarfs of, at either fixed stellar mass, or fixed [Fe/H].In addition, they find that data supports a scenario in which the disrupted dwarfs formed closer to the MW host.The higher density environment in which progenitors evolved, coupled with the proximity to the MW, are factors that may have led to the more enhanced star formation histories of these early accreted satellites compared with the surviving population.
In this study, we use the ARTEMIS suite of high-resolution cosmological hydrodynamical simulations to compare the properties of the surviving and disrupted satellites of MW-mass galaxies and test this scenario.In addition, we investigate each population in more detail, for example by analysing the scatter in various properties at  = 0 (e.g., in the [Fe/H] - * and [Mg/Fe] - * relations, in the MZRs, and in the [Mg/Fe] -[Fe/H] distributions).We also study the dependence on redshift of the MZR of disrupted population and compare the results from simulations with available observations in the MW.
The paper is organised as follows.Section 2 summarizes the main characteristics of the suite of MW-mass simulations used in this study.In Section 3, we analyse the (stellar) MZRs of the surviving and disrupted dwarf populations ( § 3.1) and the corresponding [Mg/Fe] - * and [Mg/Fe] -[Fe/H] distributions ( § 3.2); we investigate the nature of the scatter in these relations ( § 3.3); and seek to find which properties of the host progenitors (namely, the average chemical abundances, the slope of the MZR, the slope of the [Mg/Fe] - * relation) correlate with the accretion histories of their hosts ( § 3.4); we compare the stellar mass distributions of surviving and disrupted dwarfs ( § 3.5); we also follow the locations in which the disrupted and surviving satellite populations formed and evolved prior to accretion onto their hosts ( § 3.6), in order to determine the role that environment plays in their properties; and analyse the disruption times of the disrupted progenitors after accretion onto their hosts ( § 3.7).Section 4 includes a discussion of our results, and Section 5 presents the conclusions of this study.
Unless otherwise specified, all chemical abundances correspond to the stellar components of dwarf satellites, while MZR refers to the stellar mass -stellar metallicity relation.

ARTEMIS SIMULATIONS AND MERGER TREES
ARTEMIS is a suite of zoomed-in, cosmological hydrodynamical simulations of MW-mass systems run in a flat ΛCDM WMAP9 (Hinshaw et al. 2013) model.The following cosmological pa-rameters are assumed: Ω m = 0.2793, Ω b = 0.0463, ℎ = 0.70,  8 = 0.8211 and   = 0.972.Here, we use the sample of 42 MW-mass hosts simulations introduced in Font et al. (2020), with total masses ranging between 8 × 10 11 <  200 /M ⊙ < 2 × 10 12 , where  200 is the mass enclosing a mean density of 200 times the critical density of the Universe at present time.Dark matter particle masses are 1.17×10 5 M ⊙ ℎ −1 and the initial baryon particle masses are 2.23 × 10 4 M ⊙ ℎ −1 , while the force resolution (the Plummerequivalent softening) is 125 pc ℎ −1 .Below we summarise the main properties of the simulations.

Subgrid physical prescriptions
ARTEMIS simulations were carried out with the Gadget-3 code (last described in Springel et al. 2005), with an updated hydro solver and subgrid physical prescriptions developed for the EAGLE project (Schaye et al. 2015).The latter include metal-dependent radiative cooling in the presence of a photo-ionizing UV background, star formation, stellar and chemical evolution, formation of supermassive black holes, stellar feedback from supernova (SN) and stellar winds, as well as feedback from active galactic nuclei.These prescriptions are described in detail in Schaye et al. (2015) and Crain et al. (2015) (see, also, references therein).In particular, the ARTEMIS chemical enrichment model considers the yield from AGB stars, stellar winds and core collapse SN. 11 element species are followed: H, He, C, N, O, Ne, Mg, Si and Fe are tracked individually, and Ca and S track the Si abundance.The yield depends on the amount of mass reaching the end of the main sequence phase in each star particle.Nucleosynthetic yields from Marigo (2001) and Portinari et al. (1998) are applied.The yield from SNIa is implemented separately, where the rate of SN events per unit initial mass is defined as  SNIa =  −1  −/ , with  = 2 Gyr and  = 2 × 10 −3 M ⊙ −1 (see Wiersma et al. 2009 andSchaye et al. 2015).The SNIa delay time considered in this model takes into account the minimum time that is required for stars to evolve into white dwarfs 1 The SNIa yields are from Thielemann et al. (2003).
The stellar feedback scheme used in ARTEMIS is similar to that implemented in the EAGLE simulations (Schaye et al. 2015), however its associated parameters were re-calibrated (see Font et al. 2020) to obtain a better match to the stellar mass -halo mass relation for galaxies in the MW-mass range.As a result, the host galaxies in the ARTEMIS sample encompass the range of stellar masses and magnitudes of a large number of "MW analogues" from the Local Group (i.e., M31), the Local Volume and in the SAGA survey (10 -40 Mpc;Geha et al. 2017).For further details, see Font et al. (2021).The global properties of surviving satellites around MW analogues are also matched remarkably well (Font et al. 2021(Font et al. , 2022)).These include: the overall number and the radial distribution around their hosts, their luminosity functions, colours, current star formation rates and quenched fractions.The simulations match these observables in the surviving satellites in the MW, M31 and in other MW analogues (e.g., Bennet et al. 2019;Carlsten et al. 2021a,b).
In the simulations, individual galaxies (MW-mass hosts and dwarf satellites) are identified as gravitationally bound structures could model the rate of SNIa (see also Maoz et al. 2010).We note that changes to the SNIa delay time in the model can lead to variations in chemical properties of galaxies, such as their -enhancement.An analysis of this effect is beyond the scope of this paper.

Chemical abundances offsets
As discussed in Font et al. (2020), these simulations produce satellite galaxies that are more metal-rich (in both [Fe/H]) and [Mg/Fe]) than observed ones (see figure 2 of that study).A similar, albeit slightly higher, discrepancy in metallicities is seen in the EAGLE simulations (see figure 13 in Schaye et al. 2015, and also the discussion in De Rossi et al. 2017).Upon closer inspection in ARTEMIS, we find that the chemical abundance values, both [Fe/H] and [Mg/Fe], are systematically higher than the observed values across several orders in stellar mass ( * ≃ 10 6 − 10 10 M ⊙ ).This results in the MZR of surviving satellites being offset systematically upwards in [Fe/H] by about 0.4 dex, compared with the observed MZR (Kirby et al. 2013) in the MW.Similarly, the averaged [Mg/Fe] values for the simulated dwarfs at  = 0 are higher by ≃ 0.2 dex compared to observations.The cause of these discrepancies between the observed and simulated chemical abundances is not clear, however, given these are systematic differences across about the four orders of magnitude in stellar mass investigated here, it suggests that the cause may be related to the particular nucleosynthetic yields adopted in these simulations.As discussed in detail in Wiersma et al. (2009), nucleosynthetic yields are uncertain by factors of a few and the evolution of simulated nuclei abundance is sensitive to the specific choice of yield tables.On the other hand, the stellar masses and metallicities inferred from observations may also be affected by systematic uncertainties associated, for example, with the aperture of instruments, selection effects, or the different methods applied for deriving element abundances from the spectral energy distributions (SEDs) of galaxies (e.g., Maiolino & Mannucci 2019).
Because we are mainly interested in the differences between the chemical abundances of surviving and disrupted dwarf galaxies at any given redshift, and since these abundances are affected in a similar manner by the nucleosynthetic yields, the main results of this study are not affected by these systematic differences.To simplify the comparisons with the observations, we choose to subtract a fixed value of 0.4 dex from all [Fe/H] values of individual star particles in the simulations, in order to match the  = 0 MZR of observed satellites in the MW.Similarly, we substract a fixed value of 0.2 dex from all [Mg/Fe] values in the simulations.Hereafter, all [Fe/H] and [Mg/Fe] values include these offsets.As the same shifts are applied at all redshifts (), this procedure does not affect the study of chemical abundance variations across cosmic time, or the scatter, which are the main focus of this paper.

Tracking dwarf progenitors along merger trees
Merger trees for the MW-mass hosts were constructed previously, using the same methods described in McAlpine et al. (2016) for the EAGLE simulations.The 42 MW-mass systems used here have a variety of merger histories, ranging from quiescent to very active.We use these merger trees to follow the dwarf satellite galaxies back in time, prior to accretion onto their hosts, and to determine their properties over time.From the merger trees, we identify satellite dwarf galaxies as subhaloes that merge with the main branches (the main haloes, i.e., the main progenitors of the MW-mass galaxies).We define as surviving all dwarf galaxies that lie at present time ( = 0) within the virial radius of their host, namely within  200 , which is the radius where the overdensity of the main halo is 200 times the critical density of the Universe.All other dwarfs that were accreted onto their hosts, but they could not be identified anymore by SUBFIND, are considered to be disrupted.Note that SUBFIND imposes a minimum number of total particles (dark matter + baryonic), in this case assumed to be 20, in order to identify a self-bound subhalo.
Throughout the paper, we refer to two specific redshifts, which are relevant in the evolution of dwarf satellites: • For disrupted dwarfs, we refer to their redshift of accretion,  acc .Typically,  acc is defined as the redshift when the subhalo joins the FoF group of the main halo progenitor.Here, however, we define  acc as the redshift when the dwarf galaxy (subhalo) reaches its maximum  * .In general, the two definitions give similar results, as seen in Fig. 1.The almost 1:1 correspondence in this figure indicates that most dwarfs reach their maximum mass around the time they cross the  200 radius of their host (we denote  acc|  200 as the corresponding redshift), after which their star formation quickly becomes quenched.Also, many satellites begin to be tidally disrupted after they enter their host system and experience mass loss.Therefore, for most dwarf satellites, the redshift when they reach the maximum  * is also when they cross the virial radius ( acc ≈  acc|  200 ).However, for a small fraction of dwarf galaxies, this is not the case, as they continue to form stars after crossing  200 and then become disrupted.Fig. 1 shows that this occurs mostly for high mass satellites (systems for which  acc <  acc|  200 ).Since in this study we are interested in the intrinsic differences in chemical abundances between different populations of dwarf galaxies, we aim to capture the entire chemical evolution of these systems.This is more appropriately represented by the redshift when galaxies reach their maximum  * , i.e.,  acc .
• For surviving satellites, we also compute the redshift when the dwarfs reach half of their maximum  * , denoted as  50 (specifically, when  * is within 40 to 60 per cent of the maximum  * ).Here we use  50 instead of  acc because we aim to determine whether the scatter in chemical abundances of satellites depends on the accre-tion history of the host.But, as surviving dwarfs are predominantly accreted recently (usually, at  ≲ 1),  acc does not have sufficient range to correlate with any scatter in the properties of this population.Typically  50 >  acc , given that satellites grow most of their  * prior to accretion onto their host.We note that for a small number of surviving satellites (45 out of a total of 551; i.e., ≈ 8 per cent) we could not identify  50 , due to an insufficient number of simulation outputs to capture the growth between 40 − 60 per cent of the maximum  * .These systems are excluded from the plots that use  50 values (Figs. 5 and 6).

PROPERTIES OF DISRUPTED AND SURVIVING SATELLITES
In the following, we compare the properties of simulated dwarf galaxies accreted onto the MW-mass systems, both disrupted and surviving at present day, with the corresponding satellite populations in the MW.For observations, we use the data from Naidu et al. (2022), comprised of 9 disrupted satellites and 13 surviving, with properties collected from a variety of sources and listed in Tables 1 and 2 of their article (see also references therein).Specifically, the observed disrupted dwarf galaxies are Sagittarius, GES, Helmi Streams, Sequoia, Wukong/LMS-1, Cetus, Thamnos, I'itoi and Orphan/Chenab, while the surviving satellites are the Large and the Small Magellanic Clouds (LMC and SMC), Fornax, Leo I, Sculptor, Antlia 2, Leo II, Carina, Sextans, Ursa Minor, Crater 2, Draco and Canes Venatici I. We note, however, that the observational sample does not include all disrupted satellites, given the spatial limitation of the H3 survey (heliocentric distances of 3 − 50 kpc).Because of this, progenitors like Kraken, for example, are not included here.
For simulations, we choose to include all disrupted progenitors, as we intend to uncover general trends related to their host accretion histories.

Stellar mass -metallicity relations
We first compare the MZRs of the surviving and disrupted populations of all 42 simulated MW-mass systems with the corresponding observed MZRs for the MW.The left panel of Fig. 2 shows the distributions of mean [Fe/H] versus  * for all2 dwarf galaxies in the full simulated sample.For the disrupted populations in the simulations, the  * and [Fe/H] values are computed at  acc (i.e., when  * reaches its maximum value, according to our definition in Section 2.2).For the surviving populations the values correspond to  = 0, as in the observations.In agreement with the observations, simulations predict well-defined correlations between [Fe/H] and  * .This is seen in both the surviving and disrupted populations.Also, at fixed  * , the surviving satellites are typically more metalrich than the disrupted dwarfs, also in qualitative agreement with the observations.For a more quantitative comparison with the observations, we compute the best linear fits for the MZRs of the two dwarf populations in the simulations.Considering the entire dwarf populations in the 42 MW-mass systems, we obtain: [Fe/H] = (0.32 ± 0.01) log( * ) − (3.49 ± 0.10), for the surviving population (solid red line in left panel of Fig. 2), and [Fe/H] = (0.48 ± 0.02) log( * ) − (5.42 ± 0.13), for the disrupted one (solid purple line in the same left panel).Overall, the MZR fit for surviving satellites in simulations agrees with the one derived for the surviving satellites in the MW.For example, the slope of 0.32 (±0.01) in the simulations is similar to the slope corresponding to the observed sample from Naidu et al. (2022), which is ≈ 0.36, and it is also very close to the slope of ≈ 0.3 derived from a larger sample of Local Group satellites by Kirby et al. (2013) (the latter fit is shown with dashed black line in the left panel of Fig. 2).We recall that, for simulations, all [Fe/H] values were adjusted downwards by 0.4 dex, therefore the normalisation of the MZR in simulations matches the one for the observed satellites in the MW by construction.However, no additional recalibration was done to match the slope of the observed MZR.
The MZR fit for the entire disrupted population in simulations has a slope of 0.48 ± 0.02, which marginally agrees (within the errors) with the slope of 0.36 +0.12 −0.04 reported by Naidu et al. (2022) for the disrupted population in the MW.We note, however, that the latter value is derived from a sample of 7 disrupted dwarfs, which excludes Sagittarius and I'itoi (the Naidu et al. 2022 fit is shown with dashed blue line in Fig. 2).According to these authors, Sagittarius ( * = 10 8.8 M ⊙ , [Fe/H]= −0.96) is excluded because it is not yet fully disrupted and was still forming stars as recently as 2 Gyr ago, while I'itoi ( * = 10 6.3 M ⊙ , [Fe/H]= −2.39) appears to be an outlier in the relation, its metallicity being too low.However, by taking into account these two systems, the MZR for disrupted dwarfs in MW has a slope of 0.50 ± 0.08, which is in better agreement with the average value of 0.48 predicted for this population by our simulations.One determining factor for improving the comparison between simulations and observations is to better constrain the properties of I'itoi.Furthermore, discovering more debris from low-mass progenitors, similar to I'itoi, would help pin down the slope of the MZR for the disrupted population in the MW.Also, as mentioned previously, the observational sample of disrupted galaxies of Naidu et al. (2022) does not include systems whose debris are contained near the centre of the Galaxy ( < ∼ 6 kpc) or in the outskirts ( > ∼ 50 kpc), while our simulated sample includes all disrupted dwarfs.A more detailed comparison with observations, taking into account the location of merger debris and other observational selection effects, is left for a future study.
We note that eqs. 1 and 2 correspond to average fits for the entire populations of all simulated MW-mass systems, and significant variations can occur on a system-by-system basis.The right panel of Fig. 2 shows the different MZR fits for all disrupted populations in each of the simulated MW-mass hosts (purple lines).For comparison, the linear fits for the disrupted population in the MW from Naidu et al. (2022) are shown here, both excluding (dashed blue line) and including (solid blue line) Sagittarius and I'itoi.This figure shows that the two MZR slopes derived for the disrupted population in the MW (≈ 0.36 or ≈ 0.50, respectively) can be matched by several MW-mass systems.However, we can find more matches in the simulated sample of MW-mass systems to the observations, when the latter includes I'itoi and Sagittarius (solid blue line).This also indicates that the MW may have a fairly typical disrupted population for galaxies of its mass.
On the other hand, matching the MZRs of disrupted and surviving populations simultaneously in the simulations reduces the number of matches, as the MZR fits may differ in normalisation or slope than the values inferred from the observations.In addition, as it will be discussed in the next sections, fully realistic simulations should match other observables besides the MZRs, such as the [Mg/Fe] -[Fe/H] - * distributions of the two populations.For the disrupted population, we also investigate the evolution of its MZR with redshift.Since the whole population of disrupted dwarfs is composed of systems accreted at different redshifts, it is interesting to quantify the level of evolution of the MZR during different cosmic epochs, and to connect it to the scatter observed in the [Fe/H] - * distribution of disrupted dwarfs.For this analysis, we include the entire population of disrupted dwarfs in the simulations.Specifically, for each progenitor, we compute the  * and mean [Fe/H] values at accretion ( acc ).After this, the disrupted dwarfs are separated into different  acc ranges ([2.0, 3.5], [1.5, 2.0], [0.5, 1.5] and [0.0, 0.5]) and, for each range, we compute the mean [Fe/H] values in different  * bins.The mean MZRs obtained from this procedure are shown in Fig. 3. (We note that the results are similar if medians are used instead).
Fig. 3 shows that the normalization of the MZR evolves with redshift, changing by ≈ 0.25 dex between the  acc bins chosen here.As expected, the average [Fe/H] of progenitors increases with time, given that systems accreted later had more time to form stars and enrich in metals prior to their accretion onto the host.Also, there does not seem to be a marked difference between the latest accreted dwarfs that are now disrupted (i.e., the lowest  acc bin) and those that survive (grey curve), the two populations overlapping.This is because both of these populations are mainly composed of high mass dwarfs, which are more metal-rich.
This figure also shows that the slopes of the MZRs of the disrupted populations accreted at different  acc are very similar, with values consistently around 0.27 − 0.29.This suggests that the slope of the MZR of the MW progenitors does not change with redshift (only its normalization).Also, the slopes for the early progenitors are similar with the slope for surviving satellites (which is ≃ 0.32 in our simulations, see the grey line in Fig. 3; and ≃ 0.3 in the Local Group).
Overall, the slope of the combined population of disrupted dwarfs (of 0.48; see eq. 2 and the translucent purple line in Fig. 3) is steeper than the slopes for the individual  acc bins.This can be explained by the increase in the MZR normalization with decreas-ing  acc , plus the changes in the stellar mass distributions of the disrupted dwarfs with  acc .The latter are shifting progressively towards higher masses at lower redshifts, and this is because later accreted dwarfs have more time to form stars prior to accretion.
Finally, as discussed in Naidu et al. (2022), the disrupted population in the MW includes systems accreted at intermediate redshifts ( acc ≈ 0.5 − 2).In the case of ARTEMIS, Fig. 3 shows that, if we were to limit the sample of disrupted galaxies to this redshift interval, their mean MZR slope would be very similar to that derived for the simulated surviving population, in very good agreement with Naidu et al. (2022) findings.A discussion about the evolution of MZR with redshift, in the context of other studies, is presented in Section 4.

[Mg/Fe] -stellar mass and [Mg/Fe] -[Fe/H] distributions
Similar to the analysis of the [Fe/H] distributions of the disrupted and surviving dwarf populations, we investigate the different trends in terms of their alpha-enhancements.
In the left panel of Generally, the [Mg/Fe] abundances decrease with increasing stellar mass, for both surviving and disrupted satellites.This is because more massive galaxies sustain star formation for longer, and enrich more in Fe from Type Ia supernovae.Lower mass dwarfs have shorter periods of star formation, and therefore they enrich mainly in -elements (namely, Mg), which are created in Type II supernovae that occur on shorter timescales than Type Ia (Tinsley 1979).We note that the time of accretion is also an implicit factor here, as later accreted dwarfs have more extended star formation histories, and therefore have higher [Fe/H] and lower [Mg/Fe] at accretion.The dependency of chemical abundances on  acc will be investigated in more detail in the next section.
Overall, the [Mg/Fe] distributions show a plateau (or a shallow slope) at the metal-poor end ([Fe/H] ≲ −1), followed by a decrease towards higher metallicities.This trend is seen in both disrupted and surviving populations, and for both simulated dwarfs and observed ones.Also, the [Mg/Fe] -[Fe/H] trend obtained for the disrupted dwarfs is similar to the trend observed in the stellar halo of the MW (e.g., Venn et al. 2004), confirming the ability of ΛCDM models to produce realistic (accreted) stellar haloes of MW-mass galaxies.
One main result in Fig. 4 is that simulations predict disrupted populations that are more enriched in [Mg/Fe] than the surviving ones, at both fixed  * (left panel), and at fixed [Fe/H] (right panel).The trend agrees with the one seen in [Mg/Fe] between corresponding populations in the MW.Note, however, that the simulated chemical abundances show a large scatter in these offsets, compared to the MW data.This is mainly due to the inclusion of populations from all 42 simulated MW-mass systems.If, instead, we analyse the offsets on a system-by-system basis (as shown in Appendix A), their magnitudes are similar to the one measured between the MW populations.
In conclusion, the simulations are in qualitatively good agreement with the MW observations, in that:  iii) the simulations also predict that the disrupted dwarfs are more enhanced in [Mg/Fe] than the surviving ones, both at fixed  * and at fixed [Fe/H]; iv) in addition, the simulations match the observed MZRs for the two populations and retrieve the observed offset in [Fe/H] between disrupted dwarfs and surviving dwarfs at fixed  * (Fig. 2).

The scatter in the chemical abundance versus stellar mass
The nature of the scatter in the simulated distributions deserves further investigation.Figs. 2 and 4 show significant scatter in the chemical abundances of both disrupted and surviving populations.These figures also highlight two points: 1) the scatter is larger in the disrupted population than in the surviving one; and 2) the simulated chemical abundances present a larger scatter in both populations, compared with observations of disrupted and surviving satellites in the MW.
The latter point is plausibly explained by the inclusion of populations from many more MW-mass systems in the simulations (42), compared with the single system in the observations.The larger scatter in the disrupted populations compared with the surviving ones can also be explained, in part, by the larger number of disrupted systems over the lifetime of a MW-mass galaxy (Fattahi et al. 2020) compared with the number of surviving dwarfs, of which there are typically a dozen or so per host halo (see Font et al. 2021).
However, some of the scatter in the disrupted populations may be due to different accretion histories of the simulated MW-mass systems (e.g., different combinations of stellar masses and/or different accretion times in each disrupted dwarf population).In particular, the large variation in chemical abundances (both [Fe/H] and [Mg/Fe], at fixed  * ) in the disrupted populations suggests a connection with the accretion histories.For example, at fixed stellar mass satellites that were accreted relatively earlier will typically have had higher star formation rates than satellites of the same stellar mass but which were accreted later on.This will lead to different chemical abundances at the time of accretion.In addition, the star formation histories of satellites can be affected by the environmental conditions specific to their hosts; for example, satellites that have spent more time in denser environments (i.e., closer to their hosts) may have had their star formation rates curtailed which in turn will affect their chemical evolution.
To test this scenario, and determine how different accretion events can be identified in the scatter of the [Fe/H] - * and [Mg/Fe] - * distributions, we choose a few simple parameters as proxies for the accretion histories.First, we use a redshift associated with accretion or the mass growth of each dwarf, specifically,  acc for the disrupted dwarfs and  50 for the surviving ones (as discussed in Section 2.2).We also use the cold, star forming (SF) gas fraction in each dwarf satellite (  SF ), as this parameter has previously been found to correlate well with the distribution of chemical abundances (for example, with the metallicity scaling relations; see De Rossi et al. 2017Rossi et al. , 2018;;Maiolino & Mannucci 2019 and references therein).
Fig. 5 shows the [Fe/H] - * distribution for the disrupted and surviving dwarf populations in the simulations, using the same samples as in the left panel of Fig. 2. As before, [Fe/H] are the means of the stellar [Fe/H] values in these dwarfs, measured at  acc for disrupted dwarfs, and at  = 0, for the surviving ones.For the disrupted populations (top panels), the distributions are colour-coded by  acc (left) and  SF (right), while for the surviving ones (bottom panels), are colour-coded by  50 (left) and  SF (right), respectively.The cold SF gas fraction is defined as usual (see De Rossi et al. 2017):  SF =  SF /( * +  SF ), where  SF is the total gas mass in a dwarf galaxy which fulfills the required conditions to form new stars (see Schaye et al. 2015).Note that  SF are computed at  acc for the disrupted dwarfs, but at  = 0 for the surviving ones.(As it was done for Fig. 2, the values for the surviving populations are computed at  = 0 to facilitate the comparison with observations).
This figure shows that, for the disrupted dwarfs at least, the scatter is strongly dependent on both  acc and  SF .At fixed  * , dwarfs that are accreted early are more metal-poor; also, at fixed  * , dwarfs that are more metal-poor are also more gas rich.The correlation between  acc and  SF is not surprising, given that these parameters are not independent of each other.Dwarfs accreted early are expected to be more gas-rich, first because gas fractions generally were higher given that there has been less cosmic time to form stars; also, early on, dwarf galaxies can still retain their gas reservoir, whereas later, this gas gets depleted through supernovae winds and ongoing star formation.
The magnitude of the scatter in [Fe/H] also appears to vary with  * , the lower mass dwarfs showing the largest amount of scatter, and highest mass dwarfs, the least.In general, the magnitude of the scatter correlates with the range in  acc and  SF values.For example, disrupted progenitors of  * ≈ 10 6 −10 7 M ⊙ are accreted at redshifts ranging from  acc ≈ 6 to 0, whereas those of  * ≈ 10 9 − 10 10 M ⊙ are accreted from  acc ≈ 1 to 0. This is expected, as low mass dwarfs tend to form first and they could be accreted at any redshift, whereas high mass dwarfs require more time to grow and therefore are typically accreted late.A similar behaviour is seen in terms of  SF : low mass progenitors displaying a wide range of  SF values, corresponding to systems that are either gas-rich and star-forming, or that are partially quenched, whereas all LMC-mass systems have low  SF values, which indicates that they consumed most of their cold gas by the time of accretion.There is also an extreme case of progenitors without any cold gas at accretion (  SF = 0), which are indicated with empty circles in the top right panel of Fig 5 .These progenitors tend to be of very low mass, so they are prone to having their gas stripped by environmental processes, or to have their gas expelled by supernovae winds.We note, however, that towards the very low mass end, some of the scatter in dwarf properties may be due to the limited numerical resolution of the simulations.
The change in the magnitude of the scatter with  * is also seen in the surviving population (bottom panels of Fig. 5), however in this case the correlations with  50 or  SF are much weaker.Low mass surviving dwarfs still form at a wide range of redshifts ( 50 ≈ 3 − 0), and high mass ones typically form late ( 50 ≈ 1 − 0).Overall, however, the range in  50 (and, by extension, in  acc ) in the surviving dwarfs is narrower than the  acc range in the disrupted dwarfs, both at fixed  * and across all masses.This is because the surviving dwarfs are typically accreted late.Moreover, there does not seem to be a clear correlation between  50 and [Fe/H] of the surviving dwarfs at fixed  * , as it was previously the case between  acc and [Fe/H], for the disrupted ones.The dependence on  50 is more pronounced across the range of stellar masses, with progenitors that form early being typically more metal-poor.
Similarly, the present-day  SF fractions in the surviving popu- For the early accreted progenitors of disrupted dwarfs, the scatter in [Mg/Fe] values correlates, at fixed  * , with both  acc and with  SF .For progenitors of the same stellar mass, those accreted earlier have higher [Mg/Fe] abundances than those accreted late.Likewise, those progenitors with higher [Mg/Fe] at accretion, tend also to have higher gas fractions.However, the correlation of the scatter with  acc at fixed  * (left panel), appears to be stronger than the correlation with  SF values (right panel).This suggests that, although  acc and  SF are themselves correlated, they do not always track each other.The  SF values appear to be either very high, or very low, which indicates that the progenitors may experience quenching processes on a relatively short timescales.
We note that some of the low-mass, metal-poor disrupted satellites present unusually high [Mg/Fe] abundances, e.g., > 0.4.These galaxies were accreted at high redshifts ( > 4, see Section 3.3) which may indicate that the nucleosynthetic yield tables implemented in the ARTEMIS simulations are more uncertain at these ranges.The choice of the time delay between SNII and SNIa may also contribute to unusually high [Mg/Fe] abundances for the early accreted dwarfs. elements are mainly synthesized in SNII, whereas Fe are mainly produced in SNIa.Since SNII oc- cur in younger stellar populations than SNIa, a high level of enhancement in early accreted satellites is expected.Wiersma et al. (2009) found that, up to 10 8 yr, the integrated metal ejecta in their model is dominated by SNII, whereas the contribution of SNIa for Fe becomes important after 10 9 yr.
The bottom panels in Fig. 6 show the [Mg/Fe] - * distributions for the surviving dwarfs, colour-coded by  50 and by  SF , respectively.A weak correlation with  50 is seen at the low-mass end, where at fixed  * , dwarfs with higher [Mg/Fe] tend grow in mass more recently than those with lower [Mg/Fe].This behaviour is not seen at higher masses though.At low masses, surviving dwarfs with higher  50 have more time to evolve, increasing their Fe content and therefore having lower [Mg/Fe].The right panel shows that, for massive systems ( * > ∼ 10 7 M ⊙ ) and at fixed M * , dwarfs with higher [Mg/Fe] tend to have higher fractions of star-forming gas.
We have seen that, for the disrupted dwarfs,  acc and the fraction of cold gas these systems have at this redshift,  SF are correlated.Typically, progenitors accreted earlier are more gas rich.We explore this correlation in more detail in Fig. 7, where we plot  SF versus  acc for all disrupted dwarfs, with star symbols colour-coded by the mean chemical abundances of each dwarf ([Mg/Fe] and [Fe/H], respectively), measured at accretion.Generally,  SF values are high at early times and decrease sharply below  acc ≃ 4. The scatter in the relation between  SF and  acc is fairly large above redshift ≈ 4, however towards low redshift the two parameters correlate well.The galaxies lower  SF values than the average could be related to environmental effects or stochastic star formation, particularly at the low mass end.Additionally, the  SF values could be affected by numerical resolution.For example, the population of low-mass dwarfs with  SF ≃ 0 (shown with empty circles in Figs. 5 and 6) display a wide range of  acc values.This suggests that their gas has been removed/consumed at a much faster pace than the average population accreted at that redshift.Note that the apparent vertical stripes in  acc , which are most visible at high redshifts, are just due to frequency of simulation snapshots that have been output.The merger tree algorithm does not interpolate between snapshots but simply uses the redshift of the snapshot which is nearest to the accretion event.
In summary, the scatter in averaged chemical abundances ([Fe/H] and [Mg/Fe]) of disrupted dwarfs at fixed  * is larger than the corresponding scatter in the surviving dwarfs.Also, the scatter in the distributions for the disrupted dwarfs correlates with the  acc and with the cold gas fraction at the time of accretion,  SF of the progenitors, whereas in the case of surviving dwarfs, there is no correlation in the scatter with  50 (or  acc ), and their present-day  SF fractions are typically very low.Also, the chemical abundances of the disrupted dwarfs are more diverse than those of the surviving satellites, which are similar regardless of the properties of their host (e.g., total mass) or the past accretion history in that system.This suggests that the chemical abundances of disrupted population may be a good indicator of the accretion history of their host, for example using chemical abundances to extract the stellar masses and redshifts of accretion of the progenitors.This may help in the reconstruction of the formation histories of MW-mass galaxies, as chemical abundances are more readily available from the debris found in the MW today, rather than total stellar masses and accretion redshifts, which usually are inferred via modelling of the orbits of satellite progenitors.

Global properties of the disrupted dwarfs populations
Our results so far show that the chemical abundances of the disrupted dwarfs correlate with their  acc , suggesting that the accretion history of a MW-mass galaxy could be inferred from the chemical abundance properties of its disrupted population.In contrast, the chemical abundance properties of surviving populations (see bottom panels of Figs. 5 and 6), seem to have little connection with either the properties of their hosts, or with those of the disrupted population.
In this section we explore the global properties of the disrupted populations, in order to determine which diagnostics could be used to infer the accretion histories of their hosts.Specifically, we select several observables that involve the chemical abundance distributions of the disrupted populations.These are: the average slopes of the MZRs of the disrupted populations, their mean [Fe/H] and [Mg/Fe] abundances at accretion, and the average slopes of the [Mg/Fe] - * distributions.These diagnostics are motivated by the correlation with  acc in the scatter of these relations for the disrupted populations (see Figs. 3,5 and 6).The accretion histories are also measured globally, using the average redshift of accretion of the disrupted population, ⟨ acc ⟩, which is defined as the stellar massweighted mean of individual  acc of all disrupted dwarfs in a given MW-mass host.Likewise, for chemical abundances, we compute mass-weighted, global values for the entire population of disrupted dwarfs in a given MW-mass system.Therefore,

⟨𝑋⟩ = ∑︁
where the index  denotes a given disrupted dwarf galaxy,  * , is its stellar mass measure at accretion, and  is the parameter that is being computed (i.e.,  acc , [Fe/H] or [Mg/Fe]).We compute the same global parameters ⟨⟩ for the disrupted dwarfs in the MW, using the data from Naidu et al. (2022).We note that these authors use a different definition for  acc , namely as the redshift associated with the first pericentric passage of the progenitor of the disrupted dwarf inside the MW.We retain this definition in our computation of ⟨ acc ⟩ for the observed MW population, however we note that in some cases the progenitors of the disrupted dwarfs can spend a significant amount time between accretion onto the host and their first pericentric passage inside it (or the final time of disruption), as will be discussed in Section 3.7.
Fig. 8 shows the relation between the mean chemical abundance values, ⟨[Mg/Fe]⟩ and ⟨[Fe/H]⟩, for the disrupted popula-tions in all simulated MW-mass systems.Each full circle corresponds to a disrupted population of a given MW-mass system.To these, we add the corresponding values for the disrupted population in the MW.As before, we use data from Naidu et al. (2022), and consider two cases: one in which Sagittarius and I'itoi are included in the sample (blue square), and one in which they are not (black square).Hereafter, the first sample (which contains 9 dwarfs) is called "MW", whereas the second sample (containing 7 dwarfs) is called "MW*".As discussed in Section 3.1, I'itoi and Sagittarius were excluded by Naidu et al. (2022) in their estimate of the MZR for the disrupted population, however we find that by including these two systems we obtain a better agreement between the predicted MZRs of disrupted populations and the observations (i.e., the "MW full sample" in the left panel of Fig. 2).
This figure shows a strong anti-correlation between the global ⟨[Mg/Fe]⟩ and ⟨[Fe/H]⟩ values of the simulated disrupted populations.The best linear fit (shown with dashed line in both panels) has a slope of −0.40 ± 0.015, with a corresponding Spearman coefficient of   = −0.95.Moreover, both observational data points (both the MW and MW* samples) follow the predicted ⟨[Mg/Fe]⟩ -⟨[Fe/H]⟩ relation.The bottom panel of Fig. 8 shows the same relation for the simulated disrupted populations as the top panel, but now colour-coded by ⟨ acc ⟩ (computed as in eq. 3.This shows an anti-correlation between ⟨[Fe/H]⟩ and ⟨ acc ⟩ and a correlation between ⟨ acc ⟩ and ⟨[Mg/Fe]⟩.These resemble the individual relations between [Fe/H] and  acc , and between [Mg/Fe] and  acc for disrupted dwarfs seen in Figs. 5 and 6, but in this case, the correlations are much stronger. Overall, these relations between ⟨[Mg/Fe]⟩ and ⟨[Fe/H]⟩, and each in turn with ⟨ acc ⟩, suggest that these observables are good diagnostics for constraining the accretion histories of MW-mass hosts.However, we note that as these are global parameters of a given disrupted population, they are more sensitive to the most massive accretion events in this population (given that they are mass-weighted estimates of the accretion history).
For a more in-depth comparison with the disrupted population of the MW, we choose six simulated disrupted populations whose ⟨[Mg/Fe]⟩ and ⟨[Fe/H]⟩ values fall close to the measurements in the MW* sample.These populations are hereafter called the " 6 * subsample".These are located inside the green rectangle centred on the MW* data point in the top panel of Fig. 8, and are also shown with green circles.In the ARTEMIS simulations, these populations correspond to the MW-mass hosts: G1, G23, G26, G29, G38 and G42 (see Font et al. 2020).We note that four of these MW-mass hosts (G1, G23, G29 and G42) are among the disc galaxies with a GES feature identified by Dillamore et al. (2022).For a detailed characterization on a system-by-system basis in the  6 * subsample, see Appendix A. We prefer to select our subsample based in its relation to the MW* sample instead of the MW one since the analysis carried out by Naidu et al. (2022) is mainly done on the former.
As the  6 * populations match the averaged properties of the MW* sample (in ⟨[Fe/H]⟩, ⟨[Mg/Fe]⟩ and ⟨ acc ⟩), we want to know whether they also match other observables, for example the slope of the MZR in the MW* sample, or the slope of the [Mg/Fe] - * distribution.The latter are also potential useful diagnostics, as suggested by the results in previous sections.The panels in Fig. 9 show, from top to bottom, with purple circles: the relation between ⟨[Fe/H]⟩ and ⟨ acc ⟩ for all disrupted populations in the simulations, the relation between the slopes of the MZR of these disrupted populations and the corresponding ⟨ acc ⟩, and the relation between the slopes of the [Mg/Fe] - * relations of these populations and ⟨ acc ⟩, respectively.The six disrupted populations in the  6 * sub-sample are highlighted again in green.The two disrupted population in the MW are shown again with empty squares, blue for the "MW" sample and black for "MW*".
The top panel shows that ⟨[Fe/H]⟩ correlates very strongly with ⟨ acc ⟩.Our best linear regression fit to the ⟨[Fe/H]⟩ -⟨ acc ⟩ relation, shown in this panel with the black dashed line, is: while the Spearman coefficient of this relation is   = −0.95.Interestingly, the same   value was obtained for the ⟨[Mg/Fe]⟩ -⟨[Fe/H]⟩ relation in Fig. 8, which indicates that these diagnostics are equally strong for inferring the accretion history of the population.For example, one can infer the ⟨ acc ⟩ of a given disrupted population by estimating both ⟨[Fe/H]⟩ and ⟨[Mg/Fe]⟩ and using the prediction from the bottom panel of Fig. 8, or by estimating only ⟨[Fe/H]⟩ and using the fit from eq. 4.These predictions are useful as they require only chemical abundance measurements, whereas estimating ⟨ acc ⟩ from individual  acc values typically would involve additional assumptions, e.g., about the total mass of the host, or about the orbits of the infalling dwarf galaxies.
The overall trend of increasing ⟨[Fe/H]⟩ with decreasing ⟨ acc ⟩ can be understood as a result of more chemical enrichment taking place in dwarf galaxies accreted late, as well as to the iron abundance being specially sensitive to the occurrence of delayed SNIa at later times.Overall, those disrupted dwarf populations which are accreted more recently tend to include more massive dwarfs, which are more metal-rich at the time of accretion.Note that this relation includes mass-weighted values, and so it tracks mainly the evolution of most massive disrupted dwarfs.The correlation found here is also expected on the basis of the most massive progenitors obeying a MZR.
Our observational estimates for ⟨[Fe/H]⟩ and ⟨ acc ⟩ for the MW* and MW samples are consistent with the prediction from the simulations (eq.4).In addition, the data points for the six disrupted populations in the  6 * subsample (green circles) fall in the vicinity of the MW* sample in the ⟨[Fe/H]⟩ -⟨ acc ⟩ plane (note that the  6 * subsample was selected on the basis of [Mg/Fe] and [Fe/H] alone), although their corresponding ⟨ acc ⟩ values show a somewhat wide range, from 1.76 (G23) to 2.42 (G1).For the observed disrupted population (i.e., MW*), the simulated fit (eq.4) predicts an average redshift of accretion of ⟨ acc ⟩ ≃ 1.8, which agrees with our estimate based on the MW's data from Naidu et al. (2022).For the observational sample which includes both Sagittarius and I'itoi (MW sample, blue square), the simulated fit predicts a slightly lower value, of ⟨ acc ⟩ ≃ 1.6, but the estimated value from the observations of ≃ 1.2 is within the simulated scatter at that ⟨[Fe/H]⟩.The lower ⟨ acc ⟩ / higher ⟨[Fe/H]⟩ data point for the MW sample compared to MW* is due mainly to the inclusion of Sagittarius, which is accreted more recently and is also fairly massive.
The middle panel of Fig. 9 shows distribution of the slopes of the MZRs in the simulated disrupted populations and the corresponding ⟨ acc ⟩.(The individual MZRs of these populations have been shown before, in the right panel of Fig. 2).The populations from the  6 * subsample are highlighted again in green, and the estimates for the observational MW* and MW samples are added with the same symbols as in the top panel.The results do not support a correlation between the slope of the MZR and ⟨ acc ⟩, which is consistent with the results in Fig. 3.The linear regression fit to the simulated populations has Spearman coefficient of   = −0.11.The fit, shown with dashed black line in the middle panel, has a negligible slope of −0.06 ± 0.03 (the y-intercept is 0.54 ± 0.04).This result is expected, as the individual MZRs of the disrupted populations are not significantly different in slopes and normalisations, although there is some scatter in these relations.The MZR slope estimates for the disrupted population in the MW are matched, within the scatter, by the simulations.
The relation between the slope of the [Mg/Fe] - * relation of each disrupted population and the corresponding ⟨ acc ⟩ is explored in the bottom panel of Fig. 9.Here again, there is little or no correlation between the two parameters.The linear regression fit to the simulated data (shown again with black dashed line) has a slope of 0.02 ± 0.02 and a y-intercept of −0.11 ± 0.02, but a low Spearman coefficient, of   = 0.3.Compared with this fit, the two slopes derived from observational data are high, but they are still within the predicted scatter.Out of the six disrupted populations in the  6 * subsample, G26 agrees most with the estimate for the MW* sample.
In summary, simulations predict that the ⟨ acc ⟩ of a disrupted population can be very well constrained by the mass-weighted chemical abundances (⟨[Fe/H]⟩ and ⟨[Mg/Fe]⟩), but poorly constrained by other observables, such as the slope of the MZR of the population, or the slope of the [Mg/Fe] - * distribution.

The distribution of stellar masses in the disrupted and surviving populations
Our analysis indicates that the disrupted progenitors tend to be increasingly more massive at accretion towards lower redshift (see, e.g., Fig. 3).Here we investigate in more detail the distributions of stellar masses of the progenitors of disrupted and surviving dwarfs.We also track the progenitors of the six disrupted populations in the  6 * subsample separately, to assess whether the progenitors in the hosts resembling the MW (in terms of their global accretion histories) have a different stellar mass distribution from the populations in the rest of the simulated hosts.In Fig. 10 the distributions of stellar masses of populations in the  6 * subsample are shown with solid green lines, while the corresponding distributions for the population in the rest of the simulated hosts are shown with black dashed lines.This is done separately for the disrupted and surviving populations (top and bottom panels, respectively).The disrupted populations (which, again, comprise systems accreted at various redshifts) show a more even distribution in stellar masses than the surviving populations.The stellar mass distribution of surviving satellites is similar to a power-law, which is the expected mass distribution of satellites at any given time.The roughly even distribution of masses in the disrupted population is the result of a combination of factors.Extending the period of accretion to effectively a Hubble time, means that stellar mass distribution encompasses more massive satellites.The more massive satellites are also more efficiently disrupted by tidal forces; as the tidal forces scales as the square of the total mass of the satellite, this would lead to more massive satellites sinking faster and more likely to be fully disrupted (hence included in this plot).
Comparing the two types of progenitors in the  6 * subsample versus those in the rest of the simulated hosts, we find that the corresponding stellar mass distributions are generally similar.There is however a slight tendency for the  6 * hosts to have more low-mass dwarf progenitors, for both disrupted and surviving populations.The trend is more noticeable for the surviving dwarfs, however this can be attributed to low number statistics (i.e., fewer number of hosts results in even more rare high mass, recent accretions).
The finding that disrupted populations in the  6 * subsample are composed of slightly less massive dwarfs than those in the rest of the hosts, is consistent with the slightly lower than average values of ⟨[Fe/H]⟩ (and slightly higher than average ⟨[Mg/Fe]⟩) values obtained for these populations with respect to the majority of the other MW-mass systems in Fig. 8.This is corroborated by the higher ⟨ acc ⟩ of the disrupted populations in the  6 * subsample compared with the rest of the hosts in Fig. 9, which is also consistent with a more significant population of lower mass dwarfs (as higher mass dwarfs require more time to form).This suggests that the ⟨ acc ⟩ of a disrupted population is also indicative of its stellar mass spectrum, and vice versa. .Bottom: Same distributions but for the surviving dwarf galaxies.The  6 * hosts have slightly higher fractions of lower mass progenitors than the rest of the sample, for both disrupted and surviving populations.

The environments of disrupted and surviving dwarfs
The chemical abundances of the dwarf progenitors can also be affected by the environment in which they form and evolve prior to accretion.Here we investigate the scenario proposed by Naidu et al. (2022), in which the disrupted population of dwarfs originate in a denser environment, closer to the MW host, than the surviving population.To test this in the simulations, we track the location of the dwarf progenitors at different redshifts and estimate the environments in which they evolve.
To quantify the environment, we use the tidal parameter Θ 5 proposed by Karachentsev et al. (2013), and defined as: where   are the masses of the 5 nearest galaxy neighbours to dwarf galaxy  and   are the physical distances between these neighbours and the given dwarf galaxy at a given redshift .This is effectively a measurement of the tidal forces that nearby galaxies exert onto a given dwarf at that redshift.The individual Θ 5 values for each satellite dwarf galaxy/progenitor are computed until the object crosses the virial radius of its host, namely the  200 radius.
We also consider two ways to compute the masses of the neighbours: either as the cold baryonic masses,  , =  * +  SF , or as the total (dynamical) masses,  , =  * +  SF +  NSF +  DM ,  , masses) with , for both the surviving (grey) and disrupted (purple) dwarfs.Solid lines represent the median relations for each population, while the purple and grey shaded regions enclose the corresponding 25 th and 75 th percentiles of these medians, respectively.Dashed blue and black vertical lines are the median of all redshifts when the populations cross the virial radii of their hosts (i.e., the median of  acc|R 200 ), while the blue and grey vertical shades represent the 25 th and 75 th percentiles around these medians.
Higher Θ 5 values for disrupted progenitors across all  compared with Θ 5 values of surviving dwarfs, indicate that the disrupted progenitors are in denser environments.Bottom: Similar to the top panel, but using  , for computing Θ 5 .
where the last two terms indicate the mass of hot or diffuse gas (not star-forming), and of dark matter, respectively.The first method aims to reproduce the one used by Karachentsev et al. (2013) for the observations, where masses are inferred from the luminosities of galaxies, which are then multiplied by a constant factor.The second method is more physically motivated, as it tracks the intrinsic variations in mass-to-light ratios across different masses.Here we consider all galaxies as potential neighbours for the computation of Θ 5 , including the "main haloes", i.e., the progenitors of the central galaxies corresponding to the MW-mass hosts.Fig. 11 shows the evolution of Θ 5 with redshift for both the disrupted and surviving dwarf populations, and also using both methods for computing Θ 5 (using  , in the top panel, and  , in the bottom one).The solid lines represent the median Θ 5 values for the disrupted populations (purple) and surviving ones (grey), respectively.The dashed vertical lines in this figure are the median redshifts associated with the dwarf satellites crossing for the first time the  200 boundaries of their hosts,  acc|R 200 .The dashed blue line corresponds to the median  acc|R 200 of the disrupted population, while the dashed black line corresponds to the median of the surviving one.Likewise, the shaded blue and grey vertical shaded regions represent the corresponding 25 th and 75 th percentiles around these two medians, respectively.
This figure shows that the progenitors of disrupted dwarfs have higher Θ 5 values than the surviving satellites, across all redshifts considered here.This is true for both methods of computing Θ 5 (either using  , or  , ).Overall, these results suggest that the progenitors of the disrupted dwarfs form and evolve in much denser environments than the surviving dwarfs, which is in agreement with the scenario proposed by Naidu et al. (2022) for the MW populations.Moreover, our simulations indicate that this behaviour is generally valid for the populations of all MW-mass hosts, not just the MW.
Fig. 11 also shows that the progenitors of disrupted satellites cross the boundary  200 of their hosts earlier than the surviving satellites.The corresponding median  acc|R 200 redshifts for the disrupted and surviving populations are ≈ 2.2 and ≈ 0.8, respectively.Note that these populations comprise satellites from all simulated MW-mass hosts, whereas the individual median  acc|R 200 values for populations pertaining to each host may vary considerably (see the vertical shaded regions).However, for each host, the progenitors of disrupted satellites are accreted earlier than the corresponding surviving population, as expected.
The median Θ 5 values for the disrupted population show also an interesting variation with redshift, staying relatively constant until  ≈ 4 and then increasing until close to  = 0.This is mainly the result of our choice of including the main haloes as neighbours in the calculation of Θ 5 .At low redshift the progenitors of both populations are located mainly in the vecinity of their hosts.
Fig. 11 shows consistent patterns in the types of environments in which disrupted and surviving populations evolve, across all simulated MW-mass systems.However, there is a significant amount of scatter in the Θ 5 values at any given redshift.In Fig. 12 we investigate in more detail the evolution of the Θ 5 parameter for the two populations in the six simulated MW-mass systems which are more similar to the MW, namely the  6 * subsample.This allows us to assess the variations in Θ 5 versus  on a system by system basis, and also explore any differences that may occur in MWlike hosts.This figure shows that the median Θ 5 values for the disrupted dwarfs of these systems are roughly within the 25 th and 75 th percentiles around the medians of the corresponding simulated samples including all hosts, although, occasionally, the Θ 5 values reach above and below these limits.Note also that the Θ 5 parameters are not computed up to  = 0 in the case of the disrupted dwarfs in these systems, given that these populations are typically accreted at redshifts of ≈ 1.5 − 2 and tend to get fully disrupted before present time.
Overall, our findings support the scenario in which the disrupted dwarfs inhabit denser environments than surviving ones.Moreover, we find that this is the case at all redshifts, and that is a common pattern in all MW-mass hosts, regardless of whether they have similar accretion histories to the MW (i.e.,  6 * subsample) or not.
The scenario proposed by Naidu et al. (2022) envisages that the reason that environments of disrupted dwarfs are denser than those of surviving ones, is because the progenitors of the disrupted dwarfs form earlier, and closer to their hosts.To verify whether the simulated populations are located at different distances from their hosts, we trace their evolutionary paths from their time of formation until they the last time they are identified in the simulations.The top panels in Fig. 13 show the distribution of physical distances  MH  The same distributions as in the top panels, but for the populations in the six hosts in the  6 * subsample.The differences between the locations of the two satellite populations are similar to those in the full sample, being slightly more pronounced in this case.
of the progenitors of the two populations relative to the centers of potential of the corresponding progenitors of their hosts (i.e., the main haloes).The  MH distances are shown at two fixed redshifts,  = 2 and at  = 1 (left versus right panels), for both disrupted and surviving satellites (purple and grey histograms, respectively).
In general, the progenitors of disrupted dwarfs are located significantly closer to their hosts than the surviving dwarfs, at both redshifts.For example, at  = 2, about a third of the progenitors of disrupted dwarfs are already within the virial radius of their host (shown with vertical dashed line in this figure) whereas less than 10% of surviving dwarfs are that close.By  = 1, about 60% of disrupted progenitors are within the virial radius of their hosts, compared with less than 30% of surviving ones.The shapes of the distributions are also different, the distributions for disrupted satellites being narrower than those of the surviving ones, at both redshifts.The differences are more well defined at  = 2 than at  = 1, suggesting that these differences are generated early on, likely from the times of formation of these populations.
The bottom panels of Fig. 13 show the same analysis, but for the populations in the  6 * subsample.Here each population (disrupted / surviving) combines all dwarf systems of this type from the six MW-mass hosts in the subsample.The results are similar to the case of the entire sample, namely that the progenitors of the disrupted dwarfs are closer to their hosts than the surviving population, at both redshifts.The differences appear to be slightly stronger in this case (however, note the smaller sample size).Specifically, at  = 2, the distribution peaks at distance below 100 kpc from the host (which is well within the  200 at that time), whereas the corresponding distribution for surviving satellites peaks around 250 kpc (i.e., outside that radius).By  = 1, ≃ 78% of progenitors of disrupted dwarfs are within  MH = 75 kpc from their hosts, compared with ∼ 66 of surviving satellites.
In summary, the trends seen in Figs.11, 12 and 13 are consistent with the scenario proposed by Naidu et al. (2022), in which disrupted dwarfs form in denser environments, closer to their hosts.We also show that the behaviour is common for all MW-mass hosts, not just the ones matching the MW's accretion history.Additionally, we have shown that the progenitors of disrupted dwarfs continue to evolve in denser environments than the surviving satellites, although the differences between the two populations tend to decrease towards low redshift (below 1).

Disruption times
For the disrupted population we have so far focused mainly on their properties at the time of accretion ( acc ), or prior to that.In this section, we estimate the timescales in which they disrupt after they are accreted onto their hosts.As mentioned earlier, if dwarf galaxies spend long periods of time orbiting inside their hosts they may experience additional star formation and therefore more chemical evolution.We expect that mostly high mass satellites could be undergo additional star formation, as low mass satellites are effectively stripped of their cold gas and quenched at (or even before) accretion, as discussed in Section 3.3.
We therefore estimate the disruption times for the progenitors of satellites fully disrupted by  = 0, in all MW-mass hosts.These are computed as the time intervals   between the first time when the progenitors cross the virial radii (R 200 ) of their hosts,  acc|R 200 , and the last time 3 when they are identified as individual subhaloes by SUBFIND in the simulations,  f .Therefore, The redshift associated with the progenitor crossing the R 200 boundary of its host has been introduced before, and denoted  acc|R 200 .
Note that we take into account here only satellites which are fully disrupted at  = 0.According to our definition, those satellites which are still in the process of disruption at  = 0 are considered as surviving.
The left panel of Fig. 14 shows the distribution  dis versus the maximum stellar masses of the corresponding progenitors (this is  * measured at  acc , as discussed in Section 2.2).These are shown as star symbols, which are colour-coded by  acc of each progenitor.This figure shows that progenitors which are accreted very early (e.g.,  acc between 3 − 4) are typically of low mass ( * ≲ 10 8 M ⊙ ), and have very short disruption times (  ≲ 1 Gyr).This is expected given that the progenitors of the hosts are small in size at high redshift, and therefore the timescales on which satellites merge with the main halo at these epochs are very short.Satellites that infall at  acc ≲ 2 (green and blue tones) typically have higher masses, which is also expected as they have more time to form stars prior to accretion.As their hosts are also larger in size at low redshift, satellites tend to take longer to disrupt.This figure shows that massive satellites can have a wide range of disruption times (typically, between 0 − 6 Gyr for those accreted at  acc ≃ 2 and up to 9 − 10 Gyr for those accreted  acc ≲ 1).The longer disruption times for more massive satellites, coupled with them being more gas-poor at accretion (see Figs. 5 and 6), suggests the possibility of some of them undergoing additional star formation inside their hosts.However, this effect is probably not significant, because, as suggested by Fig. 1, most massive progenitors have already reached their maximum stellar mass by the time of accretion.
The right panel of Fig. 14 shows a similar analysis, this time plotting  dis versus the lookback time of accretion,  acc,L .The latter computes how long ago the progenitors have been accreted onto their host (compared to present-day), whereas   measures the time taken to disrupt since accretion time  acc|R 200 .The star symbols here are colour-coded by stellar mass at accretion,  * .This panel shows a similar result as the left panel, in that progenitors accreted a long time ago ( acc,L ≃ 10 − 12 Gyr) are typically low mass ( * ≃ 10 6 − 10 7 M ⊙ , blue tones) and are disrupted quickly (in less than 1 Gyr).Progenitors accreted between 6 − 9 Gyr ago show the largest range in   .Interestingly, the scatter in   does not appear to depend on the masses ( * ) of progenitors, suggesting that other factors are at play (e.g., the total mass of the host, or the orbital properties of the infalling dwarfs).Satellites accreted within the past 2 − 5 Gyr have short disruption times (between 0 − 5 Gyr) and are fairly massive.However the short disruption times are mainly a selection effect, as most recently accreted satellites are not yet fully disrupted, and therefore are not included in this plot.Also, no satellites accreted less than 2 Gyr ago are fully disrupted in our simulated sample, for the same reason.
3 Given the limited number of simulation outputs, both,  acc|R 200 and  f , correspond to approximated measurements of the actual values.In practice however, we consider  acc|R 200 as the time associated with the first snapshot when SUBFIND identifies the dwarf as a subhalo inside the  200 boundary of the host, while  f corresponds to the last snapshot when SUBFIND identifies the dwarf as a gravitationally bound subhalo.

DISCUSSION: THE EVOLUTION OF THE MZR WITH REDSHIFT
Our results give insight into the stellar mass -(stellar) metallicity relation of the disrupted population of satellite dwarf galaxies in MW-mass systems.The simulations indicate how this relation evolves with redshift (Fig. 3) and predict correlations of the scatter in [Fe/H] at fixed  * with other characteristics of the accreted dwarfs, namely with their redshift of accretion ( acc ) and the cold gas fraction (  SF ) at the time of accretion (Fig. 5).Here we discuss these findings in the context of previous work.
The stellar mass -stellar metallicity relation has been studied extensively in the Local Universe, where observations show it is a global scaling relation for galaxies which extends over many orders of magnitude in stellar mass (e.g., Gallazzi et al. 2005Gallazzi et al. , 2014;;Trussler et al. 2020).This scaling relation continues at the low-mass end of galaxy formation, in the regime of "classical" and even ultra-faint dwarf satellites of the MW (Kirby et al. 2013;Simon 2019).Cosmological simulations qualitatively match this relation over a wide range of stellar masses (e.g., Schaye et al. 2015;Ma et al. 2016).A similar scaling relation is observed between the stellar masses and the gas-phase metallicities of galaxies (e.g., Lequeux et al. 1979;Skillman et al. 1989;Tremonti et al. 2004;Kewley & Ellison 2008;Maiolino et al. 2008;Sanders et al. 2021), which is also broadly matched by the simulations (e.g., De Rossi et al. 2015;Ma et al. 2016;Davé et al. 2017;De Rossi et al. 2017;Torrey et al. 2019).
Several models have been proposed to explain the origin and shape of the MZR.These include: strong stellar feedback (Dekel & Silk 1986), where supernovae winds expel the metal-enriched gas from low-mass galaxies due to their shallower potential wells; an equilibrium metallicity in galaxies (Finlator & Davé 2008), which is obtained by balancing the gas inflow rate with the gas consumption rate from star formation plus the gas outflow rate;by energy-driven supernova winds (Guo et al. 2016;Davé et al. 2012), by inefficient star formation in dwarf galaxies (Brooks et al. 2007), or by a topheavy initial mass function (Köppen et al. 2007).
The validity of these models can be tested by observing the evolution of the MZR with redshift and comparing it with the model predictions.However, at higher , the MZR is much less well constrained, in both observations and simulations.While it is well established that, at fixed  * , the high  galaxies have lower metallicities than those in the Local Universe (e.g., Erb et al. 2006;Sommariva et al. 2012;Gallazzi et al. 2014;Leethochawalit et al. 2018), there are currently uncertainties in terms of the normalisation of the relation (Kewley & Ellison 2008) or the rate of increase in metallicity at different redshifts (Maiolino & Mannucci 2019).
Some comparisons of our results (Fig. 3) with the observations are possible, particularly at the higher mass end of the simulated dwarfs.For systems in the range of 8.5 < log(M * /M ⊙ ) < 10.2, we find a good agreement with the data from the VANDELS survey (Cullen et al. 2019), where a metallicity offset of ≈ 0.6 dex is found between local galaxies and those at 2.5 <  < 5.0.Our predicted MZR evolution at the high-mass end of dwarf satellites is also similar to that obtained by Gallazzi et al. (2014) for their sample of massive galaxies ( * > 10 10 M ⊙ ) at  ≈ 0 − 0.7.
Overall, we find a gradual decrease of the normalization of the MZR with redshift (see Fig. 3).This behaviour is consistent with the evolution of the MZR obtained in other simulations, for example, in the IllustrisTNG (Torrey et al. 2019) The driving force behind the evolution of the MZR is likely the fraction of cold gas in galaxies.This would also explain the correlation of the scatter in the MZR with the gas fractions at different redshifts (see Figs. 2 and 7), which is seen in other cosmological simulations (e.g., Ma et al. 2016;Davé et al. 2017;De Rossi et al. 2017;Torrey et al. 2019).The main trend obtained in our study, namely that more recently accreted dwarfs have lower gas fractions and higher metallicities, is in general agreement with previous results.Another physical process that might contribute to the evolution of the MZR is the increasing occurrence of Type Ia supernova as redshift decreases, which raises the concentration of iron in the gas that will form new stars (e.g., Leethochawalit et al. 2019).
We note that the Auriga simulations also show differences in [Fe/H] between the surviving and disrupted satellites at fixed stellar mass, and that the scatter in [Fe/H] of surviving satellites depends on the infall time, at fixed  * (see figure 5 in Fattahi et al. 2020).However, the scatter in [Fe/H] of disrupted satellites does not show a clear dependence with infall time (Fattahi, priv. comm.).One reason for this could be due to slightly different definitions of infall and accretion times.A more detailed comparison between the two simulations may shed light on these differences, however we leave this for a future study.
The evolution of the slope of the MZR is another constraint on the models.Our results (Fig. 3) suggest that the slope does not change significantly with redshift, at least in the regimes investigated here, of dwarf satellites and up to  ≈ 4.This behaviour is also in agreement with observations (e.g., Zahid et al. 2014;Sanders et al. 2021).
Our simulations also indicate a steeper MZR slope for the combined population of disrupted satellite galaxies (of ≈ 0.48) than for the MZRs of satellites accreted at any given redshift (slopes of ≈ 0.3).This prediction can be further tested with improved observations in the MW, in particular by constraining the properties of the lowest mass/ earliest accreted satellite galaxies, such as I'itoi or other similar systems.

CONCLUSIONS
We used the ARTEMIS suite of zoomed-in, cosmological hydrodynamical simulations to investigate the properties of the satellite dwarf galaxies in MW-mass hosts.The simulations comprise 42 MW-mass hosts, with a variety of accretion histories.We focused on comparing the properties of disrupted satellites versus those surviving at  = 0, and on comparisons with observations of these populations in the MW.We found several properties of disrupted satellites that connect with their redshift of accretion, and therefore can be used to constrain the accretion history of MW-type galaxies.Our main results can be summarised as follows: • Disrupted dwarf galaxies follow a stellar MZRs (Fig. 2), but with lower normalisation than the MZR for the surviving dwarfs (lower [Fe/H] at fixed  * ).Combining progenitors of disrupted dwarfs from all MW-mass hosts, and accreted at all redshifts, we find an average slope of the MZR of ≈ 0.48, which is steeper than the slope estimated for the disrupted populations of the MW (Naidu et al. 2022), of ≈ 0.3.On an individual basis, however, the MZR slopes of the disrupted dwarfs in each MW-mass system can vary significantly, and some of them can match the MW observations.We find better agreement with the observational data if Sagittarius and I'itoi are included in the observational sample (in this case the observational MZR would have a slope of ≈ 0.5).For the surviving populations, we find an excellent agreement between the predicted MZR slope, of ≈ 0.32, and the sloped measured from observations in the MW and Local Group (Kirby et al. 2013;Naidu et al. 2022), of 0.3.
• We investigated the evolution of the MZR of the disrupted population with redshift, and found that populations accreted at the similar redshifts ( acc ) follow MZRs with similar slopes as the surviving populations, of ≈ 0.3, and with lower normalizations towards higher redshift (Fig 3).The steeper slope obtained for the combined population of disrupted dwarfs (≈ 0.48), compared with the slope derived for these populations at fixed redshift (of ≈ 0.3), can be explained by the changes in the spectrum of stellar masses of dwarfs accreted towards lower redshift, and their ongoing metalenrichment.For disrupted dwarfs accreted most recently ( acc ≈ 0.5), their MZR is similar to that of surviving population.
• Disrupted dwarf satellites have significantly different chemical properties than the surviving ones: at fixed  * , they have lower [Fe/H] (Fig. 2) and higher [Mg/Fe] (Fig. 4); also, they have higher [Mg/Fe], at fixed [Fe/H].These trends are consistent with observations in the MW.
• The [Fe/H] - * and [Mg/Fe] - * distributions of disrupted dwarfs show a significant amount of scatter.We find that this scatter correlates with both  acc and the cold gas fractions at accretion,  SF (Figs. 5 and 6), of the progenitors of these dwarfs.Specifically, at fixed  * , progenitors accreted earlier have lower [Fe/H], higher [Mg/Fe] and also higher  SF than those accreted more recently.
• Generally, the cold gas fraction of the disrupted dwarfs correlates with their  acc , but not always.Typically, before  acc ≈ 4, dwarf galaxies are gas rich (with  SF roughly constant, of ≈ 80 − 100%).Below  ≈ 4, however,  SF decreases sharply (Fig. 7).A small fraction of the disrupted dwarfs in the simulations do not follow these general trends.These are low-mass systems with  SF = 0, and are likely an artifact of the limited numerical resolution.
• We investigated several global properties of the disrupted populations to test how predictive they are in terms of their host accretion history.Among these properties, we find that the mass-weighted values ⟨[Fe/H]⟩ and ⟨[Mg/Fe]⟩ of each disrupted population correlate strongly with the corresponding ⟨ acc ⟩ value.In contrast, the average slopes of the MZR or of the [Mg/Fe]- * relation in the disrupted population are not good tracers of the accretion history.Furthermore, we selected the six disrupted populations (the  6 * subsample) with similar ⟨[Fe/H]⟩ and ⟨[Mg/Fe]⟩ values as the ones in thee disrupted population in the MW (Fig. 8), and found that these were accreted, on average, at ⟨ acc ⟩ ≃ 1.8, which is in agreement with results from observations (Fig. 9).
• The disrupted dwarfs differ from the surviving counterparts also in terms of their stellar mass function.For surviving dwarfs, there are typically more low-mass dwarfs than high mass, as expected.However, for the combined disrupted population (i.e., accreted at all  acc ), the distribution is more evenly spread.This supports the differences seen in the global MZRs and other distributions between the two populations.There is also an indication that galaxies with chemical properties like the MW (i.e., the  6 * subsample) may have a higher fraction of low-mass dwarfs in their disrupted populations than other galaxies of MW-mass (Fig. 10).
• Simulations also support the scenario in which the disrupted dwarfs of the MW form and evolve in denser environments than the surviving ones (Figs. 11).This is because the progenitors of disrupted dwarfs form closer to their hosts (Fig. 13).Differences in the environments and locations of the progenitors of the two populations are more predominant at high redshift (e.g. = 2) than at lower one ( = 1).These differences are found in all MW-mass simulated galaxies.A1.Main properties associated to the MW and to the six MW-mass hosts in the  6 * subsample shown in Fig. A1): the slopes and offsets derived from the fits to the MZR of the disrupted and surviving populations, and their corresponding Spearman coefficients,   .Also included are the stellar mass weighted values ⟨ acc ⟩ and ⟨ [Fe/H] ⟩ of the disrupted populations, the latter being computed at  acc (see Section 3.4).For the MW we compute ⟨ acc ⟩ and ⟨ [Fe/H] ⟩ for the disrupted populations using data from Naidu et al. (2022), both including I'itoi and Sagittarius (MW) and excluding them (MW * ).The parameters of the fit for the MW disrupted populations are from Naidu et al. (2022), but we correct the offset to the same units used for the simulations.For the surviving populations in the MW we include the fit from Kirby et al. (2013).A1, but for the fits to the [Mg/Fe] - * relations corresponding to the six hosts in the  6 * subsample and to the observations.The ⟨ [Mg/Fe] ⟩ values for the disrupted populations in the MW are computed using the data from Naidu et al. (2022).The slopes and offsets for the relations of the disrupted and surviving populations in the MW are also computed using data from Naidu et al. (2022).

Figure 1 .
Figure 1.The redshift of accretion ( acc ) of disrupted dwarfs versus the redshift when these dwarfs cross the  200 radius of its MW-mass host ( acc|R 200 ).Values are colour-coded by the maximum  * reached by the dwarfs, i.e., measured at  acc .The blue dashed line represents the 1:1 correspondence.In general, the two redshifts are similar for most of the disrupted dwarfs.However, for the majority of the most massive satellites,  acc <  acc|R 200 , which indicates that these continue to form stars for some time while within the  200 boundary.

Figure 2 .
Figure 2. Left: The MZR of disrupted (purple stars) and surviving (grey squares) populations of dwarf galaxies in all ARTEMIS Milky Way-mass hosts.The linear fits to the distributions are shown with solid purple and red lines, for disrupted and surviving populations, respectively.The dashed black line corresponds to the MZR fit of Kirby et al. (2013) for surviving satellites in the Local Group: [Fe/H] = 0.30 log(  * /10 6 M ⊙ ) − 1.69.The dashed blue line represents the fit from Naidu et al. (2022) for the disrupted population in the MW: [Fe/H] = 0.36 log(  * /10 6 M ⊙ ) − 2.11.Right: Individual fits to the disrupted dwarf galaxies in each of the simulated MW hosts (purple lines).The two fits for the disrupted population in the MW, from Naidu et al. (2022), are indicated with blue lines: the dashed blue line (which is identical with the one in the left panel) is derived from a sample which excludes I'itoi and Sagittarius satellites (see text), while the solid blue line includes these satellites in the fit.

Figure 3 .
Figure 3. Mean MZR for surviving dwarfs at  = 0 (grey line) and disrupted dwarfs accreted at different redshifts,  =  acc (purple lines).The translucent dashed purple line corresponds to the mean relation for all disrupted dwarfs.[Fe/H] are mean values computed at different  and in different  * bins.For disrupted dwarfs, the normalization of the MZR increases towards lower  acc , by about 0.2 − 0.25 dex per  acc bin.Although disrupted dwarfs in different  acc bins trace MZRs with similar slopes (≈ 0.27 − 0.29), the slope for the entire disrupted population is significantly steeper (≈ 0.48).
Fig. 4 we show the [Mg/Fe] - * distributions in the combined disrupted and surviving populations in the simulations (purple stars and grey squares, respectively), compared with the corresponding distributions in the MW populations (data are again from Naidu et al. 2022).The right panel of the same figure shows the corresponding [Mg/Fe] -[Fe/H] distributions.
i) the observed [Fe/H] and [Mg/Fe] values for surviving and disrupted dwarfs in the MW are matched by the distribution of corresponding values in the simulations; ii) the simulations retrieve the general trends seen in the MW observations, for example, the decrease in [Mg/Fe] with  * , or the

Figure 4 .
Figure 4. Left: [Mg/Fe] - * distributions for the two populations of dwarfs in simulations, disrupted (with purple stars) and surviving (with grey squares).For comparison, the corresponding populations of the MW from Naidu et al. (2022) are shown, disrupted (with blue stars) and surviving (with black squares).At fixed  * , disrupted dwarfs tend to be more enhanced in [Mg/Fe] than the surviving ones, in both simulations and observations.Right: A similar comparison for the [Mg/Fe] -[Fe/H] distribution.At fixed [Fe/H], disrupted dwarfs tend to have higher [Mg/Fe] than surviving ones, again, in both simulations and observations.

Figure 5 .
Figure 5.The [Fe/H] - * distributions of the disrupted (top panels) and surviving (bottom panels) populations.The symbols for the disrupted dwarfs are colour-coded by the redshift of accretion,  acc (left) and by the star-forming gas fraction,  SF (right), while for the surviving ones, they are colour-coded by the redshift of reaching half of maximum stellar mass,  50 (left) and by  SF (right), respectively.Mean [Fe/H],  * and  SF values are computed at  acc for the disrupted satellites, and at  = 0 for the surviving ones, respectively.A small population of disrupted dwarfs has no star-forming gas at accretion,  SF = 0 (empty circles).For the disrupted dwarfs, a clear correlation is seen between the scatter in [Fe/H] at fixed  * and  acc or  SF .Correlations are very weak for the surviving populations.

Figure 6 .
Figure 6.Top: [Mg/Fe] - * distributions for the disrupted populations, with values colour-coded by  acc (left) and by  SF at accretion (right).Systems with  SF = 0 are shown with empty circles.A clear correlation of the scatter in [Mg/Fe] at fixed  * is seen for the disrupted dwarfs, with both  acc and  SF .Bottom:The same distributions for the surviving dwarfs, colour-coded by  50 (left) and by  SF at  = 0 (right).The correlations of the scatter with these parameters are much weaker in this case.

Figure 7 .
Figure 7. Star-forming gas fraction,  SF , versus redshift of accretion,  acc , for disrupted dwarfs.Values are colour-coded by [Mg/Fe] (top) and by [Fe/H] (bottom), respectively.The  SF values tend to decrease with redshift, with a sharp turnover around  acc ≈ 4. Before this redshift, there is a considerable scatter in the relation between these parameters, however after this redshift the two parameters correlate well.

Figure 8 .
Figure 8. Top: The ⟨ [Mg/Fe] ⟩ -⟨ [Fe/H] ⟩ relation of disrupted populations in the simulations (filled circles).Symbols represent the mass averages of each disrupted population in the 42 simulated MW-mass systems.⟨ [Mg/Fe] ⟩ anti-correlates with ⟨ [Fe/H] ⟩.The best fit is a linear relation with a slope of −0.40 ± 0.015.The blue and black squares correspond to the samples of disrupted dwarfs in the MW: including I'itoi and Sagittarius (MW sample) and excluding them (MW* sample), respectively.The six closest populations to the MW* data in this relation are shown with green circles and called the " 6 * subsample."Bottom: Same relation as above, colour-coded by the ⟨ acc ⟩ of each disrupted population.There is a strong anti-correlation between ⟨ acc ⟩ and ⟨ [Fe/H] ⟩ (recently accreted populations are more metal-rich) and a strong correlation between ⟨ acc ⟩ and ⟨ [Mg/Fe] ⟩ (recently accreted populations have lower [Mg/Fe] abundances).

Figure 9 .
Figure 9. Scaling relations between ⟨ acc ⟩ and global parameters of disrupted dwarf populations in the simulated MW-mass hosts.Top: ⟨ [Fe/H] ⟩-⟨ acc ⟩ relation.Middle: the MZR slope versus ⟨ acc ⟩.Bottom: the slope of the [Mg/Fe] - * relation versus ⟨ acc ⟩. Green circles correspond to the disrupted populations in the  6 * subsample, while purple circles correspond to those of the rest of the simulated populations.Values inferred from the MW disrupted population, both including and excluding I'itoi and Sagittarius, are shown with blue and black squares, respectively.The black dashed lines represent the linear fits for the simulated relations (from top to bottom, the corresponding   values are: −0.95, −0.11 and 0.30, respectively).

Figure 10 .
Figure10.Top: Stellar mass distribution of the progenitors of disrupted dwarfs in the  6 * subsample (solid green line) and in the remaining simulated hosts (black dashed line).Bottom: Same distributions but for the surviving dwarf galaxies.The  6 * hosts have slightly higher fractions of lower mass progenitors than the rest of the sample, for both disrupted and surviving populations.

Figure 12 .
Figure 12.Evolution of the median tidal parameter Θ 5 versus redshift, and their associated 25 and 75 percentiles around the median for the populations in all simulated hosts (purple line and shaded region for disrupted dwarfs, in the left panel; and grey line and shaded region for the surviving dwarfs, in the right panel).Superimposed are the corresponding Θ 5 medians for the two populations in the  6 * subsample.These are indicated with dashed lines for disrupted populations and full lines for the surviving ones.Lines for individual hosts in the  6 * subsample are shown with different colours, as stated in the legend.The evolution of Θ 5 is similar in the case of  6 * subsample as in the entire sample.

Figure 13 .
Figure13.Top: Distributions of distances between the progenitors of dwarf galaxies (disrupted dwarfs with purple lines and surviving ones with grey lines, respectively) and the centers of potential of the corresponding progenitors of their MW-mass hosts, at  = 2 (left) and at  = 1 (right), respectively.The black dashed lines represent the median  200 of the progenitors of MW-mass systems at the corresponding redshift, while the shaded regions are the 14 th and 86 th percentiles of these values.Dwarf satellites which are disrupted today were closer to their hosts than the surviving ones, both at  = 1 and at  = 2. Bottom: The same distributions as in the top panels, but for the populations in the six hosts in the  6 * subsample.The differences between the locations of the two satellite populations are similar to those in the full sample, being slightly more pronounced in this case.

Figure 14 .
Figure 14.Left: Disruption time,  dis , versus  * for disrupted dwarfs, with symbols colour-coded according to  acc .The most recently accreted dwarfs tend to have higher masses and greater  dis values.Right:  dis versus lookback time of accretion,  acc,L , for the disrupted dwarfs.Values are colour-coded by  * .

Figure A1 .
Figure A1.MZRs for the disrupted (purple stars) and surviving (grey squares) populations in the  6 * hosts, and the corresponding MW populations studied by Naidu et al. (2022) (blue stars and black squares, respectively).The linear fits for the simulated dwarfs are shown with purple and grey solid lines, respectively.The dashed blue and black lines correspond to similar fits to the MW dwarf populations, where the disrupted sample excludes I'itoi and Sagittarius.

Figure A2 .
Figure A2.[Mg/Fe]- * relation for the disrupted (purple stars) and surviving (grey squares) populations in the  6 * hosts and the corresponding MW populations studied by Naidu et al. (2022) (shown with blue stars and black squares, respectively).The linear fits for the simulated disrupted and surviving populations are shown with purple and grey solid lines, respectively.The dashed blue and black lines correspond to similar fits to the MW dwarf populations.For the observational disrupted sample we exclude I'itoi and Sagittarius (i.e., the MW* sample).

Table A2 .
Similar to Table