The Fundamental Plane of Massive Quiescent Galaxies at z ∼ 2

In S20, we presented a sample of MQGs at z > 2 studied with the X-shooter spectrograph (D'Odorico et al. 2006; Vernet et al. 2011) at the Very Large Telescope and Hubble Space Telescope (HST) H_F160W that are crucial to obtain both rest-frame optical stellar velocity dispersions and effective sizes. The sample is selected from the 2 square degree COSMOS field, using multi-wave band photometric fits (Muzzin et al. 2013). In summary, the sample is selected to be K-band bright and massive ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11$) UVJ quiescent galaxies at z > 2. The adopted stellar population parameters such as mass-weighted age and stellar mass were based on the COSMOS15 photometry (Laigle et al. 2016), the continuum emission modeling of the X-shooter spectra, and the choice of star formation history (see details in S20, Section 4.3). Here we consider eight of the total 15 galaxies from S20 with measured stellar velocity dispersions, which is essential to study them in the scaling relations. We find no selection bias when comparing the size, age, stellar mass, and redshift of this sample to the parent sample in S20.

In Table 1, we list the mass-weighted age, stellar mass, stellar velocity dispersions, and sizes for these eight galaxies. The adopted stellar velocity dispersion will be referred to as the effective stellar velocity dispersion in this study due to minimal correction (<5%) when following the equations based on X-shooter observations in van de Sande et al. (2013). For further details, we refer to S20. Contrary to S20, we compute the dynamical masses in this paper using the circularized sizes (${r}_{e,{circ}}=\sqrt{{ba}}$) ¹⁴ to make them consistent with the dynamical masses derived using the circularized radius from the study of local cluster scaling relations in previous work such as, e.g., Jørgensen & Chiboucas (2013). We have verified that the qualitative results from S20 remain when using circularized dynamical masses and sizes. Following local studies (e.g., Jørgensen 1999), we adopt the Bessel B-band luminosity, estimated in our case from the rest-frame fluxes obtained from the most recent COSMOS photometry (Laigle et al. 2016) using the photometric redshift code EAZY (Brammer et al. 2008). ¹⁵ Normally, the magnitudes (luminosities) from the effective radius profile fitting are used when studying the scaling relations. We adopted another approach to, instead, extract the COSMOS magnitudes to treat the high-redshift samples (S20 and the galaxies from Section 2.2) consistently. We confirm that these magnitudes, when compared to the H_F160W magnitudes from the profile fits, can be considered representative with negligible differences. As a result, the method here is similar to the standard approach. The luminosity and average effective surface brightness are estimated using the method outlined in Appendix A. Hereafter, the dynamical mass-to-light ratio in the Bessel B-band is referred to as M/L. The eight galaxies in our sample have a mean age of ∼1.5 Gyr and span a similar stellar mass and size range to the full 15 galaxy parent sample from S20.

Table 1. Sample Summary

Target ID	R.A.	Decl.	zspec	log age yr–1	$\mathrm{log}{M}_{* }/{M}_{\odot }$	$\mathrm{log}{M}_{\mathrm{dyn}}/{M}_{\odot }$	re,circ	σe	$\mathrm{log}\langle I{\rangle }_{e,{\rm{B}}}$
UV-105842	150.26265	2.0177791	2.0195	${9.19}_{-0.33}^{+0.26}$	${11.68}_{-0.17}^{+0.16}$	11.47 ± 0.19	2.91 ± 0.29	263 ± 57	4.40
UV-171687	149.88702	2.3506956	2.1020	${9.13}_{-0.32}^{+0.28}$	${11.51}_{-0.19}^{+0.18}$	11.31 ± 0.24	4.49 ± 0.45	182 ± 50	3.99
UV-90676	150.48750	2.2700379	2.4781	${9.09}_{-0.29}^{+0.29}$	${11.78}_{-0.18}^{+0.17}$	11.78 ± 0.21	4.08 ± 0.41	347 ± 82	4.46
UV-155853	149.55630	2.1672480	1.9816	${9.23}_{-0.33}^{+0.24}$	${11.62}_{-0.17}^{+0.18}$	11.57 ± 0.11	4.20 ± 0.42	247 ± 30	3.96
UV-230929	150.20842	2.7721019	2.1679	${9.10}_{-0.28}^{+0.28}$	${11.48}_{-0.16}^{+0.16}$	11.16 ± 0.07	1.48 ± 0.15	252 ± 21	5.11
CP-1243752	150.07394	2.2979755	2.0903	${9.23}_{-0.32}^{+0.24}$	${11.79}_{-0.17}^{+0.17}$	11.61 ± 0.13	2.54 ± 0.25	350 ± 53	4.59
CP-540713	150.32512	1.8185385	2.0409	${9.16}_{-0.32}^{+0.27}$	${11.26}_{-0.23}^{+0.22}$	11.53 ± 0.24	1.46 ± 0.15	353 ± 97	4.66
UDS-19627 a	34.57125	−5.3607778	2.0389	${9.08}_{-0.10}^{+0.11}$	${11.37}_{-0.10}^{+0.13}$	11.33 ± 0.14	1.43 ± 0.14	318 ± 53	5.20

Notes. Column 1: target ID from S20; Column 2: R.A. in degrees (J2000); Column 3: decl. in degrees (J2000); Column 4: spectroscopic redshift; Column 5: mass-weighted age; Column 6: stellar mass; Column 7: circularized dynamical mass calculated using β(n) from Cappellari et al. (2006); Column 8: circularized effective radius in kpc; Column 9: stellar velocity dispersion in km s⁻¹;  Column 10: average surface brightness within r_e,circ in L_B,⊙pc⁻² (see Appendix A).

^a The spectroscopic redshift, age, stellar mass, and stellar velocity dispersion are from Toft et al. (2012) and the HST/WFC3 H_F160W sizes are from S20.

Download table as:  ASCIITypeset image

In addition to the galaxies from S20, we adopt a sample of 1.5 < z < 2.5 quiescent galaxies from Belli et al. (2017). We choose 11 out of 24 galaxies with stellar velocity dispersions, which have available COSMOS photometry in Laigle et al. (2016) to ensure consistent photometry extractions similar to the S20 sample. Out of the 11 galaxies, seven are at z < 2 and four at z > 2. The 11 galaxies from Belli et al. (2017) introduced here will be referred to as B17.

From Belli et al. (2017), we adopt the effective semimajor axis, absorption line-measured stellar velocity dispersions, and dynamical masses (see their Table 2). The effective semimajor axis are derived, similar to S20, using Sérsic profile fits to the rest-frame optical HST/WFC3 H_F160W images using GALFIT (Peng et al. 2002) and converted to circularized radius. We know from the S20 sample that the COSMOS magnitudes in the H-band are a good approximation of the magnitudes in the above profile fit method and thus indirectly related to the effective radii. The stellar velocity dispersions are derived from rest-frame optical Multi-Object Spectrometer For Infra-Red Exploration (McLean et al. 2012) spectra using Penalized PiXel-Fitting (Cappellari & Emsellem 2004) to fit Bruzual & Charlot (2003) (hereafter BC03) stellar population models which are qualitative similar to S20. The dynamical masses are converted from the semimajor axis to the circularized radius using the axis ratio. As in S20, the dynamical mass is estimated using the method from Cappellari et al. (2006).

Table 2. Fundamental Plane and M/L Scaling Relations

No.	Sample	Relation	Ngal	rms	σγ
(1)	Coma a	$\mathrm{log}{r}_{e}=(1.30\pm 0.08)\,\mathrm{log}\sigma -(0.82\pm 0.03)\,\mathrm{log}\langle I{\rangle }_{e,{\rm{B}}}-0.443$	123	0.08	0.01
(2)	Combined sample (1.5 < z < 2.5) b	$\mathrm{log}{r}_{e}=(0.46\pm 0.18)\,\mathrm{log}\sigma -(0.46\pm 0.07)\,\mathrm{log}\langle I{\rangle }_{e,{\rm{B}}}+1.275$	19	0.15	0.03
(3)	Coma a	$\mathrm{log}M/L=(0.27\pm 0.03)\,\mathrm{log}{M}_{\mathrm{dyn}}-2.145$	123	0.09	0.01
(4)	Combined sample (1.5 < z < 2.5) b	$\mathrm{log}M/L=(0.51\pm 0.15)\,\mathrm{log}{M}_{\mathrm{dyn}}-6.393$	19	0.26	0.06
(5)	Coma a	$\mathrm{log}M/L=(0.96\pm 0.10)\,\mathrm{log}\sigma -1.326$	123	0.10	0.01
(6)	Combined sample (1.5 < z < 2.5) b	$\mathrm{log}M/L=(1.09\pm 0.41)\,\mathrm{log}\sigma -3.236$	19	0.28	0.06

Notes. Column 1: sample; Column 2: fitting method; Column 3: scaling relations; Column 4: number of galaxies included in fit; Column 5: rms scatter along the y direction of the scaling relation; Column 6: uncertainty on the zero-point (γ) estimated as rms/$\sqrt{{N}_{\mathrm{gal}}}$ (assuming fixed coefficients).

^a Fits are from Jørgensen et al. (2019) Table 4. ^b Samples from S20 and B17 (see Section 2)

Download table as:  ASCIITypeset image

The two z ∼ 2 samples (Section 2.1 and 2.2) used here satisfy the UVJ quiescent galaxy selection from Muzzin et al. (2013). Together they allow for exploring a larger dynamical mass range, $10.5\lt \mathrm{log}({M}_{\mathrm{dyn}}/{M}_{\odot })\lt 11.9$, and redshift, 1.5 < z < 2.5, corresponding to a cosmological time span of ∼1.7 Gyr, considered short compared to their consequent 8-10 Gyr of evolution to the present day.

The original volume-limited MASSIVE Survey sample is selected as the most massive and K-band brightest early-type galaxies within 108 Mpc of the northern hemisphere (Ma et al. 2014). Here, we use the 25 most massive, $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.7$, MASSIVE galaxies (hereafter MASSIVE(n) sample), selected at fixed cumulative number density (CND) matching our massive z ∼ 2 sample as described in S20. The CND of the z ∼ 2 sample is estimated from the massive, $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.2$, UVJ quiescent galaxies at 1.9 < z < 2.5 in the Muzzin et al. (2013) catalog. The results in S20 are shown to be robust against the choice of CND method (fixed and probabilistic; Wellons & Torrey 2017), as well as the mass–rank scatter. Utilizing this CND approach is an attempt to minimize progenitor bias by predicting the local progenitors of the high-z galaxies. A thorough discussion of the assumptions and uncertainties is covered in S20 Section 5.1.

The magnitudes for 17/25 galaxies are obtained from the SDSS DR14 catalog photometry (Blanton et al. 2017) by cross-matching the MASSIVE(n) sample using the SDSS SkyServer. ¹⁶ The de Vaucouleur photometry ("deVMag") in the u, g, r, i, z bands is extracted and converted to rest-frame Bessel B-band magnitude using the EAZY code, in the same way as in the high-redshift samples in this study. The "cmodelMag" estimate, also used for the Coma sample (Section 2.4), is compared with the de Vaucouleur fit magnitudes to establish that the latter is a good representation, which is also confirmed by the "fracDeV" parameter. The luminosity and effective surface brightness are calculated from the apparent magnitude using the methods covered in Appendix A. Instead of estimating the luminosity distance from the redshift, these galaxies are close enough that peculiar velocities have a significant impact on their distance measurement. We therefore use the distance measurement from Ma et al. (2014), who correct for this effect. The 17 galaxy subsample used here is referred to as MASSIVE(n₁₇).

The de Vaucouleur effective radii ("deVRad") corresponding to extracted photometry from the SDSS SkyServer are used. After confirming that the radii are consistent among the different bands, we choose to use the g-band effective radii as the wavelength coverage is comparable to the Bessel B-band. The radii are circularized using the axis ratio ("deVAB"). The average luminosity weighted dispersion within the effective radius is adopted (Veale et al. 2018). The dynamical masses are estimated using the method in S20 (with the prescription from Cappellari et al. 2006) using n = 4 and the circularized effective sizes.

The stellar mass, size, and stellar velocity dispersion between the MASSIVE(n) and MASSIVE(n₁₇) samples are compared in Appendix B. Here we find that the MASSIVE(n₁₇) sample is uniformly sampled from the initial distribution of stellar mass, size, and stellar velocity dispersion and is 68% complete. As a result, our MASSIVE(n₁₇) selection with available photometry is representative of the parent CND matched sample and can be considered a suitable minimal progenitor-biased reference sample.

As a local reference cluster we use Coma/A1656 complete to ${g}_{{rf}}^{{\prime} }\leqslant 16.1$ magnitudes, as used in Jørgensen et al. (2019). The scaling relations of the Coma cluster are well studied and we include them as a reference to other high-redshift comparisons. We adopt the 123 galaxies with average effective surface brightness, stellar velocity dispersion, and circularized sizes. The stellar velocity dispersion measurements (from Jørgensen et al. 2018) are derived from high signal-to-noise ratio (∼60 Å⁻¹) spectra using the recipe presented in Jørgensen et al. (2017). The surface brightness is calculated using rest-frame Bessel B magnitudes based on the SDSS cmodelmag magnitude (see Appendix A in Jørgensen et al. 2019). The circularized effective radii are constructed from SDSS parameters to be a pseudo-Sérsic effective radii with the purpose of matching the "cmodelMag" magnitudes (see the method in Jørgensen et al. 2019). Both data sets are calibrated to the Legacy data (Jørgensen 1999; Jørgensen & Chiboucas 2013) to provide a trustworthy low-redshift reference cluster. For further details on the physical properties of the Coma cluster, we refer to Section 2 in Jørgensen et al. (2018).

M_dyn/L_B reveals information about how the dynamical structure of galaxies compares to the luminosity of their stellar populations. This ratio has been found to increase in massive quiescent galaxies, believed to be driven by the non-star-forming and passive evolution of their stellar population. Here, we study the M_dyn/L_B relations for MQGs and their evolution in the past 10 billion years. We show these relations in Figure 1, both as a function of dynamical mass (1(a)) and stellar velocity dispersion (1(b)), for a sample of z ∼ 2 MQGs (S20 + B17) alongside two local reference samples (Coma, MASSIVE(n₁₇)).

Zoom In Zoom Out Reset image size

Compared to the local Coma and MASSIVE(n₁₇) galaxies, the high-redshift samples (S20 and B17) have lower M/L as expected for brighter younger systems. The majority of the galaxies from S20 have, already at z = 2, dynamical masses similar to the  ∼10% most massive galaxies in the Coma cluster (see Figure 1(a)). The MASSIVE(n₁₇) galaxies have dynamical masses similar to the 2% most dynamical massive Coma galaxies, which shows that they are among the most massive galaxies in the local universe. The distribution of the stellar velocity dispersions of S20 and B17 is consistent with the high-end measurements in Coma (see Figure 1(b)). When comparing the stellar velocity dispersions between the S20 and the MASSIVE(n₁₇) sample we find that a shallow evolution is expected between z = 2 and the present day. The MASSIVE(n₁₇) sample have high stellar velocity dispersions similar to the high end of the Coma cluster measurements.

The combined samples of S20 and B17, in Figure 1(a) and (b), are fit by minimizing the least-squares in the y-direction and the uncertainty on the slopes is estimated using a bootstrap procedure (see also Jørgensen et al. 1996; Jørgensen & Chiboucas 2013). A relation for the combined high-redshift (S20 and B17) sample is established, while this was not possible using the narrow dynamical mass range of the S20 sample alone. The fits are shown in Figure 1, together with the associated rms from the regression, and listed in Table 2. M/L with dynamical mass relation appears slightly steeper, although by less than 2σ. For both M/L versus M_dyn and M/L versus σ_e , we find best-fit slopes for the combined high-redshift sample to be consistent with the local Coma relation (see relations 4 and 6 in Table 2).

The FP is spanned by the effective size, r_e , stellar velocity dispersion, σ_e , and average effective surface brightness, 〈I〉_e,B. Its edge-on and face-on orientations are defined as

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\,{r}_{e}=\alpha \,{\mathrm{log}}_{10}\,{\sigma }_{e}+\beta \,{\mathrm{log}}_{10}\,\langle I\rangle {}_{e,{\rm{B}}}+\gamma \end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}(2.22\,{\mathrm{log}}_{10}\,{r}_{e}+\beta \,{\mathrm{log}}_{10}\,\langle I\rangle {}_{e,{\rm{B}}}+\alpha \,{\mathrm{log}}_{10}\,{\sigma }_{e})/2.7\\ \quad =\,(\alpha \,{\mathrm{log}}_{10}\,\langle I\rangle {}_{e,{\rm{B}}}-\beta \,{\mathrm{log}}_{10}\,{\sigma }_{e})/1.54,\end{array}\end{eqnarray} \tag{ 2 }$$

