X-shooter Spectroscopy and HST Imaging of 15 Massive Quiescent Galaxies at z ≳ 2

The MASSIVE Survey samples the most massive K-band-selected early-type galaxies within the local 108 Mpc northern hemisphere (Ma et al. 2014). These galaxies have central stellar ages suggesting a formation epoch at z > 2 (Greene et al. 2015). Given the similar selection for our MQGs at z > 2, stellar masses, and inferred formation epoch, this sample is adopted as the local reference sample. This sample is further motivated in Section 5.1.

The extinction-corrected absolute K-band magnitudes listed in Table 3 of Ma et al. (2014) are converted into stellar masses using Equation (1) in van de Sande et al. (2019). The NASA-Sloan Atlas semimajor axis optical effective radii, also listed in Table 3 of Ma et al. (2014), are used. These were derived from 2D Sérsic (1968) fits with Sérsic parameters varying between n = 2 and 6. For the galaxies where this is not available, the infrared Two Micron All Sky Survey (2MASS) measurements were used to convert these to semimajor axis optical effective radii using Equation (4) in Ma et al. (2014). These sizes were derived from single Sérsic and de Vaucouleurs profile fits (n = 4). The effective velocity dispersion measurements used are reported in Veale et al. (2018). They were estimated using the MILES stellar library (Falcón-Barroso et al. 2011) together with the penalized pixel-fitting algorithm (pPXF; Cappellari & Emsellem 2004). Finally, the average luminosity-weighted stellar velocity dispersion within the effective radius is adopted.

Mounted on the VLT, X-shooter is a single-object echelle spectrograph covering 3000–25000 Å with three arms: UVB (2936–5930 Å), VIS (5253–10489 Å), and NIR (9827–24807 Å). We were granted 35 and 57 service mode hours in P86 and P93, respectively (PI: Toft). The latter carried over and finished in period 96. The observations were completed using a default nodding mode to ensure a robust sky subtraction of the NIR band, probing the rest-frame optical part of the spectra for the z ∼ 2 quiescent galaxies. The majority of the P86/P93 observations (89/96%) were completed with an average airmass-corrected DIMM seeing of 08 in the NIR arm. The telluric standard stars were observed close to the science observations, both in airmass and in time, to mimic the conditions of the sky and optimize the atmospheric absorption correction. The P86/P93 observations for the NIR (VIS) frames were executed with 480/900 s (314/863 s) exposures, a 09 × 11'' slit configuration, and—for the P93 sample only—including the K-band blocking filter. We aligned the slit along the galaxy's major axis in the UltraVISTA K-band images avoiding bright nearby sources.

The data are reduced using a wrapper of the ESO X-shooter pipeline (Modigliani et al. 2010; Sparre 2015), along with customized modifications (Zabl et al. 2015). Beyond the standard pipeline processing steps for the NIR arm in nodding mode, we account for the spatial variations of the background level outside of the orders in each raw science frame by removing the median level obtained from the illuminated areas from each row of pixels in the detector. The 2D VIS and NIR individual science frames are corrected for telluric absorption with a customized and publicly available wrapper²² (Selsing et al. 2016) of the pPXF (Cappellari & Emsellem 2004), based on the PHOENIX stellar atmosphere library (Husser et al. 2013). A response function is constructed by modeling the atmosphere during the science exposures, and each individual observation block (OB) is corrected.

Finally, individual OBs are combined into an optimally weighted 2D spectrum removing flux outliers using 3σ and 5σ median clipping for the VIS and NIR, respectively. Bad pixels automatically flagged during the reduction are also excluded. Furthermore, off-trace emission is flagged and excluded in the construction of the OBs from UV-105842, UV-171687, and UV-155853 to minimize the contamination from surrounding sources. The 1D spectrum is optimally extracted (Horne 1986). Flux corrections are made, anchoring the synthetic photometry to the total magnitudes from the latest COSMOS15 catalog (Laigle et al. 2016, their Section 3.3) and accounting for point-spread function (PSF) matching in different bands and for the Galactic extinction. The H- and I-band magnitudes are used to compute independent aperture correction factors for the NIR and VIS spectra, respectively.

Eleven orbits of HST/WFC3 with HST-GO-14721 (PI: Conselice) are allocated to observe the rest-frame optical images, H_F160W, for UDS 19627 and the 10 galaxies in the P93 sample. The P86 sample is covered by the following programs: CP-1243752 (HST-GO-12440; PI: Faber) and CP-561356 (HST-HLA-14114; PI: van Dokkum) with WFC3 and CP-1291751 and CP-540713 with HST/NICMOS (HST-HLA-9999; PI: Scoville).

The WFC3/H_F160W data are reduced using the "grism redshift and line" analysis software, Grizli,²³ which is an end-to-end processing code for WFC3/IR data using ASTRODRIZZLE.²⁴ The starting point is the standard calibrated images downloaded from the MAST archive (FLT extension images). The calibrated images are 1014 × 1014 pixels with 013 pixel⁻¹. For each visit, there are four dithered exposures that are combined using Grizli. The resulting products for each visit are aligned, background-subtracted, and drizzled images with 006 pixel⁻¹. The NICMOS data for CP-1291751 and CP-540713 are reduced in a similar manner with ASTRODRIZZLE.

