The MASSIVE Survey - XII Connecting Stellar Populations of Early-Type Galaxies to Kinematics and Environment

TL;DR

This study analyzes stellar populations in 90 early-type galaxies, revealing how their properties relate to galaxy kinematics and environment, with implications for galaxy formation and evolution.

Contribution

It provides detailed measurements of stellar populations out to large radii in massive ETGs and links these to kinematic and environmental factors, highlighting new correlations.

Findings

Abundance ratios correlate with sigma and M* at large radii.

Higher radial anisotropy is associated with older, metal-poorer stars.

Environmental density influences stellar population properties.

Abstract

We measure the stellar populations as a function of radius for 90 early-type galaxies (ETGs) in the MASSIVE survey, a volume-limited integral-field spectroscopic (IFS) galaxy survey targeting all northern-sky ETGs with absolute K-band magnitude M_K < -25.3 mag, or stellar mass M* 4x10^11 M_sun, within 108 Mpc. We are able to measure reliable stellar population parameters for individual galaxies out to 10-20 kpc (1-3 R_e) depending on the galaxy. Focusing on ~R_e (~10 kpc), we find significant correlations between the abundance ratios, sigma, and M* at large radius, but we also find that the abundance ratios saturate in the highest-mass bin. We see a strong correlation between the kurtosis of the line of sight velocity distribution (h4) and the stellar population parameters beyond R_e. Galaxies with higher radial anisotropy appear to be older, with metal-poorer stars and enhanced [alpha/Fe]. We suggest that the higher radial anisotropy may derive from more accretion of small satellites. Finally, we see some evidence for correlations between environmental metrics (measured locally and on >5 Mpc scales) and the stellar populations, as expected if satellites are quenched earlier in denser environments.

Figures

Tables

Table 1. Table 1: Galaxy Gradients

Gal	$M^{\ast }$	$\sigma$	$R_e$	$\Delta {\Sigma }_{23}$	[Fe/H]c	[Fe/H]o	$\Delta$[Fe/H]/$\Delta$ log $R$	[$\alpha$/Fe]c	[$\alpha$/Fe]o	$\Delta$[$\alpha$/Fe]/$\Delta$ log $R$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
NGC0057	11.8	251	6.3	$-1.93\pm$0.10	$-0.07\pm$0.06	$-0.24\pm$0.06	$-0.27\pm$0.05	$0.38\pm$0.05	$0.36\pm$0.05	$-0.034\pm$0.05
NGC0080	11.8	222	8.4	$-2.30\pm$0.18	$-0.10\pm$0.05	$-0.31\pm$0.05	$-0.11\pm$0.11	$0.31\pm$0.03	$0.27\pm$0.04	$0.077\pm$0.07
NGC0315	12.0	341	9.2	$-2.67\pm$0.08	$-0.14\pm$0.04	$-0.16\pm$0.03	$-0.04\pm$0.06	$0.34\pm$0.03	$0.28\pm$0.03	$-0.079\pm$0.06
NGC0383	11.8	257	8.0	$-1.91\pm$0.05	$-0.13\pm$0.05	$-0.12\pm$0.05	$-0.06\pm$0.27	$0.42\pm$0.05	$0.38\pm$0.04	$-0.058\pm$0.09

Table 2. Table 2: Structural Correlations

Stellar Pop	Radius	Gal. Prop.	$N_g$	$\rho$	$P$
(1)	(2)	(3)	(4)	(5)	(6)
[$\alpha$/Fe]		$\sigma$	87	0.4	$4\times {10}^{-5}$
[$\alpha$/Fe]	$0.9R_e$	$\sigma$	74	0.5	$1\times {10}^{-5}$
[$\alpha$/Fe]	9kpc	$\sigma$	73	0.4	0.00018
[$\alpha$/Fe]	1.5$R_e$	$\sigma$	58	0.6	$4\times {10}^{-6}$
[$\alpha$/Fe]	15kpc	$\sigma$	45	0.4	0.003
[$\alpha$/Fe]		$M^{\ast }$	87	0.3	0.02
[$\alpha$/Fe]	$0.9R_e$	$M^{\ast }$	74	0.4	0.0003
[$\alpha$/Fe]	9kpc	$M^{\ast }$	73	0.4	0.0003
[$\alpha$/Fe]	15kpc	$M^{\ast }$	45	0.4	0.02
[$\alpha$/Fe]	kpc	$h4$	87	0.4	0.0002
[$\alpha$/Fe]	$0.9R_e$	$h4$	74	0.5	$2\times {10}^{-5}$
[$\alpha$/Fe]	9 kpc	$h4$	73	0.4	0.0001
[$\alpha$/Fe]	$1.5R_e$	$h4$	74	0.5	$7\times {10}^{-5}$
[$\alpha$/Fe]	15 kpc	$h4$	73	0.4	0.003
[Fe/H]	kpc	$M^{\ast }$	87	0.2	0.04
[Fe/H]	$0.9R_e$	$M^{\ast }/R$	74	0.3	0.006
[Fe/H]	$0.9R_e$	$M^{\ast }/R^2$	74	0.3	0.004
[Fe/H]	$0.9R_e$	$h4$	73	$-0.3$	0.003
[Fe/H]	$0.9R_e$	$\Delta {\Sigma }_{23}$	81	0.3	0.02

Table 3. Table 3: Environmental Correlations

Stellar Pop	Radius	Env.	$N_g$	$\rho$	$P$
(1)	(2)	(3)	(4)	(5)	(6)
[$\alpha$/Fe]	$0.9Re$	${\delta }_g$	74	0.3	0.014
[$\alpha$/Fe]	9kpc	${\delta }_g$	73	0.3	0.006
[$\alpha$/Fe]	$1.5Re$	${\delta }_g$	58	0.6	$5\times {10}^{-6}$
[$\alpha$/Fe]	$0.9Re$	$\nu$	74	0.3	0.015
[$\alpha$/Fe]	9kpc	$\nu$	73	0.3	0.012
[$\alpha$/Fe]	$1.5Re$	$\nu$	58	0.4	0.0005
d[$\alpha$/Fe]/dr	$\cdots$	${\delta }_g$	77	$0.3$	0.02
[Fe/H]	$0.9Re$	${\delta }_g$	73	$-0.3$	0.007
[Fe/H]	15kpc	${\delta }_g$	45	$-0.3$	0.036
[Fe/H]	$0.9Re$	$\nu$	73	$-0.3$	0.004

Table 4. Table 4: Mass-Environment Fits

