The Complete Redshift Distribution of Dusty Star-forming Galaxies from the SPT-SZ Survey

Introduction

Millimeter and submillimeter continuum observations have transformed our understanding of galaxy formation and evolution by demonstrating that luminous, dusty galaxies were 1000 times more abundant in the early universe than they are today (see reviews by Blain et al. 2002 andCasey et al. 2014). The most intense star formation in the universe takes place in these high-redshift (z > 1) dusty star-forming galaxies (DSFGs), which form new stars at rates of >100-1000 -M yr 1  behind dense shrouds of dust. DSFGs are thought to be the progenitors of the massive elliptical galaxies seen in the present-day universe (Blain et al. 2004). The exact details of how these galaxies form stars at such prodigious rates is still an open question (Narayanan et al. 2015), though galaxy mergers likely play a role (e.g., Tacconi et al. 2008;Engel et al. 2010;Narayanan et al. 2010;Hayward et al. 2011Hayward et al. , 2013;;Bothwell et al. 2013;Ma et al. 2016). While the first surveys of the redshift distributions of these DSFGs suggested that the population peaks at z ∼ 2 (Chapman et al. 2003(Chapman et al. , 2005)), modern redshift surveys suggest that the DSFG population peaks at higher redshifts (2.5-2.9; Simpson et al. 2014;da Cunha et al. 2015;Danielson et al. 2017;Dudzevičiūtė et al. 2020). These updated surveys are also finding objects at increasingly higher redshift, extending past z > 6 (e.g., Riechers et al. 2013;Fudamoto et al. 2017;Strandet et al. 2017;Marrone et al. 2018;Zavala et al. 2018). The presence of these extremely highredshift DSFGs challenges our understanding of the underlying distribution of these objects and their role in the cosmic star formation history.

Dust emission at high redshift (z > 1) exhibits dimming from increased cosmological distance, which is counteracted by a steep rise on the Rayleigh-Jeans side of the spectral energy distribution (SED) at a fixed observing wavelength, which leads to the so-called "negative K-correction" in the submillimeter (Blain & Longair 1993). Because of this effect, fluxes at millimeter and submillimeter wavelengths are roughly constant with redshift, facilitating the discovery of highredshift DSFGs. However, the dust-obscured nature of DSFGs suppresses emission at optical/UV wavelengths, making robust redshifts difficult to obtain, especially at high redshift. Technological advances in correlator bandwidth have enabled "blind" spectroscopic surveys to be conducted at facilities such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Plateau de Bure Interferometer (PdBI), which is now the Northern Extended Millimeter Array (NOEMA). These spectroscopic surveys search for molecular emission at millimeter wavelengths and can be conducted without prior optical/near-IR spectroscopy. Because the molecular emission can be unambiguously related to the continuum emission, these surveys provide a direct and unbiased way to derive the redshifts of DSFGs (e.g., Scott et al. 2011;Weiß et al. 2009). Carbon monoxide (CO) is the second most abundant molecule in the universe after H 2 and its rotational transitions are frequently targeted in blind searches. These rotational transitions are spaced evenly every 115 GHz and are among the brightest lines in the millimeter spectrum. The CO line brightness is due to its abundance. The low energy required to excite its rotational states, and the fact that the resulting emission lines are mostly accessible at frequencies of high atmospheric transmission make it ideal for line surveys with ALMA.

