Galaxies with Shells in the Illustris Simulation: Metallicity Signatures

Stellar shells are low surface brightness arcs of overdense stellar regions, extending to large galactocentric distances. In a companion study, we identified 39 shell galaxies in a sample of 220 massive ellipticals ($\mathrm{M}_{\mathrm{200crit}}>6\times10^{12}\,\mathrm{M}_\odot$) from the Illustris cosmological simulation. We used stellar history catalogs to trace the history of each individual star particle inside the shell substructures, and we found that shells in high-mass galaxies form through mergers with massive satellites (stellar mass ratios $\mu_{\mathrm{stars}}\gtrsim1:10$). Using the same sample of shell galaxies, the current study extends the stellar history catalogs in order to investigate the metallicity of stellar shells around massive galaxies. Our results indicate that outer shells are often times more metal-rich than the surrounding stellar material in a galaxy's halo. For a galaxy with two different satellites forming $z=0$ shells, we find a significant difference in the metallicity of the shells produced by each progenitor. We also find that shell galaxies have higher mass-weighted logarithmic metallicities ([Z/H]) at $2$-$4\,\mathrm{R}_{\mathrm{eff}}$ compared to galaxies without shells. Our results indicate that observations comparing the metallicities of stars in tidal features, such as shells, to the average metallicities in the stellar halo can provide information about the assembly histories of galaxies.


Introduction
Stellar halos are diffuse, low surface brightness regions that contain unique morphological features, such as tidal streams, stellar shells, rings, or plumes (e.g., [1,2]). Shells are a particular type of tidal debris that appear as interleaved caustics with large opening angles, which have been identified in numerous deep surveys probing the faint outskirts of galaxies [3][4][5][6][7][8][9][10]. Interest in the detailed study of stellar substructures such as shells has been stimulated by previous work indicating that tidal features are powerful tracers of galaxies' assembly histories (see e.g., [11][12][13][14][15][16]). Due to the long dynamical timescales in the outskirts of galaxies, information about the accretion histories of galaxies is well preserved at the present epoch. Previous studies indicate that stellar shells are very long-lived structures, with average lifetimes of ∼2-3 Gyr according to a recent paper [17].
Several theories have been proposed for the formation of shells, ranging from models invoking star formation in shocked galactic winds [18][19][20][21] to tidal interaction theories in which the shells are density waves induced by a passing galaxy [22,23]. Nevertheless, extensive numerical and observational studies support the merger scenario, in which shells are composed of stars stripped from accreted satellites (e.g., [24][25][26][27][28][29]). The increasing richness of stellar halo observations will give us a unique look into the past merger histories of galaxies, yet the interpretation of the data requires detailed model predictions for the formation of stellar substructures. Most of the previous studies on shell galaxies have focused on idealized mergers involving small satellites on radial orbits, (e.g., [25][26][27][28][30][31][32][33]), yet major mergers have also been shown to produce shells, (e.g., [34][35][36][37]). In a companion paper, Pop et al. ([17], hereafter P17) found that shells in high-mass galaxies are often the result of mergers with massive satellites (relative to the mass of the host galaxy). Since the efficiency of dynamical friction is higher for major mergers compared to minor mergers, massive satellites can form shells for a wide range of impact parameters; small satellites, however, need almost perfectly radial orbits in order to produce stellar shells. As a result, in P17, we show that most satellites producing shells in massive ellipticals have stellar mass ratios µ stars 1 : 10 .
Previous studies suggest that metallicity gradients can be used to characterize the relative importance of major and minor mergers in the accretion history of a galaxy (e.g., [38,39]). During the initial formation phase of massive galaxies, in-situ star formation will leave a steep metallicity gradient, with stars on the outskirts of galaxies forming out of lower density gas and metals being further removed from the outer regions of the stellar halo by stellar winds [39,40]. Overall, subsequent accretion events (z 1) have the net effect of flattening the metallicity profiles (see e.g., [41][42][43]). In a previous study of Illustris galaxies, Cook et al. [44] found that galaxies with higher fractions of accreted material tend to have less steep gradients than galaxies with a lower fraction of ex-situ stars. Based on these results, it has been suggested that metallicity gradients could provide additional information about the mass ratios of mergers that produce shells. Several groups have found that shell galaxies tend to have lower central Mg 2 values than non-shell galaxies (e.g., [45,46]), with galaxies with shells lying below the Mg 2 -σ relation from [47].
Studies of shell galaxies using deep imaging face several challenges, such as the extremely low surface brightness of shells (∼28 mag arcsec −2 or fainter in the V band, e.g., [8,48]), the varied and sometimes irregular shell morphologies, and the fact that outer shells extend far out in the outskirts of the stellar halo, at galactocentric distances ∼100 kpc. For sufficiently nearby elliptical galaxies (≤20 Mpc), direct stellar photometry has successfully provided strong constraints on the age and metallicity of stars extending to large galactocentric distances in shell galaxies (e.g., [49][50][51][52]). In this paper, we investigate the metallicity of stellar shells in the Illustris simulation, studying how different metallicity signatures can set constraints on the number and mass of shell-forming progenitors. As metallicity measurements are reaching larger radii through slit spectroscopy (e.g., [53][54][55][56][57]), integral-field spectroscopy (e.g., [58][59][60][61][62]), and deep photometric studies (e.g., [51,52,63,64]), as well as probing increasingly larger samples of galaxies [55,65], the study of stellar metallicities of tidal features will provide previously inaccessible information about the galaxies' accretion histories.

