Near-infrared spectroscopy of 5 ultra-massive galaxies at 1.7 < z < 2.7

TL;DR

This study uses near-infrared spectroscopy to analyze five ultra-massive galaxies at redshifts 1.7 to 2.7, confirming their high stellar masses and diverse star formation activities through spectroscopic and photometric data.

Contribution

First near-infrared spectroscopic analysis of ultra-massive galaxies at these redshifts, combining spectroscopy with photometry to refine stellar population properties.

Findings

Three galaxies are quiescent with low star formation rates.

One galaxy shows high star formation and dust obscuration.

Spectroscopic data confirms extreme stellar masses.

Abstract

We present the results of a pilot near-infrared (NIR) spectroscopic campaign of five very massive galaxies ($\log(\text{M}_\star/\text{M}_\odot)>11.45$) in the range of $1.7<z<2.7$. We measure an absorption feature redshift for one galaxy at $z_\text{spec}=2.000\pm0.006$. For the remaining galaxies, we combine the photometry with the continuum from the spectra to estimate continuum redshifts and stellar population properties. We define a continuum redshift ($z_{\rm cont}$ ) as one in which the redshift is estimated probabilistically using EAZY from the combination of catalog photometry and the observed spectrum. We derive the uncertainties on the stellar population synthesis properties using a Monte Carlo simulation and examine the correlations between the parameters with and without the use of the spectrum in the modeling of the spectral energy distributions (SEDs). The spectroscopic constraints confirm the extreme stellar masses of the galaxies in our sample. We find that three out of five galaxies are quiescent (star formation rate of $\lesssim 1 M_\odot~yr^{-1}$) with low levels of dust obscuration ($A_{\rm V} < 1$) , that one galaxy displays both high levels of star formation and dust obscuration (${\rm SFR} \approx 300 M_\odot~{\rm yr}^{-1}$, $A_{\rm V} \approx 1.7$~mag), and that the remaining galaxy has properties that are intermediate between the quiescent and star-forming populations.

Tables

Table 1. Table 1: Pointing and observation information for the initial sample.

id	RA	DEC	Dates	Exp. Time (sec)	Instrument
				( exposure time x exposure count )
COS-75355	10h02m28.49s	+02d02m13.70s	12/18/13	600 x 10 = 6000	GNIRSa
COS-207144	10h00m33.48s	+02d28m54.74s	03/07/13	600 x 18 = 10800	GNIRSa
COS-90676	10h01m57.00s	+02d16m12.14s	12/17/13	600 x 6 = 3600	GNIRSa
COS-189962c	10h02m14.42s	+02d35m11.92s	05/22/13, 05/23/13	600 x 18 = 10800	GNIRSa
COS-71929	10h01m40.60s	+01d58m57.47s	05/04/13	908.8 x 4 = 3635.2	FIREb
COS-37207	09h59m42.59s	+01d55m01.55s	05/03/13, 05/04/13	908.8 x 8 = 7270.4	FIREb

Table 2. Table 2: Photometry from the UltraVISTA catalog for the six galaxies in the initial sample. Total magnitudes are used. Magnitudes are given in the AB system.

		UltraVISTA DR1 ID
	75355	207144	90676	189962	71929	37207
filter
$m_{u^{\ast }}$	25.7	26.4	28.82	24.4	26.7	25.2
$m_{g^{+}}$	24.8	26.3	26.8	24.3	26.1	24.7
$m_{r^{+}}$	23.9	26.0	24.7	24.0	25.4	24.1
$m_{i^{+}}$	23.0	25.1	24.2	23.7	24.4	23.5
$m_{z^{+}}$	22.1	24.5	23.8	23.2	23.8	22.9
$m_Y$	21.5	23.3	23.1	22.9	22.9	22.5
$m_J$	20.3	22.1	22.1	21.7	21.7	21.3
$m_H$	19.8	21.0	20.6	21.0	20.6	20.5
$m_{K_s}$	19.4	20.5	20.2	20.5	20.1	20.0

Table 3. Table 3: FAST Best-fit Parameters

id	207144	37207	71929	75355	90676
${\text{A}}_{\text{V}}$ (mag)	${0.00}_{-0.00}^{+1.10}$	${1.70}_{-0.40}^{+0.10}$	${1.30}_{-0.40}^{+0.10}$	${0.30}_{-0.10}^{+0.10}$	${0.50}_{-0.50}^{+0.30}$
log (Age [yr])	${9.50}_{-0.60}^{+0.00}$	${9.00}_{-0.00}^{+0.40}$	${9.50}_{-0.10}^{+0.00}$	${9.00}_{-0.10}^{+0.10}$	${9.10}_{-0.40}^{+0.20}$
log (${\text{M}}_{\star }/{\text{M}}_{\odot }$)	${11.46}_{-0.26}^{+0.00}$	${11.67}_{-0.01}^{+0.12}$	${11.87}_{-0.12}^{+0.00}$	${11.53}_{-0.02}^{+0.05}$	${11.64}_{-0.16}^{+0.06}$
log (SFR [${\text{M}}_{\odot }{\mathrm{yr}}^{-1}$])	$-{0.18}_{-1.82}^{+0.28}$	${2.53}_{-0.46}^{+0.12}$	${0.94}_{-0.67}^{+0.02}$	$-{5.41}_{-5.80}^{+5.80}$	$-{0.54}_{-1.46}^{+1.12}$
log (sSFR [${\mathrm{yr}}^{-1}$])	$-{11.64}_{-87.36}^{+0.37}$	$-{9.14}_{-0.55}^{+0.10}$	$-{10.93}_{-0.60}^{+0.09}$	$-{16.93}_{0.00}^{+5.75}$	$-{12.18}_{-86.82}^{+1.10}$
log ($\tau$ [yr])	${8.50}_{-1.50}^{+-0.13}$	${8.50}_{0.00}^{+0.50}$	${8.60}_{-0.10}^{+0.00}$	${7.60}_{0.00}^{+0.40}$	${8.00}_{-1.00}^{+0.00}$
Z	${0.03}_{-0.02}^{+0.00}$	${0.01}_{-0.01}^{+0.00}$	${0.00}_{0.00}^{+0.00}$	${0.01}_{0.00}^{+0.02}$	${0.00}_{-0.00}^{+0.03}$
$z_{\text{peak}}$ or $z_{\text{spec}}$	${1.89}_{-0.06}^{+0.08}$	${2.11}_{-0.04}^{+0.01}$	${2.0000}_{-0.0055}^{+0.0055}$	${1.8220}_{-0.0006}^{+0.0006}$	${2.61}_{-0.08}^{+0.05}$
$q_{\text{sf}}$	$0.993\pm 0.005$	${0.276}_{-0.276}^{+0.529}$	$0.963\pm 0.028$	${1.0}_{-0.008}^{+0.0}$	${0.999}_{-0.008}^{+0.001}$