We make ample use of the 14 broadband COSMOS photometry measurements from the Laigle et al. (2016) catalog, covering the full UV-to-NIR wavelength range to model our stellar populations in Section 4.3. The total magnitudes are adopted using the method described in Appendix A.2 by the same authors. Complementary to the UV-to-NIR photometry, we check the available deep X-ray Chandra imaging (Marchesi et al. 2016) and the "superdeblended" FIR catalog (Jin et al. 2018), superseding the previous 24 μm catalog (Le Floc'h et al. 2009) used in the selection of P93. This new implementation adopts active priors from the Spitzer/MIPS 24 μm and radio observations to deblend the low-resolution imaging from Herschel/PACS and SPIRE, SCUBA2, AzTEC, and MAMBO. The sources are cross-checked with the Galaxy Evolution Explorer far- and near-UV data from Zamojski et al. (2007) and Capak et al. (2007). This search for UV or X-ray counterparts results in no detections for any of our galaxies. On the other hand, we do find hints of mid-infrared (MIR) and radio emission from part of the sample, as detailed in Section 4.4.2 and discussed in Section 6.4. Similar UV-to-NIR multiwavelength coverage is available for UDS 19627. For an in-depth discussion of the available photometric data for this object, see Toft et al. (2012).

Thirteen of 15 galaxies have HST I_F814W imaging that is part of the COSMOS publicly released data (Scoville et al. 2007; Koekemoer et al. 2007). It covers ∼2 deg² of the sky with the Advanced Camera for Surveys (ACS) in the I band and comprises 81 tiles. Each tile is observed in four dithered exposures that are combined to produce a pixel scale of 003 pixel⁻¹ and a PSF of 0095 at FWHM. COSMOS images reach a point-source limiting depth of AB(F814W) = 27.2 (5σ).

All spectra of targeted sources (P86 and P93) show prominent hydrogen absorption features, which are typical of evolved stellar populations (see Figure 2). The stellar absorption features are modeled using pPXF, and both the line-of-sight velocity centroid (i.e., the spectroscopic redshift) and the line-of-sight stellar velocity dispersion (hereafter "velocity dispersion") are measured.

The initial redshift and velocity dispersion guess is obtained from running pPXF with the Bruzual & Charlot (2003, hereafter BC03) stellar population library. The stellar population analysis is performed with complex SF history (SFH) fitting of the spectra and SED (see Section 4.3) adopting this initial estimate. The resulting best-fit model is confirmed to be stable against perturbations of Δσ = ±100 km s⁻¹. The velocity dispersion measurement is refined by rerunning pPXF with a non-velocity-broadened best-fit stellar population model.

The spectra and best-fit model are convolved to the same resolution (FWHM = 3.2 Å) and rebinned to a constant velocity scale without additional interpolation. Low-order additive (a = 2) and multiplicative (m = 2) correction polynomials are fit over the rest-frame range 3750–5950 Å. The JH band gap and the regions where emission lines might be expected²⁵ are excluded while also masking out bad pixels.

The associated systematic and statistical errors are quantified by varying the wavelength range, correction polynomials, and stellar libraries (see details in Appendix B), similar to the method used in Toft et al. (2017). In all cases (P86 and P93), we determine secure redshifts, and for 10/14 galaxies, we estimate robust velocity dispersions. The spectroscopic redshifts and velocity dispersion measurements, along with the combined systematic and statistical errors (Appendix B), are listed in Tables 2 and 3, respectively. In Table 3, we also list the velocity dispersion for UDS 19627 derived in Toft et al. (2012). In Figure 3, the derived spectroscopic redshifts are compared with the photometric estimates from Muzzin et al. (2013a) and Laigle et al. (2016). Using the normalized median absolute deviation (σ_NMAD) from Brammer et al. (2008), no catastrophic outliers are found, except for photometric redshifts being systematically below the spectroscopic redshifts for both catalogs, finding a better agreement for Muzzin et al. (2013a).

Zoom In Zoom Out Reset image size

Table 2.  The Stellar Population Model Parameters

Target ID	zspec	log(M*/M⊙)	log(Age/yr)	A(g)	SFRSSP (M⊙ yr−1)	SFRopt (M⊙ yr−1)	SFR24 (M⊙ yr−1)a
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
UV-108899	2.2312	${11.62}_{-0.18}^{+0.16}$	${9.15}_{-0.30}^{+0.27}$	${0.38}_{-0.38}^{+1.00}$	<13	6 ± 4 ([O ii])	<15
UV-250513	2.0814	${11.51}_{-0.19}^{+0.18}$	${9.16}_{-0.31}^{+0.27}$	${0.38}_{-0.38}^{+1.02}$	<12	<3 (Hα)	<13
CP-561356	2.6963	${11.62}_{-0.20}^{+0.21}$	${9.14}_{-0.32}^{+0.28}$	${0.62}_{-0.62}^{+1.17}$	<86	<19 ([O ii])	<90⊕
UV-105842	2.0195	${11.68}_{-0.17}^{+0.16}$	${9.19}_{-0.33}^{+0.26}$	${0.81}_{-0.81}^{+1.04}$	<17	<2 (Hα)	19 ± 5
UV-171687	2.1020	${11.51}_{-0.19}^{+0.18}$	${9.13}_{-0.32}^{+0.28}$	${0.64}_{-0.64}^{+1.13}$	<24	<3 (Hα)	$26\pm {6}^{\diamond \oplus }$
UV-90676	2.4781	${11.78}_{-0.18}^{+0.17}$	${9.09}_{-0.29}^{+0.29}$	${0.41}_{-0.41}^{+0.99}$	<88	<6 ([O ii])	$\lt {92}^{\diamond \oplus }$
CP-1291751	2.0253	${11.24}_{-0.22}^{+0.23}$	${9.17}_{-0.33}^{+0.27}$	${0.82}_{-0.82}^{+1.26}$	<17	<2 (Hα)	18 ± 4
UV-155853	1.9816	${11.62}_{-0.17}^{+0.18}$	${9.23}_{-0.33}^{+0.24}$	${0.88}_{-0.86}^{+1.04}$	<14	<4 (Hα)	<15
UV-171060	2.0995	${11.48}_{-0.17}^{+0.16}$	${9.16}_{-0.31}^{+0.27}$	${0.41}_{-0.41}^{+1.03}$	<14	<2 (Hα)	$\lt {15}^{\diamond }$
UV-230929	2.1679	${11.48}_{-0.16}^{+0.16}$	${9.10}_{-0.28}^{+0.28}$	${0.22}_{-0.22}^{+0.89}$	$\lt 6$	<4 ([O ii])	<7
UV-239220	2.0057	${11.57}_{-0.20}^{+0.20}$	${9.18}_{-0.33}^{+0.26}$	${0.66}_{-0.66}^{+1.14}$	<19	35 ± 15 (Hα)	$21\pm {4}^{\diamond \oplus }$
UV-773654	2.0328	${11.59}_{-0.20}^{+0.19}$	${9.20}_{-0.33}^{+0.26}$	${0.68}_{-0.68}^{+1.13}$	<12	<2 (Hα)	$13\pm {3}^{\diamond \oplus }$
CP-1243752	2.0903	${11.79}_{-0.17}^{+0.17}$	${9.23}_{-0.32}^{+0.24}$	${0.76}_{-0.76}^{+1.06}$	<11	<2 (Hα)	$\lt 12$
CP-540713	2.0409	${11.26}_{-0.23}^{+0.22}$	${9.16}_{-0.32}^{+0.27}$	${0.57}_{-0.57}^{+1.19}$	<10	<2 (Hα)	<12
UDS 19627b	2.0389	${11.37}_{-0.10}^{+0.13}$	${9.08}_{-0.10}^{+0.11}$	${0.77}_{-0.32}^{+0.36}$	⋯	<6 (Hα)	<40c

Notes. Column 1: target ID. Column 2: spectroscopic redshift. Column 3: stellar mass. Column 4: mass-weighted stellar age. Column 5: extinction in g band. Column 6: 3σ upper limit percentiles (representing 99.73% Gaussian confidence intervals) of the stellar population modeled SFR/100 Myr distribution. Column 7: 3σ SFR upper limits based on Hα or [O ii] λ3727 (Section 4.4.1). Column 8: 24 μm estimated SFR (see Section 4.4.2).

^aGalaxies with detections in 1.4 (${}^{\diamond }$) and 3 (^⊕) GHz are indicated with matching symbols. ^bThe values listed for UDS 19627 are from Toft et al. (2012). From this study, the A(v) extinction instead of the listed A(g) is quoted. ^cThe 2σ 24 μm SFR upper limit using the method from Franx et al. (2008).

Download table as:  ASCIITypeset image

Table 3.  Summary of Structural Properties

Target ID	zphot	zspec	σ (km s−1)	Rmaj (kpc)	n	q	RFluxa	${\mathrm{log}}_{10}({M}_{* ,c}/{M}_{\odot })$	log10(Mdyn/M⊙)	Classb
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
UV-108899-1c	2.19	2.2312	470 ± 82	1.36 ± 0.14	2.51	0.44	0.56	${11.38}_{-0.18}^{+0.16}$	⋯	Pd
UV-108899-2c	⋯	⋯	⋯	3.38 ± 0.34	7.15	0.56	0.44	${11.26}_{-0.17}^{+0.17}$	⋯	Pd
UV-250513-1c	2.03	2.0814	174 ± 44	3.84 ± 0.38	4.00	0.59	0.55	${11.26}_{-0.19}^{+0.17}$	⋯	Pd
UV-250513-2c	⋯	⋯	⋯	1.60 ± 0.16	4.00	0.63	0.45	${11.16}_{-0.18}^{+0.17}$	⋯	Pd
CP-561356-1c	2.43	2.6963	280 ± 128	4.14 ± 0.41	1.45	0.62	0.71	${11.47}_{-0.20}^{+0.21}$	⋯	Pd
CP-561356-2c	⋯	⋯	⋯	2.78 ± 0.28	0.90	0.39	0.29	${11.09}_{-0.21}^{+0.21}$	⋯	Pd
UV-105842-1	1.93	2.0195	263 ± 57	4.07 ± 0.41	3.51	0.51	1.00	${11.68}_{-0.17}^{+0.16}$	11.61 ± 0.19	P
UV-171687-1	2.04	2.1020	182 ± 50	5.12 ± 0.51	4.00	0.77	1.00	${11.51}_{-0.19}^{+0.18}$	11.37 ± 0.24	P
UV-90676	2.57	2.4781	347 ± 82	5.22 ± 0.51	4.98	0.61	1.00	${11.78}_{-0.18}^{+0.17}$	11.89 ± 0.21	P
CP-1291751	2.06	2.0253	⋯	3.47 ± 0.35	3.59	0.67	1.00	${11.24}_{-0.22}^{+0.23}$	⋯	P
UV-155853	1.96	1.9816	247 ± 30	4.55 ± 0.46	3.62	0.85	1.00	${11.62}_{-0.17}^{+0.18}$	11.60 ± 0.11	E
UV-171060	2.02	2.0995	⋯	1.73 ± 0.17	4.00	0.54	1.00	${11.48}_{-0.17}^{+0.16}$	⋯	E
UV-230929	2.09	2.1679	252 ± 21	1.74 ± 0.17	3.01	0.73	1.00	${11.48}_{-0.15}^{+0.16}$	11.23 ± 0.08	E
UV-239220	2.00	2.0057	⋯	5.35 ± 0.54	4.21	0.62	1.00	${11.57}_{-0.20}^{+0.20}$	⋯	E
UV-773654	1.96	2.0328	⋯	3.77 ± 0.38	3.34	0.84	1.00	${11.59}_{-0.19}^{+0.19}$	⋯	E
CP-1243752	2.01	2.0903	350 ± 53	2.85 ± 0.29	4.50	0.79	1.00	${11.79}_{-0.17}^{+0.17}$	11.66 ± 0.14	E
CP-540713	1.98	2.0409	353 ± 97	1.65 ± 0.17	0.96	0.79	1.00	${11.26}_{-0.23}^{+0.22}$	11.59 ± 0.24	E
UDS 19627	2.02	2.0389	318 ± 53	2.00 ± 0.20	3.32	0.51	1.00	${11.37}_{-0.10}^{+0.13}$	11.48 ± 0.15	E

Notes. Column 1: target ID. Column 2: photometric redshift from Muzzin et al. (2013a). Column 3: spectroscopic redshift (Section 4.1). Column 4: stellar velocity dispersion measurement (Section 4.1). Column 5: effective semimajor axis (Section 4.5). Column 6: Sérsic index (Sérsic 1968). Column 7: axis ratio q = b/a (defined by the ratio between the semiminor and major axes). Column 8: flux scaling used to estimate the corrected stellar mass, ${\mathrm{log}}_{10}{M}_{* ,c}$ (Section 4.6). Column 9: corrected stellar mass (Section 4.6). Column 10: dynamical mass (Section 5.4). Column 11: morphological classification (P: peculiar; E: elliptical) from Conselice et al. (2005).

^aRelative flux ratio = ${F}_{i}/({F}_{i}+{F}_{j})$. ^bGalaxies marked with (d) are classified as major mergers in Section 4.6. ^cDouble sources have similar photometric and spectroscopic redshifts, as well as the stellar velocity dispersions estimated from their composite spectra. ^dGalaxies classified as major mergers in Section 4.6.

Download table as:  ASCIITypeset image

No on-source nebular line emission is detected at 3σ for any object in the sample. For UV-108899 and UV-239220, we find indications of emission (∼2σ) from [O ii] 3726.2,3728.9 and Hα 6563, respectively. In Appendix C, we discuss the specifics of the fitting method and list, in Table 2, the SFR and uncertainties from [O ii] and Hα (Kennicutt 1998). Furthermore, spatially offset line emission is observed in four (UV-155853, UV-171687, UV-171060, UV-105842) 2D spectra coinciding with close-proximity sources. In 3/4 cases, this emission arises from foreground or background sources around UV-155853, UV-171687, and UV-171060 (Appendix A.2). The latter source northeast of UV-105842 shows significant [O ii] 3726.2,3728.9 Å, [O iii] 4959,5007 Å, and Hα emission with a matching redshift of z = 2.0124. This corresponds to a velocity offset of 2130 ± 120 km s⁻¹ from UV-105842. If purely due to galaxy motion, such an offset suggests that the two sources are not gravitationally bound at the time of observation. Another explanation of the asymmetric morphology might be a high-redshift analog of the locally observed offset AGN (Comerford & Greene 2014), likely caused by a recent merger event.

In order to put constraints on the physical parameters of the stellar populations, the VIS+NIR X-shooter spectra and the broadband photometry are fit with the Bayesian approach from Gallazzi et al. (2005; recently revised in Zibetti et al. 2017) using the derived spectroscopic redshift. Spectral regions of poor atmospheric transitions are not included in the calculation. Before fitting, the models are convolved by the initial velocity dispersion estimated in Section 4.1.

Models are obtained by convolving the latest revision of the BC03 simple stellar population (SSP) models using the MILES stellar libraries (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) with a large Monte Carlo library of SFHs, metal-enrichment histories, and dust attenuations. The prior distribution of models is the one described in Zibetti et al. (2017), but here it is limited to 50,000 models with formation ages younger than 5 Gyr to be consistent with the high redshift of our galaxies. A full description of the model library is given in Zibetti et al. (2017); however, the most relevant information is summarized here.

The SFHs are modeled with a continuous component parameterized à la Sandage (1986),²⁶ thus allowing for both an increasing and a decreasing SFH phase, on top of which random bursts of SF are added. Stellar metallicity evolves according to the SFH (see Zibetti et al. 2019), with initial and final values randomly generated in the range 0.02–2.5 Z_⊙. Finally, for 75% of the models, the effect of dust attenuation is included following the model of Charlot & Fall (2000) that separates the contributions of the birth clouds affecting stars younger than 10⁷ yr and the ISM affecting stars of all ages.

The Bayesian modeling approach assumes the likelihood of each model to be ∝exp(−χ²/2). The probability distribution function (PDF) of each physical parameter of interest is computed by weighing the prior distribution of the models in a given parameter by their likelihood, marginalizing over all other parameters. We additionally used the information from the MIR flux limit to restrict the sample of acceptable models to those that have an SFR consistent with 24 μm based upper limits and detections (see Section 4.4.2). The median and 16th and 84th percentiles of the PDFs are adopted as the fiducial estimates and their uncertainties for each parameter. Note that this approach allows the derivation of realistic uncertainties on the key physical parameters, accounting for both the observational errors and the intrinsic degeneracies among different parameters.

The stellar mass, mass-weighted mean stellar age, effective dust attenuation (A(g)), and SFR averaged over the last 100 Myr for our sample are reported in Table 2. In this table, the SFR limits from nebular line and 24 μm emission (see Sections 4.4.1 and 4.4.2) are also listed. Stellar masses are within the range of ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=11.23-11.79$, with a median of 11.57. Compared to Belli et al. (2017), this sample is on average more massive, which is reflected by the brighter K-band magnitudes (see Figure 1). Such MQGs have also been found over a larger area in Arcila-Osejo et al. (2019). The SFR limits and dust-corrected stellar masses, together with the mean stellar mass–weighted ages of ∼1.4 Gyr, confirm the expectations from the selection that this is, in fact, a sample of massive recently quenched galaxies. Three of the galaxies are double sources, and the stellar masses are corrected in Section 4.6.

4.4.1. Rest-frame Optical Emission Lines

In order to confirm the quiescent nature of our galaxies, upper limits on [O ii] λ3727 and Hα emission are measured. These are converted into upper limits of the unobscured SFRs following Equations (2) and (3) in Kennicutt (1998), under the assumptions of solar abundance ratio and that all massive SF is traced by ionized gas. A 3σ flux upper limit is determined by summing up the flux error density squared over a region of Δλ = 1000 km s⁻¹ (similar to 300–500 km s⁻¹ line dispersions):

$$\begin{eqnarray}&&{F}_{3\sigma \mathrm{limit}}=3\sqrt{\sum {\sigma }_{\mathrm{flux}}^{2}\delta {\lambda }^{2}}.\end{eqnarray} \tag{ 1 }$$

Here σ_flux and δλ are the flux uncertainty and bin size, respectively. Note that we do not introduce any dust extinction in this conversion, as this is largely unconstrained (see Section 4.4.3 for an estimated upper limit on the dust extinction). We find unobscured SFR upper limits that are consistent with the expectation that these galaxies are quiescent ($-10\lt {\mathrm{log}}_{10}(\mathrm{sSFR}/\mathrm{yr})\lt -11.5$). The difference between the [O ii] λ3727 and Hα SFR limits is <11 M_* yr⁻¹, and in Table 2, the lowest SFR upper limits are listed.

4.4.2. MIR Emission

4.4.3. Comparison of Different SF Tracers

Figure 4 shows the position of the sample of MQGs in the log(SFR)–log(M_*) main sequence at z = 2. For reference, the SFR main sequence at matching redshift from Speagle et al. (2014) is shown, extrapolated to the stellar mass range log₁₀(M_*/M_⊙) > 11.1 covered by our galaxies.

Zoom In Zoom Out Reset image size

The rest-frame optical SFR limits are systematically lower than the MIR estimates (both probing 10–100 Myr timescales). This suggests that the star-forming regions are either strongly obscured and/or dominated by AGN dust heating (Fumagalli et al. 2014). Under the assumption of no AGN contribution to the heating that produces the MIR emission (see also Section 6.4), the dust extinction is estimated by comparing the obscured and unobscured SFR estimates, resulting in a mean extinction of A(v) < 1–2, consistent with our SED fit derived A(g) (g-band) extinction. In order to judge if a significant contribution to the MIR heating arises from AGNs, we check whether there are any radio counterparts detected in Jin et al. (2018). Radio emission is detected in five sources at 1.4 GHz and five sources at 3 GHz (indicated with symbols in Table 2), showing that AGN heating could be responsible for the elevated MIR SFR estimates. Further treatment of the radio detections will be part of a future paper (I. Cortzen at al. 2020, in preparation).

The SFRs derived from our stellar population analysis (Section 4.3) are consistent with SFR ∼ 0 M_* yr⁻¹ for all galaxies in our sample. In Table 1, we list the 3σ upper limits on these SFR limits. However, even considering the most conservative upper limits on the SFR from the 24 μm emission, our sample of MQGs lies ∼2 dex below the SFR main sequence at their redshifts, confirming their quiescent nature (see Figure 4).

The 2D stellar light distribution traced by HST/WFC3 H_F160W imaging is modeled with the χ²-minimization fitting code GALFIT (Peng et al. 2002) in order to retrieve the structural parameters of our sample of MQGs. A first run of SExtractor (Bertin & Arnouts 1996) allows us to detect the objects in each field and obtain an initial guess for the structural parameters. A postage stamp for each target is constructed such that it encloses an ellipse with a major axis 2.5 times the Kron radius obtained by SExtractor. The local sky level in each stamp is calculated using Galapagos (Barden et al. 2012). This sky level is passed to GALFIT and kept fixed during the fitting. For the WFC3 data, a combination of the TINYTIM²⁷ simulated PSF and an empirical stacked star PSF are used. For the NICMOS data, an empirical stacked PSF is used.

Finally, GALFIT is run on each postage stamp, adopting a flexible Sérsic profile for every source (Sérsic 1968),

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{e}\exp \left\{-{\kappa }_{n}\left[{\left(\displaystyle \frac{R}{{R}_{e}}\right)}^{1/n}-1\right]\right\}.\end{eqnarray} \tag{ 2 }$$

The parameter R_e is the effective radius enclosing half of the flux from the model light profile, Σ(R_e) is the surface brightness at the effective radius, and n is the Sérsic index. The quantity κ_n is a function of the Sérsic index, which defines the global curvature of the light profile and is obtained by solving the equation Γ(2n) = 2γ(2n, κ_n), where Γ and γ are, respectively, the gamma function and incomplete gamma function.

GALFIT is run several times to ensure that the solutions correspond to a global minimum in the minimization algorithm for each image by varying the initial guesses of the total magnitude, effective radius, and Sérsic index. The parameters are constrained so as to avoid any unphysical solutions (effective radius >0.2 pixels, q > 0.1, 0.5 < n < 8). Initially, all targets are fit with n as a free parameter. In unstable cases where the maximum or minimum n is reached, the images fixing the Sérsic index at either n = 1 or n = 4 are refit, choosing the model providing the smallest χ² as the best-fit solution. These two choices represent realistic descriptions of an early-type galaxy dominated by either a disk or a bulge. Throughout the whole fitting procedure, neighboring objects are either modeled or masked, depending on their proximity to the main target. A 10% measurement uncertainty on the size (van der Wel et al. 2008; Newman et al. 2012) is shown to be a fair representation. This conservative error estimate is thus adopted. The semimajor axis, ${R}_{{\rm{e}},\mathrm{maj}}$, is adopted as the effective radius in the following sections. The best-fit parameters and their uncertainties are reported in Table 3.

In Figure 5, we present the rest-frame UV (I_F814W) and optical (H_F160W) images along with the GALFIT model and residual. The morphologies of these galaxies are classified in the H_F160W image according to Conselice et al. (2005), and they fall into the two categories for quiescent systems: ellipticals (E) and peculiars (P). When available, the spectroscopic observations are used to determine the distance in redshift space to objects that fall in the X-shooter slit (see Section 4.2). The majority of sources turn out not to be associated with the central galaxy. Nine of 15 galaxies are categorized as elliptical galaxies, while the remaining ones are categorized as peculiar galaxies with major mergers (UV-108899, UV-250513, CP-561356), minor mergers (UV-105842, CP-1291751), and/or strong tidal/post-merger features (UV-105842, UV-90676). The galaxies UV-108899, UV-250513, and CP-561356 are confirmed as ongoing major mergers in the following section. The classifications and morphological parameters are listed in Table 3.

Zoom In Zoom Out Reset image size

The RB color images in Figure 2 reveal that three galaxies (UV-108899, UV-250513, CP-561356) appear to be double systems. The spectra, shown in the same figure, are the total extraction of the combined light from the two galaxies. These objects are within close proximity, and the light in the reduced 2D frames is blended to an unknown extent (due to limited seeing). At the expense of drastically decreasing the signal-to-noise ratio (S/N), an attempt to separate the sources and determine whether their individual redshift measurements can confirm their proximity is made.

For each system, the resolved 1D HST H_F160W light profile (extracted parallel to the X-shooter slit) is overlaid on the wavelength-collapsed 2D spectrum trace. A double Gaussian profile fit allowed us to gauge the amount of blending and make a conservative extraction of each individual galaxy, minimizing cross-source contamination. In Figure 6, the individual extractions and the best fit to the composite spectrum from Section 4.3 are shown.

Zoom In Zoom Out Reset image size

Because of the low S/N of the individual conservative flux extractions, the estimation of the velocity offset is refrained, since it would be dominated by large uncertainties. However, the galaxies are within close physical proximity due to the matching absorption lines shown in the figure and can be considered ongoing quiescent (dry) major mergers. This confirmation is important as, in the following section, it can be used to correct their stellar masses, prior to presenting them in the mass–size plane (see Section 5.3).

Spectroscopic confirmation allows us to deblend the composite stellar mass of each system using the H_F160W magnitude as a proxy for tracing the bulk of the stars in the galaxies. The GALFIT-modeled H_F160W flux ratio supports the fact that these galaxies are major mergers with mass ratios of 1: 1–3. We used the flux ratio to correct the stellar masses as

$$\begin{eqnarray}&&{M}_{* ,i}={M}_{* ,\mathrm{tot}}\displaystyle \frac{{F}_{i}}{{F}_{i}+{F}_{j}}={M}_{* ,\mathrm{tot}}{R}_{\mathrm{Flux}},\end{eqnarray} \tag{ 3 }$$

where i and j refer to the two merging galaxies and F is the total flux from GALFIT. The corrected stellar masses (${M}_{* ,c}$) and the relative flux ratio scaling, R_Flux, are listed in Table 3, with sources names matching the numbering in Figure 5. Following this correction, the galaxies are still classified as MQGs with stellar masses, ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\gt 11$.

A major issue preventing us from deriving a consistent evolutionary picture connecting galaxy populations across time is the "progenitor bias" problem (e.g., van Dokkum & Franx 1996; Carollo et al. 2013). When comparing galaxies across time, the implicit assumption is that the high-redshift sample contains all progenitors of the low-redshift reference sample. However, the fraction of quenched galaxies has been found to grow over time (Buitrago et al. 2013), introducing an unknown bias when comparing samples of galaxies across different epochs.

One approach that has been suggested to minimize the progenitor bias is comparing the evolution of galaxies at fixed velocity dispersions (see, e.g., Belli et al. 2014a). Archaeological studies (Graves et al. 2009; van der Wel et al. 2009; Bezanson et al. 2012) find evidence suggesting that the velocity dispersion in quiescent galaxies remains approximately unchanged across cosmic time (z < 1.5). In such a scenario, the velocity dispersion must be weakly affected by the average merger history, which, according to the numerical study by Hilz et al. (2012), occurs for minor merger–driven evolution. A detailed discussion on fixed velocity dispersion evolution is given in Belli et al. (2014a, 2017). Another way to minimize the progenitor bias has been to study galaxy populations at a constant cumulative number density (CND) instead of a fixed velocity dispersion or stellar mass (see, e.g., Mundy et al. 2015). This approach is introduced in van Dokkum et al. (2010) and refined further in Behroozi et al. (2013) and Leja et al. (2013). In Section 2.1, a sample of massive galaxies with central stellar population ages suggesting formation at z > 2 is introduced. This sample is volume-limited and represents the most massive early-type systems observed in the local universe. In order to draw a meaningful comparison, a subgroup of the most massive galaxies at z = 0 is selected and matched with the CND at z = 2. This will hereafter be referred to as the "fixed" CND. This approach is based on the assumption that the rank of galaxies within the stellar mass function is not strongly affected across cosmic time. This occurs if the stellar mass continuously grows from z = 2 to 0, implying the availability of surrounding material to accrete (or events that trigger secondary SF, although this is not expected for the MQGs at z > 2; Brammer et al. 2011; Behroozi et al. 2013; Muzzin et al. 2013b; Marchesini et al. 2014).

First, the CND of massive ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.2$) UVJ quiescent galaxies in the redshift range 1.9 < z < 2.5 is estimated using the Muzzin et al. (2013a) catalog. The stellar mass limit represents the lower limit on the standard deviation of the mean stellar mass from the sample of galaxies studied in this paper. Our sample is 22% stellar mass complete using these selection criteria. We count 58 galaxies inside a comoving volume spanned by this redshift range, giving $n(\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.2)=9.7\,\times \,{10}^{-6}\,{\mathrm{Mpc}}^{-3}$.