Illustris
This project uses the Illustris simulation [66][67][68], a cosmological hydrodynamical simulation run using AREPO [69]-a moving-mesh code based on a quasi-Lagrangian finite volume method. We are working with the highest resolution run (Illustris-1), which corresponds to a simulation box with a periodic volume of (106.5 Mpc) 3 and mass resolution of m DM = 6.26 × 10 6 M . The Illustris simulations include a broad range of astrophysical processes (e.g., feedback from supermassive black holes, supernovae, and AGNs) in order to self-consistently reproduce galaxy formation and evolution [70,71].

The Galaxy Sample
Similar to P17, we start with an initial sample consisting of the 220 most massive central galaxies in Illustris at redshift z = 0. This corresponds to a virial mass cut M 200crit > 6 × 10 12 M , where M 200crit is defined as the total mass inside a radius enclosing a sphere with mean density 200 times greater than the critical density of the Universe. In Illustris, central galaxies are the most massive galaxies in a given halo/group, and satellites inside each halo are identified using the SUBFIND algorithm [72][73][74].

Stellar History Catalogs and Shell Galaxies Identification
The present study uses the Illustris shell galaxies previously identified in P17 using a two-step approach:

•
Step 1: Galaxies with stellar shells are visually identified using stellar surface density maps.
Each galaxy in the sample is studied using all three projections (x-y, y-z, z-x) and two different contrast levels, and each image stamp received a score ranging from 0-2 from three different members of our team. A score of 2 indicates a galaxy with two or more well-defined shells, a score of 1 corresponds to galaxies with one or two shell-like structures, and a score of 0 indicates no shell detection. We order candidate shell galaxies based on the total scores given to all their corresponding image stamps. As discussed in P17, most shell galaxies exhibit shell-like structures in at least 2/3 projections, and we use the second identification step to verify the presence of shells in each of the candidate galaxies.

