The Kinematics of Massive Quiescent Galaxies at 1.4 < z < 2.1: Dark Matter Fractions, IMF Variation, and the Relation to Local Early-type Galaxies^*

Observations of VIRIAL galaxies were carried out between 2014 and 2016 using the KMOS YJ band (1–1.36 μm). Data were taken using a standard object-sky-object pattern with individual exposure times of 300 s. Each science exposure was offset by between 01 and 06 in order to avoid bad pixels in the final extracted spectra. Along with our science targets, we assigned one IFU from each of the three KMOS spectrographs to a reference star that we used to monitor the ambient conditions (seeing, atmospheric transmission, etc.), pointing accuracy, and point-spread function (PSF) shape. Due to the relatively small angular size of the KMOS IFUs (28 × 28), sky exposures were taken nodding completely off-source. Total on-source integration times range from 440 to 740 minutes (see Table 1).

Data were reduced using a combination of the Software Package for Astronomical Reductions with KMOS pipeline tools (SPARK; Davies et al. 2013) and custom Python scripts. In the following we briefly outline the steps used to produce calibrated one-dimensional spectra. Details of the VIRIAL reduction will be described in a future paper (Mendel et al., in preparation). Calibration exposures (dark, arc, and flat) were reduced using standard SPARK routines to produce flat-field, wavelength, and spatial calibration frames. When processing science frames, we first corrected each raw image for a readout-channel-dependent bias term estimated from reference pixels around the perimeter of each detector. We then adjusted the wavelength and spatial illumination calibrations for each exposure based on the positions and relative flux of bright sky lines before subtracting the object and sky images. The brightness of atmospheric OH lines can vary significantly between object and sky exposures (∼10% on 5- to 10-minute timescales; Ramsay et al. 1992; Davies 2007), often leading to significant systematic residuals in the initial sky-subtracted frames. In order to limit the impact of these systematics on our final spectra, we performed a second-order correction to the sky for each IFU using residuals measured in other IFUs in the same detector, excluding the IFU of interest.

One-dimensional spectra were extracted directly from the flat-fielded, illumination-corrected, and sky-subtracted detector frames for each exposure separately. Since VIRIAL targets are typically undetected in individual 300 s exposures, we used the available 3D-HST/CANDELS F125W imaging (Grogin et al. 2011; Koekemoer et al. 2011; Skelton et al. 2014) to model the source flux distribution and mask neighboring objects in the optimal extraction. The HST images were convolved to match the KMOS PSF measured from the reference stars in each exposure, which were also used to adjust for changes in transmission between exposures. Individual optimally extracted spectra were then corrected for telluric absorption using synthetic atmospheric models computed with MOLECFIT (Kausch et al. 2014) and combined using inverse variance weights. Uncertainties on the output spectra were estimated using bootstrap combines of the individual 1D spectra for each object. The typical spectral resolution in the extracted 1D spectra (as measured from sky lines) ranges from R = 3000 to 3500 (σ_inst ≈ 36–42 km s⁻¹) depending on arm and detector (see also Wisnioski et al. 2019).

We estimated stellar velocity dispersions for VIRIAL galaxies using a simultaneous fit to the observed KMOS spectrum and multiband photometry from 3D-HST (Skelton et al. 2014). We generated model spectral energy distributions (SEDs) using FSPS v2.4 (Conroy et al. 2009; Conroy & Gunn 2010) assuming a lognormal SFH with