Stellar Pop	Radius	Env.	$\alpha$	$\beta$
(1)	(2)	(3)	(4)	(5)
[$\alpha$/Fe]	kpc	all	$0.265\pm 0.004$	$0.23\pm 0.03$
[$\alpha$/Fe]	kpc	low	$0.255\pm 0.005$	$0.44\pm 0.05$
[$\alpha$/Fe]	kpc	high	$0.297\pm 0.006$	$0.08\pm 0.03$
[$\alpha$/Fe]	$0.9R_e$	all	$0.247\pm 0.005$	$0.26\pm 0.03$
[$\alpha$/Fe]	$0.9R_e$	low	$0.230\pm 0.006$	$0.53\pm 0.07$
[$\alpha$/Fe]	$0.9R_e$	high	$0.300\pm 0.008$	$0.05\pm 0.04$
[Fe/H]	kpc	all	$-0.0618\pm 0.005$	$0.073\pm 0.03$
[Fe/H]	kpc	low	$-0.0417\pm 0.007$	$-0.15\pm 0.09$
[Fe/H]	kpc	high	$-0.111\pm 0.008$	$0.26\pm 0.04$
[Fe/H]	$0.9R_e$	all	$-0.197\pm 0.006$	$-0.31\pm 0.04$
[Fe/H]	$0.9R_e$	low	$-0.173\pm 0.008$	$-0.49\pm 0.10$
[Fe/H]	$0.9R_e$	high	$-0.25\pm 0.01$	$-0.08\pm 0.05$

Equations

$$f(v)\propto \frac{\exp \left[\frac{(\mathrm{v}-\mathrm{V}{)}^2}{{\sigma }^2}\right]}{\sqrt{2\pi {\sigma }^2}}\left[1+{\Sigma }_{m=3}^n{h}_m{H}_m\left(\frac{v-V}{\sigma }\right)\right]$$

$$H_m(x)=\frac{1}{\sqrt{m!}}\exp [{\mathrm{x}}^2]{\left(-\frac{1}{\sqrt{2}}\frac{\partial }{\partial \mathrm{x}}\right)}^{\mathrm{m}}\exp [-{\mathrm{x}}^2]$$

Full text

galaxies: elliptical and lenticular, cD – galaxies: evolution – galaxies: formation – galaxies: kinematics and dynamics – galaxies: stellar content

1 Introduction

The late-time assembly history of massive early-type galaxies remains a topic of ongoing interest. At early times, typical quiescent galaxies were quite compact (e.g., van Dokkum et al., 2008; van der Wel et al., 2008). At the present time, massive early-type galaxies typically have extended envelopes (e.g., Schombert, 1986; Kormendy et al., 2009; Huang et al., 2013, 2018b). To some degree, the larger sizes reflect that at fixed mass, larger galaxies join the red sequence at later times (e.g., Newman et al., 2012), but most massive galaxies also likely build up their outskirts through the accretion of smaller satellites that dissolve at large radius (e.g., Naab et al., 2009; Bezanson et al., 2009). Cosmological hydrodynamical simulations from numerous groups have shown that there are two phases to the build-up of stellar mass in massive early-type galaxies, with a gas-rich phase forming a compact core at early times ($z\approx 2$) followed by dissipationless merging at late times (e.g., Oser et al., 2010; Rodriguez-Gomez et al., 2016; Wellons et al., 2016). From an observational perspective, it is still unclear what fraction of the size growth can be explained by minor mergers as opposed to the quenching of larger galaxies (e.g., Valentinuzzi et al., 2010; Newman et al., 2012; Barro et al., 2013). Observations of the stellar populations and kinematics of local galaxies at large radius may provide complementary insights into the assembly history of massive galaxies.

In the simulations, the fraction of accreted (“ex-situ”) stars rises with both halo mass and stellar mass (e.g., Oser et al., 2010; Rodriguez-Gomez et al., 2016). Photometric observations of local massive early-type galaxies have presented some confirmation for a two-phase formation scenario (Huang et al., 2013; D’Souza et al., 2014; Oh et al., 2017), while Huang et al. (2018b, a) present empirical evidence that the ex-situ fraction correlates both with stellar and halo mass. Spectroscopic observations should provide complementary constraints on the assembly history of massive elliptical galaxies. More specifically, the radial gradients in stellar populations should encode the build-up of stellar mass, particularly if we can reach to large radius. It is very challenging to obtain high signal-to-noise (S/N) observations at large galactocentric radius, and the bulk of papers looking at radial gradients in stellar populations have worked within $R_{e}$ (e.g., Spinrad & Taylor, 1971; Mehlert et al., 2003; Annibali et al., 2007; Spolaor et al., 2010; Jimmy et al., 2013; McDermid et al., 2015; Goddard et al., 2017). In recent years, integral-field spectroscopy (IFS) has enabled stellar population measurements at large radius (e.g., Greene et al., 2013; Scott et al., 2013; Greene et al., 2015; McDermid et al., 2015; Goddard et al., 2017; Boardman et al., 2017; Barone et al., 2018; van de Sande et al., 2018).

Here we focus on the most massive galaxies in the present-day universe using the MASSIVE survey (Ma et al., 2014). We have gathered integral-field data for $90$ MASSIVE galaxies and measured their stellar kinematics including higher-order moments (Veale et al., 2017b, a), stellar velocity dispersion profiles (Veale et al., 2018), and kinemetry (Ene et al., 2018). In addition, we have HST/WFC3 imaging for 30 of the galaxies (Goullaud et al., 2018), and uniform CFHT $K-$band imaging for nearly all objects.

2 Galaxy Sample

MASSIVE is a volume-limited survey of 116 galaxies (see details in Veale et al., 2017a) in the Northern hemisphere within $D<108$ Mpc (i.e., to the distance of the Coma cluster), with $K-$band magnitudes $M_{K}<-25.3$ mag (roughly M^{*}$$>10^{11.5} $M_{\odot}$). Details of the target selection are described in Ma et al. (2014). Briefly, our targets are drawn from the Two Micron All-sky Survey (2MASS; Skrutskie et al., 2006) Extended Source Catalog (XSC; Jarrett et al., 2003) combined with distances from the high density contrast (HDC) group catalog from Crook et al. (2007), supplemented by surface brightness fluctuations when available (Blakeslee et al., 2009, 2010; Blakeslee, 2013) and with the flow model of Mould et al. (2000) when needed.

We have completed integral-field spectroscopic (IFS) observations for 90 galaxies, including complete coverage of the 75 galaxies having $M_{K}<-25.5$ mag. We use the Mitchell IFS at McDonald Observatory (Hill et al., 2008), which has a $107\times 107$ arcsec2 field-of-view, one-third of which is filled by 246 fibers of 4″ diameter. Each galaxy is observed at three dither positions with twenty minute exposures, interspersed with ten minute sky observations. The spectra span 3650-5850 Å, with a spectral resolution of $4.5$Å full-width at half-maximum. The data are reduced using the custom software Vaccine (Adams et al., 2011; Murphy et al., 2011). For more details see Greene et al. (2015).