Equations

$$q_{\text{sf}}=1-\frac{{\text{SFR}}_{\text{current}}}{⟨{\text{SFR}}_{\text{past}}⟩}=1-\frac{{\text{SFR}}_{\text{current}}}{M_{\star }/\text{age}},$$

Full text

1 Introduction

The standard $\Lambda$CDM paradigm of structure formation depicts a universe dominated by dark matter. In this paradigm, dark matter structures form hierarchically, with low-mass haloes forming before their higher mass counterparts. It is therefore expected that galaxies would also form in a hierarchical manner. However, archaeological studies of local galaxies show that today’s most massive galaxies formed the bulk of their stars rapidly in the very early universe (see, for example, Thomas et al., 2005): these observations are supported by the existence of a significant population of massive quiescent galaxies up to redshift $z\sim 3$. (Cimatti et al., 2002; Kriek et al., 2006; Mancini et al., 2009; Marchesini et al., 2010; Brammer et al., 2011; Marchesini et al., 2014; Newman et al., 2015; Hill et al., 2016).

Recent models of galaxy formation have been successful in addressing several long-standing tensions between theoretical predictions of massive galaxy growth in the early universe and observational results of such massive galaxies at high-$z$. These models predict that massive ellipticals will possess old, metal-rich stellar populations and the shortest formation timescales (De Lucia et al., 2006; De Lucia & Blaizot, 2007), but there are still persistent differences between observations and theoretical predictions of the high-mass end of the stellar mass function. These issues raise questions about the evolution of massive galaxies: Marchesini et al. (2014) propose a revised evolutionary path for the formation of today’s ultra-massive galaxies (UMGs), in which the progenitor population at $1.5<z<2.5$ is dominated by heavily dust-obscured ($A_{V}\sim 2$ mag) massive star-forming galaxies, and only $\sim 40\%$ of the population is comprised of quiescent galaxies with relatively little dust obscuration.

Though recent photometric campaigns have been highly successful in probing the evolutionary path of UMGs, little is known about the evolution of the most massive galaxies (i.e., $\log{(M_{\star}/M_{\sun})}>11.5$) in the early universe. According to measurements of the stellar mass function, the number density of such galaxies appears constant from $z\sim 4$ to $z\sim 1.5$ (Marchesini et al., 2009; Muzzin et al., 2013b). In this stellar mass regime, model predictions are, at best, marginally consistent with these observations (Vogelsberger et al., 2014; Davé et al., 2016, but see also Henriques et al., 2015).

Due to the rarity of these galaxies, however, a huge amount of data and manpower are needed to generate a large enough sample to characterize these populations. Very few high-$z$ galaxies in this mass regime have been confirmed spectroscopically, and population studies rely almost entirely on photometry for redshifts and stellar population characteristics. Recent studies have furthered the spectroscopic effort for estimating the nature of highly massive galaxies at high redshift (Onodera et al., 2010; van de Sande et al., 2011; Onodera et al., 2012; Bezanson et al., 2013; van de Sande et al., 2013; Belli et al., 2014, 2015; Marsan et al., 2015, 2016), but the total sample of spectroscopically confirmed very massive, high-$z$ galaxies known to the community at the present is still not large enough to be able to leverage satisfactory statistical power in characterizing the underlying distribution that govern the characteristics of such galaxies.

In an ongoing effort to expand the sample of spectroscopically confirmed very massive galaxies at high-$z$, we present the results of near-infrared (NIR) spectroscopy of 5 highly massive galaxies at $1.7<z<2.7$. In Section 2, we outline the sample selection process. The spectroscopic observations are detailed in Section 3. We derive robust redshifts and stellar population properties using a combination of the NIR spectroscopy and UltraVISTA photometry in Section 4 in order to confirm their status as ultra-massive, and to study the characteristics of this rare class of galaxies in the early universe. We present spectroscopic redshifts and modeling results in Section 5. We place the characteristics derived for this sample of high-$z$ ultra-massive galaxies within the context of the population of such galaxies in Section 6. In this work, we assume a Chabrier (2003) IMF, and adopt a cosmology of $\Omega_{M}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}$ = 70 km s*-1* Mpc*-1*.

2 The Sample