$$\begin{eqnarray}\mathrm{SFR}(t)=\left\{\begin{array}{ll}\displaystyle \frac{1}{t\sqrt{2\pi {\tau }^{2}}}{e}^{-\tfrac{{\left(\mathrm{ln}t-\mathrm{ln}{t}_{0}\right)}^{2}}{2{\tau }^{2}}} & {\rm{if}}\ t\leqslant {t}_{\mathrm{trunc}}\\ 0 & {\rm{if}}\ t\gt {t}_{\mathrm{trunc}},\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where t is the age of the universe, t₀ is the delay time, and τ controls the width of the distribution (see also Gladders et al. 2013). The additional parameter t_trunc allows for star formation to be abruptly truncated and provides added flexibility when modeling the SFHs of quiescent galaxies. We stress that our adoption of a lognormal SFH is motivated by its flexibility compared to more commonly used τ or delayed-τ models, rather than an assumption that galaxies' SFHs are intrinsically lognormal (e.g., Gladders et al. 2013; Abramson et al. 2016; Diemer et al. 2017). In Appendix A we show that our derived velocity dispersions are not biased by the use of a parametric SFH. We modeled the effects of dust using a two-component extinction law that includes a foreground screen and additional attenuation toward young stellar populations (<10⁷ yr), which are assumed to remain embedded within their birth clouds (see, e.g., Charlot & Fall 2000). We used the reddening curve of Calzetti et al. (2000) and, following Wuyts et al. (2013), adopted a relationship between the total V-band extinction A_V and the additional extinction toward young stellar population A_extra such that ${A}_{\mathrm{extra}}=0.9{A}_{V}-0.15{A}_{V}^{2}$. For simplicity we assume a fixed solar metallicity; in Appendix A we show that changing the metallicity by ±0.2 dex leads to systematic shifts in the derived velocity dispersion of ≲2%.

Before fitting, templates were smoothed to match the wavelength-dependent KMOS resolution measured from sky lines in extracted 1D spectra. The final matched templates include an additional (constant) offset of σ_offset = 65 km s⁻¹ to account for the resolution difference between KMOS (σ_inst ≈ 40 km s⁻¹) and the adopted MILES spectral library (σ_MILES ≈ 75 km s⁻¹; Beifiori et al. 2011). This effectively sets a floor for our velocity dispersion measurements of 65 km s⁻¹. We limited our fits to the wavelength range from 3750 to 5300 Å and included a ninth-order additive polynomial—corresponding to ∼1 order per 10,000 km s⁻¹—to minimize the effects of template mismatch on our final velocity dispersion measurements. We verified that our results are not sensitive to the adoption of an additive, as opposed to multiplicative, polynomial. In the end our model has a total of seven free parameters: redshift, z; stellar mass, M_SPS;¹² three parameters that describe the SFH, τ, t₀, and t_trunc; absolute V-band extinction, A_V; and stellar velocity dispersion, σ_*. Samples from the posterior distribution were generated using emcee (Foreman-Mackey et al. 2013), and our final estimates of velocity dispersion and stellar mass were taken as the medians of their respective marginal posterior distributions, with 1σ uncertainties estimated from the 16th and 84th percentiles. We have confirmed that the derived stellar masses do not change significantly if we refit objects using only the available photometric data (i.e., excluding spectra). The final redshifts, stellar masses, and velocity dispersions are provided in Table 2.¹³ One-dimensional spectra and the corresponding best-fit models are shown in Figure 2.

Zoom In Zoom Out Reset image size

Table 2.  Redshift, Velocity Dispersion, and Stellar Masses of KMOS Galaxies

Field	ID	z	$\mathrm{log}({M}_{\mathrm{SPS}}/{M}_{\odot })$	σea
				(km s−1)
UDS	22480	1.5288	11.08	323 ± 42
UDSb	24891	1.6031	10.99	146 ± 32
GOODS-S	39364	1.6118	11.10	203 ± 42
GOODS-S	42113	1.6140	11.20	362 ± 65
GOODS-S	43548	1.6143	10.64	169 ± 43
COSMOS	6977	1.6424	10.86	187 ± 32
UDSb	22802	1.6660	11.13	337 ± 28
UDSb	29352	1.6886	10.91	277 ± 46
UDS	10237	1.7664	11.38	233 ± 23
COSMOS	7411	1.7808	11.09	186 ± 28
UDS	35111	1.8217	10.95	228 ± 36
UDS	32892	1.8243	11.02	206 ± 27
UDS	38073	1.8245	10.94	194 ± 49
COSMOS	6396	1.8364	10.90	169 ± 33
COSMOS	9227	1.8604	10.98	273 ± 41
COSMOS	7391	1.8681	10.54	145 ± 38
COSMOS	2816	1.9230	11.26	297 ± 49

Notes. The formal statistical uncertainties on stellar masses derived from our SED fitting are of order 0.02 dex. Where relevant we include an additional 0.15 dex uncertainty on $\mathrm{log}({M}_{* }/{M}_{\odot })$ in quadrature to account for systematic uncertainties in the determination of stellar masses (see, e.g., Conroy et al. 2009; Mendel et al. 2014).

^aVelocity dispersion corrected to r_e following the procedure described by van de Sande et al. (2013). ^bThese galaxies are in common with the Belli et al. (2017) sample. See discussion in Section 2.2.

Download table as:  ASCIITypeset image

In Figure 3 we show a comparison of stellar velocity dispersions obtained with and without the inclusion of photometric data in the fit. The two estimates are generally consistent within their quoted uncertainties, though there is a clear systematic offset in the sense that spectra-only fits return velocity dispersions that are ∼5% lower on average than those that also incorporate photometric data. This stems from the fact that the photometric data generally down-weight the youngest spectral templates—as one might expect from our a priori selection of galaxies based on their U − V and V − J colors—preferring instead solutions with smaller contributions from (rapidly rotating) early stellar types.

Zoom In Zoom Out Reset image size

In our high-redshift sample 50/58 galaxies fall within the HST WFC3/F160W imaging footprint of the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011), and for these objects we used the mosaics and composite PSFs described by Skelton et al. (2014).¹⁴ The remaining eight galaxies were observed separately using HST WFC3/F160W as part of HST-GO-12167 (PI: Franx; AEGIS 17300, AEGIS 21129, COSMOS 21434, COSMOS 20866, COSMOS 07447, UDS 53937, and UDS 55531) and HST-GO-13002 (PI: Williams; UDS 19627) for a single orbit each with total exposure times of 2611 or 2411 s, respectively. Level 2 data products were retrieved from the Hubble Legacy Archive (HLA),¹⁵ and we constructed empirical PSFs for these objects by stacking the images of bright unsaturated stars in each combined frame. We generated segmentation maps for the HLA images using SExtractor (Bertin & Arnouts 1996) with parameters similar to those given by Skelton et al. (2014) for 3D-HST.

We used galfit to model the surface brightness distribution of the primary galaxy, while also including in the fit neighboring galaxies with m_F160W < 25 and projected separations r_p < 5 (r_primary + r_neighbor). Initial estimates of the galaxy sizes, i.e., r_primary and r_neighbor, were taken from the SExtractor output. The local sky background for each object was estimated using the full image by first masking all pixels within 3 Kron radii of nearby sources using the ellipse parameters produced by SExtractor. We then identified the nearest 10,000 unmasked pixels as sky. An initial estimate of the background was taken as the mode of these sky pixels, which was then iteratively refined to obtain our final estimates of the local sky background. Postage stamps for individual objects were then extracted and the local sky background removed; the background level was subsequently held fixed during fitting. The structural parameters derived in this way are consistent with those available in the literature; a direct comparison with literature values is given in Appendix B.1. An example of our photometric modeling for COSMOS 30145 is shown in Figure 5, with figures for the remaining galaxies included in Appendix C.