respectively. The best-fit Coma relation slopes (α = 1.30 ± 0.08, β = −0.82 ± 0.03) and zero-point (γ = −0.443) in the rest-frame B-band are adopted as our local reference orientation of the plane (Jørgensen et al. 2006; see also Table 2).

The best-fit Coma data, from Jørgensen et al. (2019), were fit using the same method as described below while adopting the FP parameters from Jørgensen & Chiboucas (2013) (α = 1.30 ± 0.08, β = −0.82 ± 0.03) to get the zero-point (γ = −0.443). We use this as our local reference cluster fit and list the parameters along with their uncertainties in Table 2.

In Figure 2, the FP edge-on and face-on projections, as described in Equations (1) and (2), are shown. We use this to examine how z ∼ 2 MQGs populate and evolve to z = 0 in this plane. In the edge-on FP, the dominating errors from the stellar velocity dispersion are shown on the y-axis. For the face-on plane, the errors are calculated similarly to the approximation used in Jørgensen et al. (2006).

Zoom In Zoom Out Reset image size

The COSMOS quiescent galaxies from S20 and B17 are found to be below (in the edge-on plane) and above (in the face-on plane) the local Coma FP relation. These galaxies have compact sizes and younger stellar populations (due to their high redshift and more recent quenching), effectively increasing their mean effective surface brightness.