The MASSIVE galaxy sample is trimmed starting from the most massive object of the survey and including progressively less massive systems until we reach the fixed CND of the massive UVJ quiescent galaxies at z ∼ 2. The final fixed CND-matched MASSIVE sample consists of the 25 most massive local elliptical galaxies with stellar masses of log(M_*/M_⊙) > 11.70. The fixed CND-matched MASSIVE sample is hereafter referred to as "MASSIVE(n)." The MASSIVE(n) sample is considered a minimal progenitor-biased sample and used as our local reference sample in Sections 5.2–5.4.

The CND evolution suffers from large uncertainties from individual merger histories causing scatter in the mass rank, which is the main uncertainty for the highest stellar masses (Behroozi et al. 2013; Torrey et al. 2017). In Torrey et al. (2017), they estimated the mass rank scatter for log₁₀(M_*/M_⊙) > 11 galaxies in Illustris (Genel et al. 2014; Nelson et al. 2015) by forward modeling of the CND. Their forward modeling, referred to as density distribution functions, is well described by a lognormal distribution, and the uncertainties can thus be treated as confidence intervals. For massive galaxies, the dominating uncertainty, the mass rank scatter, introduces a uncertainty of a factor of ∼2 (within 80% confidence intervals) on the CND following the evolution from z = 2 to 0. In Behroozi et al. (2013), they found a similar uncertainty for the fixed CND evolution. This uncertainty on the CND evolution from the mass rank scatter is adopted and used to repeat the selection of the local reference sample, resulting in a corresponding uncertainty on the limit of the stellar mass cut $\mathrm{log}({M}_{* }/{M}_{\odot })\gt {11.70}_{-0.10}^{+0.07}$ and thus the number of galaxies in the local reference sample.