Zoom In Zoom Out Reset image size

Although there are numerous existing catalogs of structural parameter measurements for the SDSS and GAMA (e.g., Simard et al. 2002; Kelvin et al. 2012; Meert et al. 2015), for consistency with our high-z data we chose to rederive these quantities using the methodology described above. We retrieved "corrected" r-band images from the SDSS Data Archive Server (DAS), along with their associated mask and PSF files. We then used SExtractor to generate segmentation images following the procedures described by Simard et al. (2011), and the local sky background for each source was estimated using the method described above. Individual postage stamps and PSFs¹⁶ for each galaxy were then extracted and the background removed. galfit was used to simultaneously fit the primary galaxy and any neighboring sources with m_r < 22 and r_p < 5 (r_primary + r_neighbor). All other sources were masked during the fit. A comparison of our measurements with several different literature catalogs can be found in Appendix B.2.

At both high and low redshift we derive sizes in fixed photometric bands (HST WFC3/F160W and SDSS r band, respectively), which probe different rest-frame wavelengths at different redshifts. In the presence of strong color gradients this shift in rest-frame wavelength can systematically bias our size measurements and must be taken into account. Following van der Wel et al. (2014), we define the corrected, in this case r-band, semimajor axis size ${r}_{{\rm{e}}}^{\mathrm{sma}}$ as

$$\begin{eqnarray}&&{r}_{{\rm{e}}}^{\mathrm{sma}}={r}_{{\rm{e}},\mathrm{obs}}^{\mathrm{sma}}{\left(\displaystyle \frac{1+z}{1+{z}_{p}}\right)}^{\tfrac{{\rm{\Delta }}\mathrm{log}{r}_{{\rm{e}}}}{{\rm{\Delta }}\mathrm{log}\lambda }},\end{eqnarray} \tag{ 3 }$$

where ${r}_{{\rm{e}},\mathrm{obs}}^{\mathrm{sma}}$ is the measured half-light size in either the WFC3/F160W or SDSS r-band filter and z_p is the "pivot redshift." z_p = 0 by definition for the GAMA/SDSS sample, as we are correcting to the rest-frame r-band size, while for F160W imaging z_p = 1.49. Kelvin et al. (2012) used GAMA data to show that ${\rm{\Delta }}\mathrm{log}{r}_{{\rm{e}}}/{\rm{\Delta }}\mathrm{log}\lambda =-0.3$ for early-type galaxies on average. Chan et al. (2016) and van der Wel et al. (2014) derive similar values based on their analyses of quiescent galaxies' high redshift, and we therefore adopt ${\rm{\Delta }}\mathrm{log}{r}_{{\rm{e}}}/{\rm{\Delta }}\mathrm{log}\lambda =-0.3$ for all galaxies in our sample. The typical correction derived in this way is of order 2%–3%, and we adopt these corrected r-band sizes for the remainder of this work.

While the single-component Sérsic fits described in Section 3.2.1 provide a straightforward summary of the overall surface brightness profile, Sérsic models have several drawbacks that complicate their use in constructing dynamical models. As well as providing a poor description of multicomponent profiles (e.g., bulge + disk), the coupling between inner and outer profile shapes makes the Sérsic models extremely sensitive to sky background: over-/undersubtraction of the sky level can significantly affect the inferred inner profile shape. In addition, with the exception of a few special cases, Sérsic profiles cannot be deprojected analytically, making their use for constructing dynamical models computationally expensive compared to simpler functional forms. In this context, modeling galaxies as a sum of individual Gaussian components—so-called multi-Gaussian expansion (MGE; Emsellem et al. 1994; Cappellari 2002)—provides a flexible description of surface brightness profiles that does not require any extrapolation of the profile to large radii, can accommodate multiple photometric components, and can be easily deprojected to obtain an estimate of the three-dimensional luminosity density (see Section 3.3.2).

The starting points for our MGE models were the background-subtracted postage stamps produced as described in Section 3.2.1. We used the results of our Sérsic model fits to subtract neighboring sources before identifying the primary object and producing binned two-dimensional surface brightness measurements using the find_galaxy and sectors_photometry routines described by Cappellari (2002).¹⁷ A model of the surface photometry in terms of nested Gaussians was then derived using the mge_fit_sectors method of Cappellari (2002). For high-redshift galaxies we constructed MGE-based PSF models per field using either composite PSFs provided by Skelton et al. (2014) for galaxies within the CANDELS/3D-HST footprint or else the stacked images of bright stars within the same field for stand-alone observations. In the right panel of Figure 5 we show a comparison of the observed and MGE-derived surface brightness contours for one object, COSMOS 30145.

MGE PSF models for the low-redshift SDSS/GAMA data were constructed on a galaxy-by-galaxy basis using the PSF extracted from the SDSS drField files. In all cases—that is, both high and low redshift—we tied the ellipticity of the fitted Gaussian components together to avoid large variations in the derived axis ratios for low surface brightness components; however, we confirmed that our results are not qualitatively sensitive to this assumption.

As outlined in Section 1, the tight relationship between size, stellar velocity dispersion, and mass for nearby early-type galaxies can be understood as a consequence of virial equilibrium, where for a pressure-supported system the total mass is given by

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}=\kappa (n)\displaystyle \frac{{\sigma }_{{\rm{e}}}^{2}{r}_{{\rm{e}}}^{\mathrm{sma}}}{G}.\end{eqnarray} \tag{ 4 }$$

Here κ(n) is the so-called virial coefficient, and in this case it is taken as an analytic function of Sérsic index that encapsulates the effects of structural and orbital nonhomology (e.g., Bertin et al. 2002; Cappellari et al. 2006). We adopt the relation derived by Cappellari et al. (2006) based on spherical, isotropic models,