The targets of this study were chosen from the sample of ultra-massive galaxies identified in the multiwavelength catalog constructed by Muzzin et al. (2013a) across 1.62 deg2 of the UltraVISTA/COSMOS field. The UltraVISTA survey itself observes across four broadband filters Y,J,H and Ks, as well as one narrowband filter centered on $H\alpha$ at z=0.8, which is not used in the present work (McCracken et al., 2012). The NIR coverage of UltraVISTA probes the rest-frame optical for galaxies at $1.5<z<4.0$, an important range for the estimation of stellar population properties. The catalog from which the present sample was selected was constructed using the first data release of the UltraVISTA survey; this data release spans about one season of observing time, with a $5\sigma$ depth of approximately $K_{s}<23.9$ AB in a 2*′′* aperture.

The UltraVISTA catalog spans 30 filters from $0.15\mu m$ to $24.\mu m$, including ultra-violet (UV) imaging from the GALEX satellite (Martin et al., 2005), and infrared coverage from Spitzer (Sanders et al., 2007; Frayer et al., 2009). Sources are selected from the $K_{s}$ band, which reaches 90% completeness at $K_{s,TOT}=23.4$ AB, which corresponds to a stellar mass 90% completeness limit of $\log(\text{M}_{\star}/\text{M}_{\odot})=10.43$ at $z=2$, as shown in Figure 1. This broad wavelength coverage in the UltraVISTA/COSMOS field allows us to confidently identify high-$z$ massive galaxies candidates for spectroscopic follow-up. Figure 1 shows the stellar mass ($M_{\star}$) versus redshift ($z$) and $K_{s}$ magnitude ($K_{s}$,tot) versus $M_{\star}$ diagrams with the five galaxies targeted by our spectroscopic program highlighted as colored filled circles.

We originally selected all galaxies at $1.7<z<2.7$ with $\log{(M_{\star}/M_{\odot})}>􏰀11.5$, for a total of 18 candidates spanning the range in KS-band total magnitude of $19.4<K_{\rm S,tot}<23.0$. From this sample, we selected all galaxies (eight in total) brighter than $K_{\rm S,tot}=20.6$ to observe spectroscopically. Due to scheduling constraints, two of these galaxies were not observed. The exclusion of these two galaxies does not introduce any additional biases. Therefore, the six galaxies presented in our work constitute a representative sample of the population of brightest and most massive galaxies in the redshift range $1.5\lesssim z\lesssim 3.0$.

Table 1 gives the RA and DEC of each of the sources in our sample along with observing and exposure information.

3 Data

Four of the five galaxies in the sample (COS-75355, COS-90676, COS-207144, and COS-189962) were observed with the Gemini Near-infrared Spectrograph (GNIRS, Elias et al., 2006); the remaining two galaxies were observed with the Folded-port InfraRed Echellette (FIRE, Simcoe et al., 2013) mounted on the Magellan 1 telescope. All of the galaxies in our sample were additionally imaged with HST WFC3 in the F160W band (GO-12990; PI: Muzzin).

3.1 GNIRS

Of the four galaxies observed with GNIRS, one (COS-189962) was discarded from the sample due to poor seeing during the observations.

GNIRS was set in its cross-dispersed mode with 32 line/mm grating and a 0.675” slit during the observing runs, which spanned 6 dates in March, May, and December of 2013. The observations were made in queue mode, and are taken in 600 second exposures; the details of the observation runs are found in Table 1. The telluric emission lines present in the NIR place a limit on the per-exposure observation time; we therefore opt to take many shorter exposures and subsequently stack the resulting spectra in order to boost S/N. The adopted GNIRS configuration provides a resolution of $R\approx 800$ and a wavelength range of around 0.8 to 2.4 $\micron$. A blind offset star and a B9V star were also observed for acquisition and to derive a telluric correction. Seeing averaged around 1.0*′′* across the observing run, with minimal cloud cover (photometric conditions).

3.2 FIRE Observations

FIRE observations were made in its echellete mode with a 0.6*′′* slit and practical coverage from approximately $1.0\mu{\rm m}$ to $2.0\mu{\rm m}$ with a spectral resolution of $R\approx 6000$.

Observations were taken across two nights in May of 2013. As was the case for the GNIRS observations, several short exposures were taken of each target at 908.8 seconds per exposure due to the telluric emission lines present in the observed wavelength range. An A0V star was also observed in order to perform telluric corrections. Seeing averaged between 0.7*′′* and 1.0*′′* throughout the observing run, with minimal cloud cover.

3.3 HST Imaging

The galaxies at hand span a significant variety of morphologies, as shown in Figure 2: COS-75355 appears to be consistent with a galaxy undergoing one to two minor mergers and COS-37207 is consistent with being a massive face-on red spiral, while COS-71929 and COS-90676 are relatively isolated with some diffuse emission. For a quantitative discussion of these morphologies, see Marsan et al. (2017; in prep.). COS-207144, though identified as a single object in the UltraVISTA catalog, is clearly resolved in two components in the HST image (see Figure 2). COS-207144 is the only object in our sample that falls within the 3D-HST survey (Skelton et al., 2014). In 3D-HST, COS-207144 is resolved as two objects, consistent with being an interacting pair at $z=2.02^{+0.09}_{-0.14}$ and $z=2.07\pm 0.05$. There exists a grism redshift for the latter object at $z_{\rm grism}=2.362\pm 0.015$ (Brammer et al., 2012; Momcheva et al., 2016). However, this grism redshift is based upon a low signal-to-noise, likely spurious, feature identified as $[OII]$ which, if adopted, represents a shift of $>5\sigma$ from the original photometric redshift, implying that the two resolved objects are not physically interacting. We therefore choose to compare the results of our analysis with those obtained using the original photometric redshift from Skelton et al. (2014).

3.4 Data Reduction and Extraction of Spectra

We performed the initial steps of the GNIRS data reduction (cosmic ray detection, bias subtraction, flat-fielding) using pyraf. We then performed a first-pass sky subtraction by constructing a sky frame using up to 4 dithered exposures taken at sufficiently similar times as the target expsoure. We then ran a second pass sky subtraction by using the IRAF task background on the target exposure, which fits the sky on each side of the galaxy continuum in order to remove remaining telluric emission lines present in the frame.