As an alternative approach to the fixed CND matching, the probabilistic approach from Wellons & Torrey (2017) is used to estimate the CND at z = 0. In Appendix D, the results (from Figures 7–9) for both a fixed and a probabilistic CND-matching approach is presented. The choice of CND-matching method does not affect the qualitative results of this paper presented in the following sections.

Zoom In Zoom Out Reset image size

In Figure 7, the stellar velocity dispersion–size plane that allows us to study the kinematic evolution of MQGs from z = 2 to 0 is presented. The ongoing major merger galaxies are included to show that their incorrect composite dispersion measurements increase the scatter if not properly accounted for.

The mean velocity dispersion of the sample studied in this paper is 289 ± 58 km s⁻¹ (without major mergers). This is consistent with previous z > 2 MQG literature (see studies shown in Figure 7) with a mean dispersion of 272 ± 31 km s⁻¹. Our velocity dispersion and size measurements (including other structural parameters) for CP-1242752 (indicated by the blue square in Figure 7) are consistent with previously published values (van de Sande et al. 2013; Belli et al. 2014b, 2017; Kriek et al. 2016).

Comparing the median dispersion of our study to that of the local MASSIVE(n) sample, a shallow or no kinematic evolution from z = 2 to 0 is found. In Figure 7, significant effective size evolution consistent with earlier findings is observed (Newman et al. 2012; van der Wel et al. 2014). The effect of the mass rank scatter on the fixed CND matching is shown as the purple shading around the median evolution. This shading outlines the variation on the median when using the upper and lower limits of the CND matching (based on the stellar mass cut $\mathrm{log}({M}_{* }/{M}_{\odot })\gt {11.70}_{-0.10}^{+0.07}$) from the mass rank scatter.