$$\begin{eqnarray}&&\kappa (n)=8.87-0.831n+0.0241{n}^{2},\end{eqnarray} \tag{ 5 }$$

which has been shown to provide a reliable estimate of the total mass for nearby early-type galaxies in the SAURON and ATLAS^3D samples (e.g., Cappellari et al. 2006, 2013b). Note that in Equation (4) we used the semimajor axis size, ${r}_{{\rm{e}}}^{\mathrm{sma}}$, following the discussion of Cappellari et al. (2013a, their Figure 14). The semimajor axis size is expected to be more robust to systematic changes in galaxy shapes than the harmonic mean size (e.g., $\sqrt{{ab}}$, where a and b are the semimajor- and semiminor-axis sizes, respectively), especially for (thin) disk galaxies where the observed b/a is an indicator of inclination rather than intrinsic shape.

The assumptions of spherical symmetry and isotropy discussed above appear to be reasonable at low redshift; however, high-z quiescent galaxies are known to be flatter on average—that is, have intrinsically lower b/a—than their low-redshift counterparts (e.g., van der Wel et al. 2011; Chang et al. 2013), leading to a possible bias in their derived masses when using Equation (5). We therefore consider an alternative approach to computing dynamical masses based on the Jeans Anisotropic MGE (JAM) method discussed by Cappellari (2008), which allows us to relax these assumptions. The modeling requires as input the MGE-derived surface brightness profile described in Section 3.2.2 and a measurement of the stellar velocity dispersion (see Section 2.1.2 and 2.2).

Following Cappellari (2002), the deprojected luminosity density can be computed from the best-fit MGE decomposition given assumptions about the inclination, which is related to the intrinsic axis ratio of an oblate ellipsoid q_int by

$$\begin{eqnarray}&&\cos i=\sqrt{\displaystyle \frac{{q}_{\mathrm{obs}}^{2}-{q}_{\mathrm{int}}^{2}}{1-{q}_{\mathrm{int}}^{2}}},\end{eqnarray} \tag{ 6 }$$

where i is the inclination and q_obs is the observed axis ratio. Since in this work we are concerned with the sensitivity of our dynamical mass estimates to possible changes of the intrinsic axis ratio, we computed JAM models over a grid of q_int from 0.05 ≤ q_int ≤ min(0.95, q_obs) in steps of 0.05; unless otherwise stated, our results are based on marginalizing over q_int. In our default modeling we assumed that the velocity ellipsoid is marginally anisotropic with an anisotropy parameter $\beta \equiv 1-{\sigma }_{z}^{2}/{\sigma }_{R}^{2}=0.2$ (where z and R define directions parallel and perpendicular to the symmetry axis for an axisymmetric system) based on local early-type galaxies (e.g., Cappellari et al. 2007; Thomas et al. 2009). We explored possible systematic effects over a range of anisotropies from 0 ≤ β ≤ 0.8 and found that they resulted in variations of the derived dynamical masses of at most a few percent, consistent with previous results (e.g., Wolf et al. 2010; Dutton et al. 2013); all results are therefore quoted adopting our fiducial value of β = 0.2.

We adopt two different implementations of the JAM modeling procedure distinguished by their treatment of baryonic and dark matter components. In the first instance we assume that the total mass is proportional to the light at all radii, i.e., mass-follows-light (MFL). This provides a self-consistent estimate of the dynamical mass-to-light ratio (M/L)_MFL. MFL models have been shown to reliably recover the total mass within relatively small apertures (≲r_e), even in the presence of multiple mass components (e.g., Cappellari et al. 2006; Williams et al. 2010), and provide a baseline comparison for dynamical masses computed following Equation (5). A similar approach was used by Shetty & Cappellari (2014) to study quiescent galaxies at z ∼ 0.8 in the DEEP2 survey. In the MFL case the best-fitting value of (M/L)_MFL for a given combination of q_int and β is simply given by (M/L)_MFL = (σ_e,obs/σ_e,model)², where σ_e,obs is the observed aperture velocity dispersion and σ_e,model is the model prediction assuming (M/L)_MFL = 1. Instead, our second implementation includes an explicit dark matter component described by a spherical NFW halo profile (Navarro et al. 1996). With sufficient sampling of the velocity field it is possible to independently constrain the stellar mass-to-light ratio, (M/L)_*,NFW, and properties of the dark matter halo (e.g., Cappellari et al. 2013b; Übler et al. 2018). However, the aperture velocity dispersions used here cannot be used to break the degeneracy between stellar and dark matter components, leading us to impose additional constraints on the properties of the dark matter halo. Starting from our photometric estimates of galaxy stellar mass, we assigned dark matter halo masses based on the evolving stellar-to-halo mass relation derived by Moster et al. (2013). We then used the calculations of Diemer & Kravtsov (2015) to assign a halo concentration. This halo profile was then fed back into the JAM modeling procedure along with the MGE-based stellar density profile, and a grid search was used to determine the mass-to-light ratio of the stellar component, (M/L)_*,NFW, as a function of q_int and β.

In the explicit DM halo case we obtain an estimate of the dark matter fraction within r_e, f_DM[<r_e], defined as

$$\begin{eqnarray}&&{f}_{\mathrm{DM}}[\lt {r}_{{\rm{e}}}]=\displaystyle \frac{{M}_{\mathrm{DM}}}{{M}_{* ,\mathrm{NFW}}+{M}_{\mathrm{DM}}}.\end{eqnarray} \tag{ 7 }$$