To extract the spectra to one dimension, we used the spectrum of the telluric star as a trace for any remaining non-linearity in the spectrum shape. We then performed a telluric correction using the one dimensional spectrum of the observed telluric star. The orders (GNIRS’s cross-dispersed mode yields six partial spectra corresponding to different wavelength ranges) were then concatenated using a weighted average and flux calibrated using the UltraVISTA $H$ and $K_{s}$ photometry.

The IDL-based reduction pipeline FIREHOSE (Matejek & Simcoe, 2012) was used to reduce the FIRE spectra. We adopted manual tracing for these spectra to account for the faintness of the targets, but left the pipeline otherwise intact.

In order to increase the signal-to-noise ratios of these observations, we binned each spectrum with bins of variable size. As a result, bin size is not constant across a given spectrum – areas of low signal-to-noise are binned more as compared to areas in which the signal-to-noise was relatively high in the unbinned spectrum. Areas of low transmission, strong skyline residuals, or artifacts from remaining skyline residuals were masked. The resulting bin size (observed frame) ranges from 51Å ($15^{\rm th}$ percentile) to 163Å ($85^{\rm th}$ percentile), with a median bin size of 77Å. The resulting signal-to-noise ratio ranges from 2.51 ($15^{\rm th}$ percentile) to 13.34 ($85^{\rm th}$ percentile), with a median signal-to-noise ratio of 6.99.

Figure 3 shows the final binned 1D spectra of the five sources (red filled circles), overplotted with the UltraVISTA broadband and medium-band photometry (green filled circles).

4 SED Modeling

We combined the binned spectra with the UltraVISTA photometry to model the spectral energy distributions (SEDs) of galaxies using the photometric redshift code EAZY (Brammer et al., 2008) and the stellar population synthesis code FAST (Fitting and Assessment of Stellar Templates, Kriek et al., 2009). The addition of the spectra allows us to derive continuum redshifts and stellar population properties from the two codes by taking advantage of the coverage provided by the broadband photometry and the finer observation of the 4000Å break delivered by the spectroscopy, producing robust redshift estimates in the absence of emission or absorption features.

Following Muzzin et al. (2013b), EAZY fits a linear combination of six templates from the PEGASE models (Fioc & Rocca-Volmerange, 1997), a red template from the model library of Maraston (2005), a $\sim$1 Gyr old post-starburst template, as well as a slightly dust-reddened Lyman Break template. EAZY also incorporates a template error function to account for larger uncertainties in the models in the rest-frame NIR relative to the rest-frame optical. A detailed description of this process may be found in Brammer et al. (2008). EAZY was run using the v1.0 template error function and a K-band magnitude prior.

We used FAST to estimate stellar population properties. FAST produces best-fit parameters for stellar mass, star formation rate, stellar age, characteristic timescale ($\tau$) of the assumed star formation history (i.e., delayed exponentially declining star formation history), extinction (using the extinction curve detailed in Kriek & Conroy 2013) and metallicity. In this analysis we set an allowed stellar age range of $\log{(\text{age}~{}[\text{yr}]~{})}=8-10.1$, an allowed range of $\tau$ of $\log{(\tau~{}[\text{yr}]~{})}=7-10$, and an allowed extinction range of $A_{\rm V}=0-5$ mag.” We used the high-resolution version of the flexible stellar population synthesis (FSPS) models from Conroy et al. (2009a, b) with a Chabrier (2003) IMF.

In order to estimate errors on the stellar population properties derived from FAST, we constructed a Monte Carlo (MC) simulation external to FAST and EAZY. The input photometry and spectra were perturbed according to a normal distribution described by $\sim~{}N(f_{i},\sigma_{i})$, where $f_{i}$ is the mean flux of the $i^{\text{th}}$ measurement, and $\sigma_{i}$ is the uncertainty of the $i^{\text{th}}$ measurement in the spectrum and photometry. At this juncture, we implicitly assumed that each wavelength bin is uncorrelated from its neighbors; we do expect there to be some degree of correlation between adjacent flux bins in the spectra, but the impact of a bin-to-bin correlation will be an overestimation of the MC-derived errors. We therefore chose to assume that each bin is fully independent of the state of nearby wavelength bins.

For those galaxies that have spectroscopic redshifts, the FAST redshift was held fixed for all runs. For the galaxies that do not have spectroscopic redshifts, we run EAZY using the perturbed spectrum and photometry, and run FAST with the perturbed data, holding the redshift fixed to the best-fit redshift of the corresponding EAZY output.

This method of estimation allows for the independent derivation of errors on the stellar population properties estimated by FAST and EAZY. For all galaxies, 500 simulations were run for both FAST and EAZY. A full description of the parameter correlations that result from these simulations is given in Appendix A.

5 Results

5.1 Spectroscopic Redshifts

There are spectroscopic redshifts for 2 galaxies in the sample. First, COS-75355 has a literature $z_{\text{spec}}=1.822$ from Onodera et al. (2012). Secondly, the spectrum of COS-71929 shows absorption features consistent with CaII H & K, allowing us to measure a spectroscopic redshift at $z_{\text{spec}}=2.000\pm 0.006$.

Onodera et al. (2012) finds a spectroscopic redshift for COS-75355 using the CaII H & K absorption features, in conjunction with $\text{H}_{\beta}$ in absorption. The spectrum in our sample does not have sufficient signal-to-noise to detect most of the features used to identify the spectroscopic redshift – there is a possible detection of Na in absorption at $z=1.822$, but the detection is not secure.