•
Step 2: We develop stellar history catalogs in order to identify the shell-forming progenitors and separate stars in shells from all other stars in the galaxy. In these catalogs, we trace the birth, trajectory and progenitors of all star particles inside z = 0 halos, saving information about each individual star particle at three key moments during its life: formation, accretion (when the parent satellite enters the virial radius of the host) and stripping time (when the star becomes gravitationally bound to the new host). In P17, we show that shells in Illustris are composed of ex-situ stars, which is in agreement with merger models predicting that shell-like structures correspond to overdensities of stripped stars accumulating at the apocenters of their orbits (see analytical and numerical studies by, e.g., [24][25][26][27][28][29]34,75,76]). By tracing the configuration and phase space (v r vs. r) distribution of star particles with a common parent satellite, we obtain a systematic sample of all the satellites that are responsible for z = 0 stellar shells.
Out of 220 massive ellipticals in Illustris (with M 200crit > 6 × 10 12 M ), Pop et al. [17] find that 39 galaxies exhibit stellar shells visible at z = 0, corresponding to a fraction of 18% ± 3% of all galaxies in the sample. 1 Moreover, the fraction of galaxies with shells increases with increasing mass cut, with as many as 34% ± 11% of galaxies with M 200crit > 3 × 10 13 M having shells. Despite significant uncertainties in the number of observed shell galaxies, the incidence of galaxies with shells in Illustris is in overall agreement with current observational limits (e.g., [4][5][6][7][8]).
In this paper, we take advantage of the versatility of stellar history catalogs, which can be extended to study a wide range of properties for stars inside tidal features like shells, streams or plumes (for more details about the construction of the catalogs, see P17). In particular, the current study focuses on the metallicity signatures of several shell galaxies from the Illustris simulation. We use stellar surface density maps such as those presented in Figure 1, in which we weigh each star particle in the z = 0 halo by its total metallicity (Z), normalized to that of the Sun (Z ). Each bin in the 2D histograms in Figure 1 is colored based on the median total metallicity of all star particles therein. In Illustris, the metallicity of each individual star particle is based on the metallicity of the gas cell at the star's time of birth (i.e., when the gas cell is converted into a star particle). 1 We provide Poisson errors for the fractions of galaxies with shells in Illustris.
However, the actual uncertainties in the f shells could be higher due to environment (rich groups and clusters vs. field galaxies), mass distribution of galaxies in the sample, redshift evolution f shells (z), surface brightness limits, and projection effects. For more details about the impact of these effects, see the discussion in Pop et al. [17]. Figure 1. Stellar maps in configuration space (first three columns) and phase space, i.e., radial velocity vs. radial distance (last column). The color of each bin of area A bin = (1.5 kpc) 2 in the 2D histograms above corresponds to the median total metallicity (Z) of all star particles therein. Top row: entire halo of a redshift z = 0 shell galaxy with total stellar mass M stars = 6.3 × 10 11 M , virial mass M 200crit = 1.4 × 10 13 M , and effective radius R eff = 16.3 kpc (where R eff is defined as the radius that contains half of the total light in a galaxy). Bottom row: redshift z = 0 stars that have been accreted from the same common progenitor, which is responsible for the shell structures in the top row. We find that stellar shells in this galaxy have higher total metallicities than stars in the outer regions of the stellar halo. Moreover, we find a gradient in the average metallicity of individual shells, with outer shells having lower average Z that inner shells.
In configuration space, we identify shells as interleaved stellar overdensities, often observed on both sides of the galaxy center and extending to very large galactocentric distances, in the low surface brightness regions of the stellar halo. Stellar shells also have a specific signature in phase space, i.e., in the space of galactocentric distance vs. radial velocity (r-v r ). Stars spend longer times near the apocenters of their orbits, where their radial velocities are close to zero, thus creating phase wraps peaking at v r 0, as exemplified in the right column of Figure 1. Using stellar history catalogs, we identify satellites responsible for z = 0 shells, and we investigate the metallicity distribution of stars inside these shells (see bottom rows of Figures 1 and 2). . Stellar maps in configuration space (first three columns) and phase space, i.e., radial velocity vs. radial distance (last column). The color of each bin of area A bin = (1.5 kpc) 2 in the 2D histograms above corresponds to the median total metallicity (Z) of all star particles therein. Top row: entire halo of a z = 0 shell galaxy with total stellar mass M stars = 1.3 × 10 12 M , virial mass M 200crit = 4.6 × 10 13 M , and R eff = 33.5 kpc. This galaxy had two different shell-forming progenitors responsible for the shells visible at z = 0. Middle and bottom rows: redshift z = 0 stars accreted from each of the two shell-forming satellites with stellar mass ratios µ stars = M stars satellite /M stars host = 0.15 and 0.24, measured at accretion time (i.e., when the satellite comes inside the virial radius of the host). Stellar shells in this galaxy have higher average metallicity than the other stars in the outskirts of the stellar halo. Moreover, the shells produced by the first satellite are on average more metal-rich than those produced by the second satellite. Similar to Figure 1, we find that the outer shells have lower total metallicities that inner shells.