Half of the morphologies of compact massive galaxies at z ∼ 2 have been suggested to be disk-dominated (van der Wel et al. 2011). So far, only one spatially resolved study of a rotating disk quiescent galaxy at this epoch has been discovered (Geier et al. 2013; Toft et al. 2017; Newman et al. 2018). The line-of-sight measured velocity broadening of the absorption lines could be a combination of both rotation and dispersion in the presence of a disk-dominated system (see an analytical prescription in Belli et al. 2017). Care must therefore be taken when comparing z > 2 spatially unresolved dispersion with resolved local measurements.

Wuyts et al. (2011) showed that the stellar light distribution of galaxies, measured by the Sérsic index, traces well the log(SFR)–log(M_*) relation, separating disk and spheroidal galaxies by n = 2.5 at z < 1.5. Under the assumption that this is valid at z = 2, we classify our galaxies by Sérsic index and find that 92% of our galaxies have spheroidal (n > 2.5) morphologies (when excluding the ongoing major mergers). If the Sérsic index n > 2.5 is a good tracer of dispersion-dominated systems at z > 2, it suggests that our sample of galaxy dispersion measurements is not strongly contaminated by rotation. Bezanson et al. (2018) find no strong correlation between Sérsic index and rotational support. Instead they find that galaxies with stellar masses log(M_*/M_⊙) > 11.2 (similar to our systems) appear to not be rotationally supported.

A recent study by Veale et al. (2018) presents the spatially resolved velocity dispersion measurements for the MASSIVE Survey sample. Here log₁₀(M_*/M_⊙) > 11.7 galaxies (similar to our stellar mass cut of the MASSIVE(n) sample) all have velocity dispersions in the range 200 km s⁻¹ < σ < 350 km s⁻¹ at all radii (<15–30 kpc). This rules out the possibility that the shallow dispersion evolution comparison is driven by spatial resolution. A comparison to the fixed CND-matched MASSIVE(n) sample establishes that the dispersion remains nearly unchanged.

A negligible median dispersion evolution of our MQGs across the last 10 billion yr (z = 2–0) is found in Figure 7 when comparing to the MASSIVE(n) sample. In the absence of spatially resolved spectroscopy, we make use of the morphological classification and high stellar masses, which suggests that our kinematics are unlikely to be strongly contaminated by rotation. Studying the evolution of galaxies at fixed dispersion has been suggested as a method to minimize progenitor bias (e.g., Belli et al. 2014b).

In Figure 8, the stellar mass–size plane (${\mathrm{log}}_{10}{M}_{* }\mbox{--}{R}_{{\rm{e}},\mathrm{maj}}$) is presented, which allows us to study the structural and stellar mass evolution of MQGs since z ∼ 2. The three ongoing major merger galaxies with resolved sizes of individual galaxies (Section 4.5) and their flux-corrected stellar masses (Section 4.6) are shown in the figure. The post-merger stellar masses and sizes of these are predicted using the argument of virialization from Bezanson et al. (2009). The resulting position of post-merger galaxies is consistent with the average locus of the most massive (log₁₀(M_*/M_⊙) > 11.5) individual galaxies in our sample, showing that a way to form the most MQGs in our sample could be major quiescent-to-quiescent dry galaxy mergers (Naab et al. 2006).

