Resolving Quiescent Galaxies at z ≳ 2. I. Search for Gravitationally Lensed Sources and Characterization of Their Structure, Stellar Populations, and Line Emission

Massive quiescent galaxies at high redshifts are relatively rare. The probability to find such a galaxy magnified by a given cluster is the product of the population's surface density and the source plane area magnified by the cluster above some minimum magnification of interest. The typical value of this area depends on the cluster sample, so only rough estimates are possible. An upper limit can be estimated from the Hubble Frontier Fields, which were chosen to be among the most powerful cluster lenses. For these clusters, Johnson et al. (2014) derive a typical source plane area of ≃0.3 arcmin² magnified by μ > μ_min = 3. This area approximately scales as ${\mu }_{\min }^{-2}$. These figures enable a rough estimate of the frequency of lensed sources magnified above some limiting flux. We consider a source population of quiescent galaxies at z = 1.8–3, whose magnitude distribution we estimated using the NEWFIRM Medium Band Survey (NMBS) catalogs (Whitaker et al. 2011). Here quiescence is defined using the UVJ criterion in the form presented by Whitaker et al. We find a probability of 0.02 per cluster to find such a galaxy magnified to H < 20. As we will show, this flux limit is of interest because it is brighter than the unlensed examples that are found in the widest extragalactic deep fields. Although the estimated probability is uncertain and should be considered to be an upper limit, it shows that one must search hundreds of clusters to find a few very bright, lensed quiescent galaxies at these redshifts.

Although they are rare, massive quiescent galaxies at z ≃ 2–3 can be efficiently identified with two or three filter imaging, and their high NIR surface brightnesses, imply that relatively shallow exposures are sufficient. This makes it feasible to search large numbers of clusters.

Quiescent galaxies at high redshifts can be identified on the basis of their red colors at NIR wavelengths, which are produced by the redshifted Balmer/4000 Å break (Franx et al. 2003). The addition of an optical magnitude probing the rest-frame UV further helps to eliminate dusty star-forming galaxies (Daddi et al. 2004). We have applied this basic strategy to two imaging surveys of massive clusters: (1) a dedicated campaign with the FourStar camera at the Magellan Baade telescope, and (2) a search of archival data from Wide Field Camera 3 (WFC3) onboard HST. These surveys will be described in Sections 2.3 and 2.4.

Using the NMBS catalogs, we experimented with color cuts that efficiently identify a target population of UVJ-quiescent galaxies at z = 2–3 with H < 22.5, which corresponds to the desired lensed magnitude limit of H < 20 for a maximum expected magnification of μ = 10. (The precise limit is not relevant since the colors do not vary drastically with magnitude.) The goal was to balance a high completeness for selecting this target population (orange circles in Figure 1) with minimal contamination from star-forming or lower-redshift galaxies (black points). The FourStar filter set includes the broadband J, H, and K_s filters and the medium-band J₁, J₂, J₃, H_s, and H_l filters (Persson et al. 2013). We found that the combination ${J}_{1}-H\gt 1.7$ is nearly optimal, as shown in Figure 1. A second cut of r − H > 2.8 is also plotted. For this filter combination, the optical-IR cut removes only a few galaxies, so it is not necessary to identify bright quiescent candidates in our FourStar imaging: ${J}_{1}-H$ is sufficient. For our HST archival search, we adapted these color cuts based on the available filters, as described in Section 2.4. The optical-IR cut plays a greater role for some of these filter combinations.

Zoom In Zoom Out Reset image size

The color cuts ${J}_{1}-H\gt 1.7$ and r − H > 2.8 select 85% of the target population, so our selection is reasonably complete. Furthermore, 31% of the color-selected galaxies are quiescent galaxies at z > 2, and of the remainder, an additional 33% are quiescent galaxies at slightly lower redshifts z = 1.6–2, which are still of interest. This level of purity makes it possible to pursue efficient spectroscopic follow-up to confirm the redshifts and quiescent nature of the sources. As outlined later, we have examined the fields of 232 clusters and located five color-selected galaxies that are magnified above the H < 20 flux limit. Fainter color-selected galaxies were also identified, but in this paper we confine ourselves to H < 20 sample, for which our spectroscopic follow-up is complete.

We imaged 131 clusters through the J₁ and H filters with FourStar over five semesters from 2014A to 2016A. Targets were drawn from several sources. Approximately 40% were X-ray–selected clusters from the Massive Cluster Survey (MACS; Ebeling et al. 2001; Repp & Ebeling 2018). Although these were included in the target lists of the HST snapshot programs described in Section 2.4, they had not yet been observed with WFC3-IR. Another ≃40% were selected from the first Planck catalog of Sunyaev–Zel'dovich sources (Planck Collaboration et al. 2014). To maximize the lensing efficiency, we only considered those sources that had been confirmed as clusters at redshifts $0.3\lt z\lt 0.8,$ and we gave priority to those with higher signal-to-noise ratios in the Planck maps. Finally, ≃20% of our FourStar targets were optically selected clusters that were drawn primarily from the redMaPPer DR8 catalog (Rykoff et al. 2014) and prioritized by richness. A handful were drawn from Wong et al.'s (2013) catalog of fields suggested to be powerful lenses.

The typical observing sequence was to move to the target, correct the focus and mirror figure with the facility Shack–Hartmann system and obtain 5 × 5.8 s unguided exposures in the H band at each of 11 random dither positions within a $90^{\prime\prime} \times 90^{\prime\prime}$ box, followed by 2 × 32 s exposures at 13 random dither positions through the J₁ filter. The total wall clock time per cluster was approximately 30 minutes.

The data were reduced using the automated system described by Kelson et al. (2014), which produces stacked images for each cluster and filter with astrometry tied to 2MASS. Photometric calibration was obtained using stars in the 2MASS point source catalog (Skrutskie et al. 2006). This is straightforward for the H and K_s filters, which are present in the 2MASS catalog. For the medium-band filters, we used the Pickles (1998) stellar library to derive mean relations that relate the 2MASS magnitudes of a star to its J₁, J₂, and J₃ magnitudes.

Applying the cuts ${J}_{1}-H\gt 1.7$ and H < 20 to this sample of 131 clusters yielded five sources. Two of these sources are clearly magnified by the clusters PSZ1-G295.24–21.55 and MACSJ0138.0–2155. As the right panel of Figure 1 shows, these objects occupy a region of color–magnitude space that is virtually empty in extragalactic field surveys. The other three sources are compact and are located at least 4 arcmin from the cluster center. The FourStar field of view is much larger than the high-magnification region of a cluster, so these three sources are likely to be unmagnified systems on the very bright tail of the luminosity function. Their presence is not unexpected given that the FourStar survey encompasses 4.3 deg⁻²; however, since they are not magnified, we will not discuss these sources further. The yield of $2/131\simeq 0.02$ lensed quiescent galaxies per cluster is consistent with the upper limits roughly estimated in Section 2.1.

Many galaxy clusters have been observed with HST, but a much smaller subset has been observed with WFC3-IR through two filters and with ACS or WFC3-UVIS through at least one optical filter, as required to implement our color criteria. At the time of our search (2014–2016), the vast majority of such observations had been undertaken through one of three programs: (1) a snapshot imaging campaign based on the Massive Cluster Survey (MACS), led by PIH. Ebeling in Proposal IDs 10491, 10875, 12166, and 12884; (2) multiband imaging of lensing clusters identified in the Sloan Digital Sky Survey, led by PIM. Gladders in Proposal ID 13003, and (3) the Cluster Lensing and Supernova survey with Hubble (CLASH), a multi-cycle treasury program, led by P.I. M. Postman in Proposal ID 12065.

For all of the clusters included in these programs that, at the time of the analysis, had been observed through two WFC3-IR filters and one or more ACS or WFC3-UVIS filters, we produced multiband photometric catalogs and then applied color and magnitude criteria to search for lensed quiescent galaxies. Since the surveys used a variety of filters, we adapted the color criteria shown in Figure 1 for each survey. For all of the galaxies in the NMBS catalogs that were used to construct Figure 1, we integrated the spectral energy distribution fit to produce synthetic magnitudes in the relevant HST filters. We then reproduced Figure 1 with the appropriate filter combinations and adjusted the color thresholds to match the balance of completeness and contamination shown in the figure.

Within the MACS cluster sample, we processed 46 clusters with WFC3-IR and ACS imaging. The color criteria became (${\rm{F}}110{\rm{W}}-{\rm{F}}140{\rm{W}}\gt 0.7$) and (${\rm{F}}814{\rm{W}}-{\rm{F}}140{\rm{W}}\gt 2$ or ${\rm{F}}606{\rm{W}}-{\rm{F}}140{\rm{W}}\gt 2.5$). We also imposed a flux limit of ${\rm{F}}140{\rm{W}}\lt 20.7$, which reflects our fiducial $H\lt 20$ cut adjusted by the mean ${\rm{F}}140{\rm{W}}-H$ color of the target population. We identified one lensed galaxy behind cluster MACSJ0150.3–1005.