The features used to determine a spectroscopic redshift for COS-71929 are shown in Figure 4, as labeled. The best-fit EAZY model SED is overlaid in brown. From this spectrum we constrained our spectroscopic redshift for COS-71929 to be $z_{\rm spec}=2.0000\pm 0.0055$, which is within one sigma of the continuum redshift, as detailed in the following section.

We simultaneously fit the CaII H & K and Hδ features using a triple Gaussian with the observed velocity dispersion and redshift constrained to be equal for all features. The underlying continuum was approximated as linear in this regime. From this treatment, we derived a formal error for the spectroscopic redshift of $\pm 0.0055$. We note that the spectroscopic redshift for COS-71929 is located at the secondary peak in the galaxy’s redshift probability density function, with the bimodality of the probability density function likely to be caused by a degeneracy between the locations of $\rm{CaII}~{}H$ and Hδ.

We additionally attempted to approximate a constraint on the stellar velocity dispersion of COS-71929 using the fit detailed above. The instrumental spectral resolution is given to be 50 km s*-1*; we subtracted this value in quadrature from the observed velocity dispersion to make an estimate of the stellar velocity dispersion. We estimated the uncertainty in the fit due to measurement error by running a Monte Carlo simulation to convergence in which the binned flux values were perturbed according to a normal distribution with a standard deviation equal to the measurement uncertainty. This procedure gives a velocity dispersion of $549^{+106}_{-262}$ km s*-1*, where the quoted errors are the 2-sided 95% confidence bounds. However, we stress that our measurement of the velocity dispersion is very crude. The low SNR of the spectrum does not allow for a proper analysis of the velocity dispersion, and our measurement is subject to extremely high uncertainty, with a 1-sided 99% lower confidence bound of 217 km s*-1*.

5.2 Continuum Redshifts

For all of the galaxies in the sample (including those that have spectroscopic redshifts), we estimated a continuum redshift with EAZY ($z_{\rm cont}$), in which both the photometry and spectroscopy are used in the fit. Figure 5 shows the continuum redshifts, $z_{\rm cont}$, for the galaxies in our sample against their original photometric redshifts ($z_{\rm phot}$, derived from the UltraVISTA DR1 photometry only). The quoted errors in Figure 5 are the 68% confidence intervals from the EAZY redshift probability density function.

The continuum redshifts of both COS-75355 and COS-71929 are in good agreement ($<1\sigma$ difference) with their spectroscopic redshifts. The spectroscopic redshifts of both galaxies are shown as horizontal lines in Figure 5. Although the errorbars shown on the figure are taken from the EAZY redshift probability density function (PDF), the distribution of our EAZY Monte Carlo simulation are in good agreement with the output EAZY PDF. The main peak of the EAZY redshift probability density function is somewhat higher than the spectroscopic redshift for COS-71929 (see Figure 3), but there is a secondary peak in the EAZY redshift PDF that is consistent with our spectroscopic redshift. Furthermore, the spectroscopic redshift is within $1\sigma$ of the continuum redshift for the galaxy.

The formal errors on the continuum redshift are lower than the original photometric redshift in all cases when the spectrum is added to the fit, though COS-37207 is the only galaxy whose continuum redshift is significantly different than the original photometric redshift. Although the original photometric redshift of COS-37207 lies outside of the continuum redshift 95% confidence interval, the absolute difference between the continuum and photometric redshifts is relatively small ($\sim 0.16$; i.e., $\Delta z/(1+z)\sim 0.05$”).

As the photometric redshifts for both of the galaxies that make up COS-207144 are available, it is informative to compare our results against the photometric redshifts available from (Skelton et al., 2014). Given that we expect the majority of the light in the observed spectrum to originate from the brighter galaxy (which, in the NIR, is the galaxy at $z=2.02$), it is unsurprising that our continuum redshift is in better agreement with the brighter of the two objects ($<1\sigma$ difference) than with the fainter ($<2\sigma$ difference).

5.3 Stellar Population Properties

We used the procedure outlined in Section 4 to estimate the stellar population properties of the sample. Table 3 shows the FAST output for the five galaxies, along with their 1$\sigma$ errors. The best-fit values that we quote are the original FAST best-fit values, and not average values of the Monte Carlo distributions. Thus, in some cases, the best-fit value lies on the edge of the 1$\sigma$ interval. Re-defining our best-fit values by using an average of the Monte Carlo distributions does not qualitatively change our results.

We confirm that all of the galaxies in our sample are indeed ultra-massive. All galaxies in the sample have a stellar mass exceeding $\log(\text{M}_{\star}/\text{M}_{\odot})\gtrsim 11.5$. In the case of COS-207144 (the blended object), stellar mass estimates obtained here may be compared against those obtained for the objects resolved in 3DHST by Skelton et al. (2014). The best-fit stellar mass of the brighter galaxy of the pair is in good agreement (within 68% confidence interval) with that of COS-207144. The sum of the best-fit stellar masses of the resolved objects is significantly larger ($\log(M_{*}/M_{\odot})=11.62$) than our best-fit mass for COS-207144. This is likely a result of the brighter galaxy dominating the light of the spectrum – the best-fit stellar mass of COS-207144 ($\log(M_{*}/M_{\odot})=11.46^{+0.00}_{-0.26}$) is in strong agreement with the best-fit stellar mass of this source as derived by Skelton et al. (2014) of $\log(M_{*}/M_{\odot})=11.40$.

Given that the original selection criteria from the UltraVISTA DR1 catalog was based partially on stellar mass, we plot the stellar mass FAST estimated from photometry alone against the same quantity from the photometry in conjunction with the spectroscopy in Figure 6. None of the updated best-fit stellar masses differ significantly from the stellar mass derived from the photometry alone. In the case of COS-37207, the galaxy for which the best-fit redshift changed the most significantly when the spectrum was added to the fit (see Figure 3), the formal error on the stellar mass increased when the spectrum is added, with relatively little change to the best-fit stellar mass.