A primary goal of the MASSIVE survey is to measure spatially resolved stellar kinematics for dynamical modeling (e.g., Thomas et al., 2016). To achieve a mean signal-to-noise (S/N) of at least 20 pixel*-1*, we group individual fibers into spatial bins and co-add the spectra from fibers in a given bin into a single spectrum. The binning scheme is described in detail in Veale et al. (2017b), but we summarize the procedure here briefly for completeness. Central fibers with S/N$>20$ are kept as their own individual bins. Outside of the central regions, the fibers are ‘folded’ across the major axis to boost S/N. The fibers are then grouped into annular bins, and each annulus is sub-divided into an even number of angular bins. The radial extent of each bin is chosen to achieve the target S/N of 20, subject to the constraint that the aspect ratio of each bin, [0.5($R_{\rm outer}+R_{\rm inner}\Delta\theta$)/[$R_{\rm outer}-R_{\rm inner}$], be less than or equal to 1.5. The radial extent of the bins increase outwards until it is no longer possible to achieve S/N$>20$, at which point the remaining fibers are binned into two large radial bins.

Stellar kinematics are measured in each bin using pPXF (Cappellari & Emsellem, 2004) and are presented in Veale et al. (2017b), including stellar velocity and velocity dispersion, but also higher-order moments. The kinematics as a function of environment are analyzed in Veale et al. (2017a) and Veale et al. (2018). Studies of the misalignments between kinematic and photometric axes and local velocity features such as kinematically distinct components are presented in Ene et al. (2018), where an “unfolded” binning scheme was used.

2.1 Environment Measures

We use two probes of the density field, as described in detail in Veale et al. (2017a). One useful measure of galaxy environment is the large-scale density field surrounding a galaxy on the scale of several Mpc. For this, we use the density field from Carrick et al. (2015), based on the 2M++ redshift catalog (Lavaux & Hudson, 2011). The 2M++ covers nearly the full sky to a depth of $K=12.5$ mag, and includes 69,160 galaxies from 2MRS, the Sloan Digital Sky Survey Data Release 7 (SDSS-DR7; Abazajian et al. 2009), and the 6dF galaxy redshift survey Data Release 3 (6dFGRS-DR3; Jones et al. 2009). The Carrick et al. galaxy density contrast $\delta_{g}\equiv(\rho_{g}-\overline{\rho})/\rho_{g}$ is the luminosity-weighted density contrast smoothed with a 5.7 Mpc Gaussian kernel.

We also calculate $\nu_{10}$, a local galaxy density, by taking the distance to the 10th nearest neighbour and measuring the luminosity enclosed in this region. As discussed in detail in Appendix A of Veale et al. (2017a), we adopt an absolute magnitude limit of $M_{K}=-23$ mag to identify the tenth nearest neighbor. Given the magnitude limit of the 2MASS parent sample, we start to lose neighbors for MASSIVE galaxies beyond 80 Mpc. However, Veale et al. estimate that the $\nu_{10}$ values are impacted at a low level by this incompleteness.

2.2 Photometry

Our size and surface-brightness measurements come from Canada-France-Hawai’i Telescope (CFHT) $K-$band imaging (Quenneville et al. in prep). Elliptical isophotes are fitted to each galaxy using ARCHANGEL (Schombert, 2007). Then, a curve of growth is fit to the cumulative aperture luminosities as a function of radius to yield the magnitude of the galaxy. The half-light or effective radius is determined as the radius enclosing 50% of the light determined from the curve of growth analysis. These radii are measured along the major axis and have not been circularized. We find that the radii measured in this way are systematically larger than the 2MASS radii by 20% on average (see details in Quenneville et al. in prep). The adopted effective radii for the 84 galaxies with both CFHT photometry and Mitchell IFS data are included in Table 1.

We measure the outer slope of the surface brightness profile ($\Delta\Sigma/\Delta\rm{log}R$), following Pillepich et al. (2014) and Cook et al. (2016), who show from Illustris simulations (Vogelsberger et al., 2013) that this slope tracks the ex-situ fraction (see also Huang et al., 2018b). Using the ARCHANGEL radial profiles, we fit surface brightness as a function of the log of the radius in $R_{e}$ units, following Cook et al. We do not have uniform coverage to $4R_{e}$ for all targets (Figure 1), so we explore how robust the outer slopes are as a function of radial coverage. We have 14 targets with coverage beyond $4R_{e}$. For this sample, we measure the slope from $2R_{e}$ out to $3R_{e}$, and $4R_{e}$. We find that the slopes measured within the two radial range are well correlated, although those measured within the more restricted radii are slightly smaller on average by $0.2$. We thus adopt the slope measured between $2-3R_{e}$ ($\Delta\Sigma_{23}$ hereafter), as it can be uniformly measured for 80 of the galaxies in our sample.

3 Stellar Population Measurements

In this work, Lick indices are measured to trace the stellar populations. Lick indices were originally developed as a way to extract stellar population information from spectra without flux calibration (Burstein, 1985; Faber et al., 1985; Worthey et al., 1992; Trager et al., 1998). Each Lick index is a narrow region (typically $\sim 20$Å wide) that is mostly dominated by a single element. Isolating these regions allows to study age (through the Balmer lines and particular H$\beta$), metallicity ([Fe/H] through Fe lines), and variable abundance ratios. Of course in practice no window is impacted purely by a single element; particularly at the velocity dispersion of our galaxies (200-400 km s*-1*) all indices are blends of multiple elements, as summarized nicely in Table 1 in Graves & Schiavon (2008).

We adopt the code lick_ew (Graves & Schiavon, 2008) measures the Lick indices and these are fed to EZ_Ages (Graves & Schiavon, 2008) to convert the Lick indices to physical parameters (age, [Fe/H], [$\alpha$/Fe]). The code uses pairs of indices, starting with H$\beta$ and $\langle{\rm Fe}\rangle$, to iteratively solve for the age, abundance and abundance ratios. The models use response functions from Korn et al. (2005) and synthesis models from Schiavon (2007).

We follow Schiavon (2007) and quote [Fe/H], which is directly inferred from the Fe indices, rather than quote a total metallicity. Total metallicity depends on oxygen (the most abundant heavy element) and we do not measure oxygen directly. Instead, we generally assume that [O/Fe] tracks [Mg/Fe]. In our standard runs, we utilize the $\alpha$-enhanced isochrone from Salasnich et al. (2000) and the default assumption that [O/Fe]$=0.5$ to match the $\alpha$-enhanced isochrone value. Since the development of EZ_Ages, more sophisticated modeling schemes have been developed that implicitly solve for [O/Fe] using Lick indices (Thomas et al., 2011; Johansson et al., 2012; Worthey et al., 2014). As well, full spectral modeling takes advantage of information in all of the pixels and boosts the signal-to-noise of the final determinations. This is our goal for the future. Since Conroy et al. (2014) show that Lick methods and full spectral fitting recover the same basic ages, metallicities, and abundance ratios from the same set of SDSS spectra, it is possible to intercompare our results with both literatures.

3.1 Radial Coverage

With the advent of integral-field data, it has become increasingly clear that the aperture used for stellar population measurements impacts the final result; for instance Barone et al. (2018) uncover trends with galaxy density and stellar populations within the effective radius that do not hold for “central” values. Thus, we wish to exploit our IFS data by measuring all properties within multiple standardized radii, both in physical and $R_{e}$ units. In this subsection, we first describe the final sample that we adopt, after removing a few galaxies with more limited radial coverage or S/N. We then discuss the primary measurements that we use at different radii. Ultimately, we will adopt luminosity-weighted measurements measured at $0.9R_{e}$ and 1.5 $R_{e}$, as well as fixed physical radii of 9 kpc and 15 kpc, as motivated by the radial coverage of our data.

We start with 90 MASSIVE galaxies with Mitchell IFS data. We remove five galaxies (NGC 910, NGC 1226, NGC 7052, UGC 3021, and UGC 10918) from consideration for all large-radius tests, due to their very limited radial coverage (in these cases due to poor observing conditions). We exclude them from analysis of both radial gradients and large-radius measurements. These five galaxies do not have any other properties in common (e.g., they are neither the most or least massive, or largest or smallest galaxies). We then examine the radial coverage of our IFS observations for the remaining galaxies. In Figure 1, we show the maximum radial extent of our binned data in kpc and in $R_{e}-$scaled units. In the more massive galaxies, we can reach larger physical radii, but the galaxies are also much larger. Thus, in $R_{e}$ units we reach to larger distances in the lower-mass galaxies.

From inspecting Figure 1, we choose to take a measurement at $0.8-1R_{e}$ that includes most (74) of the galaxies. We also examine trends measured at $1.5R_{e}$ for 58 galaxies, but we systematically exclude the most massive galaxies in this sample due to the very large size of these galaxies and the limited field-of-view of the instrument. The trends with galaxy mass and size are somewhat different when we consider fixed physical apertures. We determine that we can reach $\sim 10$ kpc for 73 galaxies, with no real trend in mass or $\sigma$ for the galaxies that are excluded. These do tend to be the smallest galaxies. We reach smaller physical radius in physically smaller galaxies because we require S/N$>20$ to make the stellar population measurements, and the smallest galaxies have very low surface brightness at 10 kpc. We also take a measurement at 15 kpc for 45 galaxies, although in this case we systematically exclude the lower-mass galaxies, which are both smaller on average and fainter, thus making measurements at large physical radius most challenging.

To determine the measurements at each of these radii, we simply take the measured radial profiles in [Fe/H], and [$\alpha$/Fe], and interpolate to the radius of choice by weighting each measurement by radial distance (see Appendix A). These measurements at fixed radii are the primary way that we investigate the spatial variation in stellar populations, in conjunction with the radial gradients that we describe in the next section.

3.2 Gradients in Stellar Populations

We measure the radial gradients as linear fits to each stellar population parameter as a function of the logarithm (base 10) of the radius. The fits are linear and include the error bars in the parameters. We show the radial profiles of all of the galaxies in Appendix A, with the gradient fits superposed. As emphasized above, the aperture used to measure the gradients can be very important, so we experiment with measuring the gradients over different radial ranges. We try adopting different inner and outer radial coverage by measuring gradients between 1 and 10 kpc, 2 and 10 kpc, and then just truncating the inner radii at 1 or 2 kpc but allowing the outer coverage to extend as far as possible. We find that there is a large variation in measured slopes when the inner radius is changed. For [$\alpha$/Fe], taking the same outer radius but changing the inner radius leads to a scatter in slopes of roughly 0.1. For the [Fe/H] measurements, the scatter is larger (0.15 in slope). However, if the inner radius is kept constant while the outer radius is varied, the scatter is reduced by a factor of two. Thus, we choose to maximize the radial coverage available for individual objects, and take a standard inner radius of 1 kpc, but allow the outer radius to float depending on the object. When we investigate trends with gradients, we exclude galaxies that have fewer than 10 radial points to define the gradient.

The [Fe/H] gradients are generally negative, with a median value of $-0.26$ dex per decade in radius (Figure 2). Only one galaxy (NGC 6223) has a gradient that is consistent (at the $2~{}\sigma$ level) with having a positive gradient ($\Delta$[Fe/H]/$\Delta$ log R$=0.26\pm 0.12$). Unlike the [Fe/H] gradients, the [$\alpha$/Fe] gradients are nearly all consistent with being flat (Figure 2), with a median value ($-0.03$) that is only very slightly negative. There are four galaxies that have [$\alpha$/Fe] gradients that are $>3~{}\sigma$ away from being zero. One (NGC 2256) has a positive gradient ($\Delta$[$\alpha$/Fe]/$\Delta$ log R$=0.23\pm 0.07$). Three have negative slopes (NGC 1272, UGC 02783, NGC 1453) with $\Delta$[$\alpha$/Fe]/$\Delta$ log R$=-0.19\pm 0.06,-0.23\pm 0.08,-0.19\pm 0.05$. These outliers are likely to be real; if we assume that all [$\alpha$/Fe] gradients have the same intrinsic value of $-0.03$, and then perturb each measurement by its uncertainty, we expect to detect three $>3\sigma$ outliers in the full sample only 1% of the time, and we find four outliers. On the other hand, the one positive [Fe/H] gradient could be marginally consistent with pure scatter. If we assume that all the [Fe/H] gradients are intrinsically equal to the mean $-0.3$ dex per decade, then given our error bars we expect to find a positive gradient 2% of the time.

4 Scaling Relations Between Stellar Populations and Galaxy Properties

Early-type galaxies are known to exhibit scaling relationships between their structural properties and their stellar populations (e.g., Trager et al., 2000a; Thomas et al., 2005; Graves et al., 2009). Only recently, with the advent of large IFS surveys, has it become possible to examine these scaling relationships outside of the central regions of the galaxies (e.g., Greene et al., 2015; van de Sande et al., 2018; Barone et al., 2018), and since there are gradients in metallicity with radius, different trends have different dependence on the aperture used.

In this section, we focus on possible trends between structural properties of our galaxies and their stellar populations in the center and at large radius. Specifically, we examine stellar mass as inferred from the $K-$band magnitude, using a dynamically derived mass-to-light ratio (see details in Ma et al., 2014; Cappellari, 2013). We also look at the stellar velocity dispersion, and since we are interested in our galaxies on large scales, we use the luminosity weighted average $\sigma$ over all Mitchell bins within $R_{e}$ (Table 1, col 12 in Veale et al., 2017b). We also consider proxies for the gravitational potential $\Phi\propto M^{*}/R_{e}$ and the stellar surface density $\Sigma\propto M^{*}/R_{e}^{2}$ following Barone et al. (2018), who study a galaxy sample of (typically) lower $M^{*}$ with SAMI. Motivated by the Illustris works of Pillepich et al. (2014) and Cook et al. (2016), we investigate correlations between the stellar population parameters and the slope of the surface brightness profile $\Delta\Sigma_{23}\equiv\Delta\Sigma/\Delta\rm{log}R(2-3R_{e})$ (§2.2; Figure 3).

4.1 Structural Correlations with the Galaxy Centers

By virtue of our IFS data, we have stellar population measurements out to $>10$ kpc ($\sim R_{e}$ or beyond; Fig. 1) for the majority of galaxies in our sample. Our primary interest is to investigate the stellar content in the outer parts of the galaxies, but first we show briefly that our central measurements are consistent with prior work (for more detailed comparisons, see Greene et al., 2012, 2013). We construct “central” stellar population measurements using all of the fibers contained within 2 kpc of the galaxy center. These are luminosity-weighted measurements, to mimic the process for individual SDSS fibers. We then ask whether our central measurements correlate with $\sigma$ or $M^{*}$. To test for correlation, we use the non-parametric Spearman correlation coefficient. The results of our correlation tests are presented in Table 2 and Figure 2, where we only include “significant” correlations with a probability $P<0.05$ of the null hypothesis.

In keeping with prior results, we find very strong correlations between [$\alpha$/Fe] and $\sigma$ ($\rho=0.5$, $P=1\times 10^{5}$ of a null result), significant correlations between [$\alpha$/Fe] and $M^{*}$ ($\rho=0.26$, $P=0.017$), and weak to no correlation between [Fe/H] and $M^{*}$ ($\rho=0.22$, $P=0.04$) or $\sigma$ ($\rho=0.13$, $P=0.23$; e.g., Trager et al., 2000b; Graves et al., 2009; Wake et al., 2012; Conroy et al., 2014).

We also find an interesting hint that [$\alpha$/Fe] saturates above M^{*}$$=10^{11.8} $M_{\odot}$. At the lowest masses covered in our sample, [$\alpha$/Fe] shows an increasing trend, but then seems to flatten out at the highest masses. Quantitatively, galaxies in a stellar mass bin of M^{*}$$=10^{11.8}-10^{12} $M_{\odot}$, and all galaxies more massive than this limit, both have a consistent weighted mean [$\alpha$/Fe]$=0.32\pm 0.01$. This consistency at high mass holds independently of exactly how we divide the two mass bins. This convergence may be a sign that major mergers are needed to make these most massive galaxies. Simulations also show that at the highest stellar masses, the predominance of more major mergers leads to a convergence of properties for high-mass central and satellite galaxies (e.g., Wetzel et al., 2013).

By virtue of our IFU data and high signal-to-noise ratio (SNR) spectra, we are able to look for correlations with other kinematic tracers beyond $\sigma$. We see intriguing correlations with the higher-order moments of the line-of-sight velocity distribution (LOSVD). Using a Gauss-Hermite decomposition, we model the LOSVD with velocity, velocity dispersion, and four higer-order moments ($h3-h6$) as described in detail in Veale et al. (2017b). Specifically, for spectrum $f(v)$, mean velocity $V$, velocity dispersion $\sigma$, and the number of higher orders $n=6$, we have:

[TABLE]

And the $H_{m}(x)$ are the Hermite polynomials:

[TABLE]

The third term $h3$ is the skewness, and the fourth term $h4$ is the kurtosis. To maximize SNR, we use a light-weighted average $h4$ measurement (see details in Veale et al., 2017b). In the galaxy center, we see a correlation between $h4$ and [$\alpha$/Fe] that is as strong ($\rho=0.5$, $P=2\times 10^{5}$) as the correlation between $\sigma$ and [$\alpha$/Fe] (Figure 4). We see no correlation between $h4$ and central values of [Fe/H] ($\rho=-0.0073$, $P=0.95$).

4.2 Structural Correlations with the Outskirts

We then extend the above analysis to the outer regions of the galaxies and a broader range of structural parameters. In all cases, we seek correlations between the stellar populations measured at fixed $R_{e}$-scaled radii ($0.9R_{e}$ and $1.5R_{e}$) and fixed physical radii (9 kpc and 15 kpc), and gradients measured beyond 1 kpc (§3.2). All significant correlations are included in Table 2. The correlations we saw between central [$\alpha$/Fe] and both $\sigma$ and $M^{*}$ persist when we investigate the outer parts of the galaxies ($\sigma$-[$\alpha$/Fe] at $0.9R_{e}$: $\rho=0.5$, $P=10^{-5}$; $M^{*}$-[$\alpha$/Fe] at $0.9R_{e}$: $\rho=0.4$, $P=0.0003$). This is not surprising, since we do not measure significant gradients in [$\alpha$/Fe] (Figure 2). In contrast, no compelling correlations between $\sigma$ or $M^{*}$ emerge with [Fe/H] or its gradients at large radius ($\sigma$-[Fe/H] at $0.9R_{e}$: $\rho=-0.12$, $P=0.3$; $M^{*}$-[Fe/H] at $0.9R_{e}$: $\rho=-0.13$, $P=0.26$). We will discuss the correlations with higher-order moments in the following section.

Like Barone et al. (2018), we find a significant correlation between $M^{*}/R_{e}$ ($\rho=0.3,P=0.006$; a proxy for the gravitational potential) and $M^{*}/R_{e}^{2}$ ($\rho=0.3,P=0.004$; the stellar surface density) with [Fe/H] measured at $0.9R_{e}$ (Figure 3; Table 2), and in common with that work, we find no such correlation for the central stellar populations, only those measured at larger radius. Higher metallicity in higher surface mass density and/or higher gravitational potential strongly suggests that the stars in question were formed in a similar potential to their current one, and that the metal retention rate is set by the gravitational potential.

We next examine the correlation between stellar populations and $\Delta\Sigma_{23}$, the gradient in the surface brightness profile measured between $2-3R_{e}$. We are motivated by results from Illustris (Vogelsberger et al., 2013) showing that the surface brightness slope grows shallower as galaxies accrete more of their stars from external galaxies (Pillepich et al., 2014). Cook et al. (2016) expand on this finding to show that the metallicity gradients also grow shallower as the amount of accretion increases. We do find a very similar range in surface brightness slopes in our data to that presented by Cook et al. (2016) for a similar stellar mass range. However, we do not detect a significant correlation between the gradients in stellar population parameters and the gradients in surface brightness (with $P>0.3$ in both cases). Nor do we find a correlation between $\Delta\Sigma_{23}$ and [$\alpha$/Fe] measured at any radius ($P>0.2$ in all cases). $\Delta\Sigma_{23}$ is postively correlated with [Fe/H] measured at $0.9R_{e}$ ($\rho=0.28,P=0.02$) and [Fe/H] at $1.5R_{e}$ ($\rho=0.26,P=0.05$).

Finally, we checked for correlations between the stellar population gradients and the ratio of rotational to dispersion support ($\lambda$; Emsellem et al., 2007; Veale et al., 2017b) as well as looking for any interesting stellar population properties for galaxies with counter-rotating componenets (Ene et al., 2018), but we did not find any significant trends.

4.3 The correlation between $h4$ and outer stellar populations

As described above, we are able to robustly measure not only $V$ and $\sigma$, but also deviations from Gaussian LOSVDs from our high SNR data. We are particularly interested in the kurtosis, $h4$. High kurtosis can be a sign of radial anisotropy, which in turn may point to stellar accretion events (e.g., Wu et al., 2014; Amorisco, 2017). We investigate possible correlations between $h4$ and the stellar populations in the outer parts of the galaxy (Figure 4). To maximize the SNR, we use the light-weighted average $h4$ measurement from Veale et al. (2018).

We find that the global $h4$ measurement correlates with the $r=0.9R_{e}$ measurements of both [Fe/H] ($\rho=-0.33,P=0.003$) and [$\alpha$/Fe] ($\rho=0.48$, $P=2\times 10^{-5}$). The correlation between [$\alpha$/Fe] and $h4$ is comparably strong to that between [$\alpha$/Fe] and $\sigma$. In the case of [Fe/H], we do not see a correlation with $\sigma$ at all, but detect a very significant correlation with $h4$ at $R\approx 0.9R_{e}$. This correlation persists even if we remove the four points with the lowest [Fe/H]$<-0.4$. Thus, we do not believe that the correlation with $h4$ is somehow derivative of the correlation with $\sigma$. Likewise, there is a known correlation between $h4$ and stellar mass in this sample (Veale et al., 2017b), but [Fe/H] at $0.9R_{e}$ and $M^{*}$ are not correlated either ($\rho=-0.13,P=0.26$). To double check that the correlation between stellar populations and $h4$ is not derivative of correlations with stellar mass, we also perform a linear fit to [$\alpha$/Fe] and [Fe/H] versus stellar mass, dividing into two $h4$ bins. We find significant differences in zero-point for both [$\alpha$/Fe] ($0.22\pm 0.006$ for $h4<0.029$; $0.29\pm 0.007$ for $h4>0.029$) and [Fe/H] ($-0.15\pm 0.007$ for low $h4$; $-0.28\pm 0.009$ for high $h4$) at $0.9R_{e}$. This difference confirms our claim that even at fixed mass, there are real differences in the stellar populations as a function of $h4$.

One concern is that template mismatch may cause an apparent correlation between $h4$ and stellar population properties, where $h4$ rises to compensate, for instance, if the template library does not have stars with high [$\alpha$/Fe] and low [Fe/H]. To investigate this possible bias, we look for correlation between the local values of $h4$, [Fe/H], and [$\alpha$/Fe] measured in each spatial bin of individual galaxies. Restricting our attention to systems with more than 15 radial bins, we find that only $\sim 6\%$ show correlation between local $h4$ and [$\alpha$/Fe] (with the majority showing an anticorrelation, rather than the positive correlation seen for the ensemble) and $\sim 14\%$ show correlation between [Fe/H] and $h4$, split between positive and negative correlation. We thus conclude that the correlations we see between $h4$ and the outer stellar population properties are not caused by template mismatch.

4.4 Correlations with Environment

We now look for trends between stellar populations and large-scale ($\delta_{g}$) and local ($\nu$) environmental measures (these are reviewed in §2.1). First we naively investigate correlations between [Fe/H] and [$\alpha$/Fe] (or gradients therein) with $\delta_{g}$ or $\nu$ without controlling for stellar mass. Table 3 contains all of the correlations with $P<0.05$, again using the Spearman rank correlation test. No central stellar population parameters show any correlations with environmental measures ($P>0.1$ in all cases), but we see significant correlations between both environmental measures and [Fe/H] and [$\alpha$/Fe] at larger radius (Figure 5). In general, correlations with $\delta_{g}$ (the large-scale environment measure) appear more significant than those with $\nu$, the local indicator. The sense of the trend is that [Fe/H] is lower and [$\alpha$/Fe] is higher in overdense regions. This type of trend is as we would expect if galaxies in the overdense regions formed earlier. We see a weak correlation between gradients in [$\alpha$/Fe] and $\delta_{g}$, which we were not able to detect with our earlier stacking analysis in Greene et al. (2015).

Given the strong covariance between stellar mass and environment (see review in Blanton & Moustakas, 2009), and the observed correlation between $M^{*}$ and stellar population parameters (Figure 2), we must investigate whether the correlations between stellar population parameters and environment are simply a biproduct of a dominant correlation with $M^{*}$ (e.g., Veale et al., 2017b, a). We therefore fit a relation between $M^{*}$ and [Fe/H] or [$\alpha$/Fe], in two bins of $\delta_{g}$. We divide the sample at the median value of log ($1+\delta_{g}$)$=0.6$, so that there are even numbers of galaxies in each bin. The results are shown in Figure 6 and the fits are presented in Table 4.

We find significant slope differences in the $M^{*}$-[Fe/H] and $M^{*}$-[$\alpha$/Fe] fits to the two density bins. Specifically, for both [Fe/H] and [$\alpha$/Fe], the high-density systems show a much flatter relation between $M^{*}$ and stellar population parameters than in the systems in the lower-density environments. The slope differences are significantly different in all cases, although they are larger in the $0.9R_{e}$ bin. As shown by the figure, at high mass (where we have sparse data) the two density bins are convergent, while the different slopes arise mainly because the lower-density points have higher [Fe/H] and commensurately lower [$\alpha$/Fe] than galaxies of similar mass in higher-density environments. Similar trends have been reported by a number of authors (Thomas et al., 2005; Liu et al., 2016; Gu et al., 2018, and see §5). Specifically, the ATLAS3D stellar population data from McDermid et al. (2015) show a similar split in [$\alpha$/Fe] and [Fe/H] by local overdensity in their highest mass bin.

We are looking at these trends binned in $\delta_{g}$, which looks at overdensities on $\sim 5$ Mpc scales. Much of the literature has looked at $\nu$, or local overdensity. Discussing the relative merits of different environment indicators is beyond the scope of this work, but we note here that $\nu$ in particular can be a complex environmental indicator, since at low densities it is sensitive to large-scale environment (the so-called “two-halo” term), while in dense regions it probes small scales (or the “one-halo” term; see Woo et al., 2013). Accordingly, as shown by Veale et al. (2017a), $\nu$ and $\delta_{g}$ are well correlated at low density (where $\nu$ probes $>5$ Mpc scales anyway) but diverge in high density regions, where $\nu$ can grow much more rapidly. In our case, given our limited numbers, the median $\delta_{g}$ value of log$(1+\delta_{g})\sim 1.5$ is close to the value where the two indicators diverge, and thus we see very similar trends when we play the same game with $\nu$.

5 Discussion

5.1 Ex-situ fractions, merger mass ratios, and expected stellar populations

Motivated by observations of very compact sizes for quiescent galaxies at high redshift (e.g., van Dokkum et al., 2008; van der Wel et al., 2014), theorists have examined the growth histories of massive galaxies and determined that there is often a two-phase growth. An early and rapid dissipational phase builds a compact central nugget, and is followed by a build-up of the outer parts of galaxies via minor merging (e.g., Oser et al., 2010; Lackner et al., 2012; Rodriguez-Gomez et al., 2016; Wellons et al., 2016). However, the relative roles of minor merging and progenitor bias in the aggregate growth of the galaxy population remains a matter of debate (e.g., Newman et al., 2012). Observations of light profiles of low-redshift early-type galaxies also provide some support for a two-phase picture (e.g., Huang et al., 2013; D’Souza et al., 2014; Oh et al., 2017), but adding stellar population and kinematic information at large radius may add additional insight.

A number of cosmological and cosmological zoom studies have looked at how mass is accumulated in massive galaxies with time. They find that the fraction of total stellar mass brought in via mergers rises towards more massive galaxies (e.g., Oser et al., 2010; Lackner et al., 2012; Hirschmann et al., 2015). Despite this overall trend, there is considerable scatter from system to system because of different merger histories, as emphasized by Rodriguez-Gomez et al. (2016).

In terms of the radii that are most impacted by merging, the more massive and more tightly bound satellites deposit their stars closer to the center, while the lower-mass and less-bound systems deposit mass at larger radius (see also Boylan-Kolchin & Ma, 2007; Amorisco, 2017). Accordingly, major mergers will tend to flatten the steep gradients that arise from pure in-situ formation (e.g., White, 1980; Kobayashi, 2004; Hirschmann et al., 2015). As the dominant merger type becomes more minor, the accreted systems tend to have their star formtion truncated early, when they fall into the larger halo but before they are accreted into the primary galaxy. Thus, minor mergers tend to create positive age gradients, bring in $\alpha-$enhanced stars, and steepen metallicity gradients (Hirschmann et al., 2015). This expectation is also consistent with emerging observations of the stellar populations of massive galaxies at $1<z<2$, which also show higher [$\alpha$/Fe] and lower [Fe/H] in the population that our MASSIVE galaxies may have later accreted (e.g., Onodera et al., 2015; Lonoce et al., 2015; Kriek et al., 2016).

In the context of these simulations, our goal is to ask whether various structural measurements may link to the merger history, and then whether the stellar populations show any correlations with them. For instance, at the highest $M^{*}$ in our sample, we see that the central [$\alpha$/Fe] values match those in galaxies factors of $2-3$ less massive. The most likely explanation is that major merging dominates the mass growth at the centers of these galaxies. The trend between [Fe/H] and gravitational potential also points to a significant contribution from in-situ stars $<R_{e}$, since those stars preserve a memory of the gravitational potential in which they formed. Larger radius measurements are closer to the spatial regions where more minor mergers would deposit their stars.

5.2 Observable consequences of varying ex-situ fractions and links with

stellar populations

We now consider how we might sort galaxies based on their accretion histories. Stellar population gradients alone are insufficient, since steep gradients can result from predominantly in-situ formation or a large number of minor mergers (e.g., Kobayashi, 2004; Hirschmann et al., 2015). We have explored two structural characteristics here, the average $h4$ value (Thomas et al., 2007; Wu et al., 2014; Amorisco, 2017) and the outer surface brightness slope $\Delta\Sigma_{23}$ (Pillepich et al., 2014; Cook et al., 2016). We will discuss each of these in turn.

We see strong correlations between both [Fe/H] and [$\alpha$/Fe] and $h4$. We have argued that template mismatch is unlikely to drive this correlation. We now ask whether the correlations are linked with the accretion history, or instead some other property (e.g., the gravitational potential) of the galaxy. Simulations have shown that radial anisotropy, and thus $h4$, can rise with merging (Wu et al., 2014). On the other hand, positive $h4$ can result from gradients in circular velocity as well (e.g., Gerhard, 1993; Baes et al., 2005). In the MASSIVE sample, Veale et al. (2017b, 2018) show that there are likely two causes for positive $h4$, particularly in the outer parts of our galaxies. In two thirds of the galaxies, $\sigma$ falls outward. Those galaxies with falling $\sigma$ and positive $h4$ likely have rising radial anisotropy. While major mergers can produce radial orbits at large radius (e.g., Rantala et al., 2018), N-body simulations suggest that the stronger the radial anisotropy is, the more likely it originates from minor mergers and/or accretion (Hilz et al., 2012). In the other third of galaxies (typically the most massive), we see rising $\sigma$ profiles (e.g., Dressler, 1979; Loubser et al., 2008; Veale et al., 2018, and references therein). In these galaxies, $h4$ will rise outward simply because of the gradient in the potential.

We thus revisit the correlation between $h4$ and the outer stellar populations, this time splitting the sample into those with rising and falling $\sigma$ profiles (Fig. 7). We find that the galaxies with rising $\sigma$ profiles show uniformly high [$\alpha$/Fe] and relatively low [Fe/H], while the galaxies with falling $\sigma$ profiles span the full range of [$\alpha$/Fe] and [Fe/H], and these drive the correlation with $h4$. In these galaxies, we posit that the higher $h4$ comes from radial anisotropy and that the outer stellar populations are closely linked to the anisotropy of their orbits. The higher the radial anisotropy, the more $\alpha$-enhanced (and thus older and metal-poor) are the stellar populations. The range in radial anisotropy is linked to the merger history, perhaps pointing to more minor merging in galaxies with more radial anisotropy (Hilz et al., 2012).

Based on this behavior, we propose the following picture. Galaxies with falling $\sigma$ profiles display a sequence in ex-situ fraction, whereby those with low $h4$ are also those with quiet accretion histories, leading to solar-like abundances and abundance ratios in their outskirts. Those with the highest $h4$ have high radial anisotropy due to additional merging, which brings in older, metal-poor and $\alpha-$enhanced stars. Galaxies with rising $\sigma$ profiles have high $h4$, but in these cases radial anisotropy need not be invoked. These also tend to be the galaxies in the most massive halos (Veale et al., 2018), and we are likely seeing a transition into an outer envelope or intra-cluster light component. The stellar populations that we measure in the outer parts of these galaxies (namely very high [$\alpha$/Fe] and low [Fe/H]) are consistent with the stellar populations measured in the envelopes of brightest cluster galaxies in other works (e.g., Coccato et al., 2011; Edwards et al., 2016). We will revisit these galaxy outskirts from the perspective of environment in the next section.

Pillepich et al. (2014) suggest that gradients in surface brightness profile (here we measure $\Delta\Sigma_{23}$) should be a robust way to sort galaxies of similar mass by their ex-situ fractions. Cook et al. (2016) use the same Illustris simulations to predict shallower [Fe/H] gradients in systems with the shallowest surface brightness slopes. We see a correlation between $\Delta\Sigma_{23}$ and the [Fe/H] measured at $0.9R_{e}$, but no correlation with stellar population gradients. Thus, there is a hint that the gradient in surface brightness is linked with the stellar populations, and the correlation runs in a similar way as suggested by Cook et al. However, we have a few limitations currently to making a definitive comparison. Despite our efforts, we still cannot measure the stellar populations at matching radii to the surface brightness gradients. Furthermore, we would need a wider dynamic range in $M^{*}$ to determine whether $\Delta\Sigma_{23}$ correlates with stellar population properties.

5.3 Environment and Galaxy Assembly

We observe a hint of lower [Fe/H] and enhanced [$\alpha$/Fe] at fixed mass in more dense environments, both when environment is measured as local overdensity ($\nu$) and when measured on Mpc scales ($\delta$). We can understand this result in the context of prior work focused both on the fossil record, looking at star formation histories of massive galaxies as a function of environment, and work that looks at the star formation rates in overdensities at different redshifts.

Previous work has seen clear evidence that star formation ended earlier in denser environments (e.g., Thomas et al., 2005; Scott et al., 2017), particularly at lower $M^{*}$ where “rejuvenated” spheroids are only found in low-density environments (e.g., Thomas et al., 2010; Pasquali et al., 2010). In fact, Liu et al. (2016) show that low-mass spheroids have a very large range in [$\alpha$/Fe] both in the densest and least dense environments. Liu et al. (2016) also find that [$\alpha$/Fe] in dwarfs correlates with the galaxy distance from M87, also consistent with our picture. McDermid et al. (2015) additionally show that the star-formation histories of cluster galaxies are not just truncated earlier but actually proceed more rapidly (e.g., the same stellar mass is built up in a shorter burst) in denser environments, leading potentially to lower [Fe/H] and higher [$\alpha$/Fe]. Observations of star formation in dense environments at $z>1.5$ also support the idea that cluster galaxies experience more intense bursts of star formation (e.g., Elbaz et al., 2007; Brodwin et al., 2013; Tran et al., 2010; Noirot et al., 2018).

Of particular interest to our study, Gu et al. (2018) find that all galaxies in the central parts of the Abell 3827 cluster show enhanced [$\alpha$/Fe] for their mass. They argue that all galaxies in cluster centers have experienced rapid and early truncation of their star formation, and that this “coordinated assembly” of the satellite and central galaxy leads directly to the flat [$\alpha$/Fe] gradients that they observe in the central galaxy.

Gu et al. also suggest that the radial gradient in [$\alpha$/Fe] should be flatter in overdense environments, since all the galaxies available for accreting will also contain elevated [$\alpha$/Fe] for their mass. In contrast, a relatively isolated MASSIVE galaxy might ingest a satellite with “normal” [$\alpha$/Fe] for its mass, leading to a declining [$\alpha$/Fe] gradient. Unlike in Greene et al. (2015) where we did not have the statistics, we do in fact see a weak trend between the [$\alpha$/Fe] gradients and the environment in just this sense.

6 Summary

We look at the stellar populations throughout 90 MASSIVE galaxies, with a focus on the outer parts of the galaxies. In moving beyond measurements of stellar populations weighted towards the galaxy centers, we are able to uncover new trends between stellar population parameters and the kinematics of these massive galaxies, providing new insight into their assembly history.

In galaxy centers, we see that the [$\alpha$/Fe] saturates at stellar masses $\sim 10^{11.8}$ $M_{\odot}$, suggesting that major mergers play a large role in building up the centers of the most massive galaxies. At and beyond $\sim R_{e}$, $\sigma$ and [$\alpha$/Fe] are tightly correlated. In contrast, the metallicity [Fe/H] is correlated with the gravitational potential and the surface mass density, pointing to a significant in-situ component to the stars at these radii.

We investigate two structural measurements that may correlate with the ex-situ (or accreted) fraction. We find a correlation between the outer surface brightness slope $\Delta\Sigma_{23}$ and [Fe/H] beyond $\sim R_{e}$. Our finding provides some support for the Illustris results from Pillepich et al. (2014) and Cook et al. (2016) that the outer surface brightness slope is a proxy for accreted or ex-situ fraction. To truly compare with simulations, a wider dynamic range in stellar mass is needed.

We also find a strong correlation between [$\alpha$/Fe] and $h4$, and we even see a correlation between [Fe/H] and $h4$ at $\sim R_{e}$. Galaxies with the most positive $h4$ have super-solar [$\alpha$/Fe]$\approx 0.4-0.6$ dex and low [Fe/H]$\approx-0.5$ dex. Mergers can increase the radial anisotropy and thus boost $h4$, or $h4$ can rise due to a gradient in the potential (e.g., the transition from a galaxy halo to a cluster halo). Even removing the galaxies with rising $\sigma$ profiles, we still find a strong correlation with stellar populations, pointing to a rising accreted population in galaxies with higher radial anisotropy.

The galaxies with rising $\sigma$ profiles at large radius also tend to be those in the densest environments (Veale et al., 2018), and we find evidence that the [$\alpha$/Fe] in the outer parts correlates with environment at fixed mass. Lower-density environments have lower [$\alpha$/Fe], while higher-density environments have very super-solar [$\alpha$/Fe]$\sim 0.4-0.6$. If star formation proceeds earlier and more rapidly in denser environments, then we would see older, metal-poorer, and more $\alpha-$enhanced stars in the outskirts of galaxies in richer environments.

We are interested in probing the assembly history of massive galaxies by looking at their stellar populations at large radius. Although our survey was designed to reach large radii for very massive galaxies, and does as well as any existing data set, there is some real chance that our observations still do not reach out to large enough radius to truly constrain the accreted stars. Current imaging surveys are now able to reach out to $\sim 100$ kpc for individual galaxies with photometry (e.g., Huang et al., 2018b), but we must wait for next-generation telescopes to have the collecting area to measure detailed stellar population properties out to comparable radii.

Acknowledgements

We thank the referee for a timely, supportive, but thoughtful report that improved this manuscript. The MASSIVE survey is supported in part by NSF AST-1411945, AST-1411642, AST-1815417, AST-1817100, HST GO-14210, GO-15265 and AR-14573. We are grateful for useful discussions with M. Kriek.