For the Gladders sample, we processed 30 clusters that had been imaged with WFC3-IR and WFC3-UVIS. The color criteria became (${\rm{F}}125{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 0.7$ or ${\rm{F}}110{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 1.0$ or ${\rm{F}}105{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 1.4$) and (${\rm{F}}606{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 3$) and ${\rm{F}}160{\rm{W}}\lt 20.3$. We identified one source magnified by the cluster SDSSJ1522 + 2535.

For the CLASH sample, 16 bands of imaging are available for 25 clusters. The color criteria became ${\rm{F}}110{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 0.9$ and ${\rm{F}}814{\rm{W}}-{\rm{F}}160{\rm{W}}\gt 2.5$ and ${\rm{F}}160{\rm{W}}\lt 20.3$. We identified one source magnified by MACSJ2129.4–0741.

The coordinates and basic properties of the five lensed color-selected galaxies that were located in our survey are listed in Table 1. We name these galaxies MRG-M0138, MRG-M0150, MRG-P0918, MRG-S1522, and MRG-M2129 based on the name of the lensing cluster, where the prefix MRG denotes Magnified Red Galaxy. Of these objects, MRG-M2129 was previously identified by Geier et al. (2013) and studied by Toft et al. (2017), while the other four were discovered in this survey. A pilot study of MRG-M0150 was presented by Newman et al. (2015a).

Table 1.  Sample of Lensed Quiescent Galaxies

Lensed Galaxy	R.A. (hr)	Decl. (deg)	${z}_{\mathrm{spec}}$	H (mag)	Lensing Cluster	zlens	Source of zlens
MRG-M0138	01:38:03.9	−21:55:49	1.95	17.1	MACSJ0138.0-2155	0.338	This paper
MRG-M0150	01:50:21.0	−10:05:14	2.64	19.6	MACSJ0150.3-1005	0.365	SDSS DR12 (Alam et al. 2015)
MRG-P0918	09:18:34.1	−81:03:08	2.36	19.4	PSZ1 G295.24-21.55	0.61	Planck Collaboration et al. (2014)
MRG-S1522	15:22:53.6	+25:35:49	2.45	19.8	SDSSJ1522 + 2535	0.58	SDSS DR12 (photometric)
MRG-M2129	21:29:22.3	−07:41:31	2.15	20.0	MACSJ2129.4-0741	0.589	Ebeling et al. (2007)

Note. In cases of multiple images, the coordinates and magnitudes listed are those of the brightest image.

Download table as:  ASCIITypeset image

Figure 2 shows HST images of the five lensed galaxies. We will now briefly describe each of these galaxies, their lensing configurations, and the lensing clusters. We note that the distinctive colors and high surface brightnesses of the lensed quiescent galaxies makes it straightforward to identify multiple images, which are confirmed by lens modeling in Section 5.

Zoom In Zoom Out Reset image size

MRG-M0138 is a remarkable system that presents five multiple images: two merging images forming a giant arc southward of the cluster center, a single image on the western side, a radial arc near the brightest cluster galaxy (BCG), and a central image (not visible in Figure 2). The brightest image (Image 1) has H = 17.1, which is astoundingly bright for this redshift: Figure 3 shows that this image is ∼3 mag brighter than the brightest unlensed quiescent galaxies at similar redshifts found in the 1.6 deg⁻² UltraVISTA field (Muzzin et al. 2013a)! Summing the two merging images 1 and 2, the total AB magnitudes of the giant arc are H =  16.5 and K_s =16.1. We believe this is the brightest giant arc known at NIR wavelengths, at least among z ≳ 1 sources. It is 1.1 mag brighter at K_s than the giant arc discovered by Dahle et al. (2016) and 2.3 mag brighter than the source described by Wuyts et al. (2010). For the lens MACSJ0138.0-2155, we measured a redshift of 0.338 based on a Magellan/LDSS-3 spectrum of the brightest cluster galaxy.

Zoom In Zoom Out Reset image size

MRG-M0150 is also a multiple image system, presenting three multiple images in a "naked cusp" configuration analyzed by Newman et al. (2015a). MRG-P0918, MRG-S1522, and MRG-M2129 are all singly imaged. The images of MRG-P0918 and MRG-M2129 are both highly elongated, whereas that of MRG-S1522 is bright but only modestly elliptical. The clusters lensing MRG-P0918 and MRG-S1522 are not well-known and, as we will discuss below, we do not have the constraints needed to determine the magnification factor or reconstruct the source in these cases. MACSJ2129.4-0741, on the other hand, is one of the most X-ray-luminous clusters (Ebeling et al. 2007; z = 0.589) and is a well-studied lens with multiple public lens models.

Figure 3 compares the observed (lensed) magnitudes of this sample to unlensed quiescent galaxies in the UltraVISTA field (Muzzin et al. 2013a). All of the sources in our lensed sample are substantially brighter than even the brightest quiescent galaxies located at comparable redshifts, even within the largest deep fields. This makes them premier targets for absorption line spectroscopy in the near-infrared and enables a number of new observations, including their spatially resolved stellar kinematics (Newman et al. 2018, hereafter Paper II) and populations. In the following sections, we will describe follow-up observations to verify the quiescent nature and redshifts z = 1.95–2.64 of these sources (Sections 3 and 4) before deriving the lensing magnifications (Section 5) and the source intrinsic properties (Sections 6–7).

A variety of archival data exists for the lensed galaxies that we located in our HST search. The sources are listed in Section 2.4. For MRG-M0150, we used archival WFC3-IR/F140W, WFC3-IR/F110W, ACS/F814W, and WFPC2/F606W images, as described in Newman et al. (2015a). In order to obtain deeper images with improved sub-pixel sampling and a filter combination that better straddles the Balmer break, we obtained additional WFC3-IR images through the F160W and F125W filters, each at half-orbit depth (Proposal ID 14205, PI: A. Newman). For MRG-S1522, we used WFC3-IR/F160W, WFC3-IR/F105W, and WFC3-UVIS/F606W archival images from the Gladders program (see Section 2.4). For MRG-M2129, we used archival CLASH imaging through 13 ACS and WFC3 filters (see Table 2).

No archival data existed for the two lensed galaxies that were located in our FourStar survey. Therefore, for MRG-P0918 we obtained WFC3-IR/F160W and F105W images, each with half-orbit depth, and ACS/F814W and F555W images, each with one-orbit depth, through the aforementioned program. For MRG-M0138, observations were undertaken through a mid-cycle program (Proposal ID 14496, PI: A. Newman). Two orbits were split between WFC3-IR/F160W (1.6 ks) and F105W (3.6 ks) exposures, and two further orbits were devoted to an F555W exposure.

All of the new HST observations employed standard sub-pixel dither patterns and were reduced using the MultiDrizzle package (Koekemoer et al. 2003). For the source reconstructions that we will perform in Section 5, we require an estimate of the WFC3 F160W point-spread function (PSF). For MRG-M0150 and MRG-M2129, we constructed PSFs from stars in the mosaic. Suitably bright and isolated stars were not present for MRG-M0138, so we instead generated a synthetic TinyTim (Krist et al. 2011) PSF and convolved it with a circular Gaussian to best match a (slightly broadened) star in the mosaic. As expected, the curves of growth for all three PSFs are very similar.

We analyzed archival images from the Infrared Array Camera (IRAC) onboard the Spitzer Space Telescope for three galaxies: MRG-P0918 was observed in Program ID 90233 by PI: C. Lawrence, MRG-S1522 was observed in Program ID 70154 by PI M. Gladders, and MRG-M2129 was observed in Program ID 90009 by PI M. Bradač. We obtained new IRAC images of MRG-M0138 via a Director's Discretionary Time program (Program ID 12127) and of MRG-M0150 via a joint HST program in cycle 23 (Program ID 12003, PI: A. Newman).

The IRAC images cover the 3.6 μm and 4.5 μm channels, with the exception of MRG-P0918. This galaxy falls near the edge of the 4.5 μm mosaic and suffers from various artifacts. Therefore, we excluded this image from our analysis. We performed photometry using the standard calibrated mosaics produced by the IRAC pipeline with a $0\buildrel{\prime\prime}\over{.} 6$ pixel scale.

We imaged the entire sample through various filters using FourStar at the Magellan Baade telescope. For MRG-M0138, MRG-P0918, and MRG-S1522, we obtained images through the J₁, J₂, J₃, H, and K_s filters, using the medium bands to improve the sampling around the Balmer/4000 Å break. For MRG-M0150 and MRG-M2129, due to the greater number of HST/WFC3-IR observations, we obtained FourStar images in J, H, and K_s for the former and only in K_s for the latter.

At optical wavelengths, we observed MRG-M0138 through the g, r, and z filters using the LDSS-3 imaging spectrograph at the Magellan Clay telescope. We also observed MRG-S1522 through the z filter. Photometric calibration was obtained using images of stellar fields in the Sloan Digital Sky Survey.

In each field, we produced pixel registered images with matching PSFs following the procedures described by Newman et al. (2012). Briefly, the PSF in each image was constructed by combining suitable bright stars. The ground-based and HST images were convolved to match the image with the worst seeing, which was always $\leqslant 1^{\prime\prime}$. Since the IRAC PSFs are considerably broader, we did not convolve all images to match the IRAC resolution. The arc flux measured in the IRAC images was instead scaled by the fraction of light lost from the photometric aperture when the F160W image was convolved to match the IRAC PSF (see Newman et al. 2012).

We measured colors in rectangular apertures aligned with the direction of magnification that had the following dimensions: 4'' × 15 on MRG-M0138 Image 3, 16 × 12 on MRG-M0150 Image 1, 2'' × 15 on MRG-P0918, 15 × 15 on MRG-S1522, and 32 × 15 on MRG-M2129. The lengths of these apertures (in the direction of maximum magnification) approximately match the spectroscopic apertures described in Section 6.1 and so are appropriate for a joint analysis of the spectra and photometry.⁵ The rectangular apertures are wider (15) than the spectrograph slit (06–075), because we preferred to avoid selecting an aperture much smaller than the IRAC PSF. However, this does not introduce a significant mismatch between the spectroscopic and photometric apertures because the images are narrow in the direction of minimum magnification.

For each galaxy, the fluxes were then uniformly scaled to match the total F160W flux within a large aperture and corrected for Galactic extinction following Schlafly & Finkbeiner (2011). Uncertainties of 5% (3% for HST measurements) were added in quadrature to account for uncertainties in the PSF matching and photometric calibration. The resulting photometric measurements are listed in Table 2.

For MRG-M0138, colors were measured on Image 3, which is the most isolated and affords the cleanest photometry. In Table 2, we have scaled these fluxes to match the total observed F160W flux of the brightest image and our spectroscopic target, which is Image 1. (Although there is some ambiguity in separating the merging Images 1 and 2, this only affects the total fluxes in Table 2 and has no consequence for any of the inferred source properties, which will ultimately be scaled based on a model of the source.)

We generally used the 075-wide slit, which provides a spectral resolution of σ_inst = 33 km s⁻¹. For a portion of the MRG-M2129 observations that were conducted in excellent seeing, we used the 060-wide slit. To minimize read noise, we operated the detector in the up-the-ramp sampling mode for integration times of 20–30 minutes.

During twilight, we obtained sky exposures to measure the illumination of the slit. Although some of the target arcs could be acquired directly on the slit-viewing acquisition camera, we usually acquired a nearby star and then offset the telescope to the target. After the offset, we compared the position of the offset star and other sources to their expected pixel coordinates in images taken with the acquisition camera. In some cases, these differed by up to 04, and we offset the telescope to place the sources at the expected positions. We monitored these positions throughout the exposure sequence and corrected gradual drift when it occurred.

The 6''-long FIRE slit was oriented in the direction of elongation of the target lensed galaxy (see P.A. in Table 3). Observations were made in an AB pattern with short dithers of 08–25 depending on the angular size of the target. A simple A-B subtraction is not feasible due to the resulting overlap of the extended target in the two dither positions and also to the long exposures needed to minimize the read noise. The spectra of the sky and target must instead be modeled in each exposure, as we describe below. Exposures of the internal quartz and ThAr lamps were interspersed throughout the FIRE observations and used for flat fielding and wavelength calibrations. To remove telluric absorption, A0V stars were observed both before and after each target and usually in the middle of longer exposure sequences.

Table 3.  Spectroscopic Observing Log

Target	Instrument	Dates	Exposure Time (hr)	Seeing	Slit PA (deg)
MRG-M0138 Image 1	MOSFIRE J	2015 Nov 3	1.0	073	82
MRG-M0138 Image 1	MOSFIRE H	2015 Nov 3, 6	2.0	079	82
MRG-M0138 Image 2	FIRE	2016 Sep 8–9	5.5	048	−60
MRG-M0150 Image 1	MOSFIRE H	2014 Nov 26–27	4.3	067	−7.5
MRG-M0150 Image 1	FIRE	2014 Nov 1, 3	6.5	055	−7.5
MRG-P0918	FIRE	2014 Apr 14–15	7.0	042	176
MRG-S1522	FIRE	2014 Feb 28–29,	9.0	057	136
		2014 Apr 13–15
MRG-M2129	FIRE	2015 Sep 27–28,	16.7	049	−13
		2016 Sep 8–9

Download table as:  ASCIITypeset image

The seeing is an important ingredient in our dynamical modeling. To measure the seeing, we monitored stars on the acquisition camera and on the facility guide camera. The seeing was also estimated from the science spectra themselves through a comparison with HST images (see Section 4.2). By comparing these methods, we estimate that the uncertainty in the seeing is ≲01. The mean seeing during the observations of each galaxy is listed in Table 3 and ranges from 042 to 057.

The FIREHOSE pipeline⁶ was used to flat field the data and to provide an initial wavelength solution and initial rectification of each spectral order. Since FIREHOSE was designed primarily for the reduction of point sources, in the subsequent steps we relied on custom IDL routines that were based on the FIREHOSE code.

First, traces of bright stars were used to make small corrections to the rectifications in each order. For each science exposure, we then masked all orders to isolate the inter-order background. We fit and subtracted a smooth variation with column within each of the four amplifier regions, which was necessary to remove discontinuities at the boundaries. We then modeled and subtracted a smooth scattered light background. Low-order corrections to the initial FIREHOSE wavelength solution, which was derived from ThAr lamp exposures, were then made using the OH lines. Once these low-order corrections were derived for one science exposure in a series, cross-correlations were used to correct for instrumental flexure.

Each spectral order of each exposure was then modeled as the sum of sky and galaxy emission using iterative bspline techniques as employed by FIREHOSE (see, e.g., Kelson 2003). This method requires that the spatial distribution of the emission be specified. For the sky background, we found that the residual intensity variation along the slit can be modeled as a quadratic polynomial. For the targeted galaxy, we first measured its spatial profile in the H band using the initial FIREHOSE sky model. However, this estimate is imprecise due to the difficulty of separating an extended source from the sky background over a short slit with no prior information on the source structure. Therefore, we used the HST WFC3/F160W image to measure the expected galaxy flux profile along the slit, taking into account its width and orientation. This HST-based flux profile was then shifted and convolved by a Gaussian PSF to best match the profile in the spectrum. This procedure provides an estimate of the seeing and the position of the target in each exposure, and it produces a galaxy profile suitable for accurate sky modeling.

Observations of A0V stars were analyzed with the xtellcorr package (Vacca et al. 2003), as implemented within FIREHOSE, to provide flux calibration and removal of telluric absorption. While the A0V observations are needed to track the temporal variation of telluric features, we found that the relative flux calibration could be improved through observations of white dwarf standards. These were reduced in the same way as the galaxy spectra and were then used to derive low-order corrections to the flux calibration in each order. Additionally, in a few cases, we made low-order corrections to the continuum shape based on comparisons to the fitted stellar population synthesis models (Section 6.2).

For each science exposure, two-dimensional rectified spectra were then produced from each order. The exposures were then normalized to a common flux level, measured in H band. All of the exposures in each order were then spatially registered and averaged using inverse variance weighting. Residual outlier pixels were then identified and interpolated over. One-dimensional spectra were then extracted in each order in a specified aperture. Small multiplicative offsets were applied to ensure that spectra extracted in adjacent orders have consistent fluxes within the wavelength regions of overlap. Finally, the orders were combined into a single one-dimensional spectrum with a scale of 12.5 km s⁻¹ pixel⁻¹.

We observed two targets with MOSFIRE: Image 1 of MRG-M0150 in H band and Image 1 of MRG-M0138 in both J and H bands. We formed a 07-wide long slit on each target and on several stars in the field, which were used to align the mask and to measure the seeing. The observations were conducted with an AB dither pattern. Since subtraction of consecutive dithered exposures is fundamental to the MOSFIRE Data Reduction Pipeline (DRP), we took care to ensure that the dither distance ($5\buildrel{\prime\prime}\over{.} 4$ for MRG-M0150 and $9\buildrel{\prime\prime}\over{.} 2$ for MRG-M0138) was sufficiently wide to avoid self-subtraction of the extended arcs.

The observations and reduction of MRG-M0150 were described by Newman et al. (2015a). To reduce the MRG-M0138 observations, we used the DRP to produce coadded 2D spectra for each night of observations (see Table 3). Relative flux calibration was then performed using twilight observations of the white dwarf GD71. Given the very high signal-to-noise ratio of the MRG-M0138 spectrum, we performed telluric absorption corrections differently from the other observations. In the H band, we iteratively modeled the galaxy stellar continuum and the telluric absorption using the radiative transfer code molecfit (Smette et al. 2015). A similar procedure was described by Newman et al. (2017). In the J band, where the signal-to-noise ratio is lower, we instead used molecfit to model the telluric absorption in observations of GD71. We then divided this synthetic absorption spectrum from the galaxy observations. Finally, observations made on different nights were registered, scaled to a common flux level, and averaged with inverse variance weighting.

5.1.1. Lens and Source Modeling Methods

We used the ray tracing code introduced by Newman et al. (2015b) and employed by Newman et al. (2015a) to fit a simply parameterized model of the lens mass distribution and the source light distribution to the pixel-level data. The regions of the WFC3/F160W image that were used to constrain the model are outlined in blue in Figure 4(a). The BCG light has been modeled and subtracted in this image.

Zoom In Zoom Out Reset image size

The lensing cluster was modeled as a superposition of dual pseudo-isothermal elliptical (dPIE) mass distributions (see Elíasdóttir et al.'s 2007 appendix). Each dPIE was described by a center, PA, ellipticity, two characteristic radii r_core and r_cut, and a normalization σ₀. Since the lensing cluster appears to be a simple, relaxed system, we found that a single dPIE component is adequate to model the cluster dark matter halo. We left all of its parameters free except for r_cut, which lies well beyond the strong lensing zone and can be fixed to 1 Mpc.

We modeled the stellar mass in the BCG by fitting a dPIE profile to the surface brightness distribution to set r_core and r_cut. We used Gaussian priors of PA = 46° ± 10° and b/a = 0.42 ± 0.05 that were informed by the photometry. The normalization σ_0,BCG was allowed to vary freely.

Other cluster galaxies were generally included in the mass model using scaling relations with luminosity: ${\sigma }_{0}={\sigma }_{0}^{* }{(L/{L}_{* })}^{1/4}$ and ${r}_{\mathrm{cut}}={r}_{\mathrm{cut}}^{* }{(L/{L}_{* })}^{1/2}$, where ${\sigma }_{0}^{* }$ and ${r}_{\mathrm{cut}}^{* }$ are free parameters with observationally motivated priors (see Newman et al. 2013). This approach ties the center, PA, and ellipticity of the mass distribution to the galaxy light and imposes a constant mass-to-light ratio. The radial distribution of the mass, encoded by r_cut, can differ from the light. Similar approaches have been widely used in other parameterized lens models (e.g., Jullo et al. 2007; Richard et al. 2010; Newman et al. 2013).

Three galaxies either significantly perturb the critical lines or are likely to deviate from these general scaling relations and, therefore, were modeled independently of them. The galaxy labeled P1 in Figure 4(a) is located very close to Image 1, and it significantly alters the critical line and increases the magnification. This perturbing galaxy is faint and not easily visible in Figure 4, but it can be seen in Figure 2. We allowed all of the parameters describing the P1 mass distribution to vary freely. Galaxies P2 and P3 are blue cluster galaxies that show spectacular examples of ram pressure stripping in Figure 2. Given their blue colors, we do not expect them to share a common mass-to-light ratio with the other cluster members. We allowed σ₀ to vary freely for both P2 and P3, as well as r_cut for P2. The final ingredient in the lens model is an external shear that is uniform across the image, which we ultimately found to be small, $| \gamma | =0.05$.

We initially modeled the source light distribution using a single elliptical Sérsic profile (although, as described below, we ultimately adopted a two-component model). We denote the effective radius as R_e,maj to emphasize that it is the semimajor axis of the ellipse containing half of the light. For a given set of parameters that describe the lens model and source, we first cast the centroids of Images 1–4 back to the source plane and then computed the dispersion in their positions. Models with a large dispersion were immediately discarded. For the remaining models, we traced pixels in the image plane back to the source to compute the surface brightness. Care was taken to integrate the surface brightness within the pixels near the center of the source. The model image plane was then convolved by the PSF and compared to the data. A Markov Chain Monte Carlo (MCMC) sampler, MultiNest (Feroz et al. 2009), was used to explore the 29-dimensional parameter space.

Although this model was able to reproduce the overall lens configuration, we found that it could not match the internal structure of all of the images in detail. The clearest example of this can be found in Image 2, which we show in Figure 5. This image presents a flattened inner component that is misaligned from the overall direction of the arc (upper panel). A model with a single-component Sérsic source (lower panel) is not able to reproduce this. However, this deficiency can be eliminated by introducing a two-component source: this model matches the morphology of the image much better in the inner parts (middle panel). Specifically, we modeled the source as the sum of two Sérsic components that share a common center and PA, but which have different axis ratios b/a, effective radii R_e,maj, Sérsic indices n, and magnitudes.

Zoom In Zoom Out Reset image size

To assess the fit quality, we compare the observed image to the model image plane in panels (c)–(f) of Figure 4. Although some areas of mismatch are visible, the model contours (blue) generally follow the data (orange) well. (The contours near the "foreground P1" label include flux from the foreground galaxy P1, which was masked during the fit.) We note that the radial images 4 and 5 were not used to constrain the model except via their approximate positions, so a close match to their detailed structure is not expected.

The magnification factors⁷ for images 1, 2 and 3 are μ₁ = 12.5 ± 5.4, μ₂ = 10.3 ± 3.1, and μ₃ = 4.9 ± 1.6, whose uncertainties are described in the next section.

5.1.2. Magnification Uncertainties

5.1.3. Reconstructing the Source Plane

Our reconstruction of MRG-M0138 in the source plane is shown in Figure 6. Given the nonlinear nature of the lens mapping, convolution by the PSF in the image plane can have very complex effects in the source plane. This makes it difficult to compare the consistency of reconstructions derived from different images. To address this issue, we do not directly cast the observed pixels back to the source plane. Instead, we take the unconvolved image plane model, add the residuals, and then cast this image (i.e., Data − Model_convolved + Model_unconvolved) back to the source plane. This technique effectively deconvolves the inner regions of the images under the assumption of a particular source model, while still allowing for deviations from the model at larger radii (see Szomoru et al. 2012 for an application in a non-lensing context). All three images consistently show a highly flattened, disk-dominated source. We note that since Images 1 and 2 merge into a giant arc, the southwestern portion of the source is not present in these images.

Zoom In Zoom Out Reset image size

In Figure 6, we also compare the surface brightness profiles derived from these three source plane reconstructions. The three multiple images reproduce a single surface brightness profile (solid colored lines) with impressive consistency. This profile is well fit by the two-component Sérsic model (solid black line), leaving azimuthally averaged residuals of ≲10% out to R = 10 kpc. The leftward arrow indicates the radius at which the effects of PSF convolution are significant; at smaller radii, the shape of the plotted profile is dominated by the assumption of the double Sérsic form.

5.1.4. Source Structure

In the best-fitting two-component model, 74% of the stellar light is in an extended (R_e,maj = 7.1 kpc), flattened (b/a = 0.19), disk-like component with a Sérsic index of n = 1.3. A fainter (26% of total flux), rounder (b/a = 0.67), and much more compact (R_e,maj = 0.8 kpc) component is also present. It has a Sérsic index n = 1, the smallest value allowed in the fits, which could indicate a structure analogous to a pseudobulge. However, we caution that the images are significantly affected by the PSF within the effective radius of the bulge-like component (Figure 6).

Table 4 lists the Sérsic parameters and uncertainties. The uncertainties were derived from the image-to-image scatter: specifically, we fixed the lens mapping, fit the source to each of the three images individually, and then measured the standard deviation of each parameter. For the total flux and R_e, we additionally list the systematic lens model uncertainties discussed in Section 5.1.2.

Table 4.  Sérsic Model Parameters

Quantity	Units	Single Sérsic Model	Double Sérsic Model
MRG-M0138
mF160W	mag	20.0 ± 0.1 ± 0.4	20.4 ± 0.1 ± 0.4	21.5 ± 0.1 ± 0.4
Re,maj	arcsec	057 ± 007 ± 012	085 ± 007 ± 018	010 ± 002 ± 002
Re,maj	kpc	4.8 ± 0.6 ± 1.0	7.1 ± 0.6 ± 1.5	0.8 ± 0.2 ± 0.2
n	...	2.9 ± 0.7	1.3 ± 0.3	1a
b/a	...	0.26 ± 0.04	0.19 ± 0.03	0.67 ± 0.16
PA	deg	36 ± 2	36 ± 1
MRG-M0150
mF160W	mag	21.3 ± 0.1 ± 0.4
Re,maj	arcsec	021 ± 003 ± 004
Re,maj	kpc	1.7 ± 0.2 ± 0.3
n	...	3.5 ± 0.4
b/a	...	0.87 ± 0.05
PA	deg	−6 ± 10
MRG-M2129
mF160W	mag	21.8 ± 0.2	22.0 ± 0.2	23.4 ± 0.1
Re,maj	arcsec	029 ± 002	027 ± 0.02	047 ± 0.05
Re,maj	kpc	2.4 ± 0.2	2.2 ± 0.2	3.9 ± 0.4
n	...	1a	1a	1a
b/a	...	0.29 ± 0.03	0.24 ± 0.03	0.80 ± 0.07
PA	deg	−39 ± 5	−39 ± 5
Gaussian mF160W	mag	25.9 ± 0.5	26.1 ± 0.3
Gaussian Re	arcsec	<0007 (95%)	<0008 (95%)

Notes. For the extrinsic parameters (magnitudes and R_e) of the multiply imaged systems, we first list the uncertainty derived from the scatter among the multiple images (Section 5.1.4) and then list the systematic uncertainty arising from the overall magnification uncertainties (Section 5.1.2).

^aEdge of prior.

Download table as:  ASCIITypeset image

In Newman et al. (2015a) we constructed a parametrized model of the lensing cluster and source that was constrained by the archival HST WFC3/F140W image of MACSJ0150.03–1005 that was then available. This technique was very similar to that described in the previous subsection for MRG-M0138, and its application to MRG-M0150 was detailed by Newman et al. (2015a). Subsequently, we obtained a deeper WFC3 image at the slightly redder wavelength of F160W with improved sub-pixel sampling (Section 3.1). Here we briefly describe our improved analysis using this new image.

Figure 7 shows the three multiple images 1–3 (panel (a)), the single-component Sérsic source traced through the lensing potential and convolved by the PSF (panel (b)), and zooms on the image regions to compare the model contours and residuals (panels (c)–(f)). The model successfully reproduces the structure of the images in detail. The magnifications of the images are μ₁ = 4.4 ± 1.1, μ₂ = 2.6 ± 1.0, and μ₃ = 4.6 ± 1.3. These uncertainties are estimated by varying the parameterization of the mass model and then comparing them to the results from an independent set of lens models constructed with the Lenstool code, as described in Section 5.1.2 and Newman et al. (2015a).

Zoom In Zoom Out Reset image size

Source reconstructions and surface brightness profiles from each multiple image are shown in Figure 6.⁸ The lens model produces three consistent reconstructions of the source. Likewise, the surface brightness profiles of the source from the observed multiple images are mutually consistent at the ≃5% level out to 4R_e.

The Sérsic parameters are listed in Table 4. The source is a compact galaxy (R_e,maj = 1.7 kpc) with a nearly de Vaucoulers profile (n = 3.5) and nearly round isophotes (b/a = 0.86). Unlike MRG-M0138, we find that a single Sérsic component is adequate to describe the source, but we note that it would be much more difficult to discern the presence of multiple components in this case since MRG-M0150 is nearly round in projection. These parameters are consistent with those measured using the shallower F140W image by Newman et al. (2015a) but the uncertainties are reduced using the new deeper data.

MRG-M2129 is singly imaged, so unlike MRG-M0138 and MRG-M0150, we require other multiple images to constrain the source structure. The lensing cluster MACSJ2129.4–0741 presents a large number of multiple images that are available to constrain the mass distribution, and several authors have published models. We used the deflection angle maps produced by several published lens models to trace a model source through the lens mapping. The model image plane can then be convolved by the PSF and fit to the HST/WFC3 F160W image. Since this system was discovered by Geier et al. (2013) and studied by Toft et al. (2017), we also compare our results to theirs in the Appendix.

We considered three lens models. Monna et al. (2017) kindly supplied maps of the deflection angle for their parametric model. We also obtained two mass models constructed by A. Zitrin (Zitrin et al. 2009, 2013) as high-level science products of the CLASH program (Postman et al. 2012). We found that the Monna et al. (2017) lens model produced the lowest ${\chi }^{2}$. Therefore, we used the parameters derived from it and estimated the uncertainties in the galaxy structural parameters by comparing to those obtained using the two Zitrin et al. lens models.

Figure 8 shows the region of the HST/WFC3 F160W image that we fit. We first modeled and subtracted a foreground cluster member that is located to the west of the arc. To avoid any residual contamination from this galaxy, the fit region is not centered on the arc. We then subtracted the sky background measured in an empty region to the east of the arc. In this case, the lensed galaxy is sufficiently distant from the BCG so that it was unnecessary to model the latter. A second source is visible in the upper right corner, which is disconnected from the main galaxy (see also Geier et al. 2013). We masked this object, which is potentially a satellite galaxy.

Zoom In Zoom Out Reset image size

We found that a single Sérsic model is inadequate to model MRG-M2129, both in the central and outer regions. First, there is a centrally located, point-like source visible in the top row of Figure 8. We modeled this by adding a second component to source model: a circular Gaussian at the galaxy center with a free effective radius R_e and flux. The inclusion of the central source improved the fit by Δχ² = 759 and, therefore, is clearly justified. Its magnitude is m_F160W = 25.9 ± 0.5 (demagnified), or ∼1% of the galaxy flux, and it is very compact: ${R}_{e}\lt 0\buildrel{\prime\prime}\over{.} 007=60\,\mathrm{pc}$ (95% confidence). The compact central source could be continuum emission from a Seyfert nucleus, which is supported by the emission line ratios (Section 7.2). The properties of the single Sérsic+Gaussian source model are listed in Table 4.

Second, we found that a single-component Sérsic model does not reproduce the shape of the outer isophotes well, which is evident in the upper-left and lower-right regions of the images in Figure 8. The addition of a second Sérsic component with a center and PA tied to those of the first component improved the fit. This can be seen qualitatively in the right panels of Figure 8 and is supported quantitatively by an improvement of Δχ² = 1287. Like MRG-M0138, MRG-M2129 shows evidence for two components with different ellipticities.

The source model consists of: (1) the central source described above, (2) a flattened (q = 0.24) component containing 79% of the flux in a compact (R_e,maj = 027 = 2.2 kpc) exponential disk (n = 1), and (3) a rounder (q = 0.80) component containing 21% of the flux that is more extended (R_e,maj = 047 = 3.9 kpc) and also nearly exponential (n = 1). The parameters are summarized in Table 4. We restricted the Sérsic index to values $n\gt 1$ in the fits. Since both components hit this limit, their Sérsic indices should be regarded with caution. However, the presence of two components—one compact and flattened, one fainter, rounder, and more extended—seems robust and holds for all three lens models. The uncertainties in Table 4 are derived from the standard deviation of values obtained using the three lens models.

The structures of MRG-M0138 and MRG-M2129 are both dominated by a highly flattened exponential component, i.e., a nearly edge-on disk. Both galaxies also contain additional components. In MRG-M0138, we find a compact rounder component emitting 26% of the luminosity, which is potentially a nascent bulge that is embedded in a very massive and extended (R_e = 7 kpc) disk. In MRG-M2129 we find an extended rounder component and, remarkably, an unresolved central point source that could be Seyfert nucleus. This is particularly interesting given the evidence for an active galactic nucleus (AGN) from the emission line ratios, which we will discuss in Section 7.2. MRG-M0150 is instead well described by a single compact component with a nearly de Vaucouleurs profile (n = 3.5) and round isophotes. With imaging alone, we cannot tell whether MRG-M0150 is intrinsically rounder than the other two systems or is merely less inclined, but we will address this question using stellar kinematics in Paper II.

The disk-dominated structures of MRG-M0138 and MRG-M2129 are quite distinct from massive early-type galaxies in the local universe. In Section 8.2, we will compare our sample to local systems and will discuss the implications for their future evolution.

To define the extraction region for each target, we examined the flux distribution along the slit. The extraction box was chosen to cover the region where the flux profile was ≳0.15 times the peak intensity.⁹ As discussed in Section 3.4, these approximately match the apertures in which the colors were measured. For our analysis, we consider the wavelength range from λ_rest = 3600 Å to ${\lambda }_{\mathrm{obs}}\simeq 2.3$ μm. The spectrum outside the range is generally too noisy to be useful. Only for the highest-redshift source, MRG-M0150, did we reduce the spectrum out to λ_obs = 2.45 μm to include Hα and [N ii].

The spectra are shown in Figure 9, where they have been rebinned for display purposes, as indicated in the caption. The quality is remarkable considering the redshifts of these sources. The median signal-to-noise ratio in the H band generally ranges from 21–32 per 300 km s⁻¹ bin, which is approximately one velocity dispersion element. For the ultra-bright MRG-M0138, this figure reaches 137 and 77 for the MOSFIRE and FIRE spectra of Images  1 and 2, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We jointly modeled the spectra and photometry using the pyspecfit code described by Newman et al. (2014). We used the Bruzual & Charlot (2003, BC03) population synthesis models and adopted an exponentially declining star formation history $\mathrm{SFR}\propto \exp (-t/\tau )$, the Calzetti et al. (2000) dust attenuation curve, and the Chabrier (2003) initial mass function. For each galaxy we fit the redshift, velocity dispersion, age,¹⁰ τ, metallicity Z, A_V, and the emission line parameters that are described below. As described in Newman et al. (2014), when comparing the data to a given model, we warped the spectrum continuum shape by the polynomial that minimizes the ${\chi }^{2}$. This accounts for flux calibration errors in the spectrum and anchors its continuum shape to the photometry. We allowed a separate polynomial of order ≃Δλ_rest/200 Å in each of the J, H, and K bands, where Δλ_rest is the unmasked rest-frame wavelength interval.

Most of the galaxies in the sample exhibit weak emission lines whose measurement requires accurate modeling of the stellar continuum. Therefore, we fit both simultaneously. Emission lines from [O ii], [O iii], [N ii], [S ii] and the Balmer series were modeled as Gaussians with a velocity and velocity dispersion ${\sigma }_{\mathrm{em}}$ that are common to all emission lines but are distinct from those of the stellar component. For each galaxy, we modeled only the lines falling in the atmospheric transparency windows. The intensity ratios [O iii] $\lambda 5008/\lambda 4960$ and [N ii] $\lambda 6585$/$\lambda 6550$ were fixed to 2.98 and 3.05 (Storey & Zeippen 2000), while for the Balmer series we fixed the relative intensities assuming Case B recombination. Since Hβ is always in net absorption, we cannot separately constrain the emission line attenuation via the Balmer decrement. Therefore, once the emission line spectrum was constructed based on the aforementioned ratios, we attenuated it at each wavelength by the same factor as the starlight. We will consider the effects of possible differential extinction in our interpretation; these are expected to be small because the continuum attenuation is mild.

The posterior distributions were sampled using MultiNest. In addition to the model parameters, we also computed derived quantities, such as the specific star formation rate (sSFR) and emission line ratios. The priors were broad and uninformative, with two exceptions: we restricted the stellar metallicity to Z = 0.01–0.05 (the solar value is 0.02 in the BC03 models) and placed a Gaussian prior on ${v}_{\mathrm{gas}}-{v}_{\mathrm{stars}}$ with a mean of 0 and a dispersion of 200 km s⁻¹.¹¹

The posterior constraints on the stellar populations and emission lines are listed in Table 5. As Figure 9 shows, these models generally fit the photometry and spectra well. For the photometry, the reduced ${\chi }_{\mathrm{phot}}^{2}$ is in the range 0.4–1.5. Since deep NIR spectra often exhibit noise somewhat above the formal estimate, we rescaled the spectral uncertainties by a constant factor in the range 1.0–2.5 so that the reduced ${\chi }_{\mathrm{spec}}^{2}\simeq 1$. We masked Mg b since the galaxies may have non-solar abundance ratios. For MRG-M0138, we also masked Na D, which clearly has a non-solar abundance or is affected by interstellar absorption. For the MOSFIRE spectrum of MRG-M0138, we also masked the region around the G band and Hγ due to imperfect correction of telluric absorption in this region.

Table 5.  Spectroscopic Measurements

Quantity	Units	MRG-M0150	MRG-P0918	MRG-S1522	MRG-M2129	MRG-M0138	MRG-M0138
						Image 1	Image 2
Stellar Population Properties
z	...	2.6355	2.3559	2.4503	2.1487	1.9486	1.9469
σ	km s−1	261 ± 30	223 ± 16	241 ± 18	266 ± 21	298 ± 7	409 ± 11
Age	Gyr	0.76 ± 0.08	0.51 ± 0.02	0.61 ± 0.06	0.80 ± 0.10	1.35 ± 0.08	1.39 ± 0.16
τ	Myr	95 ± 35	<43	<71	103 ± 24	181 ± 17	178 ± 27
AV	mag	0.61 ± 0.09	0.18 ± 0.06	0.34 ± 0.07	0.33 ± 0.09	0.35 ± 0.05	0.11 ± 0.05
[Z/H]	...	<0.33	0.02 ± 0.03	−0.03 ± 0.13	0.16 ± 0.13	0.01 ± 0.04	0.25 ± 0.09
log sSFR	yr−1	<−10.72	<−12.80	<−11.77	−11.18 ± 0.54	−11.28 ± 0.11	−11.44 ± 0.22
$\mathrm{log}\mu {M}_{* }$	M⊙	12.06 ± 0.04	11.72 ± 0.02	11.74 ± 0.03	11.62 ± 0.05	12.77 ± 0.03	12.56 ± 0.04
μ SFR	M⊙ yr−1	<22	<1	<1	3 ± 2	31 ± 9	13 ± 3
$\mathrm{log}{M}_{* }$	M⊙	11.50 ± 0.17	...	...	10.96 ± 0.10	11.69 ± 0.19	11.68 ± 0.19
SFR	M⊙ yr−1	<6.1	...	...	0.6 ± 0.4	2.6 ± 1.4	1.8 ± 0.9
Emission Line Properties
σem	km s−1	213 ± 20	190 ± 49	345 ± 55	364 ± 22	...	...
vem − vstars	km s−1	−20 ± 24	27 ± 37	40 ± 41	−12 ± 30	...	...
Hα EW	Å	4.3 ± 1.2	0.3 ± 0.5	4.2 ± 0.9	2.6 ± 0.9	...	...
[N ii] λ6585 EW	Å	24.5 ± 2.3	3.7 ± 0.6	11.1 ± 1.8	17.1 ± 1.1	...	...
[O iii] λ5008 EW	Å	...	0.1 ± 0.2	2.0 ± 0.3	5.3 ± 0.5	...	...
[O ii] EW	Å	...	<3.5a	5.1 ± 0.9	2.5 ± 1.3	...	−0.1 ± 0.6
[S ii] EW	Å	...	...	...	<3.6a	...	...
log [N ii]/Hα	...	0.74 ± 0.12	>0.33	0.40 ± 0.11	0.79 ± 0.16	...	...
log [O iii]/Hβ	...	...	...	${0.28}_{-0.11}^{+0.32}$	${0.86}_{-0.15}^{+0.34}$	...	...
Line Fluxes Uncorrected for Magnification or Extinction:
Hα flux	10−18 cgs	73 ± 20	5 ± 7	51 ± 11	30 ± 10	...	...
[N ii] λ6585 flux	10−18 cgs	391 ± 38	54 ± 9	129 ± 21	183 ± 13	...	...
[O iii] λ5008 flux	10−18 cgs	...	2 ± 4	31 ± 5	68 ± 6	...	...
[O ii] flux	10−18 cgs	...	<44a	42 ± 7	17 ± 9	...	−3 ± 27
[S ii] flux	10−18 cgs	...	...	...	<39a	...	...

Notes. BC03 models and a Chabrier (2003) IMF are assumed. Uncertainties are statistical and do not include model systematics (see text). All upper limits are 2σ. The flux magnification factor is μ; hence, $\mathrm{log}{M}_{* }$ is demagnified and $\mathrm{log}\mu {M}_{* }$ is not. EWs are reported in the rest frame. Line flux units are 10⁻¹⁸ erg cm⁻² s⁻¹ Å⁻¹. The [O ii]λλ3727,3730 and [S ii]λλ6718,6733 fluxes and EWs refer to the sum of the two doublet components. The spectra and derived extrinsic quantities, including line fluxes, are scaled to match the total photometric magnitudes.

^aThese lines were formally detected at marginal significance, but since inspection of the spectra showed these may be spurious, we have chosen to quote only upper limits.

Download table as:  ASCIITypeset image

Consistent with the goal of our color–magnitude selection, our modeling of the stellar continuum indicates high stellar masses and low sSFRs for all five lensed galaxies in the sample. Their stellar masses, uncorrected for magnification, span the range $\mu {M}_{* }={10}^{11.6-12.8}\,{M}_{\odot }$. For the three systems with estimated magnification factors μ, we find ${M}_{* }\,={10}^{11.0-11.7}\,{M}_{\odot }$. The specific star formation rates (sSFRs), which are independent of magnification, are ${10}^{-10.7}$ yr⁻¹ or smaller. As we will discuss in Section 7, emission lines (when present) do not have ratios indicative of star formation.

Inspection of the spectra immediately shows a diversity of ages from 0.5–1.4 Gyr. Coupled with the low current sSFRs, this implies a substantial decline in star formation rate over the past Gyr with e-folding times of τ ≲ 100–200 Myr. A mild amount of reddening, corresponding to A_V = 0.1–0.6, is inferred for the full sample. We find stellar abundances near solar. For the most massive and oldest galaxy, MRG-M0138, the FIRE spectrum indicates $Z\simeq 2{Z}_{\odot }$. We will analyze the detailed chemical abundance pattern of this galaxy in a forthcoming paper.

Although the stellar continuum suggests only mild reddening, mid- and far-infrared observations are needed to address the possibility of highly extinguished star-forming regions. Stacking analyses have shown that quiescent galaxies at z < 2.5 identified by the UVJ diagram or similar criteria usually do not harbor much obscured star formation (Fumagalli et al. 2014; Man et al. 2016). For the galaxies in our sample, the only existing observation with adequate depth to address this question is a Spitzer MIPS 24  μm image that covers MRG-M2129 (PI: M. Yun, Program ID 50610). MRG-M2129 is not detected with a 2σ upper limit of 86 μJy, corresponding to μL_IR < 8 × 10¹¹ L_⊙ assuming the Wuyts et al. (2008) template. In the limit where the dust is heated only by star formation, this limit implies $\mathrm{SFR}\lt 18\,{M}_{\odot }$ yr⁻¹ (corrected for magnification).¹² This would still place the galaxy below the main sequence of star-forming galaxies, given its high mass; furthermore, this limit is conservative because the emission lines in MRG-M2129 show evidence for an active galactic nucleus (AGN; see Section 7.2) that could contribute to ${L}_{\mathrm{IR}}$.

The formal uncertainties in Table 5 are usually very small, which is expected given the high quality of the data and the simplicity of the star formation history and other aspects of the model. We emphasize that the listed uncertainties are purely statistical. Systematic errors in the models certainly dominate (e.g., Muzzin et al. 2009). As a rough estimate of these systematic errors, we also analyzed the data using the FSPS v3.0 (Conroy et al. 2009; Conroy & Gunn 2010) stellar population synthesis models using the MIST isochrones (Choi et al. 2016). We find systematic differences of ${t}_{\mathrm{FSPS}}-{t}_{\mathrm{BC}03}\,=+0.2\,\mathrm{Gyr}$ in age, +50 Myr in τ, $+0.1$ mag in A_V, $+0.2$ dex in sSFR, and $+0.14$ dex in M_*, and $+0.3$ Å in emission line EWs. In [Z/H] the differences are not systematic, but there is a scatter of 0.2 dex, indicating that the metallicity constraints are the least robust.

For MRG-M0138, we have analyzed spectra of both Image  1 and 2, which provides a consistency test. The same colors were used in both fits, but we rescaled the flux level to match Image 2 when fitting the FIRE spectrum. Because MRG-M0138 is rotating (Paper II) and the slit cuts through the source at very different angles with respect to the major axis for Images  1 and 2, we do not expect to measure the same velocity dispersion σ in the two spectra. However, the stellar population parameters are reasonably consistent.

It is clear that the stellar absorption correction significantly affects the inferred Hα emission. Are the high inferred [N ii]/Hα ratios reliable? To test this, we fit only the [N ii]+Hα region of the spectrum with the continuum fixed to a simple stellar population with an age of 400 Myr. Since the Hα absorption is maximum around this age, the test supplies an upper limit to the Hα flux and thus a lower limit to [N ii]/Hα. Due to the weak sensitivity of the Hα absorption to age in the relevant range of ages, [N ii]/Hα never shifted by more than half of its random uncertainty in this test. Systematic shifts in the line ratios that arise when using the FSPS versus BC03 stellar models are even smaller (≃0.02 dex), with the exception that the lower limit on log [N ii]/Hα is weakened from >0.54 to >0.33 for MRG-P0918 when using the FSPS models. We conservatively adopted this weaker constraint in Table 5. In all other cases, the [N ii]/Hα ratio is robust to the stellar continuum modeling details.

The [O iii]/Hβ ratio is inferred less directly. Since the Hβ and higher-order emission is much weaker than the stellar absorption, the intensity of the Balmer emission lines in the fit is driven by Hα infilling, and we cannot constrain the Balmer decrement Hβ/Hα. Our spectrum model assumes that the nebular lines and the stars are equally attenuated. If the line emission is more extinguished, then this would push the inferred [O iii]/Hβ upward. Reddy et al. (2015) found that the attenuation of the Hα photons can exceed that of the stellar continuum by up to 1.5 mag in high-z star-forming galaxies. We think that such a large difference is unlikely to hold in quiescent galaxies where attenuation is mild, but assuming 1.5 mag of differential attenuation as a limiting case, we would infer [O iii]/Hβ typically 0.3 dex higher. The asymmetric error bars in Figure 11 and Table 5 include this uncertainty added in quadrature.

Determining the excitation mechanisms in the Seyfert/LIER¹⁴ region of the BPT diagram is difficult, even in local galaxies with much more detailed information. Although we cannot expect to definitively identify the excitation mechanisms in every case, we can consider which scenarios are consistent with the available evidence. The ratio [N ii]/Hα ≃2–6 is high throughout our sample, but the strength of the oxygen lines and the spatial distribution and kinematics of the gas relative to the stars provide additional constraints and show more diversity.

The spatial variation of emission line EWs is shown in Figure 12. These measurements were made from spectra extracted in bins along each image. (These bins are the same that are used to extract stellar kinematics in Paper II, where further details can be found.) We then modeled each spectrum as described in Section 6.2, except that we did not include photometric constraints, which are not necessary to measure the emission line EWs.

Zoom In Zoom Out Reset image size

Photoionization by an AGN is one way to produce a high [N ii]/Hα ratio. In this scenario, we expect high [O iii]/Hβ and [O iii]/[O ii] ratios and a negative gradient in the emission line EWs (i.e., centrally concentrated emission). MRG-M2129 is the only galaxy that fits these criteria. Its high values of log [O iii]/Hβ = 0.86 ±  0.15 and log [O iii]/[O ii] = 0.5 ±0.2 are characteristic of Seyfert rather than LIER-type emission (Veilleux & Osterbrock 1987; Kewley et al. 2006).¹⁵ The emitting gas extends out to radii of at least 2 kpc, with a declining [O iii] EW that is qualitatively consistent with photoionization by a central source. No other emission line in any of the galaxies in our sample shows >2σ evidence for radial variation in Figure 12. A compact central source is detected in the HST images of MRG-M2129 (Section 5.3), which could be a Seyfert nucleus. All of these lines of evidence point toward photoionization by an AGN as one component of the excitation. Toft et al. (2017) came to similar conclusions on MRG-M2129 with some quantitative differences, which we discuss in the Appendix.

Although present in MRG-M2129, AGN photoionization is likely to not be the dominant source of excitation in MRG-S1522 or MRG-P0918. The lower [O iii]/Hβ ratio of MRG-S1522 is more consistent with a LIER than a Seyfert classification in the BPT diagram, although the measurement uncertainties make this ambiguous. Confirmation comes from its stronger [O ii] emission than [O iii], as seen in LIER-type spectra, with log [O iii]/[O ii] = −0.2 ± 0.1. In MRG-P0918, we do not detect [O iii] to sensitive limits. Although we cannot constrain the ratio [O iii]/Hβ and place this galaxy on the BPT diagram, the low EW of [O iii] and its weakness compared to [N ii] (with a ratio ≲0.1) are not typical of Seyfert-type spectra. For MRG-M0150, the only strong emission lines that are observable from the ground are [N ii] and Hα; so, for more information on this galaxy, we will have to await improved spectral coverage with JWST.

Photoionization of diffuse interstellar gas by hot evolved stars, such as post-asymptotic giant branch (post-AGB) stars, is often invoked to explain LIER emission in local early-type galaxies (e.g., Yan & Blanton 2012; Singh et al. 2013; Belfiore et al. 2016). At the ages relevant to quiescent galaxies at z = 2, i.e., ≃0.5–3 Gyr, stellar population synthesis models predict Hα EWs of ≃0.1–0.8 Å (Cid Fernandes et al. 2011, Figure 2). For MRG-M0150, MRG-S1522, and MRG-2129, we observe Hα EWs in the range 2.6–4.3 Å, which means there is insufficient ionizing flux from post-AGB stars to explain most of the emission. In MRG-P0918, however, the lower Hα EW of 0.3 ± 0.5 Å is consistent with expectations for post-AGB stars. Furthermore, in this scenario the gas should share the distribution and kinematics of the stars, which we observe. The gaseous and stellar velocity dispersions are consistent (Table 5) and the emission line EWs are essentially constant (Figure 12). The main difficulty with this interpretation is that photoionization models of 3–13 Gyr old populations produce log [N ii]/Hα ≲ 0.1 (Binette et al. 1994; Byler et al. 2017), which is smaller than observed in MRG-P0918. Since models at the relevant age of 0.5 Gyr have not been explored and probably have uncertainties in the shape of the ionizing spectrum, we consider that post-AGB stars may still be a viable explanation for most of the line emission in MRG-P0918 but not for the other galaxies in our sample.

Shocked gas is another possible source of high-[N ii]/Hα emission with low ionization. Shocks are consistent with the lower values of [O iii]/Hβ and [O iii]/[O ii] seen in MRG-S1522. Supporting evidence comes from the high line width of the ionized gas (σ = 345 km s⁻¹) relative to the stars (241 km s⁻¹), which indicates outflowing or turbulent ionized gas. We think that shocks are likely to contribute to the line emission in MRG-M2129 as well, since the line width of the gas is also elevated well above that of the stars (364 versus 266 km s⁻¹). Such a situation could in general be explained by a differing spatial distribution of the gas and stars. However, in the case of MRG-M2129, the EW of the [N ii] emission does not change much across the image (Figure 12), and the elevated ${\sigma }_{\mathrm{gas}}$ is not confined to the nucleus but extends to R ≃ 2 kpc. This suggests that both AGN photoionization and shocks likely contribute to the line emission in MRG-M2129.

We have explored whether the Allen et al. (2008) shock models can reproduce the line ratios in Figure 11. For MRG-S1522, we find that a model with twice solar metallicity, n = 1 cm⁻³, and B = 2 μG can reproduce the [N ii]/Hα and [O iii]/Hβ ratios for shock velocities v ≳ 400 km, comparable to ${\sigma }_{\mathrm{gas}}$.¹⁶ However, the agreement only holds if we consider the shock emission and not the radiative precursor, which may be an unphysical scenario. The emitting area required to produce the observed Hα luminosity is ≲100 (μ/5)⁻¹ kpc⁻², which corresponds to a spherical radius ≲3 (μ/5)^−1/2 kpc. Since this is comparable to the galaxy size, it requires a large-scale shock, but the energetics are feasible. On the other hand, the high [N ii]/Hα ratios of MRG-M2129 and MRG-0150 cannot be reproduced by any of the Allen et al. (2008) shock models or the Groves et al. (2004a, 2004b) AGN models. This is demonstrated by the grids shown in Figure 11, where we have chosen parameters that maximize the [N ii]/Hα ratio. Extracting more detailed physical information about these intriguing systems may require further development in these models.

In summary, four of the five lensed galaxies in our sample show line emission with very high [N ii]/Hα ratios that are inconsistent with a star formation origin. (For the fifth, MRG-M0138, all of the strong nebular lines except [O ii] are inaccessible from the ground.) MRG-M2129 shows clear evidence of AGN photoionization. Based on emission line ratios and the distribution and kinematics of the ionized gas, we argue that shocked gas is present in MRG-S1522 and MRG-M2129. Shocks may well power the line emission in MRG-M0150, but we lack the information needed to distinguish AGN photoionization. They could also contribute in MRG-P0918, but post-AGB stars might instead be the main ionizing source in that system. The high incidence of line emission in galaxies seen soon (<800 Myr) after quenching, the indications of shocked gas in many (and potentially all) cases, and the lack of ongoing star formation are all consistent with the idea that AGN-driven outflows or turbulence could play an important role in maintaining the quiescence of these galaxies. We will discuss these implications further in Section 8.3.

Further papers will present spatially resolved measurements of the stellar kinematics and populations. In order to interpret these measurements, we must be able to place this sample in the context of the full galaxy population. Here we compare our lensed galaxies to coeval samples selected from deep field surveys.

Figure 13 (left panel) compares the sample to the star formation "main sequence" derived in the 3D-HST and MOSDEF surveys. The lensed sample falls ≳1.5 dex below the main sequence. The right panel of Figure 13 shows that all galaxies in our sample fall in the quiescent region of the UVJ diagram. Consistent with the wide range of ages spanning 0.5–1.4 Gyr that we measured from the spectra and photometry (Section 6), our sample spans nearly the full range of the UVJ quiescent sequence. Only the reddest and presumably oldest galaxies are not represented.

Zoom In Zoom Out Reset image size

In Figure 14, the structural properties of the three galaxies with a lens model are compared to coeval quiescent galaxies in the 3D-HST survey. The left panel shows that the lensed galaxies scatter around the mass–size relation and so can be considered typically "compact." The right panel compares the axis ratios and Sérsic indices of our sample to those of similarly massive, coeval quiescent galaxies. The lensed sample displays Sérsic indices of n ≃ 1–4 and is consistent with being drawn from the 3D-HST distribution. On the other hand, the axis ratios of the lensed galaxies are at the round and flat extremes of the 3D-HST distribution.

Zoom In Zoom Out Reset image size

Two aspects of the sample that are perhaps surprising are the high ellipticities of two systems (MRG-M0138 and MRG-M2129) and the high stellar masses ($\mathrm{log}{M}_{* }=11.69\pm 0.19$ for MRG-M0138 and $\mathrm{log}{M}_{* }=11.50\pm 0.17$ for MRG-M0150) of two others. Considering first the ellipticities, we could suppose that errors in the lens model could lead to spuriously high values. However, we found that the ellipticity is robustly recovered across three lens models for MRG-M2129. In Paper II, we will show that their stellar kinematics imply that both MRG-M0138 and MRG-M2129 are intrinsically flat. These comparisons suggest that the high ellipticities are genuine. One might suppose that a selection effect could enhance their probability of inclusion in our sample. This might be the case if images were selected for follow-up based on their length-to-width ratio, for example, but our follow-up observations were based on a complete color–magnitude–selected sample in the central regions of the galaxy clusters in our survey. Given the small number of objects, it may be that high ellipticities are simply randomly over-represented in the three lensed galaxies with lens models. Regardless of its origin, this over-representation is important to bear in mind for our kinematic analysis in Paper II.

Turning to the remarkably high stellar masses inferred for MRG-M0150 and MRG-M0138, we note that since we are studying galaxies on the bright tail of the luminosity function, we expect the lensed galaxies to pile up near the flux limit of our selection box. This is generally the case for the H band flux (Figure 1), with the exception of the ultra-bright MRG-M0138. In Paper II, we will show that its dynamical mass agrees with the stellar mass inferred in this paper. Since the stellar and dynamical masses have different dependences on the lens mapping, this consistency provides some reassurance that the magnification is not very far in error, i.e., well beyond the ∼0.2 dex uncertainty we estimated in Section 5.1.2. For MRG-M0150, we will show in Paper II that the dynamical mass is ≃0.2 dex smaller than the stellar mass. This could indicate that the magnification is underestimated and the stellar mass is overestimated, but the difference is within the estimated uncertainties and, in any case, still implies a very massive galaxy.

We conclude that the sample is broadly representative of the colors and sizes of massive quiescent galaxies at z = 2–2.6. The stellar masses and ellipticities are likely to be accurate within our estimated uncertainties, but the subset of three galaxies with lens models has a high proportion of high-ellipticity galaxies compared to unlensed samples and so is not fully representative in that property. Future lens models with additional constraints from any newly identified multiple image systems would be useful to help validate our magnification estimates.

The structures of the galaxies in our sample are very different from typical early-type galaxies in the local universe with a similar mass. Their sizes are comparable to other z = 2 quiescent galaxies (Figure 14) and are smaller than z = 0 galaxies of equal mass. Even more striking, however, is the presence of a dominant disk component in two of the three galaxies for which we have reconstructed the source, MRG-M0138 and MRG-M2129. Both galaxies have ellipticities $e\gt 0.7$. To place these shapes in the context of low-redshift galaxies, we selected galaxies with $\mathrm{log}{M}_{* }\gt 11$ and z = 0.05–0.25 from the GAMA survey (Driver et al. 2009; Taylor et al. 2011; Kelvin et al. 2012). Only 4% of the local galaxies have ellipticities $e\gt 0.7$. Even more rare is the apparent lack of a bulge in MRG-M2129, which is well-described by a pure exponential disk (n = 1), as also remarked by Toft et al. (2017). In our local comparison sample, only 0.9% of galaxies have $n\leqslant 1$.

Therefore, these galaxies need to evolve in size, shape, and concentration to resemble their z ∼ 0 descendants. The extreme differences in structure imply that relatively little of the total transformation from star-forming galaxies into local early-type systems occurred when these galaxies were quenched. Instead, these changes must have occurred later, likely through a series of major and minor mergers. The structures of the galaxies in our sample also have implications for their earlier evolution. In particular, whatever processes quenched star formation in these galaxies did not destroy the stellar disk or, in the case of MRG-M2129, even produce a significant bulge. In Paper II, we will discuss the past and future evolution of these galaxies in greater detail by bringing additional kinematic evidence to bear.

The nearly ubiquitous presence of line emission with high [N ii]/Hα ratios in our sample of massive, quiescent galaxies is striking. While massive (${M}_{* }\gtrsim {10}^{10.9}\,{M}_{\odot }$) star-forming galaxies at z = 1–3 commonly host nuclear outflows that are thought to be driven by AGN (Genzel et al. 2014), less is known about quiescent galaxies in this redshift range. Belli et al. (2017a) studied the incidence of line emission in z = 0.7–2.7 galaxies in the KMOS^3D survey. They detected unambiguous line emission in 17% of the quiescent galaxies (as classified by UVJ colors) and in about half of these attributed the emission to shocks. It is interesting that in our sample the incidence is much higher—emission lines are detected in all four galaxies for which the strong lines are observable from the ground—and that our sample reaches [N ii]/Hα ≃ 6, whereas values above ∼1 are absent from the Belli et al. sample.

These differences might be explained in part by differences in the sensitivities of the observations and in the stellar masses characterizing each sample. Belli et al. found a trend for galaxies with lower Hα EWs and higher stellar masses to have higher [N ii]/Hα ratios. Since the lensed galaxies in our sample have, on average, lower Hα EWs and higher masses than the KMOS^3D quiescent galaxies, it is possible that future observations will show that the KMOS^3D and lensed galaxies sample different parts of a common sequence. We emphasize that we can characterize the weak emission lines in our sample of quiescent galaxies only because of the galaxies' brightness (due to their high masses combined with lensing magnification) and our accurate modeling of the stellar continuum. With shallower data, we could not measure the infilling of the stellar Hα absorption. If a stellar redshift could not be measured, then one could easily mistake the stronger [N ii] $\lambda 6585$ line for Hα, leading to a completely different interpretation.

Shocks are most clearly present in two cases (MRG-S1522 and MRG-M2129), as judged by the line ratios, gas kinematics, and in the case of MRG-M2129, the extended distribution of the [N ii] emission. Although shocks could also be present in MRG-P0918 and MRG-M0150, other explanations remain viable in these cases because of the fewer constraints on MRG-M0150 and the low Hα EW seen in MRG-P0918 (see Section 7.2). Given that the current rates of star formation are very low in MRG-S1522 and MRG-M2129 ($\lesssim 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ ; Table 5), sources of energy are needed other than massive stars and core-collapse supernovae. Type Ia supernovae occurring at the expected rate (Maoz et al. 2012), which is much higher than in local massive ellipticals, can supply kinetic power comparable only to the Hα luminosities seen in our sample, which is expected to be ≲2% of the shock luminosity (Allen et al. 2008). Therefore, shocks powered by type Ia supernovae ejecta cannot generate the observed line emission. Mergers could produce shocks, but we do not find evidence for ongoing mergers in our sample except for a possible faint companion of MRG-M2129, and this explanation seems inconsistent with the spatial uniformity of the [N ii] emission. Therefore, AGNs provide the most obvious energy source.

Observations of some low-redshift quiescent galaxies have revealed significant reservoirs of turbulent molecular gas with ∼kpc (or larger) sizes and suppressed star formation relative to the Kennicutt–Schmidt relation, often combined with warm shocked gas (e.g., Alatalo et al. 2015; Guillard et al. 2015; Lanz et al. 2016). A significant fraction of intermediate-mass quiescent galaxies in the local universe appear to have ionized gas outflows (Cheung et al. 2016). The injection of turbulence into the interstellar medium via an AGN jet has been suggested to suppress star formation in these galaxies. Our observations are at least consistent with a similar mechanism acting in MRG-S1522 and MRG-M2129, and potentially in all of the galaxies in our sample that have ages ≲1 Gyr.

The high fraction of massive systems in which an AGN seems to be affecting the gas—both in high-mass star-forming galaxies at z ∼  2 and in recently quenched galaxies—suggests than AGN activity might play an important role in quenching star formation and maintaining low star formation rates in these galaxies, although proving causality is very difficult. Our sample affords detailed observations but is small. More deep spectra of z ≳ 2 quiescent galaxies are needed to firmly establish the prevalence of emission lines and the variation of their properties with the time since quenching. Observations of the distribution and kinematics of any molecular gas in such galaxies may also be revealing.

In four of the lensed quiescent galaxies discussed in this paper (i.e., all but MRG-S1522), we are able to spatially resolve the stellar continuum in our NIR spectra. In the companion Paper II, we present the resolved stellar kinematics of these galaxies. In a future paper, we plan to study their resolved stellar populations in order to dissect the star formation histories of early quiescent galaxies. We will also use observations of the extraordinarily bright MRG-M0138 to measure its multi-element stellar abundance pattern. Many of these measurements are currently only possible for quiescent galaxies that are gravitationally lensed. Observations with JWST will only marginally resolve a typical compact quiescent galaxy at z ∼ 2. While such data will be very valuable, lensed quiescent galaxies will continue to offer the highest-resolution views of these galaxies and will provide unique insights into their formation, even in the JWST era.