An edge-on FP cannot be clearly established using the S20 sample alone. However, when fitting the S20 and B17 samples together, an FP is in place at 1.5 < z < 2.5 (hereafter referred to as FP_z∼2). The FP is fitted using the least-squares method, minimizing the least-squares in the y-direction, with uncertainties from a bootstrapping method (see relation 2 in Table 2). In van de Sande et al. (2014), an FP was indicated for a similar epoch (1.5 < z < 2.5). However, in this study, the sample of z > 2 galaxies is three times  more numerous, robustly confirming the existence of a plane at 1.5 < z < 2.5.

An FP_z∼2 is established using 19 massive quiescent galaxies that as a result must have been a relatively homogeneous population already at this epoch, signaling even earlier formation and significant evolution from z = 2 to the present day.

The expected passive evolution of the stellar population, from the redshift of formation to the present day, is estimated based on the evolution in M_*/L of a BC03 simple stellar population model with Chabrier IMF and solar metallicity. The formation redshift is estimated from the mass-weighted age and the redshift of observation. The passive evolution takes into account mass loss during stellar evolution, and the M_* is thus the mass locked into stars at a given age. The mass-weighted mean stellar age is indicative of the epoch by which the bulk of the stellar mass was formed. The M_*/L uncertainties are dominated by the 1σ standard deviation of the age (∼0.2 dex) when assuming solar metallicity. The definition of passive evolution assumes that no new stellar mass is added and thus M_*/L directly relates to the ${\rm{\Delta }}{\mathrm{log}}_{10}\,{L}_{\mathrm{passive}}$ from Equation (8).