As a probe of the evolutionary stage of the galaxies in our sample, we computed the quenching factor $q_{\text{sf}}$ (Kriek et al., 2009) , which is $q_{\rm sf}=1-b$, where $b$ is the birthrate parameter (Kennicutt, 1983; Scalo, 1986). The quenching factor is therefore

[TABLE]

such that a fully quenched galaxy will have a quenching factor of 1. With the exception of COS-37207, which has a quenching factor of $0.276^{+0.529}_{-0.276}$, all of the galaxies in our sample have quenching factors consistent with $q_{\text{sf}}>0.90$.

COS-37207 is the only galaxy in the sample to show high levels of star formation, with a best-fit SED star formation rate (SFR) of $339^{+108}_{-221}M_{\odot}~{}{\rm yr}^{-1}$. The error bars quoted here correspond to the 68% confidence interval. The 99% confidence interval bounds are [31,788] $M_{\odot}~{}{\rm yr}^{-1}$, the lower bound of which is greater than the best-fit SFR for any of the other galaxies in the sample. We therefore conclude that this galaxy is securely identified as star forming. In light of the fact that this galaxy has a high estimated SFR, but is also characterized by a large dust extinction of $A_{V}=1.7^{+0.1}_{-0.4}$, it is of interest to note that the spectrum of COS-37207 does not show any secure emission lines. Using the best-fit values of SFR and $A_{V}$, we estimated the expected observed flux for $[OII]$, $H_{\beta}$, and $[OIII]$ using the relations from Calzetti et al. (2000) and Price et al. (2014). From these estimates, we conclude that neither $H_{\beta}$ nor $[OII]$ would be observable given the noise properties of the spectrum. $[OIII]$ may have been visible, but coincides with a strong skyline at the best-fit redshift of $z=2.11$. It is therefore expected that strong emission lines are absent from the spectrum of COS-37207.

Though COS-71929 also shows evidence of some star formation (SFR $=8.7^{+0.4}_{-6.9}M_{\odot}~{}{\rm yr}^{-1}$), a measurement of the quenching factor of this galaxy indicates that the SFR was significantly higher in the past ($q_{\text{sf}}\approx 0.95$).

Additionally, the estimated dust extinction for COS-37207 is the highest in the sample, with $A_{V}=1.70^{+0.10}_{-0.40}$ mag. COS-71929 displays the second-highest level of dust extinction at $A_{V}=1.30^{+0.10}_{-0.40}$ mag. We therefore characterize COS-37207 as a dusty star-forming galaxy, COS-71929 as a galaxy with mostly suppressed star-formation intermediate between star-froming and quiescent galaxies, and the rest of the galaxies in our sample as quiescent.

We plot our stellar mass and SFR estimates against the mass-complete sample of galaxies at $1.7<z<2.7$ from the UltraVISTA catalog in Figure 7, along with the star-forming main sequences measured by Speagle et al. (2014), Whitaker et al. (2014), and Tomczak et al. (2016). For the UltraVISTA catalog, we preferentially use UV+IR SFR measurements. In the abscence of such a measurement (or a non-detection), we use the SED SFR. The FAST-derived stellar masses and SFRs for our sample place the three quiescent galaxies (COS-75355, COS-207144, and COS-90676) on the high-mass end of the quiescent population. COS-37207 lies along the high-mass end of the star-forming main sequence, making it the most massive star-forming galaxy in the targeted redshift range. Finally, COS-71929, with a SFR$\approx$10 M*⊙* yr*-1* is seen to lie well below the main sequence of star-forming galaxies, indicating that the star formation activity of COS-71929 is mostly quenched. We note that this picture qualitatively holds if the SED-modeling SFRs are replaced with the SFR from the combination of ultra-violet (i.e., unobscured) and infrared (dust reprocessed) estimates of the SFRs (see Section 5.5).

All of the galaxies in our sample have nominal best-fit stellar ages greater than or equal to 1 Gyr. Here, stellar ages are given as the age since the onset of star formation. COS-207144 and COS-71929 have best-fit ages comparable to the age the universe at their redshift, whereas COS-75355 and COS-90676 have best-fit stellar ages a factor of 2-3 younger than the age of the universe at their redshift. The star-forming galaxy COS-37207 is 0.5 dex younger that the age of the universe at z=2.1. We note however that the stellar age is typically one of the most uncertain stellar population properties (Muzzin et al., 2009).

Given the analysis presented in Appendix A and summarized in Section 5.4 the metallicity is not constrained due to degeneracies in the fitted stellar population synthesis parameters.

5.4 SPS Parameter Degeneracies

Because much of this analysis is based upon conclusions reached through the use of SPS modeling, it is imperative to examine the degeneracies between the derived parameters of our SPS models.

In Appendix A, we plot the correlations between EAZY estimates of redshift, FAST estimates of stellar mass, SFR, dust attenuation, and population age for the Monte Carlo simulations of both the photometry only and the photometry + spectrum. We use a Gaussian kernel density estimate in order to visualize the distribution of the output parameters.

As expected, the FAST estimates increase in precision the most when a spectroscopic redshift is present. Though COS-37207 does not have a spectroscopic redshift, its best-fit redshift shifted signficantly with the addition of the spectrum (see Figure 3); the error on the FAST parameter estimates decreased slightly with the addition of the spectrum, but the best-fit FAST parameter estimates remain relatively unchanged.