Results and Implications
We begin by investigating the distribution of stellar metallicities for two shell galaxies, one for which the z = 0 shells are composed of stars accreted from a single satellite, while the shells in the second galaxy are the result of mergers with two different satellites. In Figure 1, we study the configuration and phase space distribution of all star particles inside the halo of a galaxy with total stellar mass M stars (z = 0) = 6.3 × 10 11 M and virial mass M 200crit (z = 0) = 1.4 × 10 13 M . The shell-forming satellite depicted in the bottom row corresponds to a close to 1 : 1 merger (in total mass) and a stellar mass ratio of µ stars = M stars satellite /M stars host = 0.22 . The satellite was accreted ∼4 Gyr ago on a low angular momentum orbit, with a radial velocity ratio 2 at accretion time of v r /|v| = −0.97. 2 The radial velocity ratio (v r /|v|) is defined as the radial component of the relative velocity of the satellite (v r ), normalized to the modulus of the relative velocity (v). This shell-forming satellite has properties that agree with the overall trends for shell formation found in P17: most of the shell galaxies in the current sample form through mergers with stellar mass ratios µ stars 1 : 10, and the shell-forming satellites are accreted on low angular momentum orbits, at intermediate accretion times (∼4-8 Gyr ago). We determine the color of each bin in the stellar 2D histograms in Figure 1 by computing the median of the total metallicity (Z) for all star particles in that bin. We find metal-rich shells visible in all three projections and that extend to large galactocentric distances (>50 kpc). While in the top panels of Figure 1 we determine the color of each bin by considering both shell stars and other halo stars that are inside the same spatial bin, in the bottom panels, we show only those stars accreted from the shell-forming satellite (identified as described in Section 2.3). Thus, since in the bottom row we are isolating the stellar shells, it becomes more evident that, for the galaxy in Figure 1, stars in shells are more metal-rich than the overall stellar population of the host galaxy.
For the second study case, we consider a shell galaxy that exhibits z = 0 shells composed of stars accreted from two different satellites. The host galaxy in Figure 2 has a total stellar mass M stars = 1.3 × 10 12 M and virial mass M 200crit = 4.6 × 10 13 M . As before, we color each bin based on the median total metallicity of all star particles therein, and we find significantly more stellar shells at large radii, compared to the example in Figure 1. However, similarly to the previous study case, stars in shells have higher total metallicities than the other halo stars at similarly large radii. The two shell-forming satellites in Figure 2 were accreted between ∼6-7 Gyr ago, and they were stripped ∼3 and ∼2 Gyr ago, for satellites 1 and 2, respectively. Both of these merger events correspond to relatively major mergers with µ stars(t acc ) = 0.15 and 0.24, and the satellites had radial infall trajectories, with v r /|v| = −0.99 and −0.84 at t acc . The middle and bottom rows of Figure 2 depict only those stars accreted from each of these two satellites. On average, both satellites contribute stars with higher total metallicities than the average metallicity of stars in the host galaxy. Moreover, we find that the shells created by the first satellite are more metal-rich than those created by the second satellite. This opens an interesting avenue for future observational studies-measurements of the metallicities of stars in shells could be used to identify shell galaxies with more than one shell-forming progenitor.
In both Figures 1 and 2, we find a gradient in the metallicity of stars in shells, with outer shells having lower average Z than inner shells. Studies of the satellites' infall trajectories show that shells in the current sample are often composed of stars stripped while the parent satellite passes close to the center of the host several times (see P17). Stars stripped during the satellite's first pericenter passage form the first generation of shells, with subsequent stripping events forming several generations of shells that are located at decreasing galactocentric distances compared to previous generations, (e.g., [31,32,77,78]). As a result, we expect the outer-most shells to be composed of stars that were initially located on the outskirts of the parent satellite, and, thus, they had lower average metallicities than stars closer to the satellite's core. This explains the trends observed in Figures 1 and 2, with the first generation of shells having lower metallicities than shells located closer to the galaxy center.
In Figure 3, we show total metallicity Z-weighted stellar surface density maps for nine galaxies with shells in our sample. The virial masses of the host galaxies at z = 0 span one order of magnitude: M 200crit ∈ 6.5 × 10 12 M , 6.3 × 10 13 M , and total stellar masses are in the range M stars ∈ 3.1 × 10 11 M , 2.2 × 10 12 M . The examples presented in Figure 3 represent a selection of shell galaxies in Illustris that have outer shells more metal-rich than the average metallicity of surrounding stars at similar galactocentric distances in the host galaxy. Nonetheless, we do not exclude that some galaxies could have shells composed of stars with similar or even lower metallicities than halo stars located at similar radii. Figure 3. Examples of nine different massive elliptical galaxies with shells from the Illustris simulation. The color of each bin of area A bin = (1.5 kpc) 2 in the 2D histograms above corresponds to the median total metallicity (Z) of all star particles therein. The total stellar masses of these galaxies range from M stars = 3.1 × 10 11 M to M stars = 2.2 × 10 12 M , while the virial masses are in the range M 200crit ∈ 6.5 × 10 12 M , 6.3 × 10 13 M . For several of the shell galaxies in our sample, stellar shells are metal-rich compared to the total metallicity (Z) of stars in the outer stellar halo.
In order to estimate the difference in total metallicity Z between outer shells and stars at similar galactocentric distances, and the resolution levels necessary to detect metal-rich shells such as those presented in Figure 3, we next study a series of Z-weighted stellar maps at increasingly lower resolutions. In the top row of Figure 4, we show the (y-z) projection of the same galaxy previously presented in Figure 2, and we vary the area of the 2D histogram bins from A bin = (1 kpc) 2 to A bin = (5 kpc) 2 , (10 kpc) 2 , and finally, (15 kpc) 2 . This galaxy has an effective radius R eff = 33.5 kpc, so the right-most stellar map in Figure 4 corresponds to a low-resolution image, with the length of each bin being L bin ≈ R eff /2 . Some of the fine structure associated with tidal debris is completely erased for A bin > (5 kpc) 2 , and the number of detectable shells drops significantly as we lower the resolution. Nonetheless, there is sufficiently high contrast between the metallicities of some of the outer shells and the metallicities of stars located at similar radii in the halo, such that several outer shells can still be observed even for A bin (15 kpc) 2 . In the bottom row of Figure 4, we present a larger version of the image stamp with A bin = (10 kpc) 2 , on which we mark six different targets, each having an area A target = (50 kpc) 2 . In the table accompanying the figure, we provide galactocentric distances for the centers of each of the six targets in units of kpc and R eff , as well as the median log 10 (Z/Z ) computed for all star particles located inside each target area 3 . The median log 10 (Z/Z ) is highest for target area 1, which is centered on the core of the galaxy. Areas 2-4 are targeting three different shells located at galactocentric distances: R 2 = 4.2 R eff , R 3 = 8.1 R eff , and R 4 = 7.1 R eff , respectively. The median log 10 (Z/Z ) of all stars in these target areas varies between log 10 (Z/Z ) = 0.18 to −0.021. Despite the median total metallicities of stars in outer shells being lower than the median Z of stars near the galaxy center (where log 10 (Z/Z ) = 0.33), there is still a significant contrast in metallicity between targets at similar radial distances depending on whether they are centered on a shell or not. For example, for target 2, we measure a median log 10 (Z/Z ) = 0.18, while target 5 is positioned at the same galactocentric distance (R = 4.2 R eff ), in a region of the stellar halo where we detect no tidal features, and it has a median log 10 (Z/Z ) = −0.22. This corresponds to a metallicity contrast ∆log 10 (Z/Z ) = 0.4 between targets 2 and 5. We can similarly compare targets 3 and 6, with the former being centered on a shell in the upper left quadrant and having median log 10 (Z/Z ) = 0.036, while target 6 is positioned at a similar radial distance (R = 8.1 R eff ) but in a low stellar density region of the halo, where the median log 10 (Z/Z ) = −0.64, leading to a contrast ∆log 10 (Z/Z ) ≈ 0.7 between targets 3 and 6. In addition, we find that the median total metallicities of outer shells (log 10 (Z/Z ) = 0.036 and −0.021 for targets 3 and 4, respectively) are lower than the median metallicity for the shell located at a smaller radius, near target number 2 (log 10 (Z/Z ) = 0.18). These results are in agreement with the metallicity gradients in shells observed in Figures 1 and 2: outer shells are less metal-rich than inner shells. Due to the extremely low surface brightness of shells, reaching the precision levels above is challenging for current observations, though not technically impossible.
Next, we quantify the difference in the median stellar metallicities of shell vs. non-shell galaxies by measuring the mass-weighted logarithmic metallicities [Z/H] of stars around all 220 galaxies in our sample (see Figure 5). Following the definitions adopted in [44], we measure [Z/H] within three radial ranges: the inner galaxy (0.1-1 R eff ), outer galaxy (1-2 R eff ), and stellar halo (2-4 R eff ). The average effective radius of the entire sample is 20.9 kpc, with shell galaxies having a higher average (R eff = 24.5 kpc) than galaxies without shells (R eff = 20.1 kpc). We do not find a significant difference in the average metallicities of stars within ∼2 R eff of the galaxies' centers. However, in the outer regions of the stellar halo (2-4 R eff ), shell galaxies have higher average metallicities than galaxies without shells: the median values for the two distributions are [Z/H] shells = −0.28 and [Z/H] non-shells = −0.36 . In turn, if stars on the outskirts of shell galaxies are more metal-rich, this implies that the metallicitiy gradients of galaxies with shells could be shallower than those of non-shell galaxies. Despite the higher average [Z/H] for shell galaxies, we find that there is significant scatter in the metallicity distributions for galaxies in our sample. Most galaxies do not show azimuthal symmetry, and, therefore, by computing [Z/H] averaged over radial profiles, we are ignoring the relative influence of different stellar substructures present at similar radii in the halo. We thus expect measurements of stellar metallicities targeted towards regions with significant stellar overdensities or clumpy substructures (see, for example, Figure 4) to be better suited at identifying new shells and characterizing the mergers producing them than measurements of overall metallicity gradients or [Z/H] averaged over large radial bins.  Figure 2). The resolution of the image stamps decreases from left to right, with the areas of 2D histogram bins given by: A bin ∈ (1 kpc) 2 , (5 kpc) 2 , (10 kpc) 2 , (15 kpc) 2 . Bottom row: Zoom-in of the third image stamp in the top row, corresponding to a histogram with A bin = (10 kpc) 2 . We mark six different target regions, each having a cross-section A target = (50 kpc) 2 . Table: We provide the distance R from the center of each of the target areas to the center of the galaxy in units of kpc and R eff , respectively, as well as the median log 10 (Z/Z ) measured over all star particles inside each target area. Targets centered around outer shells (e.g., 2 and 3) have significantly higher median log 10 (Z/Z ) than targets at the same galactocentric distances but centered on regions without tidal features (e.g., targets 5 and 6).
In P17, we show that most merger events responsible for forming z = 0 shells in massive elliptical Illustris galaxies can be described by an order-zero recipe, involving only three parameters: stellar mass ratio (µ stars ), radial velocity ratio (v r /|v|), and accretion time (t acc ). Satellites that successfully form z = 0 shells correspond to mergers with stellar mass ratios µ stars 1 : 10 and are accreted on very low angular momentum orbits, about ∼4-8 Gyr ago. While the mean total mass ratio (µ total = M total satellite /M total host ) of shell-forming merger events in P17 is lower than the mean stellar mass ratio (µ stars ), both µ total and µ stars are biased towards close-to-major mergers. Due to dynamical friction being more efficient at radializing the orbits of satellites involved in close to 1 : 1 mergers (see, e.g., [79,80]), massive satellites are allowed to probe a wider range of impact parameters and still be successful at forming shells (see discussion in P17). On the other hand, in order to produce shells in massive ellipticals, small satellites need very fine-tuned, almost perfectly radial orbits. Thus, P17 find that shells in massive galaxies are produced more frequently by major or close-to-major mergers than by minor mergers. Figure 5. Histograms of the average logarithmic metallicities of stars in a sample of 220 massive elliptical galaxies. The shell galaxies identified in P17 are shown in red, with the other galaxies in the sample shown in blue. We measure [Z/H] within three radial ranges, based on the effective radius of each galaxy, R eff , defined as the radius that contains half of the total light. The median [Z/H] for all shell/non-shell galaxies is marked by a continuous line, while dashed lines mark the 10% and 90% percentiles. The average total stellar metallicities are similar for shell and non-shell galaxies at galactocentric distances < 2 R eff . In the extended stellar halo, we find that shell galaxies in our sample have higher average metallicities (shallower metallicity gradients) than galaxies without shells.
Previous studies have shown that major mergers can deposit material closer to the center of the host, unlike minor mergers which mostly contribute to the ex-situ stellar fraction at large galactocentric radii [80,81]. Major mergers with shell-forming satellites containing metal-rich stars will thus deposit high-metallicity material over a wide range of radii-from the inner regions of the halo out to distances 100 kpc. As a result, there are two competing effects that need to be considered when estimating the metallicity of outer shells relative to the metallicity of halo stars at similar radii. On one hand, the outer-most shells (or the first generation of shells) are composed of stars stripped from the outskirts of the satellite's stellar halo during its first pericenter approach, while subsequent shells at lower radii correspond to more metal-rich stars that were located closer to the center of the shell-forming progenitor. This has the effect of generating a metallicity gradient in the shells, with outer shells being more metal-poor than inner shells (as observed in Figures 1 and 2). On the other hand, shells extend to extremely large distances in the low surface brightness and metal-poor outskirts of the host's stellar halo. Despite the fact that the very first ∼1-2 shells at large radii are generally not very metal-rich, the next shells will start to contain relatively high-Z stars stripped from closer to the center of the satellite. After the merger event, these stars are located in shells extending to much larger effective radii in the host galaxy compared to their initial galactocentric distances with respect to the center of the shell-forming progenitor. Therefore, in the case of major mergers, many of the outer shells can be on average more metal-rich than halo stars located at similar radii in the host galaxy, and we observe this effect in several of the shell galaxies in our sample (see Figure 3).
When averaging stellar metallicities over large radial bins, we find that shell galaxies have higher median [Z/H] than non-shell galaxies in the outer regions of the stellar halo (2-4 R eff ). Hirschmann et al. [38] used zoom hydrodynamical simulations of 10 galaxies with masses above 6 × 10 12 M , probing a similar mass range as our current study. They found that systems that have undergone major mergers have significantly shallower metallicity gradients at galactocentric distances > 2 R eff compared to systems dominated by minor mergers. For Illustris galaxies, Cook et al. [44] showed that galaxies with high ex-situ fractions tend to have less steep metallicity gradients (in the radial range 2-4 R eff ) than those with smaller accreted fractions, yet the metallicity gradients are not strongly correlated with the mean merger mass-ratio. As a result, it can be difficult to draw conclusions based solely on the distribution of overall metallicity gradients in shell galaxies, and precision metallicity measurements targeted around individual shells can provide more information about the shell-forming progenitor.
Based on the results presented above, we expect even outer shells (formed from material stripped from the outskirts of satellites) to be relatively metal-rich if the shells are produced by a massive parent satellite. Therefore, observations of metal-rich shells could indicate that a high fraction of the shells observed in massive early-type galaxies have a major merger origin. In this way, precision metallicity measurements targeting stars in and around the shells can provide additional information about the formation processes of shell galaxies.