We compute f_DM[<r_e] within a volume defined by $4\pi {r}_{{\rm{e}}}^{3}/3$, where for consistency with the literature r_e is the circularized half-light radius ($\equiv {r}_{{\rm{e}}}^{\mathrm{sma}}\times \sqrt{{q}_{\mathrm{obs}}}$), and the relevant masses are computed using the derived M/L values and deprojected MGE luminosity densities. For all galaxies we use the rest-frame r-band sizes computed following Equation (3).

The tendency for high-redshift quiescent galaxies to have lower dynamical-to-stellar mass ratios compared to low redshift has been reported in a number of previous studies (e.g., Toft et al. 2012; van de Sande et al. 2013; Belli et al. 2017) and is generally interpreted as reflecting a systematic decrease in the central dark matter fraction, f_DM[<r_e]. This decline in dynamical-to-stellar mass ratio appears to occur relatively smoothly with increasing redshift, as shown by a number of studies based on large spectroscopic surveys at z < 1 (e.g., Beifiori et al. 2014; Tortora et al. 2014, 2018). Figure 7 shows the cumulative distribution of dynamical-to-stellar mass ratio in the samples considered here. We find an offset in the mean dynamical-to-stellar mass ratio of −0.20 dex when moving to high redshift—$\mathrm{log}({M}_{\mathrm{MFL}}/{M}_{* })=0.29\pm 0.01$ at z = 0 compared to 0.09 ± 0.05 at 1.4 ≤ z ≤ 2.1—which is consistent with the results of previous studies. The magnitude of this offset is independent of the dynamical mass estimator used (e.g., M_MFL vs. M_vir) and does not change when considering only the oldest galaxies at z = 0 (shown as contours in Figure 6 and light-gray lines in Figure 7). The right panel of Figure 7 shows the same distribution of dynamical-to-stellar mass ratios as the left panel, but with individual measurements normalized by their uncertainties. These can be compared to the dotted–dashed (black) line, which shows the prediction for a standard normal distribution. Although nearly 40% of galaxies in the high-redshift sample have photometrically derived stellar masses that exceed their dynamical masses, given measurement uncertainties the overall distribution is consistent with a positive (albeit small) dynamical-to-stellar mass ratio on average.

Zoom In Zoom Out Reset image size

We can examine the evolution of f_DM[<r_e] more directly using our dynamical models that include an explicit dark matter component, where dark matter fractions are computed following Equation (7). In Figure 8 we show f_DM[<r_e] as a function of M_*,NFW, the dynamical mass of the stellar component. While there is significant uncertainty in the individual measurements of f_DM[<r_e] at high redshift, the overall trends support our interpretation of Figures 6 and 7 in terms of an evolution in the central dark matter fraction: galaxies at 1.4 < z < 2.1 have a mean f_DM[<r_e] = 6.6 ± 1.0%, a factor of >2 lower than galaxies of a similar mass in our SDSS/GAMA sample at z = 0 (f_DM[<r_e] ≈ 16.3 ± 0.3%; cf. 17% from Cappellari et al. 2013b). Furthermore, our low-redshift f_DM[<r_e] measurements are consistent with the values derived by Thomas et al. (2011) and Cappellari et al. (2013b) based on more detailed dynamical modeling of low-redshift galaxies, suggesting that the observed offset in f_DM[<r_e] between different redshifts is unlikely to be due to differences in the modeling approach. The above comparison between high- and low-redshift galaxies at fixed mass must nevertheless be made with some caution, as individual galaxies are expected to evolve from z = 2 to 0; we will revisit the evolution of f_DM[<r_e] using more carefully matched progenitor and descendant samples in Section 5.

Zoom In Zoom Out Reset image size

Using data from the SINS survey, Förster Schreiber et al. (2009) found that star-forming galaxies at z ∼ 2 are strongly baryon dominated, even for a Chabrier (2003) IMF, suggesting little room for either a bottom-heavy Salpeter IMF or significant dark matter. These results have recently been supported by kinematic data for hundreds of early star-forming disks in the KMOS^3D (Wisnioski et al. 2015, 2019) and MOSDEF (Kriek et al. 2015) surveys (e.g., Price et al. 2016; Wuyts et al. 2016; Lang et al. 2017; Übler et al. 2017), as well as the detailed analysis of outer rotation curves for individual high-redshift disks (Genzel et al. 2017; Genzel et al. 2020; but see also Tiley et al. 2019). For comparison, in Figure 8 we show the dark matter fractions derived by Genzel et al. (2017, shown as squares and upper limits for a subsample of high-redshift, star-forming galaxies), which are in good agreement with the f_DM[<r_e] measurements derived here.

One of our goals in comparing multiple dynamical mass estimators is to assess the impact and importance of different modeling assumptions on the inferred properties of high-redshift galaxies. In that regard, the main distinctions between M_vir and M_MFL are the assumptions of spherical symmetry and isotropy inherited through the application of Equations (4) and (5).

In practice, the dynamical masses derived here are only weakly dependent on changes of the anisotropy, β, at fixed q_int for a range of values consistent with local early-type galaxies (0 ≤ β ≤ 0.6; Cappellari et al. 2007). It is therefore unlikely that the assumption of isotropy has a significant impact on the results presented in Figures 6–8, particularly given the typical uncertainties on measurements of σ_e (20%–30%). On the other hand, changes in assumed galaxy structure—for example, from spherically symmetric to oblate and axisymmetric—can systematically bias dynamical mass estimates depending on the degree of intrinsic flattening and relative importance of rotation versus pressure support.