As the only vigorously star-forming galaxy in the sample, it is of interest to note that there exists a strong correlation between dust attenuation and SFR in the Monte Carlo distribution of COS-37207. This correlation is not found in any of the other galaxies’ best-fit parameter distributions. This finding does not change the classification of COS-37207 as dusty and star-forming, i.e., even at the lowest dust attenuation and SFR produced by the simulation, COS-37207 would be classified as dusty and star-forming. However, it reaffirms that low-resolution spectra are not sufficient to break the degeneracy between SFR and dust attenuation (Förster Schreiber et al., 2004; Muzzin et al., 2009).

As expected from other similar studies (e.g., Worthey, 1994; Muzzin et al., 2009; López Fernández et al., 2016), we find several other correlations between derived stellar population properties. Stellar population age, dust attenuation, and metallicity all redden the overall galaxy SED, and in the absence of prior knowledge that can break degeneracies (Conroy & van Dokkum, 2012), it is expected that there will be anti-correlations between these parameters as derived from SED modeling. For all of galaxies in our sample, we see the expected negative correlation between stellar population age and dust attenuation to some degree. The age-extinction correlation is especially prominent in COS-207144 and COS-75355. This is likely because the addition of the spectrum (and spectroscopic redshift, in the case of COS-75355) eliminated several competing parameter sets that produced discrete populations within the Monte Carlo distribution, as can be seen in the third panel of the left column of Figure 11.

5.5 $\text{M}_{\star}$ and SFR: Comparison with Independent Estimates

In order to assess the accuracy of our previously-derived stellar population properties, we compared our results to quiescent/star-forming classifications using the rest-frame U-V versus V-J color-color diagram (hereafter UVJ diagram; following, e.g. Whitaker et al., 2011; Muzzin et al., 2013b; Marchesini et al., 2014; Whitaker et al., 2015; Martis et al., 2016), and SFRs derived from the UV and infrared (IR), as provided by Muzzin et al. (2013a).

Figure 8 shows the targeted sample in the UVJ diagram using their original photometric redshifts from Muzzin et al. (2013a). Also plotted is the overall galaxy population in a stellar mass complete sample at the same redshift range as obtained using the UltraVISTA DR1 catalog of Muzzin et al. (2013a). Figure 8 shows that the location of the targeted galaxies in the UVJ diagram is consistent with the characterization of their stellar population properties as derived from FAST. In particular, we note that COS-37207 is in the region of the UVJ diagram typically populated by dusty star-forming galaxies (e.g., Martis et al., 2016), while COS-71929 is in an intermediate zone between the dusty star-forming and quiescent galaxies.

We checked the SED star formation rates for our sample against star formation rates derived using LIR + LUV measurements based on the UltraVISTA catalog. We determine the rest-frame UV flux following Muzzin et al. (2013a), estimating $L_{2800}$ by integrating the best-fit template generated by EAZY (Brammer et al., 2008) from 2600–2950 Å.

Infrared luminosity ($L_{\text{IR}}$) was estimated using the $24~{}\mu m$ emission alone; we caution that this approach can produce individual measurements of $L_{\text{IR}}$ which are uncertain up to several factors. A full explanation of the method to determine $L_{\text{IR}}$ for the sample at hand can be found in Wuyts et al. (2008); briefly, $L_{\text{IR}}$ is computed for the Dale & Helou (2002) infrared SEDs of star-forming galaxies at several heating levels of the interstellar environment. The best estimate of $L_{\text{IR}}$ is then taken to be the log-average of the set of $L_{\text{IR}}$ estimates for the template across the range of heating levels sampled.

We calculate IR+UV SFRs following Kennicutt (1998) using the calibration of Whitaker et al. (2014), i.e. $\text{SFR}_{\rm IR+UV}=1.09\times 10^{-10}({\rm L}_{\rm IR}+2.2{\rm L}_{\rm UV})L_{\odot}$. Bell et al. (2005) presents a near-identical calibration of $\text{SFR}_{\rm IR+UV}$, which would yield the identical values of $\text{SFR}_{\rm IR+UV}$, within errors. The upper limits in the IR+UV star formation rates indicate sources for which the MIPS $24\mu m$ flux has a signal-to-noise ratio ${\rm S/N}<3$, and the $3\sigma$ upper limit is adopted to estimate $L_{\rm IR}$.

As shown in Figure 9, for four out of the five galaxies in our sample, the SFR estimated from the SED and IR+UV are consistent with each other. Given that this includes three upper limits, however, this statement demonstrates only that the SED SFR estimates are not significantly higher than expected. COS-207144 is the only galaxy in our sample with a SFR derived from UV+IR much larger than the SFR from SED modeling. This result may be due to a difference in the UV and IR properties of the two galaxies which make up the source.

6 Discussion and Conclusions

In this work, we have presented NIR spectra of 5 ultra-massive galaxies in the redshift range of $1.7<z<2.7$. We determined a spectroscopic redshift for one of these galaxies based upon the presence of the CaII H & K absorption features in its spectrum, and estimated continuum redshifts for the remaining galaxies in our sample using information from both broad- and medium-band UV-to-8 $\mu m$ photometry and the aforementioned spectra. The stellar population parameters of the galaxies were also estimated, confirming the ultra-massive status of these galaxies.