The average estimate of the passive evolution is shown by the red arrow in Figure 3. This is what is expected as Equation (8) reduces to ${\rm{\Delta }}{\mathrm{log}}_{10}\,M/L={\rm{\Delta }}\gamma /\beta =-{\rm{\Delta }}{\mathrm{log}}_{10}{L}_{\mathrm{passive}}$ for no structural (${\rm{\Delta }}{\mathrm{log}}_{10}{r}_{e}={\rm{\Delta }}{\mathrm{log}}_{10}{\sigma }_{e}=0$) and no luminosity change (from a stellar mass increase). The passive evolution cannot fully explain the direct evolution, of MQGs from z = 2 to 0, to Coma galaxies at fixed dynamical mass (or size). Furthermore, it is inadequate in alone explaining the expected evolution to the MASSIVE(n₁₇) sample.

We learned from S20 that the size growth with shallow stellar velocity dispersion change was driven by minor merger structural evolution. Using this information we can simplify Equation (8):

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{log}}_{10}M/L=\displaystyle \frac{{\rm{\Delta }}\gamma }{\beta }=\displaystyle \frac{1+2\beta }{\beta }\,{\rm{\Delta }}{\mathrm{log}}_{10}{r}_{e}.\end{eqnarray} \tag{ 9 }$$