One of the first millimeter-wave redshift searches with ALMA identified the redshifts of 26 gravitationally lensed DSFGs selected with the South Pole Telescope (SPT; Vieira et al. 2013;Weiß et al. 2013). Because gravitationally lensed sources (with magnifications μ ∼ 10) are apparently brighter than unlensed sources, they require significantly less on-source time to survey ( m µ - t obs 2 ). Blind CO surveys can be conducted with comparative ease and enable larger spectroscopic surveys (Weiß et al. 2013;Strandet et al. 2016;Neri et al. 2020). While some unlensed blind CO surveys have been conducted (Chapman et al. 2015), the majority of unlensed redshifts were obtained through optical spectroscopic searches (Chapman et al. 2003(Chapman et al. , 2005;;Casey 2012;Koprowski et al. 2014;Brisbin et al. 2017;Danielson et al. 2017) or photometry (Simpson et al. 2014(Simpson et al. , 2017;;Michałowski et al. 2017). The discrepancies between the peaks of unlensed and lensed distributions can be explained by a combination of the selection wavelength, survey depth, and the redshift-dependent probability of strong lensing (e.g., Béthermin et al. 2015b;Strandet et al. 2016).

Wide-area surveys conducted with the Atacama Cosmology Telescope (ATCA; Marsden et al. 2014), Herschel (Eales et al. 2010;Oliver et al. 2012), Planck (Planck Collaboration et al. 2015), and the SPT (Vieira et al. 2010(Vieira et al. , 2013;;Carlstrom et al. 2011) span hundreds to thousands of square degrees and have enabled the discovery of hundreds of gravitationally lensed DSFGs at millimeter and submillimeter wavelengths. For DSFGs at high-redshift, the rest-wavelength peak in the SED at ∼ m 100 m is shifted into the observing bands of these millimeter and submillimeter instruments. Most of the brighter sources are gravitationally lensed as well (e.g., Blain 1996;Negrello et al. 2010;Wardlow et al. 2013;Spilker et al. 2016). The few intrinsically bright unlensed sources in these samples are typically major mergers of DSFGs (e.g., Fu et al. 2013;Ivison et al. 2013) or protoclusters (e.g., Casey 2016;Overzier 2016;Lewis et al. 2018;Miller et al. 2018;Oteo et al. 2018;Hill et al. 2020). The magnification due to gravitational lensing enables us to obtain spectroscopic redshifts for a complete sample of DSFGs, measure the redshift distribution of DSFGs, and ascertain the prevalence of the highest-redshift DSFGs. Samples of lensed DSFGs also afford the opportunity to study fainter observational diagnostics, in greater detail, than would otherwise be possible(e.g., Bothwell et al. 2017;Béthermin et al. 2018;Spilker et al. 2018;Zhang et al. 2018;Dong et al. 2019;Jarugula et al. 2019;Litke et al. 2019;Cunningham et al. 2020).

In this paper, we finalize the SPT-selected ALMA redshift survey, which started in Weiß et al. (2013) and continued in Strandet et al. (2016). The final catalog contains spectroscopic redshifts for all sources, making it the largest and most complete catalog of its kind to date. In Section 2, we present the 3 mm line scans obtained from ALMA (Section 2.2), including 40 new sources. We also present the photometry from SPT, ALMA, Herschel and Atacama Pathfinder Experiment (APEX; Section 2.3.1). In Section 2.3.2, we present a methodology for fitting SEDs and deriving intrinsic source properties. Section 3 is divided into two parts: spectroscopic results (Section 3.1) and photometric results (Section 3.2). In Section 3.1, we present the 3 mm spectra and the resulting spectroscopic redshifts. In Section 3.2, we describe the fitted SEDs and the resulting intrinsic properties. In Section 4, we discuss the redshift distribution of the complete SPT sample (Section 4.1), the possibility of temperature evolution (Section 4.2), the extreme nature of the SPT sources (Section 4.3), and the resulting high-redshift tail of this distribution (Section 4.4).

For this paper, we adopt a flat Lambda cold dark matter cosmology, with W =

ALMA 3 mm Blind Scans

In order to obtain redshifts for the SPT sample, a blind spectroscopic redshift search was started in ALMA Cycle 0 (project ID: 2011.0.00957.S). This program resulted in a 90% line detection rate for the 26 sources surveyed (Weiß et al. 2013). An updated distribution was presented in Strandet et al. (2016) with an additional 15 sources observed in ALMA Cycle 1 (project ID: 2012.1.00844.S). This work represents the conclusion of the SPT blind redshift survey, and presents spectroscopic scans for the remaining 40 sources from ALMA Cycles 3, 4, and 7, and 41 new spectroscopic redshifts. The individual 3 mm scans for all sources can be found in Appendix B.

The blind spectroscopic search was conducted using ALMA's Band 3 receiver, which operates between 84 and 116 GHz. The correlator has a total bandwidth of 7.5 GHz, which is split across two sidebands. The ALMA Band 3 3 mm atmospheric transmission window can be covered in five tunings, as shown in the left panel of Figure 2. This configuration results in overlapping coverage in the 96.2-102.8 GHz region. ALMA's primary beam ranges from 45″-61″over the entire scanned frequency range.

Figure 2 also demonstrates this search's sensitivity to CO lines between the CO 1 0 ( -) and CO 8 7 ( -) transitions. Scanning this region results in redshift coverage of 0.0 < z < 0.4 and 1.0 < z < 8.6 with a narrow redshift desert at 1.74 < z < 2.00. As previously shown in Weiß et al. (2013) and Strandet et al. (2016), the SPT sources have a median redshift of z∼4, enabling more than one emission line to be observed in each spectrum. Notably, the [C I](1-0) line is predicted in the range 3.3 < z < 4.8 and has been shown in Weiß et al. (2013) and Strandet et al. (2016) to be bright enough to provide ). At a fixed T dust , this color can be used as a crude proxy for redshift. secure spectroscopic redshifts at the observed sensitivities. The H 2 O(2 0,2 -1 1,1 ) and [C I](2-1) lines are accessible at the highest redshifts (Figures 2 and3).

The Cycle 3 ALMA observations were conducted from 2015 December to 2016 August (project ID: 2015.1.00504.S). During Cycle 3, between 34 and 41 antennas were employed, and resulted in typical synthesized beams of 3.9″×4.7″ to 3.3″×4.1″ (FWHM) from the low to high frequency ends of the band. Each target was observed for roughly 6minutes onsource. Sources were grouped such that at least four targets were observed in a single execution block, so as to minimize overheads. This resulted in total observation times of 25-45 minutes per source. Typical system temperatures were measured to be = T 55 84 K sys -. Flux calibration was performed on Uranus, Neptune, Ganymede, J0519-4546, and J0538-4405. Bandpass and phase calibration were determined from nearby quasars.

The data were processed using the Common Astronomy Software Application (CASA; McMullin et al. 2007;Petry et al. 2012). Calibrated data cubes were constructed using CASAʼs TCLEAN package. The cubes have a channel width of 62.5 MHz (∼ -220 km s 1 ). Observations from Cycle 3 had a typical noise per channel of 0.5-0.8 mJy beam -1 over the 37 GHz bandwidth. The TCLEANed continuum images have typical noise levels of 50 μJy beam -1 .

The Cycle 4 observations were conducted from 2016 November to 2017 May (project ID: 2016.1.00672.S). In Cycle 4, each scan utilized between 38 and 46 antennas, resulting in typical synthesized beams of 4.1″×5.0″ to 3.5″×4.3″ (FWHM) from the low to high frequency ends of the band. Each target was observed for ∼12-15 minutes onsource for each source, not including overheads. Sources were again grouped such that each execution block contained multiple targets. However, because groupings were not possible for every source, the total observation time ranged between 38 and 90 minutes per source, including overheads. Typical system temperatures were measured to be = T sys 60 89 K -. Flux calibration was performed on Mars, Uranus, Neptune, J0334-4008, J0538-4405, J2056-4714, and J0519-4546. Bandpass and phase calibration were determined from nearby quasars. Because more antennas were available in Cycle 4, the typical noise per channel decreased to 0.4-0.6 mJy beam -1 . The TCLEANed continuum images from Cycle 4 had typical noise levels of 40 μJy beam -1 .

Because a total of three sources did not exhibit lines in their initial 3 mm line scans (SPT0112-55, SPT0457-49, and SPT2340-59; more detail in Section 3.1), additional deeper 3 mm scans for these sources were conducted in Cycle 7. These observations were conducted from 2019 November to 2020 January (project ID: 2019.1.00486.S) with the aim of observing Greve et al. (2012) and Strandet et al. (2016). By performing a 1.4 mm flux-weighted average of the observed continuum-subtracted rest frame spectra, we obtain the composite spectrum presented in the bottom panel. This was first done in Spilker et al. (2014), but has been updated to include the final SPT-selected sample. Emission lines used for redshift confirmation are shown in dark gray, while other important emission lines are shown in light gray.

possible line features with fluxes of 1-2 mJy identified in the earlier scans. Each scan utilized between 42 and 49 antennas, resulting in minimum and maximum angular resolutions of 2.0″-4.0″ from the high to low frequency ends of the observed band. Each target was observed for 45-91 minutes on-source, for a total time of 130-300 minutes when including overheads. Typical system temperatures were = T 62 72 K sys resulting in a noise per channel of 0.2-0.3 mJy beam -1 . Flux calibration was performed on J0006-0623, J0519-4546, and J0238 +1636. Bandpass and phase calibration were determined from nearby quasars. The continua for the three sources reobserved in Cycle 7 had typical noise levels of 10-70 μJy beam -1 . The decreased noise in Cycle 7 is to due a combination of longer integration times and more antennas being available.

Additional Spectroscopic Observations

Because many observations from the blind 3 mm line scans described in Section 2.2 contain single CO lines, additional observations were required in order to break degeneracies between redshift solutions and obtain unambiguous spectroscopic redshifts. In addition to the 3 mm ALMA survey conducted in this work, we have also conducted surveys of the [C II] emission line and low-J CO lines. The [C II] emission line survey was conducted using the First Light APEX Submillimetre Heterodyne receiver (FLASH; Heyminck et al. 2006). A subset of these observations were published in Gullberg et al. (2015) and used to confirm three redshifts from Strandet et al. (2016). More recently conducted [C II] observations were used to confirm four spectroscopic redshifts in this work. The ATCA was used to conduct a survey of CO(1-0) and CO(2-1) and confirmed a subset of observations (Aravena et al. 2016). One of these observations was obtained after the publication of Aravena et al. (2016) is used to confirm one additional redshift in this work. The details for all of these observations can be found in Appendix C.

An additional 13 redshifts were confirmed through targeted CO line searches in the ALMA 2 mm Band. CO emission lines were targeted using the possible redshift solutions from the observed 3 mm emission line. To further guide which of the degenerate solutions should be targeted, the photometric redshift obtained from SED fitting was used to guide the target selection. The ALMA correlator sidebands were configured such that at least one CO emission line would be observed for a given redshift. The expected line strength was determined using the SPT-DSFG CO spectral line energy distribution from Spilker et al. (2014). These observations were carried out in ALMA Cycles 6 and 7 (project IDs: 2018.1.01254.S and 2019.1.00486.S) and all sources targeted yielded CO emission detections at sufficient significance levels to confirm the redshifts. Details on these observations and the obtained spectra can be found in Appendix C.1.

Of the spectroscopic redshifts presented, only two are based on a single line and are still awaiting additional observations to confirm the redshift (see Section 3.1.2 for a detailed discussion).

Photometry

The SPT DSFGs have superb alternatively, far-infrared (FIR) through millimeter photometric coverage, with flux densities measured at 3 mm (ALMA), 2 mm, 1.4 mm (SPT), m 870 m (APEX/LABOCA), and 500 μm, m 350 m, and m 250 m (Herschel/Spectral and Photometric Imaging Receiver (SPIRE)) for all sources. Additional Herschel/Photodetector Array Camera and Spectrometer (PACS) m 160 m and m 100 m observations were obtained for a subset of 65 sources. Despite the large range in redshifts (1.9 < z < 6.9), the photometry is complete between m 71 m < l rest < m 380 m, and the peak of the FIR SED at ∼100 μm is always well constrained. The flux densities for all photometric points can be found in Appendix D. The absolute calibration uncertainties of 10% for Herschel/PACS, 7% for Herschel/SPIRE data, 12% for APEX/LABOCA, 7% for SPT, and 10% for ALMA data added in quadrature to the errors quoted in Appendix D.

ALMA-The ALMA 3 mm continuum maps were obtained as a result of the observations described in Section 2.2.1. The continuum images were created using CASAʼs TCLEAN procedure with the full observed bandwidth (84.2-114.9 GHz) with natural weighting in order to optimize sensitivity. In cases where the source is unresolved, meaning >90% of the total flux detected was contained within one beam, we extract flux from the brightest pixel detected from the continuum to obtain a spectrum. The error on flux density was calculated using the rms of the residual map produced by the TCLEAN procedure. However, half of the SPT-selected sources are marginally resolved (>80% of the source's flux is contained within one beam). In order to obtain spectra for these sources, CASAʼs imfit routine is used to fit a 2D Gaussian to the source in each datacube slice and extract the flux and associated error.

SPT-The SPT 1.4 and 2.0 mm flux densities were extracted from cosmic microwave background (CMB) maps acquired from the first survey, SPT-SZ. This survey was completed in 2011 November and covered 2500 deg 2 of the southern sky in three frequency bands, 95, 150, and 220 GHz (corresponding to 3.2, 2.0, and 1.4 mm, respectively) with arcminute angular resolution. Absolute calibration for both the 1.4 and 2.0 mm bands is derived from the CMB and the calibration uncertainty is 10%. The data were extracted and deboosted according to the procedure described in Everett et al. (2020).

APEX-The sources were observed at m 870 m with LABOCA at APEX and the flux densities were extracted. LABOCA is a 295element bolometer array with an 11.4′ field of view and a measured angular resolution of 19.7″ (FWHM). The center frequency of LABOCA is 345 GHz ( m 870 m) with a passband FWHM of ~60 GHz. The measured noise performance for these observations was 60 mJy s 1 2 . These observations were performed between 2010 September and 2013 October (project IDs: M-085. F-0008-2010, M-087.F-0015-2011, E-087.A-0968B-2011, M-089. F-0009-2012, E-089.A-0906A-2012, M-091.F-0031-2013, E-091. A-0835B-2013, M-092.F-0021-2013). For more details on the observations, see Strandet et al. (2016).

APEX/LABOCA maps were created for each source using the Bolometer Array Analysis Software (BoA; Schuller 2012). The resulting time-ordered data undergo various calibration, noise removal, and flagging procedures detailed fully in Greve et al. (2012). The data is then gridded and individual maps are co-added with inverse variance weighting. The flux densities were either extracted from the peak flux density, in case of pointlike sources or by integrating over the emission region, in cases where LABOCA resolves the emission.

Herschel/SPIRE-Flux densities at 250, 350, and 500 μm for all sources were measured by the SPIRE (Griffin et al. 2010) onboard the Herschel Space Observatory. The data were observed in two programs (project IDs: OT1_jvieira_4 and OT2_jvieira_5) conducted between 2012 August and 2013 March. The Herschel/SPIRE data consists of triple repetition maps, with coverage complete to a radius of 5′from the nominal SPT position. The maps were produced using the standard reduction pipeline HIPE v9.0 (Ott 2011;Balm 2012). Flux densities were extracted by fitting a Gaussian profile to the SPIRE counterpart of the SPT detection and the noise was estimated by taking the rms in the central 5′of the map.

Herschel/PACS-Additional data were obtained for a subsample of 65 sources at 100 and m 160 m using the PACS onboard Herschel (project IDs: OT1_jvieira_4, OT1_dmar-rone_1, OT2_jvieira_5 and DDT_mstrande_1). The additional data ensures that the thermal peak is well-sampled and has been constrained for all SPT sources with z < 2.5. The data were acquired using approximately orthogonal scans centered on the target at medium speed (i.e., with the telescope tracking at ¢ - 20 s 1 ), spending a total of 180 s on-source per program. Each scan was composed of 10 separate 3′strips, each offset orthogonally by 4″ and both wavelengths were observed simultaneously. The scans were co-added and weighted by coverage. The data were then handled by a variant of the reduction pipeline presented in Ibar et al. (2010). The resulting noise levels were calculated using random aperture photometry and were found to be σ≈4 and 7 mJy at 100 and m 160 m, respectively.

SED Fitting

We fit each source in each sample with a modified blackbody law (e.g., Blain et al. 2003;Casey 2012) given by

where B ν is the Planck function for a value of T dust or the CMB temperature, T CMB . In order to reduce the number of free parameters and mitigate the degeneracies between redshift and T dust , we fix the Rayleigh-Jeans spectral slope, β. Empirically, the value β=2 was well matched to the data, so we fix this parameter for our modified blackbody fits in a similar fashion to what was done in Greve et al. (2012). However, rather than fix n » 3000 GHz 0 (λ 0 ≈100 μm), we use the empirical relationship between λ 0 and T dust given in Equation (2) of Spilker et al. (2016) to constrain λ 0 . Using this relation provides a better alternative to assuming a single value for λ 0 when an independent estimate of the size of the emission region is not available. We find that the introduction of this dependency between T dust and λ 0 improves both the reduced χ 2 value and photometric redshift. It should be noted, however, that this procedure tends to increase the value of the dust temperature by ∼20%. The only free parameters in this SED fit are the overall SED normalization, dust temperature, and redshift.

Because a modified blackbody fit alone does not typically describe the mid-IR excess found in the Wien side ( l < rest = m 50 m) of the thermal emission peak, we perform another fit including an additional power law component (Blain et al. 2003). The power law component introduces another free parameter, α, which is the power law slope. The combined modified blackbody and power law fit empirically describes all of the available photometry, including the data on the Wien side of the thermal emission peak. In this work, we use this SED fit to define the frequency at which thermal emission peaks, l peak , obtain a best fit to β for the M d calculation, and to determine total L IR .

In order to fit the data, we employ a Markov Chain Monte Carlo (MCMC) algorithm using the emcee package (Foreman-Mackey et al. 2013) to sample the posterior probability function. To ensure uniform photometric coverage in the fit region, we mask data shortward of l rest = m 50 m (Greve et al. 2012) for the modified blackbody fits. Fitting done with an additional power law included all available photometry points. The results of this fitting procedure are described in Section 3.2.

Calculating Intrinsic Source Properties

In this section, we describe the intrinsic source properties, which are calculated using the SED fits. Because we constrain λ 0 as a function of T dust (as described in Section 2.3.2), we are able to better understand the T dust distribution of the sample. The apparent FIR luminosity (L FIR ) is calculated by integrating the fitted SED over the wavelength range m 42.5 122.5 m - (Helou et al. 1988). In order to obtain the IR luminosity (L IR ), we integrate the modified blackbody function with an additional power law over the m 8 1000 m range. With FIR luminosity and T dust values, we derive star formation rates (SFRs) and dust masses for each source. The dust masses are calculated according to

where n S is the flux density at 345 GHz in the rest frame, determined from our SED fit, D L is defined as the luminosity distance, T z CMB ( ) is the CMB temperature at redshift z, and μ is the magnification factor. We adopt

) (Weingartner & Draine 2001;Dunne et al. 2003), where β=2.0 is the dust emissivity index.

To derive total SFRs, we use the following conversion from Murphy et al. (2011):

where the infrared luminosity, L IR , is calculated from the m 8 1000 m range using the modified blackbody fit with an additional power law to describe the mid-IR excess. This conversion was calculated using Starburst99 (Leitherer et al. 1999) for a Kroupa initial mass function (Kroupa 2001).

In order to calculate the intrinsic properties of these sources, the magnifications presented in Spilker et al. (2016) are used when available. Spilker et al. (2016) surveyed a sample of 47 SPT DSFGs and constructed gravitational lens models. However, because of the final flux density cut ( m S 870 m > 25 mJy) used to define the sample presented in this work, not all sources modeled in Spilker et al. (2016) were retained. As a result, 39 of the sources presented in this work have detailed lens modeling. For sources with multiple components, a flux-weighted average is used as the magnification. For the remaining sources without lens models, the median magnification ( m á ñ m 870 m = 5.5) of the modeled sources is adopted.

Spectroscopy Results

Building on the work of Weiß et al. (2013) and Strandet et al. (2016), the final SPT-selected DSFG sample is composed of 81 sources. We have obtained spectroscopic redshifts for the complete SPT-selected sample, making our catalog the largest and most complete redshift survey of high-redshift DSFGs to date. We begin by presenting the final blind 3 mm CO line scans obtained by ALMA and the resulting spectroscopic redshifts. For the cases where only a single CO line was detected, we discuss any ancillary spectroscopic data used to confirm the redshift in Appendix C.4. We also discuss the three 3 mm spectra where no spectroscopic lines were present in the initial 3 mm scans. These sources were reobserved and deeper scans enabled secure redshifts to be obtained. Finally, we discuss the two single line spectra where no ancillary data has yet been obtained and discuss the most probable redshift.

All of the SPT-selected DSFG spectra, including those originally published in Weiß et al. (2013) and Strandet et al. (2016), are summarized in Figure 3. A complete summary of the spectroscopic lines detected for each source can be found in Appendix E. In this work, we detect 62 strong line features from 12 CO and [C I] with an integrated signal-to-noise ratio > 5σ. We detect an additional 29 weaker features (>3σ), which include HCN, HCO + , H 2 O, 13 CO, and CN.

We detect 3 mm continuum emission for all the 40 previously unpublished SPT-selected DSFGs presented in this work. The positions for these sources were obtained by fitting Gaussian profiles to the ALMA3 mm data and are listed in Appendix A. The continuum flux densities were also obtained in these fits and are given with the other photometric observations in Appendix D.

SED Fits, Dust Temperature, and Photometric Redshifts

More information can be obtained by fitting the photometry using the process outlined in Section 2.3.2. We first fit a modified blackbody to the FIR thermal emission peak for rest wavelengths >50 μm. We also fix the redshift parameter to the spectroscopic values obtained in Section 3.1, shown in Figure 4 in red. This model describes the FIR thermal emission peak well, with a median reduced χ 2 value of 1.8. We also perform another fit with an additional power law component in order to fit the complete L IR wavelength range (shown in Figure 4 in blue), which gives a median reduced χ 2 of 0.7.

Sources with large χ 2 values generally exhibit a discrepancy between the 3 mm flux and the value predicted by the modified blackbody. Because of the discrepancy in beam sizes between ALMA and SPT, objects that are broken into multiple components are sometimes below the ALMA 3 mm detection threshold, leading to an underestimation of the total 3 mm flux. Using ALMA m 870 m imaging (Spilker et al. 2016), we verify that the source is split into multiple components. In cases of multiple components, the 3 mm point is then masked when fitting the thermal emission peak. One of these sources, SPT2349-56, has already been identified as a protocluster (Miller et al. 2018;Hill et al. 2020). This discrepancy in 3 mm flux from over-resolving the emission may be a good way to , such that all sources have roughly uniform photometric coverage and excess emission on the Wien side of the blackbody does not artificially drive the dust temperature toward higher values. Sources that exhibit multiple components in 3 mm are fitted by masking the 3 mm point. Any photometry points that were masked for the modified blackbody fit are represented by the gray squares. In order to account for the mid-IR excess, a power law can be added on to the modified blackbody function (Blain et al. 2003) to better describe all available data (blue).

separate unlensed protoclusters from lensed DSFGs. The remaining sources are being investigated as potential protocluster candidates. As previously stated, any photometry <50 μm in the rest frame was masked for the modified blackbody fit. Any masked photometry points are shown in Figure 4 in gray.

By extracting a dust temperature from the SED fit, a dust temperature probability distribution was created by sampling each source's dust temperatures 10 3 times using a Monte Carlo procedure, and is shown in Figure 5. Though the distribution peaks near T dust =40 K, the median of the distribution is significantly higher at = T 52.4 dust ±2.3K, with a tail extending past ∼100 K. Given the relationship between λ 0 and T dust discussed in Section 2.3.2, the implied median of the λ 0 distribution is m  155 7 m. These warm dust temperatures suggest that there could be a correlation between redshift and dust temperature, which we discuss in Section 4.2.

Although spectroscopic redshifts have been obtained for all sources, the photometric redshift can still be used in order to break degenerate redshift solutions in the two cases where a single CO transition was detected, as well as inform future spectroscopic and photometric surveys. In order to find the photometric redshift, we use the dust temperature distribution given in Figure 5 as a prior and simultaneously fit the T dust and redshift parameters. While this method exhibits good agreement with the spectroscopic redshift, as seen in Figure 6, the associated errors are large. However, this method can be used to ascertain a rough estimate of redshift in the absence of observationally expensive spectroscopic data. (See, e.g., Casey 2020 for a similar method and comparison to template fitting.)

This photometric redshift fitting can also be used to break degeneracies for the two single 3 mm line sources (SPT0150-59 and SPT0314-44). For a given degenerate redshift solution, we fix the redshift parameter to that solution and use the SED fitting procedure to find the associated dust temperature. Using the dust temperature distribution given in Figure 5 as a prior, we assign a likelihood of the source being at each possible redshift solution (see Appendix C.4 and Figure C5 for details). For both single line sources, CO(3-2) is the most probable identification. If the detected line was in fact the higher-J transition (e.g., CO(4-3)), an additional CO line should have been detected. Taken together, these two independent pieces of evidence indicate that there is no ambiguity in the redshift of these sources.

Though we are complete in spectroscopic redshifts, photometry provides insights that inform future surveys. Because the SPT DSFGs have excellent FIR-millimeter coverage, it is possible to use single flux density ratios to obtain information about the redshift. For instance, the m m S S 350 m 870 m ratio can be used as a rough indicator of redshift. We find relatively good agreement between our catalog, Herschel/SPIRE sources (Negrello et al. 2010;Conley et al. 2011;Omont et al. 2011;Harris et al. 2012;Wardlow et al. 2013;Bussmann et al. 2013;Gladders et al. 2013;Riechers et al. 2013;Messias et al. 2014) and the redshifted SED of Arp220 (Silva et al. 1998), shown in Figure 7. We perform a maximum likelihood estimation assuming the 350-870 μm flux density ratios for the Herschel/SPIRE and SPT sources can be described by an exponential model with an extra Gaussian variance term to account for intrinsic scatter. We find the following exponential fit best describes the available data:

which is shown with a 1σlimit in Figure 7. In performing this fit, we find an intrinsic logarithmic scatter of 4.2±1.1. While a redshift determination using this method would have large uncertainties, it is nevertheless useful in instances where limited IR or millimeter photometry is available.

Figure 5. The probability distribution of the dust temperature for all sources in the SPT-DSFG sample. Though this probability distribution peaks at 40 K, the distribution is skewed toward warmer sources. The median of this distribution is at 52.4 K (solid line) and the inner quartiles of the distribution are given by the dashed lines. The adopted median for each fitted T dust (Table F1) is shown for comparison in red, and is described by the right axis. Figure 6. Photometric redshift compared with spectrosocpic redshift. We find good agreement between photometric and spectroscopic redshifts, though photometric redshift has large associated errors. Unity is shown as the dashed line.

Individual Source Properties

To demonstrate the extreme star-forming nature of the SPT DSFGs, we compare our sources to a flux-limited blank-field sample throughout this work. The original LABOCA ECDFS submillimetre survey (LESS) mapped the full Extended Chandra Deep Field South (ECDFS; Weiß et al. 2009). The survey area encompassed 0.5×0.5 deg 2 and detected a total of 126 DSFGs above a significance level of 3.7σ. The individual source detections were followed up with ALMA in Cycle 0 and the sample was thus called ALESS (Hodge et al. 2013;Swinbank et al. 2014). Photometry was obtained at 1.4 GHz, 870, 500, 350, 250 160, 100, and m 70 m in order to obtain the photometric redshifts (Swinbank et al. 2014). The MAGPHYS SED modeling code was used to obtain photometric properties (da Cunha et al. 2015). Spectroscopic redshifts were also obtained for a subset of these sources (Danielson et al. 2017). Because the ALESS sample is complete for sufficiently bright sources, we use it as an example of the classic unlensed 850 μm selected DSFGs with a flux-limited selection. Unless otherwise stated, the same formalism used for the SPT analysis was also applied to the ALESS sources throughout this work.

The distributions of T dust , apparent FIR luminosity, and redshift are shown for both the SPT sample and the ALESS sample (da Cunha et al. 2015) in Figure 8. The FIR luminosity and T dust values for both samples are derived from the SED fitting procedure described in Section 2.3.2. Because many of the ALESS sources do not have robust FIR detections due to their relative faintness and the Herschel/SPIRE confusion limit (Nguyen et al. 2010), their fitted SED values have large associated uncertainties when using the same SED fitting routine described in Section 2.3.2. In the left panel of Figure 8, the lensed SPT sources are offset from the largely unlensed ALESS sample, as expected from gravitational lensing. Gravitational lensing randomly samples sources from the background and increases the solid angle they subtend. This effect makes the lensed sources appear more luminous at a given T dust than an unlensed source. It would take a median magnification of ∼18 to make the median apparent L FIR for the SPT sources without lens models (3.88× L 10 13  ) on average intrinsically identical to the median L FIR for the ALESS sources

On the right panel of Figure 8, both samples exhibit a roughly constant FIR luminosity as a function of redshift, demonstrating the negative K-correction inherent in submillimeter/millimeter surveys of DSFGs. The approximate detection limits for both surveys and the Herschel/SPIRE confusion limit are also shown, assuming a standard template from Chary & Elbaz (2001). These limits illustrate the importance of selection wavelength and survey depth on the resultant redshift distribution of the sample. Though SPT is only sensitive to the most luminous sources, the detection threshold at the selection wavelength of 1.4 mm corresponds to a decreasing (apparent) luminosity at higher redshift. This is in contrast to the selection curve for the Herschel 500 μm selection, which corresponds to an increasing luminosity at higher redshift because 500 μm corresponds to rest wavelengths near or beyond the peak of the dust SED at higher redshifts.

In order to determine the intrinsic properties for each source, as in Figure 9, lens models are needed to calculate the magnification factor. Using the 39 SPT DSFGs with lens models (Spilker et al. 2016), the magnification factors for individual sources are used to calculate the intrinsic luminosities, dust mass, and SFR values, as described in Section 2.3.3. All of the derived properties can be found in Appendix F. It should be noted that for the remaining sources without lens models, the median magnification ( m á ñ m 870 m = 5.5) is adopted. This is a reasonable assumption, given that SPT sources with lens models from high resolution ALMA imaging presented in Spilker et al. (2016) were effectively drawn at random from the larger SPT sample. However, given that the median magnifications are merely an estimate, the intrinsic properties for these sources should be also treated as estimates. At the time of this publication, ALMA m 870 m imaging has been obtained for all sources and the construction of a complete lens model catalog is currently underway.

The intrinsic dust masses and FIR luminosities are shown in Figure 9 compared with T dust and redshift. We find that the SPT sources have a median intrinsic dust mass of 1.4(1)× M 10 9  , while ALESS sources have a median of 7.4(9)× M 10 8  . The median SPT intrinsic FIR luminosity is 7.1(5)× L 10 12  while ALESS sources exhibit a lower median = Ĺ 2.37 2 FIR ( ) L 10 12  . Thus, the SPT sources exhibit both higher luminosities as well as dust masses compared to the ALESS sources. These results are discussed in detail in Section 4.3.

Sources with Multiplicity

Due to the large (1′) beam size of the initial SPT selection, the SPT-DSFG catalog contains multiple sources that break up into multiple individual galaxies at the same redshift when observed at higher resolution (1′). As previously discussed in Section 3.2, the differences between the ALMA and SPT beam sizes can result in over-resolving the emission, causing the 3 mm flux to be underestimated in SED fitting. Though we treat these systems as individual DSFGs in this work, at least two of these sources are known to break into multiple components. SPT0311-58 (Strandet et al. 2016) was shown in Marrone et al. (2018) to be composed of two individual galaxies at z=6.900, a western source with an intrinsic = Ĺ L 33 10 IR 12  and an eastern source with an intrinsic = L 4.6 IR × L 10 12  . SPT2349-56 is among the most actively star-forming high-redshift protoclusters known (Miller et al. 2018) at z=4.3. This system is composed of at least 29 individual galaxies, six of which are at L IR > 3 × L 10 12  (Hill et al. 2020). Additional APEX/LABOCA observations show that SPT0348-62, SPT0457-49, SPT0553-50, and SPT2335-53 also have significant overdensities and lie at high redshift (3<z<7; G. Wang et al. 2020, in preparation). The multiplicities in these systems would potentially change the overall shape of the observed SPT-DSFG redshift distribution, favoring even higher redshifts. For the purpose of simplicity and homogeneity, we treat each of these systems as a single source in this analysis.

Selection Effects

The two major selection effects that impact our distribution are due to the long wavelength selection (1.4 mm) and the high m S 870 m flux cut, which selects primarily gravitationally lensed sources. Both of these effects are known to bias a redshift distribution toward higher redshifts. In this section, we compare the model of Béthermin et al. (2015b) with our observed distribution. This model was selected because it best described the distribution in Weiß et al. (2013) and takes into account both selection wavelength and lensing. These effects were previously shown to account for the redshift distribution for the subset of the SPT-DSFG sample published in Strandet et al. (2016) and Béthermin et al. (2015b). Some works (Blain et al. 2002;Greve et al. 2008;da Cunha et al. 2013;Staguhn et al. 2014) suggest that the CMB could make cold DSFGs at high redshifts difficult to detect. However, since the SPTselected DSFGs are warm, with a median dust temperature of T dust =52.4 K and minimum of T dust =20.9 K, this effect only becomes relevant at very high redshifts (z > 6-18).

In Figure 10, we examine the effect of a long wavelength selection on the observed redshift distribution using the Béthermin et al. (2015b) model. By applying a different selection wavelength on a realistic population of galaxies, we see the median shift from z∼2.6 at m 850 m to z∼2.8 at 1.4 mm. A Kolmogorov-Smirnov (K-S) test reveals that the observed SPT-DSFG distribution is distinct from the modeled data populations at levels of 7.7σ and 5.4σ for the 850 μm and 1.4 mm wavelength selections, respectively. Long wavelength selection alone does not explain the difference between the Béthermin et al. (2015b) model and the observed SPT data.

Based on models of the high-redshift DSFG population (e.g., Baugh et al. 2005;Lacey et al. 2010;Béthermin et al. 2012bBéthermin et al. , 2017;;Hayward et al. 2013;Lagos et al. 2019;Lovell et al. 2020), very few sources would be intrinsically bright The offset in apparent L FIR between the lensed (SPT) and unlensed (ALESS) sources is due to gravitational magnification. It should be noted that for many of the ALESS sources, the photometry is not well constrained, which is reflected in the size of the associated error bars. Right: the FIR luminosity as a function of redshift for the same sources. The teal, purple, and red curves represent the limiting FIR luminosity as a function of redshift corresponding to 3×the LESS limit (rms∼1 mJy), 3×the SPIRE 500 μm confusion limit (30 mJy), and the 3σ SPT 2.0 mm survey limit (3.9 mJy), given the SED model of this paper and assuming a 35 K dust temperature. Many of the sources discovered in the LESS survey split into multiple components, and were treated as individual sources in the ALESS survey and consequently fall below the plotted LESS survey limit. The SPT survey, with its longer wavelength selection, is more sensitive to sources at the highest redshifts (z>5) than Herschel.

enough to exceed our adopted flux density threshold at m 870 m (>25 mJy). Because of the high apparent luminosities, the SPT-DSFG sample was expected to consist almost solely of gravitationally lensed sources (Blain 1996;Negrello et al. 2007). Indeed, the lens models provided in Spilker et al. (2016) demonstrate that at least ~74% of the flux-limited SPT sample is strongly lensed, while a small fraction of sources are unlensed (∼5%) and the remaining sources are weakly lensed. Strong lensing preferentially selects sources at high redshift, biasing an intrinsic redshift distribution to higher redshifts. The probability of sources at z<1.5 undergoing strong lensing is heavily suppressed relative to sources at higher redshifts (z>4) (e.g., Hezaveh & Holder 2011). The 1.4 mm lensed model in Figure 10 peaked at lower redshift. The original Béthermin et al. (2012b) model is based on the stellar mass function and evolution of the main sequence and relies on SED template fitting. As these parameters become better constrained with more observational data from the SPT catalog and others, more realistic models will be possible in the future. The lensing analysis in Béthermin et al. (2015b) assumes that the sources do not undergo significant size evolution over cosmic time. Compact sources have a higher probability of being highly magnified than more diffuse sources (Hezaveh et al. 2012). Thus, size evolution versus redshift coupled with gravitational lensing could potentially bias a redshift distribution toward higher redshift, as discussed in Weiß et al. (2013). The fitted values of L FIR and T dust can then be used to determine the effective blackbody radius of the source via the modified blackbody function (e.g., Equation (1) of Spilker et al. 2016). The median effective blackbody radius for SPT-selected DSFGs is 1.0±0.1 kpc, which is consistent with the median radius found using lens modeling (1.04 ± 0.07 kpc ; Spilker et al. 2016). The sizes of the SPT DSFGs are also compatible with the sizes (0.3-3 kpc) observed in unlensed DSFGs (Ikarashi et al. 2015(Ikarashi et al. , 2017;;Simpson et al. 2015;Hodge et al. 2016). None of these observations show evidence of size evolution with redshift above z > 1.

Dust Temperature Evolution

With the first spectroscopically complete sample of DSFGs in hand and excellent FIR photometric coverage for all sources, we are now in a position to investigate whether the typical FIR properties of DSFGs evolve with redshift. For example, observations of "normal" star-forming galaxies (Magdis et al. 2012;Magnelli et al. 2014) and backwards evolution modeling (Béthermin et al. 2012b;Schreiber et al. 2018) imply that higher redshift sources should exhibit progressively warmer dust temperatures. One possible explanation for this effect would be that T dust is proportional to the L M IR dust ratio, which is proportional to the specific SFR (Narayanan et al. 2018;Liang et al. 2019;Ma et al. 2019). While these studies focused on dust temperature evolution in lower-SFR objects than our sample, determining whether DSFGs show temperature evolution has been difficult due to the highly uncertain photometric redshifts, confusion limits, and possible selection effects toward more luminous (and therefore, potentially warmer) galaxies in flux-limited surveys.

We examine this relationship by considering the dust temperatures from Section 3.2.1 versus the spectroscopic redshifts measured in Section 3.1, shown in Figure 11. In order to fit the data, we perform a maximum likelihood estimation assuming a linear model with an extra Gaussian variance term to account for intrinsic scatter. Using this method, the best fit line is given in Figure 11 as the dashed line, with a slope of 4.7±1.5K/z. However, the fitted intrinsic scatter is substantial (8.4 ±1.3 K). While there is evidence of temperature evolution, the large amount of intrinsic scatter implies the fit is also consistent with very shallow to no evolution with redshift.

Performing a least squares fit without the intrinsic scatter term will also give a nonzero linear dependence. However, such a fit is largely driven by the presence of significant outliers (e.g., SPT0452-50 with T dust =21 K and SPT0346-52 with T dust =79 K). One way to mitigate the presence of outliers in this approach is to consider the binned temperature distribution. Sources were binned with a bin width of Δz=1 and the median dust temperature was adopted, with bootstrapped errors representing the error on T dust . Though fitting these points again presents a positive slope of  1.9 0.5 K z, a one-sided reduced χ 2 test yields a p-value of 0.275, which is not statistically significant (< s 1 ). Moreover, hypothesis testing reveals that a line with no slope (p-value=0.184) is favored.

Since the dust temperature is highly sensitive to the exact SED fit function, we also examine the rest frame l peak . The wide sampling of the SED makes the l peak constraint insensitive to the SED fitting function. The l peak values used were extracted by fitting a modified blackbody fit with an added power law (individual fits shown in Figure 4 in blue). Repeating the linear fit with maximum likelihood estimation, we find that l peak also exhibits a nonzero slope of -5.9± m 1.3 m z. However, like dust temperature, any evidence of evolution with redshift is obscured by the 9.5± m 1.1 m intrinsic scatter. Like dust temperature, the presence of outliers (e.g., SPT0245-63, SPT0529-54, and SPT0452-50) could affect the significance of a slope, so we performed the same binned analysis. While the binned medians (2012a, 2015b, 2017). Top: the observed redshift distribution compared with modeled source catalogs with selections of 850 μm (blue), 1.4 mm (purple), and 1.4 mm (light blue) with gravitational lensing taken into account (cyan).

The filled triangles at the top of the panel represent the median redshifts of their respective samples. Bottom: the cumulative source density as a function of redshift. yielded a negative slope of -3± m 2 m z, it again had a low significance (<1σ).

Taken together, SPT sample's dust temperature and peak wavelengths show evidence of the thermal emission peak evolving with redshift, proportional to (1 + z), though it is not strong evidence. The intensely star-forming nature of our sample implies that our sources are starburst galaxies. Unlike the "normal" main sequence galaxies, starbursts have been shown to have almost constant radiation field intensities from z=0 to z ∼ 2 (Béthermin et al. 2015a;Schreiber et al. 2018). However, recent observations of GN20 (Cortzen et al. 2020) suggest that optically thick dust could obscure temperature evolution in [C I] at highredshift. While it is possible that luminous starburst galaxies could evolve in dust temperature versus redshift, strong evidence for temperature evolution is obscured by the large scatter on T dust .

SPT Sources are Extreme Ultraluminous IR galaxies/ Hyper-luminous IR galaxies

A consequence of modifying our SED fitting function to have the optically thick transition wavelength (λ 0 ) vary as a function of dust temperature is that the fitted dust temperatures obtained are systematically higher than if λ 0 was fixed to 100 μm, as in Greve et al. (2012). Fitting the ALESS photometry with this same SED fit gives a median dust temperature of T dust =54±4 K. While this dust temperature is significantly higher than what is commonly reported for DSFGs, which typically range between 30 and 40 K (Chapman et al. 2005;Greve et al. 2012;Swinbank et al. 2014;da Cunha et al. 2015), this is consistent with the values we obtain for the SPT DSFGs (T dust =52 ± 2 K).

We use the procedure outlined in Section 2.3.3 to calculate dust mass for the SPT sources and show the demagnified, intrinsic values throughout this section. Molecular gas mass was obtained by assuming a constant gas-to-dust ratio (taken to be 100; e.g., Casey et al. 2014). Because the FIR peak for the ALESS sample is not always well constrained, fitting the full IR range as described in Section 2.3.2 is not always possible. We instead adopt the values derived from MAGPHYS for comparison purposes (da Cunha et al. 2015), which takes into account the available photometry, as detailed in Swinbank et al. (2014). When a direct SED fit is possible, the two methods produce results that agree within the stated errors, with no bias. We derive a median SPT dust mass of 1.4× M 10 9  , compared with a median ALESS dust mass of 7.3× M 10 8  . This implies that the SPT-selected sources have a factor of 2 higher dust reservoirs than the ALESS sample, despite being earlier in cosmic history (Figure 9).

The SPT sources exhibit higher FIR luminosities than their unlensed analogs. Because of the relatively shallow SPT survey limit, shown in Figure 8, the SPT is sensitive to all but the brightest sources. Generally, these sources are lensed and have apparent luminosities that are an order of magnitude larger than unlensed sources. However, even after correcting for the gravitational magnification, the intrinsic SPT FIR luminosities range from 1.0-37× L 10 12  with a median of 7.1( 5 to ∼18, these discrepancies would also be resolved. We find this possibility unlikely, however, as the SPT sources with lens models from high resolution ALMA imaging presented in Spilker et al. (2016) were effectively drawn at random from the larger SPT sample of 81 sources presented here. Thus, SPT sources, in addition to being strongly lensed, are also intrinsically more luminous than standard DSFGs found in blank-field surveys.

Taking into account the full IR range of the SED, we find a median L IR =1.5(1)× L 10 13  , which corresponds to a median SFR of

´-M 2.3 2 10 yr 3 1 ( )

for the SPT sample. Figure 12 shows the gas mass (derived from M dust ) versus SFR (derived from L IR ). The SPT sample exhibits a correlation between the two quantities. By taking the ratio of gas mass to SFR, we calculate the gas depletion times (t depl ) shown in Figure 13. Because the magnification correction enters into both the SFR and gas mass calculations in the same way, it cancels out, assuming no differential magnification (Hezaveh et al. 2012). While (t depl ) in "normal" main sequence galaxies are known to evolve with redshift (Saintonge et al. 2013;Genzel et al. 2015;Tacconi et al. 2018), both the ALESS and SPT samples show comparable depletion times below the main sequence. We also find an evolution in t depl with redshift, which has already been observed for main sequence galaxies. Alternatively, this effect could also be the result of a strong Outliers (>2σ) are shown in gray, and removing them does not appreciably change the fitted slope. Binning the dust temperatures with bin widths z=1 also yields a linear dependence (yellow). Bottom: peak wavelength (l peak ) vs. redshift. The wide sampling of the SED makes the l peak constraint insensitive to the SED fitting procedure. The same analysis was repeated and the evolution with redshift is weaker, but still produces a nonzero slope. redshift evolution of metallicity, with the highest-redshift galaxies being the most metal poor. This would cause an underestimation of gas mass, which would also explain the lower depletion times. In order to address this properly, however, we would need stellar masses for our systems. Because of the high-redshift dust-obscured nature of these sources, this requires the use of future facilities such as the James Webb Space Telescope.

Though the dust temperatures of the SPT and ALESS sources are similarly distributed, with medians of T dust =52(2) K and T dust =54(4) K, respectively, SPT sources are on average more luminous than their unlensed analogs, even after accounting for the gravitational lensing. The simple blackbody relationship would imply that these sources should exhibit similar luminosities, if they were the same sizes. The SPT sources have a median radius of 1.0 1 kpc ( ) , while ALESS sources have a median of 0.46(4) kpc. The higher intrinsic luminosities could therefore stem from larger sizes, as discussed in Section 4.1.2. This fact, coupled with the larger dust masses, implies that the SPT sources could have larger gas reservoirs, which is supported by their systematically higher dust masses.

High-z Tail

As previously discussed, 850 μm selected samples of DSFG from blank fields (e.g., Chapman et al. 2005;Swinbank et al. 2014) peak at a lower redshift than the SPT sample, typically between z∼2.3 and 2.9. This discrepancy is attributed to selection effects, as discussed in the previous section (Section 4.1.2). While the SPT sample includes some sources equivalent to those selected in the blank-field surveys, it also includes the previously undiscovered high redshift and extreme tail of DSFGs. In addition to discovering the highest-redshift DSFG to date (SPT0311-58 at z=6.9; Marrone et al. 2018;Strandet et al. 2017), and with 12 sources at z > 5, the SPT-DSFG catalog contains roughly half of the DSFGs discovered at z > 5 (Walter et al. 2012;Riechers et al. 2013Riechers et al. , 2014Riechers et al. , 2017Riechers et al. , 2020;;Dowell et al. 2014;Rawle et al. 2014;Asboth et al. 2016;Fudamoto et al. 2017;Oteo et al. 2017;Strandet et al. 2017;Pavesi et al. 2018;Zavala et al. 2018;Jin et al. 2019).

A way to assess the relative rarity of SPT-selected sources is to consider the space density. However, because the majority of SPT-selected sources are lensed, the volume probed changes as a function of lensing and the effective volume is difficult to calculate. In order to estimate the survey volume, we consider the ALESS survey, which has a well-defined survey area and can be considered complete above a luminosity threshold of >3× L 10 12  , as shown in Figure 14. This particular luminosity threshold was selected because it retains most of the SPT sources, while ensuring there are sufficiently many ALESS sources to constitute an accurate comparison sample. We assume the two surveys have equal source density in the region 2.5<z<4.5 for >3× L 10 12  . We use this region to estimate the volume for the SPT survey and compare the two redshift distributions directly in Figure 15. The redshift region Figure 12. Gas mass vs. SFR for the demagnified SPT-selected sources and the ALESS sample. Though the SPT-selected sources have comparable amounts of gas mass to the ALESS sample, the SPT sources also experience higher SFRs. (2016) where available. The dashed horizontal line represents the L FIR above which the SPT and ALESS samples are complete, while the vertical lines correspond to 2.5 < z < 4.5, where the surveys have equal source density.

2.5<z<4.5 was chosen because it is the largest region covered by both surveys, which helps mitigate the effects of cosmic variance. The error bars for SPT are taken to the statistical errors from SPT added in quadrature with the ALESS errors in order to account for the scaling. Error estimation using bootstrapping produced similar results to the errors shown. Various combinations of luminosity threshold and redshift scaling region were tested, and the one chosen proved to be the best compromise between catalog completeness and constraining power.

Because the SPT distribution is both selected at long wavelength and spectroscopically complete, we have measured the high-redshift tail of the overall DSFG distribution. By summing the contributions in Figure 15, we infer that the expected number of DSFGs above z>4 is 2.6 7 ( ) × --10 Mpc 6 3 above 3 × L 10 12  . We can also predict the number of z>6 sources above 3× L 10 12  , which we take to be 0.95 1.0 ( ) × - -10 Mpc 7

3 , but this number should be used with caution, as it is based upon just one source (Marrone et al. 2018). To place this number in context, we compare with other observations (Coppin et al. 2009;Ivison et al. 2016;Cooke et al. 2018) and theoretical predictions (Hayward et al. 2013;Béthermin et al. 2017;Lagos et al. 2019;Lovell et al. 2020) in Table 2 and Figure 16. The stated number densities are taken to have an upper limit of z=8 and the errors are again calculated using the statistical errors from both the SPT and ALESS surveys in quadrature.

The Hayward et al. (2013) model combines a semi-empirical model with 3D hydrodynamical simulations and a 3D dust radiative transfer. Strong lensing is not included in the modeling and the model predicted dn/dz is determined using sources with S 1 mm > 1 mJy, consistent with the expected intrinsic flux densities of our sample. While the Hayward et al. (2013) model overpredicts the number of z>4 and z>5 galaxies, it does not predict any galaxies at z>6, which is contrary to our observations.

As previously discussed, the Béthermin et al. (2015b) model is an empirical model that starts from the observed FIR number counts. This model also includes the effects of magnification by strong lensing, so it can directly predict the dn/dz for the SPT sample. We find that the updated Béthermin et al. (2017) model is, in general, in good agreement with the SPT measurement.

The Shark semi-analytic model Lagos et al. (2018) is also compared in Table 2 andFigure 16. The number densities of galaxies with L IR > ´L 3 10 12  were computed using the simulated lightcone with an area of 107.9 deg 2 and covers a redshift range 0z8, and includes all galaxies in the model with stellar masses >10 M 8  and with a (dummy) AB magnitude (of mass-to-light-ratio=1) 32 mag. Note that here, L IR corresponds to the total luminosity that is reemitted by dust in Shark at wavelengths 1000 μm. While this model underpredicts the number of sources at z>4 and z>5, it exactly matches at z∼6. With SPT, we have measured the "high-redshift tail" of the distribution of luminous DSFGs. The SPT measurement for z>4 luminous sources lies within the range of model predictions and previous measurements and can already constrain models of galaxy evolution. The next-generation millimeter-wave surveys, combined with spectroscopic followup, will further improve this measurement and extend the redshift constraints out to z∼8.

Table 1

Empirical Counts in Each Redshift Bin for the Final SPT dN/dz Distribution Redshift Bin dN/dz

improved the SED fitting and provided better photometric redshifts than the original method presented in Greve et al. (2012). In the absence of observationally expensive spectroscopic data, photometric redshifts can give a rough estimate of redshift. They can also serve to identify the correct redshift when multiple spectroscopic solutions are available. In instances of limited photometric data, we also demonstrated that the m m S S 350 m 870 m ratio can also serve as a rough redshift indicator. Taken together, these results demonstrate the utility of photometric data in determining redshift and can serve to guide prioritization of targets for follow-up, even when limited IR/millimeter data is available. 3. The dust temperatures obtained in Section 3.2.1 suggest a mild correlation between redshift and dust temperature (e.g., Cooke et al. 2018), but this result is not yet statistically significant. 4. Due to gravitational lensing, the apparent luminosity of SPT sources is an order of magnitude greater than the unlensed DSFG population. However, the intrinsic luminosities are still significantly higher than the unlensed 850 μm selected population, with higher dust masses and lower depletion times for the SPT sources. Thus, even after accounting for magnification due to strong gravitational lensing, SPT-selected DSFGs are more extreme than those typically found in blank-field surveys. 5. Even though the SPT survey encompassed a contiguous 2500 deg 2 area and the sample is flux limited, the gravitational lensing intrinsic to the SPT sample makes it difficult to estimate the total volume probed by the survey. Scaling the SPT sample to that of the ALESS sample allows us to estimate the source density of luminous DSFGs at high redshift. From our data and analysis, we estimate 4.7±5.1 deg 2 and 0.95±1× 10 -7 Mpc -3 DSFGs with L FIR >3×10 12 L  at z>6. 6. While the 3 mm data was obtained in this paper, it will serve as a starting point for several papers currently in preparation. The 3 mm line fluxes from high-J CO lines will be combined with the line fluxes from low-J CO observations with ATCA. Spectral line energy distributions will be fitted for all sources, which is currently in preparation. The composite 3 mm spectrum, shown in the bottom panel of Figure 3, will be further examined in C. Reuter et al. (2020, in preparation). Finally, ALMA m 870 m imaging has been obtained for the complete SPT sample, enabling lens models to be constructed for the complete catalog of SPT DSFGs. The fitted lens models and their associated properties, including source size and magnifications, will be studied in C. Reuter et al. (2020, in preparation). 7. Though this SPT-selected redshift catalog is now complete, the next generation of CMB experiments are already in operation. SPT-3G (Benson et al. 2014) has already begun the process of surveying a 1500 deg 2 field with a sensitivity of ∼5 times the original SPT-SZ survey and will detect two orders of magnitude more sources, extending all the way to z > 8. CMB-S4 (Abazajian et al. 2019), still in the design phase, is projected to find an order of magnitude more DSFGs than SPT-3G. The next generation of millimeter-wave surveys with spectroscopic follow-up will be key to uncovering the the earliest stages of the dust-obscured universe. Future upgrades (e.g., Mroczkowski et al. 2019) to significantly extend the instantaneous bandwidth of the ALMA receivers will be crucial for ascertaining the redshifts for the next generation of millimeter surveys.

We thank the anonymous referee whose careful reading and insightful suggestions greatly improved this paper. We thank A. M. Swinbank, C. Lagos, and C. Lovell for useful

C.2. ALMA Targeted Line Scans

A total of 17 sources in the ALMA 3 mm spectral scans that contained a single CO line were confirmed through targeted line scans. Though a photometric redshift could be used to further identify which CO transition was likely detected, the large error bars associated with photometric redshifts allow for some ambiguity. Using the existing CO line detections and the photometric redshift as a guide, we are able to specifically target CO lines that would be visible for degenerate redshift solution. The sidebands were configured such that an observation would yield at least one CO line detection for each source at a likely redshift option. The first 13 sources were observed between 2018 December and 2019 January (project ID: 2018.1.01254.S), using ALMA Band 4 (125-163 GHz). The line scans are shown in Figure C2. An additional four sources were confirmed between 2019 December and 2020 January (project ID: 2019.1.00486.S), using ALMA Bands 4 and 5 (163-211 GHz). The line scans are shown in Figure C3.

For both sets of observations, the flux density and bandpass calibrations were based on observations of J0519-4546, J0538 -4405, J1924-2914, J2056-4714, J2258-2758, and J2357 -5311. The phase calibration was determined using nearby quasars. Between 42 and 48 antennas were used to complete the observations, resulting in synthesized beam sizes between 3.1″ and 0.8″. The observing time for each science block ranged from 5-13 minutes on-source, excluding overheads, with average single-sideband system temperatures of T sys = 56-86 K. The data were processed using the CASA package. The images were created using natural weighting and the subsequent spectra were created with a channel width of 16 MHz. The average noise per channel in the resulting spectra is -0.4 0.6 mJy beam 1 -.

Appendix E All ALMA 3 mm Spectroscopic Line Observations

A summary of all spectroscopic data obtained for the final SPT-selected DSFG sample is given below in Table E1. Blind 3 mm scans were conducted for all 81 sources across ALMA Cycles 0, 1, 3, and 4. The Cycle 0 data was published in Weiß et al. (2013) and Cycle 1 data was published in Strandet et al. (2016). The remaining cycles are presented for the first time in this work. While the spectra are given in Figure B1, the lines detected can be found in the "Lines from 3 mm Scans" column. Any ancillary spectroscopic data obtained is presented in the "New lines and comments" column. All of the combined observational efforts have yielded secure redshifts for the complete flux-limited sample of 81 sources from the 2500 deg 2 SPT survey. Of these, 79 sources had multiple spectroscopic lines, detected either solely from the 3 mm window or with ancillary spectroscopic observations. Only two of the redshifts are based on a single line, but with a robust line identification based on our analysis of the dust temperature distribution. All together, this is the largest and most complete collection of spectroscopic redshifts for high-redshift DSFGs obtained so far at millimeter wavelengths. 2.9(1.9) SPT2353-50 5.578 5.7(0.9) 66( 16) 81(10) 5.5 a,d 0.5(0.1) Notes. Estimates of dust temperatures and FIR luminosities are derived from the modified blackbody fits described in Section 2.3.2. Using these parameters, the dust masses and SFRs are then derived using Equation (2) and Equation (3), respectively. While not explicitly given in the table, the total IR luminosity can be obtained using the SFR, as shown in Equation (3). The FIR luminosities, SFRs and dust masses were corrected for gravitational amplification (μ) according to Table D1 of Spilker et al. (2016) where such models exist. Where multiple source components exist, we use a flux-weighted average to estimate μ. Though the intrinsic values are presented above, the apparent values can be obtained by multiplying by the magnification factor. This table is available in machine-readable format at https://github. com/spt-smg/publicdata. a The sources that have not been modeled have parameters that have been corrected by median gravitational amplification m á ñ =