Out of the five galaxies, only one – COS-37207 – is found to be star-forming, with SFR $\approx 340~{}\text{M}_{\odot}~{}yr^{-1}$. This galaxy is also found to have the highest level of dust attenuation, with $A_{V}$ $\approx 1.70$ mag. A second galaxy, COS-71929, possesses properties that are intermediate with respect to the quiescent and star-forming populations at this redshift. The galaxy shows residual star formation activity, i.e., SFR$\approx$9 M⊙ yr*-1* and elevated dust obscuration, i.e., A${}_{\rm V}\approx$1.3 mag, and sits well below the main sequence of star forming galaxies at the redshift range of interest. The remaining galaxies in the sample are classified as quiescent with substantially less dust extinction by both SED modeling and UVJ diagram classification. This finding is in agreement with refined evolutionary path for the formation of local ultra-massive galaxies recently presented by Marchesini et al. (2014). Though the progenitors of today’s ultra-massive galaxies (i.e., $\text{log}~{}(\text{M}_{\star}/\text{M}_{\odot})~{}11.8$ at $z\sim 0$) were massive, highly obscured, dusty star-forming galaxies at $2.5<z<3$, we expect the sample at hand to quench at earlier times on average, as the galaxies in the present sample have significantly higher stellar masses than the majority of the galaxies from Marchesini et al. (2014). Though the present sample is too small to estimate population characteristics, if we consider COS-37207 to be the only star-forming galaxy (grouping COS-71929, given its suppressed star formation, with the quiescent galaxies in the sample) we obtain a quiescent fraction of 0.80 for such galaxies ($11.45<\text{log}~{}(\text{M}_{\star}/\text{M}_{\odot})<11.9$) at $1.7<z<2.7$.

There is not significant evidence that this quiescent fraction is systematically larger than the fraction of $\sim 0.40$ estimated by Marchesini et al. (2014) for the population of the progenitors of today’s ultra-massive galaxies at $1.5<z<2.5$. Marchesini et al. (2014) found a quiescent fraction of $\approx 0.40$ for their sample at $1.5<z<2.5$, which covered a range of stellar masses of $11.2\lesssim\text{log}~{}(\text{M}_{\star}/\text{M}_{\odot})\lesssim 11.5$. Similarly, one out of the three galaxies with stellar masses $\text{log}~{}(\text{M}_{\star}/\text{M}_{\odot})\gtrsim 11.4$ in Belli et al. (2014) was found to be star-forming at $2.<z<2.5$, and Muzzin et al. (2013b) estimates that the quiescent fraction should meet or exceed 0.50 at $\text{log}~{}(\text{M}_{\star}/\text{M}_{\odot})\approx 11.3$ for the redshift range considered in this work.

Assuming a true quiescent fraction of 0.40, there is a significant ($\approx 8\%$) chance that a sample of five galaxies will contain at least four quiescent galaxies if the sample were drawn without bias for the same distribution. However, given previous evidence that more massive galaxies enter quiescence at earlier times relative to lower mass galaxies (see, e.g., Brammer et al., 2011; Muzzin et al., 2013b), we find that our data are consistent with these claims. However, in order to probe the true mass dependence of the quiescent fraction at a given redshift of such ultra-massive galaxies a larger sample is needed in order to properly marginalize over confounding variables, and to suppress the uncertainty due to low number statistics.

In order to better characterize the population of ultra-massive galaxies in the early universe, we must increase the number of identified ultra-massive galaxies and construct a larger number of well-sampled SEDs from which to derive stellar population property estimates and photometric redshifts. The fact that the inclusion of our NIR spectra, in the abscence of a spectroscopic redshift, do not greatly affect the stellar masses and redshifts of the galaxies in our sample implies that such galaxies may be reliably identified using photometry alone, in agreement with Muzzin et al. (2009). The NEWFIRM Medium-Band Survey II (NMBS-II) will identify these ultra-massive galaxies across 5.2 ${\rm deg}^{2}$ of the sky (the UltraVISTA DR1 catalog, upon which this work is based, spans 1.62 square degrees), and will provide medium-band photometry to produce well-sampled SEDs. Moreover, the VISTA VIDEO survey (Jarvis et al., 2013) includes oft-visited areas of the sky and covers nearly 12 ${\rm deg}^{2}$, giving it both the wavelength range and spatial coverage to robustly identify many new ultra-massive galaxies.

Deeper spectroscopy will also greatly aid the study of this sample – the galaxies at hand are bright enough that deep spectroscopy is feasible from the ground, given an 8-10m class telescope. Though our detection of absorption features in COS-71929 was not high enough S/N to perform a detailed measurement of the velocity dispersion, preliminary estimates indicate a substantial velocity dispersion. Such measurements would also be greatly helpful in constraining the systematics of SED modeling by providing independent measurements of output parameters.

Acknowledgments. EKF acknowledges the support of the Tufts University Summer Scholar Program. DM acknowledges the support of the Research Corporation for Science Advancement’s Cottrell Scholarship and from Tufts University Mellon Research Fellowship in Arts and Sciences. ZCM gratefully acknowledges support from the John F. Burlingame and the Kathryn McCarthy Graduate Fellowships in Physics at Tufts University. DM and ZCM acknowledge support from the program HST-GO-12990, provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under the NASA contract NAS5-26555, and support from the National Aeronautics and Space Administration under Grant NNX13AH38G issued through the 12-ADAP12-0020 program. KEW gratefully acknowledge support by NASA through Hubble Fellowship grant#HF2-51368 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA. GHR acknowledges the support of NASA grant HST-GO-12590.011-A, NSF grants 1211358 and 1517815, the support of an ESO visiting fellowship, and the hospitality of the Max Planck Institute for Astronomy, the Max Planck Institue for Extraterrestrial Physics, and the Hamburg Observatory. GHR also acknowledges the support of a Alexander von Humboldt Foundation Fellowship for experienced researchers. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile. This work made used of the products from the 3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

Appendix A Monte Carlo Parameter Correlations

Here we present the results of the Monte Carlo simulations run to constrain the error on our stellar population property estimates. For each run, fluxes for both the spectroscopy and photometry are drawn from a normal distribution characterized by the error on the original measurement. Each instance is run with a fixed redshift drawn according to the best-fit EAZY redshift PDF, as shown in Figure 3.