The effect of the structural evolution adopted from S20 is based on the size increase from the stellar mass–size plane, their Figure 8 (${\rm{\Delta }}\mathrm{log}{r}_{e,{circ}}\sim 0.6$), alongside insignificant stellar velocity dispersion evolution (${\rm{\Delta }}\mathrm{log}{\sigma }_{e}\sim 0$). This effect is shown in the edge-on FP in Figure 3 as the sloped blue arrow. The combined effect of the mean passive evolution of the stellar population (red arrow) and mean the structural evolution from minor merger-driven size increase (blue arrow) brings the S20 MQGs in the vicinity of the MASSIVE(n₁₇) sample.

S20 found that the median stellar-to-dynamical mass ratio ${\rm{x}}={\rm{\Delta }}{\mathrm{log}}_{10}{M}_{* }/{\rm{\Delta }}{\mathrm{log}}_{10}{M}_{\mathrm{dyn}}\sim 0.5$. Specifically the stellar mass doubles (${\rm{\Delta }}\mathrm{log}{M}_{* }\sim 0.3$ dex) from z = 2 to 0. Following Equation (8), we show how this ratio affects the evolution of the scaling relations when not unity:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{log}}_{10}\,M/L=\displaystyle \frac{{\rm{\Delta }}\gamma }{\beta }={\rm{x}}\,{\rm{\Delta }}{\mathrm{log}}_{10}{M}_{\mathrm{dyn}}.\end{eqnarray} \tag{ 10 }$$