Crucially, there is growing observational evidence that quiescent galaxies at high redshift may indeed be rotationally flattened, violating the assumption of spherical symmetry inherent in Equation (5). Bezanson et al. (2018) showed that passive galaxies at z ∼ 1 have on average a higher proportion of rotational support (higher V/σ) than galaxies of the same mass at low redshift (see also van der Wel & van der Marel 2008). These results are consistent with the observed evolution of photometrically derived axis ratios over the same redshift range, which favor a significant portion of the quiescent galaxy population having 0.2 ≤ q_int ≤ 0.3 (van der Wel et al. 2011; Chang et al. 2013; Hill et al. 2019). Belli et al. (2017) argued that the dynamical masses of quiescent galaxies at z > 1.5 are statistically consistent with a factor of ∼2 increase in V/σ compared to z = 0 based on their correlation with observed axis ratios. More directly, a handful of strongly lensed passive galaxies at z > 2 have resolved kinematic profiles that are consistent with being rotationally flattened disks (e.g., Newman et al. 2015; Toft et al. 2017; Newman et al. 2018).

In the case of integrated (as opposed to resolved) absorption-line kinematics, rotation is expected to manifest as a dependence of the measured velocity dispersion—and, by extension, dynamical mass—on galaxy inclination. For an oblate model observed at inclination i with no azimuthal variation of the velocity ellipsoid (i.e., σ_ϕ/σ_r = 1, with σ_ϕ and σ_r describing velocity dispersion in the azimuthal and radial directions, respectively), the second moment of the velocity distribution σ_obs can be written as

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{obs}}^{2}={\sigma }^{2}(1-\beta {\cos }^{2}i)\\ \qquad \ \times \ \left[1+{\left(\displaystyle \frac{V}{\sigma }\right)}_{e}^{2}\displaystyle \frac{{\sin }^{2}i}{{\left(1-\beta {\cos }^{2}i\right)}^{2}}\right],\end{array}\end{eqnarray} \tag{ 8 }$$

where V and σ are the flux-weighted mean circular velocity and velocity dispersion within the effective radius for the edge-on case (i = 90°) with $(V/\sigma {)}_{e}^{2}\equiv {V}_{{\rm{e}}}^{2}/{\sigma }_{{\rm{e}}}^{2}$, and β is the anisotropy parameter as defined in Section 3.3.2. Following Belli et al. (2017), substituting Equation (8) into Equation (4) and normalizing by the dynamical mass predicted for the face-on case (i.e., i = 0°) gives

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{M}_{\mathrm{vir}}(i)}{{M}_{\mathrm{vir}}(i=0^\circ )}=\displaystyle \frac{1-\beta {\cos }^{2}i}{1-\beta }\\ \qquad \ \times \ \left[1+{\left(\displaystyle \frac{V}{\sigma }\right)}_{e}^{2}\displaystyle \frac{{\sin }^{2}i}{{\left(1-\beta {\cos }^{2}i\right)}^{2}}\right],\end{array}\end{eqnarray} \tag{ 9 }$$

with the relationship between q_int, q_obs, and i given by Equation (6). In the isotropic case where β = 0, Equation (9) reduces to Equation (5) of Belli et al. (2017) modulo a factor γ = 1.51.¹⁸

In Figure 9 we show the dynamical-to-stellar mass ratio as a function of q_obs for both M_vir and M_MFL estimates. In the virial theorem case we find evidence for a weak negative correlation between M_vir and q_obs, with a Spearman rank coefficient ρ = −0.23 ± 0.10 (p = 0.007), while for the MFL models the correlation is not significant (ρ = −0.14 ± 0.10; p = 0.068). Individual galaxies are color-coded according to their Sérsic indices as derived from the profile fits described in Section 3.2.1. In contrast to Belli et al. (2017), we find no significant dependence of the dynamical-to-stellar mass ratio on Sérsic index in either case with ρ = −0.20 ± 0.10 (p = 0.04) and −0.08 ± 0.10 (p = 0.44), which appears to preclude a simple exclusion of disk-dominated systems based on their structure and motivates a more detailed examination of the correlation between q_obs and dynamical-to-stellar mass ratio.

Zoom In Zoom Out Reset image size

Lines in the top panel of Figure 9 show predicted behavior of the dynamical-to-stellar mass ratio for a galaxy with q_int = 0.25 observed at different inclinations as given by Equation (9). We set β = 0.7(1−q_int) based on the results of Cappellari et al. (2007) and Emsellem et al. (2011) for nearby fast-rotating early-type galaxies. In the case of an oblate system, V/σ and anisotropy are related by (Binney & Tremaine 1987; Binney 2005)

$$\begin{eqnarray*}&&\beta =1-\displaystyle \frac{1+{\left(V/\sigma \right)}^{2}}{1-\alpha {\left(V/\sigma \right)}^{2}}\left(\displaystyle \frac{{W}_{{zz}}}{{W}_{{xx}}}\right),\end{eqnarray*}$$

where α is a dimensionless number that quantifies the contribution of streaming motions to the line-of-sight velocity dispersion and (W_zz/W_xx) is a shape parameter related to the intrinsic axis ratio q as

$$\begin{eqnarray*}&&\left(\displaystyle \frac{{W}_{{zz}}}{{W}_{{xx}}}\right)=\displaystyle \frac{2\left(q\sqrt{1-{q}^{2}}-{q}^{2}\arccos \,q\right)}{\arccos q-q\sqrt{1-{q}^{2}}}.\end{eqnarray*}$$

We adopt a value of α = 0.15, which provides a good description for nearby galaxies (Cappellari et al. 2007). Models are offset in M_vir/M_SPS to reflect a range of dynamical-to-stellar mass ratios. The predicted trends qualitatively reproduce the observed correlation between mass ratio and q_obs, supporting previous statistical evidence of rotational support among a fraction of high-redshift quiescent galaxies (e.g., Belli et al. 2017).