Zoom In Zoom Out Reset image size

A best-fit relation to the galaxies in this study, including the major merger–separated galaxies, reveals a shallower slope than what is found in the van der Wel et al. (2014) z = 2.25 mass–size relation but in better agreement with Mowla et al. (2018). The best-fit parameters, using a similar parameterization ($r/\mathrm{kpc}=A{\left({M}_{* }/(5\times {10}^{10}\right)}^{\alpha }$), are log(A) = 0.19 and α = 0.42. The stellar mass for CP-1243752 (blue square in Figure 8) is consistent within a 1σ standard deviation with van de Sande et al. (2013) and Belli et al. (2017) and within 1.1σ for the stellar mass published in Kriek et al. (2016).

The distribution of our sample shows that z > 2 MQGs are ∼two times more compact than objects with the same stellar mass in the local universe (Shen et al. 2003), which is a well-established result in previous works (van de Sande et al. 2013; Belli et al. 2017). The median stellar mass and size for our (MASSIVE(n)) sample, log(M_*/M_⊙) = 11.48 (11.77) and ${R}_{{\rm{e}},\mathrm{maj}}/\mathrm{kpc}=3.42\ (13.55)$, show that a doubling (∼0.3 dex) in stellar mass and a factor of 4 in size evolution are required to bring the two samples into qualitative agreement.

Using the method from Bezanson et al. (2009) for predicting stellar mass and size growth, minor and major merger tracks are shown in the mass–size plane. The median mass–size evolution between our z > 2 MQGs and the local MASSIVE(n) sample could be explained by minor merger–predicted size and stellar mass growth. The tracks start at the median size and stellar mass of our sample (only red symbols). The qualitative conclusions remain the same when using a mean instead of a median or changing the choice of reference (with/without the major merger galaxies).

The median logarithmic mass–size slope is $\alpha ={1.78}_{-0.29}^{+0.37}$ ($r\propto {M}_{* }^{\alpha }$). The uncertainties are determined based on the CND mass rank scatter, shown as the purple shaded area in Figure 8. This confirms the suggestion that minor mergers (α = 2), compared to major mergers (α = 1), are the preferred evolutionary path in the mass–size plane.

In line with earlier studies (van de Sande et al. 2013; van der Wel et al. 2014; Belli et al. 2017; Mowla et al. 2018), we find that our sample of z > 2 MQGs is compact in the stellar mass–size plane and further suggests that minor merger–driven size evolution (Bluck et al. 2012; Hilz et al. 2012, 2013; Newman et al. 2012; Oogi & Habe 2013; Fagioli et al. 2016) is preferred when comparing to the fixed CND-matched MASSIVE(n) sample.

In Figure 9, the dynamical-to-stellar mass relation for MQGs is plotted in order to study the interplay between the stellar and total (dynamical) mass potential over time. The dynamical mass derived from the Jeans equation (Jeans 1902) for symmetrical systems is as follows:

$$\begin{eqnarray}&&M(r)=\beta \displaystyle \frac{{R}_{{\rm{e}},\mathrm{maj}}{\sigma }^{2}}{G}.\end{eqnarray} \tag{ 4 }$$

Here ${R}_{{\rm{e}},\mathrm{maj}}$ is the effective semimajor axis, σ is the stellar velocity dispersion, G is the gravitational constant, and β is a parameter incorporating the full complexity of collisionless systems with radial-dependent parameters of density, dispersion, and velocity anisotropy. Following Cappellari et al. (2006), β(n) = 8.87–0.831n + 0.0241n² is adopted, where n is the Sérsic index (Sérsic 1968). The representation of β is a good approximation for symmetric systems, such as an elliptical galaxy that is well represented by a de Vaucouleurs profile. Taylor et al. (2010a) and Cappellari et al. (2013) showed that using such a parameterization of β yields dynamical masses in better agreement with the stellar masses when the sizes are estimated using a 2D Sérsic fitting method, rather than a fixed value of β.

Zoom In Zoom Out Reset image size

The galaxies of this study are consistent with the stellar-to-dynamical mass ratio, M_*/M_dyn < 1, within the large uncertainties. A ratio >1 is referred to as a nonphysical (forbidden) region, where the total mass is smaller than the mass of the stars. The galaxy, UV-230929, is located in this region at 1.1σ standard deviation from the M_*/M_dyn = 1 relation. Unfortunately, our large uncertainties prohibit trustworthy estimates of the total dust+gas mass for our sample. In Belli et al. (2017), it is suggested that dispersion-dominated systems with n > 2.5 lie closer to the M_*/M_dyn = 1 relation at z ∼ 2.

Compared to previous z > 2 MQG studies (see legend in Figure 9), our sample occupies a similar dynamical mass range but has larger stellar masses. This is further discussed in Section 6.3. The dynamical mass for CP-1243752 (indicated by a blue square) is consistent with the previous measurements in van de Sande et al. (2013) and Belli et al. (2017).

A comparison between our study and the MASSIVE(n) sample is made to learn about the fixed CND evolution in the dynamical–stellar mass plane. The median evolution in Figure 9 illustrates that the dynamical mass evolves 2× faster than the stellar mass within the effective radii. This means that the galaxies evolve such that the M_*/M_dyn ratio decreases from z = 2 to 0.

The minor and major merger evolution are shown for constant velocity dispersion evolution (${\rm{\Delta }}r\propto {M}_{* }^{\alpha }$), with α = 1 for major merger and 2 for minor merger evolution. This is motivated by the shallow/constant dispersion evolution found in Section 5.2 when also comparing to the MASSIVE(n) sample. The median evolution from z = 2 to the present day prefers the minor merger–predicted evolution when comparing our study to the MASSIVE(n) sample in the dynamical–stellar mass plane.

The median evolution from our study to the MASSIVE(n) sample at the present day in the dynamical–stellar mass plane is consistent with minor merger evolution that is similar to what is found in Figures 7 and 8.

Three galaxies in our sample, initially unresolved in ground-based imaging, are found in HST images to be double sources and confirmed with X-shooter to be ongoing major merger systems (see Section 4.6). In this section, we discuss how the high major merger fraction (6/18) affects the prevalence of minor merger structural evolution of MQGs (found in Section 5) and if the high fraction could be caused by a selection bias.

Following the definition in Man et al. (2012), we find a pair fraction of 20% ± 12% when assuming no projected sources ($\langle {N}_{\mathrm{projected}}\rangle =0$) and using the Poisson error estimate. In COSMOS and UDS, a pair fraction of 10% major mergers is found for massive (log(M_*/M_⊙) > 11) galaxies at z = 2 (Mundy et al. 2017). In the case that the observed major mergers are representative of the complete sample of MQGs, we can estimate the number of major mergers each galaxy undergoes (N_merger) following the prescription in Man et al. (2016b). Under the assumption that the merger rate is constant from z = 2 to 0, Equation (3) in Man et al. (2016b) can be written as ${N}_{\mathrm{merger}}={\rm{\Delta }}{{tf}}_{\mathrm{pair}}/{t}_{\mathrm{obs}}$. The pair fraction f_pair = 0.2, observation time t_obs = 0.8 Gyr (from z = 1.9 to 2.5), and time of evolution Δt ∼ 10 Gyr is used to estimate the number of mergers from z = 2 to 0. These numbers reproduce a major merger rate of ∼1 for a pair fraction of 10%, similar to what was suggested in Man et al. (2012). For a 20% pair fraction, we find that, on average, each galaxy undergoes N_merger = 2.5 mergers in the range 0 < z < 2. This number of 1:1 major mergers corresponds to a stellar mass increase of 0.5 dex, which is too high compared to the median stellar mass of the MASSIVE(n) sample (see Figure 8). In the case of a 10% pair fraction, the stellar mass increase is consistent with the average stellar mass of the MASSIVE(n) sample; however, in this case, another mechanism must then be in place to produce the large size growth observed between our sample and the MASSIVE(n) galaxies.

Our sample was selected to be UVJ quiescent and K-band bright, which could have introduced a bias for ground-based unresolved bright red systems like quiescent-to-quiescent galaxy major mergers (see also Section 6.5). See also Mowla et al. (2018) and Marsan et al. (2019), who addressed the issue of close pairs of MQGs at z ∼ 2. If this selection bias is responsible for the high pair fractions, this could explain why we observe that the stellar mass–size evolution from z = 2 to 0 is dominated by minor mergers (see Figure 8). The majority of our major merger targets are in the low stellar mass end of our sample (log(M_*/M_⊙) < 11.5). This could indicate that a possible way to produce ultramassive (log(M_*/M_⊙) > 11.5) QGs could be via quiescent-to-quiescent galaxy major mergers at z > 2. A scenario involving early-time major and late-time minor merger evolution will be testable with larger samples of MQGs at z > 2.5.

In Figure 8, a slope of $\alpha ={1.78}_{-0.29}^{+0.37}$ is found for the mass–size evolution of our MQGs from z = 2 to 0. Such an evolution can be interpreted using the analytical framework from Bezanson et al. (2009) and Naab et al. (2009), who found that minor merger–driven growth is needed to produce a mass–size slope of α = 2. An extended numerical treatment from Hilz et al. (2012) finds that when including the effect of escaping particles (a process arising from virialization following merger interaction), they recovered a steeper mass–size slope (α ≈ 2.4) alongside a constant dispersion evolution for minor merger–driven growth. Such a scenario could explain the observed size growth and shallow dispersion evolution observed.

The scenario presented in Hilz et al. (2012) occurs for two-component (stellar+halo) systems when they undergo 1:10 minor merger evolution. They reproduce the structural evolution found in Bezanson et al. (2009) and Naab et al. (2009) when simulating minor merger evolution of stellar-only systems. According to Hilz et al. (2012), this suggests that the growth of the dark matter halo is an important ingredient necessary to cause the shallow dispersion evolution together with the expected size growth evolution we find in this study. Moreover, Hilz et al. (2012) showed that major mergers increase the dispersion and size proportional to the stellar mass. This is not what is found when comparing the size and dispersion evolution with the MASSIVE(n) sample (see Figures 7 and 8). In the minor merger scenario, the velocity dispersion would be maintained in the inner region of the galaxy, as additional stellar mass is accreted in the outer parts from tidally stripped satellite systems. Over time, this would change the stellar light distribution on the outskirts of the galaxy, causing continuous growth of the half-light radius (van Dokkum et al. 2010; Hill et al. 2017).

In UV-105842, we may be observing a direct example of the minor merger–driven size increase. A small satellite system within close (spectroscopically confirmed) proximity of the central galaxy is found. Based on the flux ratio estimated from the GALFIT modeling, we estimate a stellar mass ratio of 1:${12}_{-3}^{+6}$ for this minor merger, consistent with the average 1:16 ratio estimated by Newman et al. (2012). To double its stellar mass (as suggested by the median ∼0.3 dex increase derived for our sample), the galaxy would need to go through ∼12 such minor mergers between z = 2 and 0. Other minor merger stellar mass ratios of 1:5, 1:10, and 1:20, suggested by Hilz et al. (2013) and Bédorf & Portegies Zwart (2013), would correspond to 5, 10, and 20 minor mergers between z = 2 and 0 for a similar stellar mass increase. In Man et al. (2016b), issues related to the translation of the H-band flux ratio to a stellar mass ratio (e.g., due to M/L ratio variation in galaxies) directly affecting the above argument are discussed.

Many observational (Bluck et al. 2012; McLure et al. 2013; Fagioli et al. 2016; Matharu et al. 2019; Zahid et al. 2019) and numerical (Naab et al. 2009, 2014; Oser et al. 2012; Oogi & Habe 2013; Tapia et al. 2014; Remus et al. 2017) studies find that minor mergers could be a dominant process for the size growth of massive galaxies, but they may not be able to explain the full size evolution (Cimatti et al. 2012; Newman et al. 2012). Feedback processes have been shown to also affect the size growth (e.g., Lackner et al. 2012; Hirschmann et al. 2013). Specifically, AGN feedback is shown by modern simulations to be necessary to reproduce the observed size evolution (see Dubois et al. 2013; Choi et al. 2018).

We found that the dynamical-to-stellar mass ratio shown in Figure 9 increases by a factor of 2 within MQGs from z = 2 to 0. This could be attributed to either IMF changes of the stellar population (Cappellari et al. 2012) affecting the stellar mass estimates or an increase in the dark matter fraction within the effective half-light radius.

Numerical simulations find that minor merger–driven evolution alters the distribution of stars over time from a core to a core–envelope system by accretion of particles in the outskirts of the galaxy (Hopkins et al. 2009; Hilz et al. 2012, 2013; Frigo & Balcells 2017; Lagos et al. 2018). A consequence of this is that the central dispersion remains constant while the half-light radius grows, encompassing a larger part of the dark matter halo and effectively increasing the dark matter fraction over time (Hilz et al. 2012).

A mass–size evolution similar to what we find is, according to Hilz et al. (2013), caused by a massive dark matter halo that drives the accretion of dry (collisionless) minor mergers at large radii through tidal stripping. This inside-out growth increases the effective half-mass radius to encompass dark matter–dominated regions, which might explain the increase of the dynamical-to-stellar mass fraction within the half-light radius that we observe.

Care must be taken when interpreting the observations in terms of idealized numerical simulations. However, Remus et al. (2017) also found that the central dark matter fraction increases with decreasing redshift when comparing different cosmological simulations. Furthermore, observational evidence for inside-out growth in massive galaxies is presented in Szomoru et al. (2012).

In Figure 9, we find that our sample is consistent with the dynamical-to-stellar mass ratio of one suggesting low dark matter fractions at z ∼ 2. For a stellar mass increase of 0.3 dex (similar to our median evolution), Hilz et al. (2012) predicted a dark matter fraction increase of ∼70% within the effective radius. If we assume that the mass of the galaxy consists only of dark matter and stars, we can estimate the dark matter mass fractions (${M}_{\mathrm{DM}}/{M}_{\mathrm{dyn}}=1-{M}_{* }/{M}_{\mathrm{dyn}}$) from the dynamical-to-stellar median ratio at z = 2 and 0 to be ${7}_{-7}^{+24} \%$ and 56% ± 8%, respectively. This suggests an increase of the dark matter fraction within the effective radius of 17%–64%. Note, however, that this increase cannot purely be associated with the dark matter from the minor mergers, as the growing half-light radius similarly encompasses more of the central dark matter halo that therefore also contributes to this increase.

According to Remus et al. (2017), the mass growth of massive galaxies can be explained by two stages: (1) high-redshift in situ mass growth resulting in a dense stellar component in the center of the potential where the dark matter fraction is low and (2) dry merger events dominating the mass growth at lower redshift (with major mergers being rare) resulting in the buildup of a stellar envelope increasing the half-light radius and thus the dark matter fraction (similar to the interpretation above).

The 24 μm SFR limit, used to restrict the stellar population models, results in sSFRs for our galaxies of log₁₀(sSFR/yr) < −10. Nonetheless, stronger limits on the sSFR can be obtained if the source of dust heating is not caused by recent SF. In Section 4.4.3, the information from optical nebular emission and MIR is combined to set stringent limits on the SFR of our sample (see also Figure 4). This information reveals that our sample lies 1.5 dex below the z = 2, SFR-stellar mass relation mass relation of Speagle et al. (2014; extrapolated to ${\mathrm{log}}_{10}({M}_{\odot }/{M}_{* })\sim 11.5$).

Low-luminosity AGNs are shown to be common in massive, log₁₀(M_*/M_⊙) > 11, quiescent galaxies at z < 1.5 through excess radio emission in stacked samples (Man et al. 2016a; Gobat et al. 2018). Six galaxies in our sample have direct radio detections, three of them with matching MIR detections (see Table 2). This could be evidence in line with the results from Olsen et al. (2013), who found a high fraction of AGNs in MQGs at 1.5 < z < 2.5 using X-ray stacking. Low-luminosity AGN activity has, in Schawinski et al. (2009) and Best & Heckman (2012), been associated with the suppression of SF, which is an important effect in maintaining quiescent galaxies. Low levels of dust heating have also been associated with evolved stellar populations as a significant source to emit at wavelengths beyond >160 μm (Salim et al. 2009; Bendo et al. 2012; Fumagalli et al. 2014; Utomo et al. 2014). However, with no detections in the Herschel/PACS bands, we cannot rule this scenario out. In the case where AGNs are indeed the dominant dust heating source in the galaxies, we can expect that the 24 μm flux does not arise from residual SF. This is consistent with Whitaker et al. (2017), who found no strongly obscured SF in MQGs at z > 2. Assuming the 24 μm emission is not due to obscured SF, we find an sSFR, log₁₀(sSFR/yr) < −11, based purely on the optical emission limits/detections. The MIR-to-radio emission of the sample will, in a future publication, be investigated in detail (I. Cortzen et al. 2020, in preparation).

The main limitations of the results are presented as follows.

• 1.  
  Overestimated stellar masses would lead to a shallower mass–size evolution and dynamical-to-stellar mass ratio evolution. Nonetheless, substantially overestimated stellar masses are ruled out by our dynamical masses being in agreement with previous kinematic studies of MQGs at z > 2 (Toft et al. 2012; Bezanson et al. 2013; van de Sande et al. 2013; Belli et al. 2014b, 2017).
• 2.  
  If rotation is significant in MQGs at z > 2, the measured velocity dispersion, depending on the inclination, could have an unknown contribution from rotation resulting in heightened dispersion measurements. On the other hand, dispersion measurements from face-on rotation-dominated galaxies could result in low values. This would further drive the dynamical mass artificially down. Such issues should be addressed by spatially resolved spectroscopy where the V_rot/σ can be estimated.
• 3.  
  Previous studies (Mancini et al. 2010) have suggested that sizes might be underestimated due to nondetection of low-luminosity profile wings. However, ultradeep imaging out to many effective radii does not find that this is the case (Szomoru et al. 2010, 2011).
• 4.  
  Dynamical-to-stellar mass evolution is sensitive to the determination of β(n). The prescription from Cappellari et al. (2006) is used, yet this relation is determined from local galaxies and assumed to be representative of dynamical systems at z ∼ 2. When comparing with the MASSIVE(n) sample, we assume a Sérsic index of n = 4 to be a fair representation of a spheroidal system. When changing the choice of β = 2–6 for the MASSIVE(n) sample, the conclusion that the ratio must evolve from z = 2 to 0 remains.
• 5.  
  The sample is 60% mass complete for the massive (log₁₀(M_*/M_⊙) > 11) and K-band-brightest (K < 20.5) UVJ quiescent galaxies at 1.9 < z < 2.5. This selection depends strongly on the performance of the photometric redshift estimate. In Figure 3, we show that this works well for our sample using the catalog from Muzzin et al. (2013a). This suggests that the sample studied in this paper is representative of the selection we presented in Section 2. However, the photometry is used to select red systems and, consequently, introduce a selection bias toward mergers between red galaxies. An unresolved merger of a quiescent galaxy with a star-forming galaxy would produce a resulting bluer system that might be excluded from the selection (see Section 6.1).

The H_F160W images from our program and the ancillary COSMOS F814W images employed in this work do not share the same world coordinate system. We need to guarantee that the astrometry is common and accurate in both bands. Therefore, we choose to align the images to the COSMOS ACS F814W image as the reference frame, which is registered to the fundamental astrometric frame of the COSMOS field, ensuring an absolute astrometric accuracy of 005–01 or better. Following Gómez-Guijarro et al. (2018), we use TweakReg along with SExtractor (Bertin & Arnouts 1996) catalogs of the two bands with the F814W catalog and frame as references to register the images. After this, the images in both bands are resampled to a common grid and a pixel scale of 006 pixel⁻¹ using SWarp (Bertin et al. 2002). In addition, the spatial resolution of the two HST bands is also different. Following Gómez-Guijarro et al. (2018), we degrade the F814W to the resolution of the F160W data (018 FWHM). We calculate the kernel to match the ACS F814W to the PSF in the F160W images employing the task PSFMATCH in IRAF, including a cosine bell function tapered in frequency space to avoid introducing artifacts in the resulting kernel from the highest frequencies. Then, we convolve this kernel to the F814W image to achieve a common spatial resolution.

Based on the spatially offset emission in the 2D X-shooter spectra, we determine whether candidate sources are within close proximity to the central galaxy. In Figure 5, we show the galaxies together with their spatially offset companions. The object UV-171687 shows offset Hα and [N ii] emission arising from a southwestern source that we establish to be a foreground galaxy at z = 1.51. We find another foreground galaxy northeast of UV-171060 at z = 1.37 based on assuming that the single emission line detection is Hα. Northeast of UV-155853, we find a background galaxy at z = 2.36 (best visible in the GALFIT modeling residuals of Figure 5) determined from the [O iii] doublet at 4959,5007 Å. For UV-105842, we find two spatially offset sources (1) ∼3'' northeast and (2) ∼1'' northeast. Source (1) is a foreground galaxy at z = 0.44 based on the detected strong O[iii] doublet at 4959,5007 Å and Hα emission. For source (2), we find the [O ii] doublet at 3726.2,3728.9 Å, O[iii] doublet at 4959,5007 Å, and Hα corresponding to a redshift z = 2.0124. The latter redshift corresponds to a velocity offset of 2130 ± 120 km s⁻¹ (uncertainty is calculated based on the spread of the individual redshift measurements), suggesting that it is not gravitationally bound to the central galaxy. Another option could be an offset AGN with high peculiar velocity following a merger.

To estimate the statistical error, we measure the spread of the velocity dispersion distribution obtained from running pPXF on 1000 data realizations. The data realizations are made by perturbing the pPXF best-fit model with the pipeline-estimated error spectrum by linearly drawing values from a Gaussian with a mean of zero and spread of the initial errors. The X-shooter pipeline-estimated noise map is subjected to a wavelength-dependent correlation of the pixels. We take this effect into account by scaling our noise spectrum to a reduced χ_red² = 1 (assuming the errors are Gaussian). We follow the method used in Toft et al. (2017) and fit a second-order polynomial to a 50 pixel running reduced χ² that we use to make a correction noise map, ${\sigma }_{{\chi }_{\mathrm{corr}}^{2}}={\sigma }_{{\chi }_{\mathrm{original}}^{2}}\sqrt{{\chi }_{\mathrm{fit}}^{2}}$.

We estimate the systematic error by testing how the dispersion is changing with the correction polynomial and implemented wavelength range. We construct a grid of correction polynomials up to 24th order of both additive and multiplicative polynomials, where we find an average of 20% variation from the fiducial dispersion, except for UV-239220, UV-773654, UV-171060, and CP-1291751. When varying the starting wavelength range ([λ_start, λ_end]) within the interval [3750–4050, 5950] and the ending wavelength within the interval [3750, 4050–5950], we find that overall, the dispersions are stable. In a few cases, the velocity dispersion increases well above the median dispersion (with varying wavelengths) with 50%–100% when excluding the higher-order Balmer and Ca H+K lines, highlighting their importance. When including the ending wavelength λ > 4500, we find more stable dispersion measurements; this is not surprising, as otherwise only half of the spectrum is included. The low-S/N cases have more unstable dispersion values when excluding wavelength areas, highlighting the importance of understanding the systematic uncertainties. We sum up the wavelength and polynomial test by confirming that our fiducial velocity dispersions are robust (except for UV-239220, UV-773654, UV-171060, and CP-1291751). The systematic error is primarily due to template mismatch, and as a result, we estimate the systematic error from the minimum and maximum values of the dispersion when using the full wavelength range and varying the additive and multiplicative correction polynomials, ${\sigma }_{\mathrm{sys}}=2/3\cdot ({\sigma }_{\max }-{\sigma }_{\min })/2$. This method is subjected to catastrophic outliers, and prior to the systematic error estimate, we exclude dispersion values more than 5σ outside of a Gaussian mean. We find that the systematic errors are on the order of the statistical uncertainties.

B.1.1. Additional Tests