Here the effect of a changing stellar-to-dynamical mass ratio can be seen to directly affect the change in the scaling relations. The effect for x = 0.5 is shown in Figure 3. When taking into account the stellar-to-dynamical mass ratio, which affects the conversion between the stellar and dynamical M/L, we obtain a predicted evolution closer to the high-mass end of the Coma relation and the median position of the MASSIVE(n₁₇) sample.

From S20 we learned that z = 2 MQGs grow their size and stellar mass through minor mergers in their evolution to z = 0. In addition to the decrease in the B-band luminosity, due to passive evolution, it is expected that such merger events could add to the B-band luminosity if these are star-forming galaxies at the time of merging. From now on, this type of merger is referred to as "wet," contrary to the "dry" minor mergers that are passive before merging. In Equation (8) this reduces to the effect of the stellar M/L change from new stars accreted from minor mergers:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{log}}_{10}\,M/L={\rm{\Delta }}{\mathrm{log}}_{10}{\left(\displaystyle \frac{{M}_{* }}{L}\right)}_{{M}_{* }}.\end{eqnarray} \tag{ 11 }$$

The B-band luminosity increase from merging galaxies between z = 2−0 is modeled with composite stellar population models from the BC03 library with solar metallicity to obtain the stellar M/L (similar to Section 4.1). The star formation history follows the evolution of the main sequence (Speagle et al. 2014). It is assumed that, after merging, the galaxies stop forming stars and follow a passive evolution.

The median stellar mass increase (${\rm{\Delta }}\mathrm{log}{M}_{* }=0.3$ dex) from minor mergers, predicted in S20, is used assuming a 1:20 merger ratio. Note that in our simplistic model, the correct mass ratio does not play a significant role. We also investigate a more realistic scenario with mergers distributed across redshift (z = 1.8–0.1) together with two extreme cases of all the mass added at z = 1.8 or 0.1. A B-band luminosity increase of 0.4–0.45 dex is found in all cases (for more details see Appendix D).

In Figure 3 the effect of increasing luminosity, due to wet minor merger stellar populations, is shown by the green arrow (0.4 dex). This effect appears to be in disagreement with the location of the Coma relation and the MASSIVE(n₁₇) sample. This suggests that dry minor merger galaxies, with no additional luminosity increase, are a preferred evolution scenario for MQGs.