In the bottom panel of Figure 9 the correlation between q_obs and dynamical-to-stellar mass ratio for MFL models is notably weaker than in the virial theorem case, both visually and as measured by Spearman ρ, though galaxies with higher q_obs still tend toward lower dynamical-to-stellar mass ratios. Unlike the virial theorem case, we can explicitly test the impact of intrinsic structure on our MFL mass estimates through application of a prior on q_int in our modeling.¹⁹ The vertical lines in the bottom panel of Figure 9 show the effect of assuming that galaxies are intrinsically flat, with q_int = 0.25, as opposed to the default case where we adopt a uniform prior on q_int. For galaxies with low q_obs, the effect of assuming a different intrinsic structure is minimal, but for galaxies with q_obs > 0.6 (i < 55°) the estimated dynamical mass can increase by as much as ∼65% (0.22 dex), with a median increase of ∼15% (0.06 dex). The resulting correlation between dynamical-to-stellar mass ratio and q_obs is also flatter, with ρ = 0.03 ± 0.06 (p = 0.95). Furthermore, assuming an intrinsically flat structure for these objects reduces the number of galaxies with dynamical masses significantly lower than their photometrically derived stellar mass M_SPS.

In summary, the data considered here support the conclusions of previous studies suggesting that rotational support is prevalent among quiescent galaxies at high redshift (e.g., Chang et al. 2013; Newman et al. 2015; Belli et al. 2017; Toft et al. 2017; Newman et al. 2018; Hill et al. 2019). While we expect rotational flattening to have a minimal impact on dynamical mass estimates for galaxies with q_obs < 0.6, galaxies with high q_obs can have their dynamical masses underestimated by 0.2 dex or more depending on their intrinsic structure (i.e., if they are intrinsically spherical vs. flattened systems viewed face-on). As mentioned in Section 4.1, such a discrepancy between intrinsic and assumed structure can at least partially explain those galaxies in our sample that have dynamical masses formally less than their photometrically derived stellar mass, though it may not be the only factor affecting this comparison. Finally, we note that enforcing an intrinsically flat structure for all galaxies in our sample (e.g., q_int = 0.25) shifts the results presented in Section 4.1.1 toward lower central dark matter fractions and cannot explain the apparent evolution of f_DM[<r_e] without also appealing to significant changes in the stellar initial mass function (see Section 4.2).

One of the key assumptions in computing M_*,NFW is that the dark matter halo component is well represented by an NFW profile, with no accounting for the possible influence of baryons on the dark matter profile shape. However, if the timescale for galaxy formation is long compared to the halo dynamical time, then the halo is expected to contract adiabatically as a result of baryonic collapse (e.g., Blumenthal et al. 1986; Gnedin et al. 2004). Dutton et al. (2016) argue that the dark matter halo can contract or expand depending on the relative balance of inflows, outflows, and feedback (see also Lovell et al. 2018), suggesting that our assumption of a static halo may bias the derived values of M_*,NFW and, by extension, α. In this section we therefore explore a broader set of dynamical models that explicitly probe the effect of a variable dark matter halo response on our results.

In the case of spherical symmetry and circular dark matter particle orbits, the adiabatic invariant is given by rM_tot(r)—where M_tot(r) is the total (baryonic plus dark matter) mass within radius r—so that r_f/r_i = M_tot,i(r_i)/M_tot,f(r_f). Therefore, given an initial mass distribution M_tot,i(r) and a final baryonic mass profile M_bar,f(r), we can derive the final dark matter profile M_DM,f(r). Here we assume that the initial dark matter distribution is described by an NFW profile with mass and concentration parameter set by the scaling relations adopted in Section 3.3, and that the baryonic mass is distributed in the same way, i.e., M_bar,i(r) = f_barM_tot,i(r), with f_bar set by the stellar-to-halo mass relation. The final baryonic profile M_bar,f(r) is given by the deprojected MGE luminosity density scaled to match the galaxy stellar mass. We note that this assumes that star formation is distributed throughout the halo and is therefore likely an upper limit on the expected contraction.

If we assume no shell crossing of the dark matter, then M_DM,i(r_i) = M_DM,f(r_f), and the final mapping between r_f and r_i is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}={M}_{\mathrm{DM},{\rm{i}}}({r}_{{\rm{i}}})/\left[{M}_{\mathrm{bar},{\rm{f}}}({r}_{{\rm{f}}})+(1-{f}_{\mathrm{bar}}){M}_{\mathrm{DM},{\rm{i}}}({r}_{{\rm{i}}})\right],\end{eqnarray} \tag{ 10 }$$

with ${\rm{\Gamma }}={({r}_{f}/{r}_{i})}^{\nu }$ following the generalized contraction formula suggested by Dutton et al. (2007). In this framework ν = 1 reproduces the standard adiabatic contraction derived by Blumenthal et al. (1986), while ν = 0.8 reproduces the modified contraction scenario described by Gnedin et al. (2004). ν = 0 is equivalent to an unmodified NFW profile. We also include a milder model for the halo response derived from the NIHAO simulations discussed by Dutton et al. (2015), such that

$$\begin{eqnarray}&&{r}_{{\rm{f}}}/{r}_{{\rm{i}}}=0.5+0.5{\left({M}_{\mathrm{tot},{\rm{i}}}/{M}_{\mathrm{tot},{\rm{f}}}\right)}^{2}.\end{eqnarray} \tag{ 11 }$$