Summary
In the current study, we investigate the average total metallicity of stellar shells around massive elliptical galaxies from the Illustris simulation. We present several shell galaxies for which we identify shells in both configuration space (x-y, y-z, z-x) and phase space (r-v r ). Stars in the outer shells of these galaxies have higher total metallicities than other stars at similar radii in the host galaxies. This result suggests that the shells studied in this paper are the by-product of major mergers. Moreover, in the particular case of a galaxy with two different shell-forming progenitors, we find a significant difference in the metallicity of the shells produced by each of the two satellites. This indicates that high-precision metallicity measurements could potentially identify shell galaxies with multiple progenitors based on their different metallicity signatures.
Previous analytical and numerical studies (e.g., [78,82]) have indicated that several generations of shells are necessary to explain the wide range of shell radii measured in observations of early-type galaxies [83]. In agreement with these results, we find that shells in Illustris are often formed by stripping stars while the parent satellite performs several pericenter passages. Stars in the first generation of shells (the outer-most shells) correspond to the least bound stars on the outskirts of the progenitor galaxy and, thus, they are the first stars stripped from the satellite. We find a gradient in the average metallicity of shells, with outer shells in our two study cases having lower metallicities than inner shells. Since stars on the outskirts of galaxies tend to have lower metallicities, first generation shells would be expected to be more metal-poor than subsequent generations, in agreement with our findings.
Finally, we investigate the average mass-weighted metallicities [Z/H] in three different radial bins corresponding to the inner galaxy, outer galaxy and stellar halo. Both galaxies with and without shells in our sample seem to have a similar average [Z/H] within 2 R eff of the galaxy center. We find, however, a bias towards a higher average [Z/H] for stars in the outer stellar halos of shell galaxies.
The results presented in this paper indicate that several of the shell galaxies in Illustris have outer shells more metal-rich than the average metallicity in the stellar halo at similar radii. We predict that multiple shell-forming progenitors, as well as different generations of shells, could leave specific signatures in the metallicity distribution of stars in shells. In addition to having the potential to find previously undetected shells, future metallicity observations could unveil information about the formation processes of shell galaxies, such as distinguishing between minor/major mergers producing shells or the number of different satellites responsible for the observed structure.
Author Contributions: A.R.P., A.P., N.C.A., and L.H. designed the research; A.R.P, A.P., and N.C.A. visually identified shell galaxies using stellar surface density maps; A.R.P. developed stellar history catalogs and post-processing tools, identified shell-forming progenitors and performed the numerical analysis; A.R.P. wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.