Studies of passive galaxies in 0.8 < z < 1.8 clusters (among others Jørgensen et al. 2006, 2019; van der Marel & van Dokkum 2007; Jørgensen & Chiboucas 2013; Beifiori et al. 2017 ) find that the change in M_*/L can be explained by passive evolution to z = 0. Below we explore if a similar analysis can account for the evolution of the scaling relations at z ∼ 2.

Jørgensen & Chiboucas (2013) derived the M_dyn /L_B evolution (based on models from Maraston 2005), as a function of stellar age and metallicity, to be $\mathrm{log}{M}_{{dyn}}/{L}_{{\rm{B}}}\,=0.935\ \mathrm{log}\,\mathrm{age}+0.337[{\rm{M}}/{\rm{H}}]-0.053$. Assuming passive evolution from z = 2 to the best-fit Coma relation, at $\mathrm{log}({M}_{\mathrm{dyn}}/{M}_{\odot })=11.5$, we find a formation redshift of ${z}_{\mathrm{form},\mathrm{Coma}}={2.01}_{-0.04}^{+0.1}$ (for details see Appendix E). The formation redshift is similar to the redshift of observation, which leaves too short a time to form the S20 MQGs at this epoch. The formation redshift, derived from our stellar population mass-weighted ages (assuming the median age), is ${z}_{\mathrm{form}}={3.41}_{-0.91}^{+4.92}$. The uncertainties are estimated using the 1σ age uncertainties. Based on this, we conclude that the S20 MQGs at z > 2 cannot evolve to the Coma relation by passive evolution alone.

The fixed CND-matched MASSIVE(n₁₇) sample allows us to study the evolution of massive galaxy scaling relations from z = 2 to 0 with minimal progenitor bias. Evidence against purely passive evolution to z = 0 is present in both the scaling relations (Figures 3 and C1) as explicitly shown in the previous section. The S20 MQGs at z > 2 are consistent with evolving into the most massive Coma galaxies, and the MASSIVE(n₁₇) sample, through passive, structural, and stellar-to-dynamical mass ratio evolution.