For each contraction model we solve for the mapping between r_f and r_i, and we use this modified dark matter profile as input to the JAM modeling procedure. While we do not explicitly include any models for halo expansion (i.e., ν < 0), our default MFL models can be interpreted as the extreme case where dark matter is completely evacuated within r_e, setting an upper limit for the dynamical effects of an expanded halo.

In Figure 11 we show the distributions of α derived for these different models of halo response, along with the standard MFL and NFW cases described in the previous section. The expected trend of a decreasing stellar contribution—that is, a lighter overall IMF normalization—with increasingly strong halo contraction is clearly visible, with pure adiabatic contraction (e.g., Blumenthal et al. 1986) predicting stellar masses that are a factor of >3 lighter than those obtained assuming a Chabrier (2003) IMF. Even the mildest model for halo response, Dutton et al. (2015), predicts values of α that are ∼60% lighter than Chabrier, and all of the contraction models considered here predict IMF normalizations that are lighter than observed or inferred for nearby stellar systems (e.g., Chabrier 2003; Bastian et al. 2010). Taken together, the results in Figure 11 suggest that any contraction of the dark matter halo due to gas inflow should be counterbalanced by equally violent outflows during the formation process.

Zoom In Zoom Out Reset image size

We showed in Section 4.1.1 and Figure 8 that, at fixed stellar mass, f_DM[<r_e] increases by a factor of ∼2 from high redshift to the present day. However, based on the matched progenitor and descendant samples described in Section 5, we find that the growth of f_DM[<r_e] for individual galaxies is likely even larger, reaching 0.64 ± 0.05 dex (0.52 ± 0.07 dex) on average for progenitors and descendants matched on evolving stellar mass (fixed σ_e).

Just as major and minor mergers are expected to have different effects on a given galaxy's evolution in size, stellar mass, and stellar velocity dispersion (e.g., Figure 13), they also have a distinct influence on the evolution of f_DM. Hilz et al. (2012, 2013) show that minor mergers can dramatically increase f_DM[<r_e], so that a factor of 2 growth in stellar mass can nearly double the central dark matter fraction. In contrast, a single equal-mass (major) merger may only increase f_DM[<r_e] by 50%. This difference is driven by the relatively efficient growth of sizes in minor compared to major mergers; for an NFW-like halo the dark matter mass within a given (small) radius r scales roughly with the virial mass of the halo as ${M}_{\mathrm{DM}}(r)\propto {M}_{\mathrm{vir}}{r}^{2}$ (Boylan-Kolchin et al. 2005). In an equal-mass merger where both r_e and the stellar mass within r_e double, the dark matter mass within r_e can increase by up to a factor of 8 (assuming that the stellar mass and halo mass increase by a similar amount). By comparison, a similar doubling of stellar mass through multiple minor mergers can increase r_e by a factor of 4 and the enclosed dark matter mass by a factor of ≳30.

These differences are particularly apparent when viewed in terms of stellar mass surface density, Σ_e ($\equiv {M}_{\mathrm{SPS}}/2\pi {r}_{{\rm{e}}}^{2}$), and f_DM[<r_e], which we show in Figure 14. We find that Σ_e and f_DM[<r_e] are anticorrelated, in the sense that those galaxies with the highest stellar mass surface densities have the lowest f_DM[<r_e]; a similar anticorrelation has been demonstrated previously for both late- and early-type galaxies (e.g., McGaugh 2005; Sonnenfeld et al. 2015). Figure 14 additionally shows that Σ_e decreases by ∼1 dex on average from z > 1.5 to z ≈ 0, which is a direct consequence of the size and stellar mass evolution discussed in Section 5.1 (see also Table 5 and Figure 12). This can be compared to the results of Remus et al. (2017), who show that simulated galaxies follow well-defined tracks in Σ_e–f_DM regardless of redshift (shown as dashed and dotted–dashed lines in Figure 14). Nevertheless, the simulation results roughly reproduce the observational trends at any given epoch, with the primary difference between the two simulations discussed by Remus et al. (2017) being their treatment of black hole feedback.

Zoom In Zoom Out Reset image size

We can use simple scaling relations to better understand galaxies' expected evolution in Figure 14 given various assumptions. Following Figure 12, we adopt a simple model for size growth during mergers such that ${\rm{\Delta }}{r}_{{\rm{e}}}\propto {\rm{\Delta }}{M}_{* }^{\alpha }$, with α = 1 or 2 for major and minor mergers, respectively. We additionally assume that the enclosed dark matter mass scales with the total (virial) mass of the halo and radius as M_DM(r) ∝ M_virr², and that the stellar and halo mass grow at the same rate (${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{vir}}={\rm{\Delta }}\mathrm{log}{M}_{* }$). Arrows in the right panel of Figure 14 show the predicted evolution for a galaxy doubling its stellar mass through either a single major merger or successive 10:1 (minor) mergers. As expected, the efficient size growth associated with minor mergers in our toy model drives rapid evolution in both Σ_e and f_DM[<r_e]. However, while minor mergers can explain the observed decrease in Σ_e from high to low redshift, our toy model overpredicts the increase of f_DM[<r_e]. Tortora et al. (2018) report a similar result for galaxies at z ≈ 0.7 and suggest that allowing for variation in the stellar-to-halo mass ratio of accreted galaxies (i.e., ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{vir}}\ne {\rm{\Delta }}\mathrm{log}{M}_{* }$) can help to lessen the tension between predicted and observed evolution in f_DM. This is particularly likely for the massive quiescent galaxies in our high-redshift sample, which are expected to trace the most massive halos at their respective redshifts (e.g., Lin et al. 2019).