The size increase of MQGs in cosmological simulations could be explained by adiabatic expansion due to AGNs, decreasing the central mass density and puffing up the galaxies (Dubois et al. 2013; Choi et al. 2018). Major mergers, as the dominant mechanism for size growth, have become less popular as they make the galaxies too massive to be consistent with massive nearby galaxies. In S20, the structural evolution is interpreted to be from minor mergers in line with the scenario presented in the idealized simulations from Hopkins et al. (2009), Naab et al. (2009), and Hilz et al. (2012, 2013). Here, the effective half-light radius grows by adding stars to the outskirts of the galaxy from tidally stripped minor mergers. This scenario is shown to cause inside-out growth, starting from a compact elliptical galaxy (core) that causes, through minor mergers, a build-up of the surface density profile wings, a present-day analog of a giant elliptical galaxy (core–envelope). A consequence of the inside-out minor merger growth scenario from Hilz et al. (2012) is an increasing dark matter fraction which has been suggested to cause a tilt in the FP over time (Boylan-Kolchin et al. 2005; D'Onofrio et al. 2013). Essentially this is a consequence of using the effective radius as a foundation for analyzing the evolution, as this results in the stellar velocity dispersion mainly tracing the stars at z = 2, but at z = 0 a higher fraction of dark matter to stars. This effectively introduces a systematic difference between comparing stellar velocity dispersions across epochs that only trace stars with a similar sample where the stellar-to-dynamical mass ratio evolves. This systematic difference, included in Section 4.3, must be taking into account and can be seen as the black arrow in Figures 3 and C1.

In Section 4.4, the predicted luminosity increase from wet minor mergers is modeled under the assumption that they are the primary drivers of the size growth. For the realistic scenario of adding wet minor mergers continuously from z = 2 to 0, we find that the luminosity increases by roughly 0.4 dex.

In Figures 3 and C1, the predicted position of the S20 MQGs (following passive, structural, and stellar-to-dynamical mass evolution) is indicated alongside the effect of the luminosity from wet minor mergers. The green arrow indicates how the predictions would move compared to the local best-fit relation of Coma and around the locus of the MASSIVE(n₁₇) sample, strongly favoring dry minor mergers.

In the inside-out growth scenario, the rest-frame B-band luminosity increase takes place in the outer parts of the galaxy. The luminosity from the MASSIVE(n₁₇) sample is measured using SDSS deVMag which represents the luminosity of the galaxy out to 8r_e. Thus, an underestimation of the luminosity, by only sampling the central part of the galaxy and missing the outskirts, is unlikely.

The wet minor merger luminosity increase offsets the predictions from the local relation and the MASSIVE(n₁₇) sample and thus appears to not be a favored way to grow MQGs at z ∼ 2. Another possibility is that the minor mergers already have quenched stellar populations (before their merger) with low rest-frame B-band luminosity (Oogi & Habe 2013; Naab et al. 2014; Tapia et al. 2014). The evolution from z = 2 to 0 of the FP and M/L scaling relations is consistent with such a scenario, caused primarily by dry minor mergers, and passive and stellar-to-dynamical mass ratio evolution for MQGs.

Data from S20 and B17 are combined to establish the FP_z∼2 at 1.5 < z < 2.5 with more than half of the sources at z > 2. A large fraction of the quiescent galaxies from B17 is found to be disk-like (based on Sérsic index, n < 2.5; Belli et al. 2017), which could mean that an unknown contribution from rotation is included in the measured stellar velocity dispersion. For spherical dispersion-dominated systems, the circularized radius and semimajor axis are comparable methods of size measurement. However, for more disk-like systems the difference grows between the two size measuring methods, further causing a bias between dispersion- and rotation-dominated galaxies. We estimate, based on the axis ratios, that the circularized sizes differ by 7%–30% compared to the semimajor axes. This is well within the quoted uncertainties of the predicted position of the S20 MQGs at z = 0. This bias could potentially affect the zero point and coefficients of the best-fit in Figures 1 and 2. This issue could be solved by spatially resolved spectroscopy disentangling the contribution from rotation and dispersion.

The SDSS modelMag luminosities of the most luminous galaxies have been underestimated (Bernardi et al. 2017). In our study, the deVMag and cmodelMag methods have been used to estimate the luminosity of the MASSIVE and Coma sample, respectively. The median offset for the brightest galaxies (M_r  ∼ −24), in the r-band, are  ∼15% (see Figure 7 in Bernardi et al. 2017). This translates to a difference in ${\rm{\Delta }}{\mathrm{log}}_{10}L\propto {\mathrm{log}}_{10}({10}^{0.4{\rm{\Delta }}{M}_{r}})=0.4\,\cdot \,0.15\,=\,0.06$, assuming that the r-band magnitudes are representative for the Bessel B-band used in our study. In the case that this assumption is valid, the M/L and average effective surface brightness would change by 0.06 dex, thus moving the local galaxies in the positive y-direction, by the same amount, in Figures 1, 2, and 3. Such an effect is minimal and would not affect the general trends, results, and conclusions made in this paper.

The dominating uncertainty of the mass–rank scatter from the CND-matching of the local MASSIVE(n₁₇) sample does not affect the conclusions of this study. Furthermore, if using a probabilistic CND-matching approach (see Wellons & Torrey 2017), this would increase the number of galaxies in the MASSIVE(n₁₇) sample from 17 to 30. In Figures 1 and 2 this would correspond to a greater number of white hexagons becoming purple, which causes no noticeable effects on the trends in the figures. On the other hand, if not all massive galaxies at z > 2 have similar merger histories, descendants that become quiescent systems at late times would have been missed (see, e.g., Naab et al. 2014). A study of the average stellar population age within the effective radius of the most massive (σ_* > 220 km s⁻¹) MASSIVE Survey galaxies (Greene et al. 2015) find their ages to be  >10 Gyr. This suggest that our sample of MASSIVE(n₁₇) galaxies already were quenched at z > 2, and likely have similar merger histories from z = 2 to 0. Studies also show that the most massive end of the stellar mass function ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.5$) evolves very little, if at all, between 0 < z < 2 (McLeod et al. 2020, their Figure 5).