The frequencies and spacings of these lines unambiguously place the galaxy at a redshift of z = 6.900 (2), which corresponds to a cosmic age of 780 Myr (using cosmological parameters 10 of Hubble constant H 0 = 67.7 km s -1 Mpc -1 , matter density Ω m = 0.309 and vacuum energy density Ω Λ = 0.691). An elongated faint object is seen at optical and near-infrared wavelengths, consistent with a nearly edge-on spiral galaxy at z = 1.4 ± 0.4 that acts as a gravitational lens for the background source (see Methods section 'Modelling the SED'; here and elsewhere the error range quoted corresponds to a 1σ uncertainty). Together, these observations indicate that SPT0311-58 is the most distant known member of the population of massive, infrared-bright but optically dim, dusty galaxies that were identified from ground-and space-based wide-field surveys 11 .

The far-infrared emission from SPT0311-58 provides an opportunity to study its structure with little confusion from the foreground galaxy. We conducted ALMA observations at about 0.3″ resolution at three different frequencies (see Methods): 240 GHz, 350 GHz and 420 GHz, corresponding to rest-frame wavelengths of 160 µm, 110 µm and 90 µm. The observations at 240 GHz include the 158-µm fine-structure line of ionized carbon ([C ii]) and those at 420 GHz the 88-µm fine-structure line of doubly ionized oxygen ([O iii]). The 160-µm continuum and the [C ii] and [O iii] line emission maps of the source are shown in Fig. 1. Two emissive structures are visible in the map, denoted SPT0311-58 E and SPT0311-58 W, which are separated by less than 2″ on the sky before correction for gravitational deflection. Although the morphology of SPT0311-58 E and SPT0311-58 W is reminiscent of a lensing arc (SPT0311-58 W) and counter-image (SPT0311-58 E), the [C ii] line clarifies the physical situation: SPT0311-58 E is separated from the brighter source SPT0311-58 W by 700 km s -1 and is therefore a distinct galaxy.

Lens modelling of the 160-µm, 110-µm and 90-µm continuum emission from SPT0311-58 was performed using a pixelated reconstruction technique 12 (Fig. 1c, Extended Data Fig. 5, Methods section 'Gravitational lens modelling'). Its structure and lensing geometry is consistent between the observations, and indicates that the two galaxies are separated by a projected (proper) distance of 8 kpc in the source plane. SPT0311-58 E has an effective radius of 1.1 kpc, whereas SPT0311-58 W has a clumpy, elongated structure that is 7.5 kpc across. The (flux-weighted) source-averaged magnifications of each galaxy and of the system as a whole are quite low (µ E = 1.3, µ W = 2.2, µ tot = 2.0) because SPT0311-58 W is extended relative to the lensing caustic and SPT0311-58 E is far from the region of high magnification. The same lensing model applied to the channelized [C ii] data reveals a clear velocity gradient across SPT0311-58 W, which could be due to either rotational motions or a more complicated source structure coalescing at the end of a merger.

Having characterized the lensing geometry, it is clear that the two galaxies that comprise SPT0311-58 are extremely luminous. Their intrinsic infrared (8-1,000 µm) luminosities have been determined from observations of rest-frame ultraviolet-to-submillimetre emission (see Methods section 'Modelling the SED') to be L IR = (4.6 ± 1.2) × 10 12 L  and L IR = (33 ± 7) × 10 12 L  for SPT0311-58 E and SPT0311-58 W, respectively, where L  is the luminosity of the Sun. Assuming that these sources are powered by star formation, as suggested by their extended far-infrared emission, these luminosities are unprecedented at z > 6. The implied (magnification-corrected) star-formation rates are correspondingly enormous-(540 ± 175)M  yr -1 and (2,900 ± 1,800)M  yr -1 , where M  is the mass of the Sun-probably owing to the increased instability associated with the tidal forces experienced by merging galaxies 13 . The components of SPT0311-58 have luminosities and star-formation rates similar to the other massive, z > 6 galaxies identified by their dust emission, including HFLS3 (z = 6.34), which has a star-formation rate of 1,300M  yr -1 after correcting for a magnification factor 14 of 2.2, and a close quasargalaxy pair 15 at z = 6.59, the components of which are forming stars at rates of 1,900M  yr -1 and 800M  yr -1 , respectively. However, unlike the latter case, there is no evidence of a black hole in either source in SPT0311-58.

Unlike any other massive dusty source at z > 6, the rest-frame ultraviolet emission of SPT0311-58 E is clearly detectable with modest integration by the Hubble Space Telescope. The detected ultraviolet luminosity (L UV = (7.4 ± 0.7) × 10 10 L  ) suggests a star-formation rate of only 13M  yr -1 , 2% of the rate derived from the far-infrared emission, consistent with SPT0311-58 E forming most of its stars behind an obscuring veil of dust. The inferred stellar mass for this galaxy (see Methods section 'Modelling the SED') is (3.5 ± 1.5) × 10 10 M  . Although no stellar light is convincingly seen from SPT0311-58 W, the absence of rest-frame ultraviolet emission is probably explained by heavy dust obscuration and is not unusual 16 . Although SPT0311-58 E is the less massive of the two components, even it is rare among ultravioletdetected galaxies at z ≈ 7. Such galaxies are found in blank-field surveys to have a sky density of just one per 30 square arcminutes 17 . The far-infrared continuum and line emission of SPT0311-58 E and SPT0311-58 W (Fig. 1d-f) imply substantial differences in the physical conditions in these objects. Compared to SPT0311-58 W, SPT0311-58 E has a higher ratio of [C ii] line emission to 160-µm continuum emission and a much larger luminosity ratio between [O iii] and [C ii]. The [O iii] emission is much more luminous in SPT0311-58 E, with most of SPT0311-58 W (excluding the southern end) showing no emission at all. Because the formation of O ++ ions requires photons with energies of more than 35.1 eV, this line arises only in ionized regions around the hottest stars and near active galactic nuclei 18 . It is unlikely that active galactic nuclei are the origin of the [O iii] line in SPT0311-58 E, because the continuum and line emission both extend across most of the galaxy rather than being concentrated in a putative nuclear region. Observations of [O iii] 88-µm emission in actively star-forming galaxies at low 7,19 and high 20 redshift have found that the line luminosity ratio between [O iii] and [C ii] increases as gas metallicity decreases. The ultraviolet photons capable of forming O ++ have a longer mean free path in a lowermetallicity interstellar medium than in a higher-metallicity one, and the electron temperature remains higher for the same ionizing flux, both of which favour increased [O iii] emission 21 . The difference in the [C ii] line-to-continuum ratio may result from multiple effects: the known suppression [22][23][24] of the [C ii]-to-L IR ratio in regions of increased star-formation surface density (higher in SPT0311-58 W), and the increased [C ii]-to-L IR ratio in star-forming galaxies of lower metallicity 7 . Whether SPT0311-58 E (or the southern edge of SPT0311-58 W, which is similar to SPT0311-58 E in these properties) has a more primordial interstellar medium than does the bulk of SPT0311-58 W can be tested with future observations. The masses of the components of SPT0311-58 are remarkable for a time only 780 Myr after the Big Bang. In Fig. 2 we compare SPT0311-58 to objects at z > 5 for which we have estimates of dust mass (M dust ) or total gas mass (M gas ). For SPT0311-58, the best constraints on both of these quantities come from the joint analysis 6 of its far-infrared continuum and line emission, specifically the rotational transitions of carbon monoxide and neutral carbon. Here we have divided these masses between the two galaxies according to the lensing-corrected ratio of dust continuum emission (6.7) that we determined from our three high-resolution ALMA continuum observations because the dust continuum luminosity is roughly proportional to the dust mass. The corresponding dust and gas masses for SPT0311-58 W are M gas = (2.7 ± 1.7) × 10 11 M  and M dust = (2.5 ± 1.6) × 10 9 M  , and for SPT0311-58 E are M gas = (0.4 ± 0.2) × 10 11 M  and M dust = (0.4 ± 0.2) × 10 9 M  . The gas mass can also be estimated using the carbon monoxide luminosity, although the conversion between luminosity and gas mass in this optically thick line is known to vary substantially depending on many factors, including star-formation intensity and metallicity 25 . Taking the observed 6 luminosity in the J = 3-2 line of carbon monoxide, converting it to J = 1-0 under the conservative assumption of thermalized emission, and connecting luminosity to mass using a standard value of α CO = 1.0M  (K km s -1 pc 2 ) -1 , we derive M gas = (6.6 ± 1.7) × 10 10 M  for SPT0311-58 W and M gas = (1.0 ± 0.3) × 10 10 M  for SPT0311-58 E. The gas mass of SPT0311-58 W is well above those of all of the known galaxies at z > 6, that is, during the first approximately 900 Myr of cosmic history. SPT0311-58 highlights an early and extreme peak in the cosmic density field and presents an opportunity to test the predictions for the growth of structure in the current cosmological model. The mass of the dark-matter halo that hosts SPT0311-58 is uncertain, but can be estimated in several ways. For most massive star-forming galaxies 26,27 the gas mass represents the dominant component of baryons that have cooled and assembled at the centre of the dark-matter halo. In this case, for the lower (α CO -based) estimate of gas mass, the cosmic baryon fraction 10 f b = 0.19 places a hard lower bound on the total halo mass of 4 × 10 11 M  . A less conservative assumption incorporates the knowledge, based on observations across a wide range of redshifts, that only a fraction of the baryons in a dark-matter halo (less than one-quarter, M b /M halo = 0.05; see figure 15 of ref. 3) are destined to accrue to the stellar mass of the central galaxy 3 . In this case, a total halo mass of (1.4-7.0) × 10 12 M  is implied, depending on which estimate of gas mass is adopted. To understand the rareness of the dark-matter halo that hosts SPT0311-58, we calculate curves that describe the rarest haloes that should exist in the Universe at any redshift 28 . In Fig. 3, we show the halo masses that are inferred for many high-redshift galaxies, using the same methods for converting gas mass to halo mass as described above. We find that SPT0311-58 is indeed closest to the exclusion curves and therefore marks an exceptional peak in the cosmic density field at this time in cosmic history.

We have found a system of massive, rapidly star-forming, dusty galaxies at z = 6.900, the most distant galaxies of this type discovered so far. Two compact and infrared-luminous galaxies are seen, separated by less than 8 kpc in projection and 700 km s -1 in velocity, probably in the process of forming one of the most massive galaxies of the era. Even before coalescence, the larger galaxy in the pair is more massive than any other known galaxy at z > 6. Although the discovery of such M gas (M )

Figure 2 | Mass measurements for high-redshift galaxies. Dust masses (M dust ) are taken from the literature, as described in Methods. Gas masses (M gas ) are primarily derived from observations of various rotational transitions of carbon monoxide, with previously reported line luminosities converted to molecular gas masses under standardized assumptions (see Methods). The comparison sample (small filled circles) includes three classes of object: dusty star-forming galaxies (DSFGs; red), quasars (QSOs; blue) and Lyman-break galaxies (LBGs; green). These objects are typically selected by far-infrared emission (DSFGs) or by optical or infrared emission (QSOs and LBGs). Three additional DSFGs-SPT0311-58 E (yellow pentagons), SPT0311-58 W (yellow hexagons) and HFLS3 (red squares)-have extensive photometry and line measurements, which enable more sophisticated estimates of their dust and gas masses 6,29 from a combined analysis of the dust and carbon monoxide line emission. For these objects we also show masses derived under a simpler assumption as open symbols (for SPT0311-58 the methods give very similar answers for M dust ). Error bars represent 1σ uncertainties.

a system at this high redshift and in a survey that covered less than 10% of the sky is unprecedented, its existence is not precluded by the current cosmological paradigm. Figure 3 | Halo masses for rare, high-redshift, massive galaxies. The mass of the dark-matter halo (M halo ; defined at an density of 200 times the mean density of matter in the Universe) is inferred for galaxies in the first 2 Gyr after the Big Bang (see Methods). These masses present a range of lower limits, from the most conservative assumption (lower bars) that all baryons in the initial halo have been accounted for in the molecular gas mass to the observationally motivated assumption (upper triangles) that the baryonic mass (M b ) in gas is a fixed ratio of the halo mass M b /M halo = 0.05, calibrated through a comparison 3 of simulations and observations spanning z = 0-8. The most massive haloes that are expected to be observable 28 within the whole sky (dotted line), within the 2,500-deg 2 area of the South Pole Telescope (SPT) survey (dashed line) and within the subset of that area that is magnified by a factor of two or more (solid line) are also plotted as a function of redshift. As SPT0311-58 E and SPT0311-58 W reside within the same halo, they are combined for this analysis. As in Fig. The redshift of SPT0311-58 and the 3-mm continuum flux density were determined from an 84.2-114.9-GHz spectrum assembled from five separate tunings in ALMA band 3 under ALMA Cycle 3 project 2015.1.00504.S. The observing strategy has been used to discover the redshifts of more than 50 SPT dusty sources, and further details on the redshift coverage are provided in previous works 30,31 . Data were taken on 2015 December 28 and 2016 January 2 in ALMA configuration C36-1 (baseline lengths of 15-310 m) using 34 and 41 antennas, respectively. The resulting image has a resolution of 3.3″ × 3.5″, although there is spatial information on finer scales that allows us to estimate flux densities separately for the E and W sources, which are separated by about 2″. Further details of the analysis are provided elsewhere 6 .

ALMA observed SPT0311-58 a second time under project 2015.1.00504.S in band 7 (LO = 343.48 GHz) to produce a continuum image suitable for gravitational lens modelling. Similar observations were used to produce lens models of SPT sources in previous cycles 24,32 . The observations were performed with 41 antennas in the C40-4 configuration, providing 15-770-m baselines. The resulting image has an angular resolution of 0.3″ × 0.5″, although, because it lacks any spectral lines, it was found to be insufficient to provide an unambiguous determination of the lensing configuration.

The ALMA Cycle 4 project 2016.1.01293.S was intended to follow up on the discovery of this very distant source through spectroscopic observations. The 158-µm line of [C ii] was observed on 2016 November 3 in ALMA configuration C40-5, which provided baseline lengths of 18-1,120 m. This provides the primary imaging for this work, because it yielded an extremely sensitive detection of the [C ii] line and continuum structure at high resolution.

A final observation was obtained in ALMA band 8 (LO = 423.63 GHz), in configuration C40-4 (baselines 15-920 m). The observations were repeated in four segments to yield the required integration time. The resulting data have 0.2″ × 0.3″ resolution. These data provide a final spatially resolved continuum observation, at 90-µm rest-frame wavelength, along with spectroscopic images of the 88-µm line of [O iii]. The ALMA continuum images are shown in Extended Data Fig. 1. Spitzer infrared imaging. Infrared observations of SPT0311-58 were acquired with the Infrared Array Camera (IRAC) instrument 33 on the Spitzer Space Telescope as a part of Cycle 24 Hubble Space Telescope (HST) programme 14740. The observations consisted of 95 dithered 100-s exposures on-source in both operable IRAC arrays at 3.6 µm and 4.5 µm. A large dither throw was used. The dataset thus has sufficiently high redundancy to support our standard reduction procedure, which involves constructing an object-masked median stack of all 95 exposures in each band and then subtracting the median stack from the raw frames to compensate for bad pixels not automatically masked by the pipeline and to remove gradients in the background. After these initial preparatory steps, the background-subtracted exposures were combined in the standard way 34 with IRACproc 35 and MOPEX to create mosaics with 0.6″ pixels. The mosaics achieved an effective total integration time of about 9,000 s after masking cosmic rays and other artefacts. Two flanking fields were covered to the same depth but separately, each in one IRAC passband.

Photometry was performed on the mosaics using Source Extractor 36 in dual-image mode after trimming to exclude the flanking fields and unexposed areas. The lens galaxy associated with SPT0311-58 was well detected with no evidence for saturation or even nonlinear detector behaviour. During this process background and object images were generated and inspected to verify that Source Extractor performed as expected and generated valid photometry. HST imaging. SPT0311-58 was observed for five orbits of HST imaging with ACS and WFC3/IR in Cycle 24 (PID 14740) to determine the morphology of the foreground lens and to better constrain the spectral energy distribution (SED) of both the lens and source. All observations were acquired on 2017 April 30. The ACS imaging consists of a single orbit divided between the F606W and F775W filters. Exposure times are 844 s and 1.5 ks, respectively. Four orbits of WFC3/IR observing was split evenly between the F125W and F160W filters. Although the nominal exposure times are 5.6 ks, a subset of the data in both filters was compromised by substantial contamination from scattered earthlight. We reprocessed the imaging to remove contaminated data, resulting in final exposure times of 4.9 ks in each band. Gemini optical and infrared imaging and spectroscopy. With the Gemini Multi-Object Spectrograph 37 (GMOS) of Gemini-South, we obtained deep i and z images of SPT0311-58 (PID GS-2015B-Q-51) on 2016 January 29 and 31. The instrument consists of three 2,048 × 4,176 pixel CCDs, separated by two 6.46″ (80 pixel) gaps, with a scale of 0.0807″ per pixel. The field of view of the GMOS camera is 5.5′ × 5.5′. Our images were taken under photometric conditions and using a 2 × 2 binning, which gives a scale of 0.161″ per pixel. The total integration times were 3,600 s for the i band and 6,600 s for the z band, with average seeing conditions of 1.3″ and 1.0″ in the i and z bands, respectively. The resulting 5σ point source depths were i AB = 25.2 and z AB = 25.0. SPT0311-58 was observed using the Facility Near-Infrared Wide-Field Imager and Multi-Object Spectrograph for Gemini (FLAMINGOS-2) 38 at the Gemini-South Observatory on the nights of ut 2016 September 23 and 2017 February 06, under PID GS-2016B-Q-68. The instrument was used in imaging mode, with 0.181″ pixels, and yielded an unvignetted circular field of view of approximately 5.5′ diameter. Our observing sequence for the survey consisted of a randomly ordered dither pattern, with 15″ offsets about the pointing centre. This pattern was repeated until the required total exposure time was achieved. The individual K s -band exposure time was set at 15 s in the first observation and 10 s in the second observation, yielding a typical background sky level in the K s band of 10,000-12,000 counts (detector nonlinearity can be corrected to better than 1% up to 45,000 counts). These counts ensure that 2MASS stars with K s > 13 do not saturate and can be used for photometric calibration. The data were reduced using the Python-based FLAMINGOS-2 Data Pipeline, FATBOY 39,40 . In brief, a calibration dark was subtracted from the dataset, a flat field image and a bad pixel map were created, and the flat field was divided through the data. Sky subtraction was performed to remove small-scale structure, with a subsequent low-order correction for the large-scale structure. Finally, the data were aligned and stacked. The seeing conditions averaged 0.7″ in the final image comprising 44 min of integration, reaching K s,AB = 23.6 at 5σ. Spectroscopy was obtained with the GMOS-S instrument on the nights of ut 2016 February 1 and 2 (PID GS-2016B-Q-68) using the 1″-wide long slit at a position angle -10° east of north and the instrument configured with the R400 grating and 2 × 2 detector binning. For a source that fills the 1″ slit this set-up results in a spectral resolution of about 7 Å. The observations were spectrally dithered, using two central wavelength settings (8,300 Å and 8,400 Å) to cover the chip gaps. The data comprise a series of individual 900-s exposures, dithering the source spatially between two positions (' A' and 'B') along the slit in an ABBA pattern, repeated four times, two at each central wavelength setting. The total integration time is 4 h. A bright foreground object was positioned along the slit midway between the acquisition star and SPT0311-58, providing an additional reference point for locating traces along the slit.

The spectra were reduced, beginning with bias subtraction and bad pixel masking using the IRAF GMOS package provided by Gemini. The individual chips were combined into a single mosaic for each exposure and the mosaicked frames were then sky-subtracted by differencing neighbouring A-B exposure pairs; this method resulted in nearly Poisson noise, even under the numerous bright sky lines. A flatfield slit illumination correction was applied and a wavelength calibration derived for each mosaic. The two-dimensional spectrum was created by median-combining the individual exposure frames.

The spectrum shows a faint continuum beginning above 9,000 Å at the location of SPT0311-58. A one-dimensional extraction of the faint trace yields no reliable redshift measurement, but is consistent with the redshifted 4,000-Å break that is expected for the foreground galaxy at z ≈ 1.4. Calibrated against the nearby R = 16.4 star spectrum we find no flux at the expected location of Lyα redshifted to z = 6.900 (about 9,600 Å) down to a 3σ flux limit of 3.0 × 10 -17 erg s -1 cm -2 for a emission line 500 km s -1 wide. Image de-blending. At the position of SPT0311-58, our optical and infrared images (Extended Data Fig. 2) show a prominent lower-redshift galaxy that is responsible for lensing the W source, and the HST images, which have the highest resolution, show direct stellar emission from the E source (Extended Data Fig. 3). To extract reliable photometry for SPT0311-58 E, particularly in the lowresolution Spitzer images that cover the rest-frame optical, and to search for emission from the W source underneath the lens galaxy, we must model and remove the lens emission. We follow procedures similar to those used previously 41 , using the HST/WFC3 images as the source of the lens galaxy model to de-blend the IRAC image. The foreground lens can be fitted with a single Sérsic profile with an index n = 1.77. As seen in Extended Data Fig. 4, there is no clear rest-frame ultraviolet emission from SPT0311-58 W in the HST bands after removal of the lens model. To remove the lens from the IRAC image, the WFC3 model is convolved with the IRAC point spread function and then subtracted from the 3.6-µm and 4.5-µm images. Residual emission is seen near the positions of the E and W sources. Unfortunately, because SPT0311-58 W lies right on top of the lens, the residuals are extremely susceptible to image de-convolution errors and we do not believe

the Spitzer/IRAC fluxes to be reliable. By contrast, SPT0311-58 E is one full IRAC resolution element, 1.7″, from the lens centroid, and we consider the residual emission at this position to be usable in our subsequent analyses. Images of the model and residuals are provided in Extended Data Fig. 4 and the resulting photometry is provided in Extended Data Table 2. Gravitational lens modelling. Gravitational lens modelling of SPT0311-58 was performed using two different codes which model the source-plane emission in different ways. Both codes fit to the visibilities measured by ALMA or other interferometers directly to avoid the correlated noise between pixels in inverted images. In each, the lens galaxy is modelled as a singular isothermal ellipsoid, and posterior parameter distributions are sampled using a Markov chain Monte Carlo technique, marginalizing over several sources of residual calibration uncertainty (such as antenna-based phase errors).

Initial lens models were created using the visilens code, which is described in detail elsewhere 24 . The source plane is modelled as one or more elliptical Sérsic profiles. Because of the simplicity of this source-plane representation, the code is able to sample large and complex parameter spaces quickly. The continuum emission at 160 µm, 110 µm and 90 µm was modelled with four Sérsic components, one for SPT0311-58 E and three for SPT0311-58 W. These models leave approximately 8σ peak residuals in the 160-µm and 90-µm data, which both reaching peak signalto-noise ratios of more than 150.

After determining the lens parameters using visilens, we used the best-fitting values as initial input to a pixelated reconstruction code 12 . This code represents the source plane as an array of pixels, rather than an analytic model, and determines the most probable pixel intensity values for each trial lens model while imposing a gradient-type regularization 42 to avoid over-fitting the data. For each dataset, we fit for the strength of this regularization. At 160 µm and 90 µm we re-fit for the lens model parameters and compare to the visilens models as a test of the robustness of the lens modelling. Within each code, the best-fitting lens parameters at the two independent wavelengths are consistent to within 10%. Further, the lens parameters and the source structure are consistent between the two independent codes, with intrinsic source flux densities, sizes and magnifications that agree to within 15%. The increased freedom in the source plane afforded by the pixelated reconstruction means that the lens parameters are not independently well constrained by the 110-µm data, which have lower signal-to-noise ratio and spatial resolution. For these data, we apply the lensing deflections determined from the other two datasets to reconstruct the source-plane emission. The pixelated reconstructions of the three continuum wavelengths are shown in Extended Data Fig. 5.

The channelized [C ii] line is modelled using the same pixelated reconstruction technique, using 39 consecutive channels of 40 km s -1 width, each with a peak signal-to-noise ratio ranging from 9 to 34. For each channel, we apply the lensing deflections from the best-fitting model of the 160-µm data, which were observed simultaneously. We fit for the strength of the source-plane regularization 12,42 at each channel, which varies across the line profile as some velocities (such those multiply imaged from -280 km s -1 to +80 km s -1 ) experience higher magnification than others (such as the entire eastern source at >+560 km s -1 ). The models of each [C ii] channel are represented in Extended Data Fig. 6.

We determine the source magnifications using the 90-µm pixelated model, in which the E source is detected at the highest signal-to-noise ratio and so the effects of varying the aperture used to measure the intrinsic flux density are minimized. Because the source-plane morphology is very similar between the three continuum wavelengths, the magnification is also essentially identical between them. We find flux-weighted, source-averaged magnifications for the E source, the W source and the system as a whole of µ E = 1.3, µ W = 2.2 and µ tot = 2.0, respectively. These magnifications are substantially lower than the median magnification of 5.5 within the sample of 47 SPT-discovered dusty galaxies 24 for which we have data adequate to construct lens models or to conclude that sources are unlensed. In this case the low magnification is a consequence of the low mass of the lensing halo, which is typically expressed as an 'Einstein' radius θ E . The lens model for this source indicates θ E = 0.29″, which is around the 10th percentile for SPT lensed sources 24 , and the background source is both much larger than and offset from the regions of highest magnification. A large portion of the source is therefore only weakly magnified and the source-averaged values are low.

Finally, we also construct a lens model of the 95-GHz ALMA data (rest-frame 380 µm; Extended Data Table 1). Because the spatial resolution of these data are low (3.5″), we model them using only the visilens code, which is more suited to low-resolution data. We allow only the lens parameters and source structural parameters (such as position and radius) to vary within the ranges determined from the higher-resolution 160-µm, 110-µm and 90-µm continuum data, leaving only the flux densities of the E and W sources as free parameters. This modelling indicates that essentially all of the observed 380-µm emission can be ascribed to the W source, with the E source 'detected' at about 1σ.

In addition to the ALMA data, we use Herschel photometry 6 to constrain the SED of SPT0311-58 E and SPT0311-58 W to rest-frame 30 µm (250 µm observed). The resolution of Herschel SPIRE is not adequate to separate the two components, so we divide the total flux density observed in the three SPIRE bands between the E and W sources according to the ratios observed in the ALMA bands. These photometric points are then corrected for the continuum magnification derived from the ALMA data and used in the SED modelling described below. The total and intrinsic flux densities are reported in Extended Data Table 3. Modelling the SED. In Extended Data Fig. 7 we present the SEDs of SPT0311-58 E, SPT0311-58 W and the foreground lens galaxy.

A photometric redshift for the lens is calculated with EAZY 43 using the data in Extended Data Table 2. The resulting redshift is 1.43, with a 1σ confidence interval of 1.08-1.85. The lens SED fitting is performed with the Code Investigating GALaxy Emission (CIGALE 44,45 ) assuming z = 1.43.

The multiple rest-frame ultraviolet to rest-frame optical detections of SPT0311-58 E allow us to constrain the stellar mass using reasonable assumptions about the star-formation history at this early point in cosmic history. The SED is fitted by varying the e-folding time and age of a previously reported stellar population model 46 under single-and two-component formation histories, assuming solar metallicity and previously reported 47 initial mass function. The minimum radiation field, power-law slope and gamma, the fraction of dust mass exposed to radiation intensities above the minimum, from one dust model 48 , and the colour excess and attenuation slope from other dust models 49,50 are kept free in the SED fitting. The AGN contribution is set to zero because there are no photometric points to constrain the spectral range that is most affected by AGN power (mid-infrared) and thus any fraction between 0% and 60% of the dust luminosity is attributable to AGNs with nearly equal probability. However, this ignores the spatial distributions of the dust and line emission, which are not strongly peaked as is usually observed in AGN-dominated galaxies, so we deem this wide range to be unphysical. The inferred stellar mass and star formation rates are (3.5 ± 1.5) × 10 10 M  and (540 ± 175)M  yr -1 , respectively, for the two-component star-formation history. These values agree within the uncertainties for a single-component star-formation history. The infrared luminosity (L IR ; integrated over 8-1,000 µm) is (4.6 ± 1.2) × 10 12 L  and the extinction is

For the W source, we have only upper limits and the potentially contaminated IRAC detections to constrain the rest-frame optical and ultraviolet emission. Accordingly, we use the IRAC photometry as upper limits, along with the HST limits and far-infrared data in Extended Data Table 3, and model the SED with CIGALE. We find a luminosity of L IR = (33 ± 7) × 10 12 L  , seven times larger than for the E source. A consistent luminosity is obtained by fitting the far-infrared SED with a modified blackbody 51 . The inferred star-formation rate, which is closely connected to L IR , is (2,900 ± 1,800)M  yr -1 . As for the E source, the SED allows the AGN fraction to fall between 0% and 60% with roughly equal probability, so we take the absence of a dominant infrared emission region (see Fig. 1c and Extended Data Fig. 5) as an indication that the AGN contribution is unlikely to be important and fix the AGN fraction to zero. The dust luminosity due to star formation could therefore in principle be up to a factor of two smaller if the spatial distribution of the emission is ignored. Given that the photometry reaches to only the rest-frame V band, it is possible to hide a very large stellar mass behind dust obscuration for plausible values of the visual extinction (A V ≤ 6, as seen in other massive dusty galaxies 14,16,52,53 ). Considering the IRAC flux densities alone, we can calculate rest-frame mass-to-light ratios for the observed bands to see what masses could exist without relying on the poorly constrained CIGALE SED modelling. We use a stellar population synthesis code 54,55 to compute a stellar mass-to-light ratio under a range of assumptions: stellar ages of 0.1-0.8 Gyr (from a reasonably 'young' population to the approximate age of the Universe at the time) and metallicity of 0.1-1 times that of the Sun, with no dust attenuation. The mass then ranges from (2-10) × 10 10 M  per µJy of measured flux density. Taking the measured and de-magnified flux density (averaged between the two wavelengths) of 0.5 µJy, we find a stellar mass of (1-5) × 10 10 M  before correcting for extinction. If the extinction is as large as 5 mag, the true stellar mass could be unphysically large (>10 12 M  ), demonstrating that we have no useful constraint without greater certainty about the reliability of the IRAC flux densities or more photometric data points. Galaxy and halo masses. In Figs 2 and 3 we compile mass measurements for high-redshift galaxies discovered by various techniques. The galaxy sample comprises primarily galaxies identified through their luminous dust emission (DSFGs) and optically identified quasars (QSOs), which are typically the objects with the largest gas, dust or stellar masses at these redshifts. At the very highest redshifts, where very few galaxies have been found, objects selected on the basis of their ultraviolet emission are also included. The subsets of galaxies included in each LETTER RESEARCH figure overlap considerably, but are not identical because not all of the requisite information is available for each source. Dust mass. Mass estimates are unmodified from literature values 4,15,[56][57][58][59][60][61][62] , owing to the heterogeneity of the data available across the sample. The dust masses are generally derived from the far-infrared continuum emission, using one to several wavelengths. Differences between the cosmology assumed here and previously result in unimportant corrections and are ignored. Gas mass. Following standard observational practice, the primary source for the gas masses 4,27,[56][57][58]61,[63][64][65][66][67] shown in Fig. 2b is measurement of the luminosity of rotational transitions of CO. The lowest available rotational transition is typically used; any translation between the observed transition and the J = 1-0 line, which is most commonly used as a molecular gas indicator, is taken from the original source. Rather than accepting the varying coefficients for the conversion of CO luminosity to gas mass, we re-calculate all masses using a common value of α CO = 1.0M  (K km s -1 pc 2 ) -1 , which is a typical value for actively star-forming galaxies 68,69 . For one source the gas mass is estimated through the star-formation surface density 59 . Halo mass. The halo masses of Fig. 3 are derived from the gas mass sample above. Each halo mass is represented using a range of values, starting with a conservative and hard lower limit found by dividing the measured gas mass by the universal baryon fraction 10 f b = 0.19. This lower limit ignores any baryonic mass that has been converted into stars or hot or cool atomic gas phases, which would increase the inferred halo mass. A more realistic, but still conservative, lower limit is represented by the top of the plotted symbols in Fig. 3. Here we assume that the ratio of baryonic mass to halo mass is M b /M halo = 0.05. This value is a factor of about four less than the universal baryon fraction but still higher than the typical stellar-tohalo mass ratio inferred for haloes of any mass and redshift via subhalo abundance matching 3 . Given that we do not expect high-mass galaxies such as SPT0311-58 to expel a large fraction of their molecular gas content 70 or to later accrete dark matter without also accreting gas in proportion to the universal baryon fraction, it is reasonable to expect that the baryon-to-halo mass ratio should be less than this inferred upper limit on the stellar-to-halo mass ratio across all masses and redshifts. HFLS3 and SPT0311-58 masses. For the two most distant DSFGs, HFLS3 4 and SPT0311-58, which have extensive far-infrared photometry and atomic and molecular line measurements, we also compute the gas mass using a joint continuum-line radiative transfer model 6,29 . The mass for SPT0311-58 has been computed previously 6 without spatially resolved (CO and [C i]) line emission. For Fig. 3, only the total gas mass of the two SPT0311-58 sources is important for estimating the halo mass. For Fig. 2, the dust mass is divided between the two sources according to the ratio of dust continuum emission in our resolved observations. The gas mass is similarly divided, although the velocity profile of the CO lines provides weak evidence that the molecular gas is concentrated in SPT0311-58 W, which would increase the gas mass for this source by 15%. Calculation of halo rareness. Figure 3 demonstrates the 'rareness' of SPT0311-58 by considering its position in the dark-matter halo mass-redshift plane compared with other extreme high-redshift objects (DSFGs, QSOs and an LBG) that are believed to be hosted by massive dark-matter haloes. To quantify the rareness of these extreme objects we use a previously reported method 28 , including a MATLAB script (https://bitbucket.org/itrharrison/hh13-cluster-rareness) that we modified slightly to extend the calculation to z = 10. This method enables us to compute (z, M halo ) contours ('exclusion curves') above which the Poisson probability of such an object being detected in the standard ΛCDM cosmology is less than α < 1; the existence of a single object above such an exclusion curve is sufficient to rule out ΛCDM at the 100(1α)% confidence level. In Fig. 3, we plot 1σ exclusion curves (α = 0.32). Of the three different statistical measures of rareness proposed 28 , we use the '>ν' measure, which quantifies the rareness according to the minimum height of the primordial density perturbation from which a halo of mass M halo and redshift z could have formed: ν(M halo , z) ∝ [D + (z)σ(M halo )] -1 , where D + (z) is the normalized linear growth function and σ 2 (M halo ) is the variance of the matter power spectrum smoothed on the co-moving spatial scale that corresponds to the mass M halo . This statistic is sensitive to changes in the ΛCDM initial conditions, such as primordial non-Gaussianity (which would lead to more high-mass dark-matter haloes at a given redshift than expected in the standard ΛCDM cosmology). For the purposes of this calculation, we assume a ΛCDM cosmology with parameters 10 Ω m = 0.309, Ω Λ = 0.691, h 0 = 0.677 and σ 8 = 0.816 and use a previously reported halo mass function 71 .

The >ν rareness statistic (and the corresponding exclusion curves) depends on the region of the M halo -z plane to which the survey is sensitive. We assume that the SPT sample of lensed DSFGs is complete for z > 1.5. At lower redshift, the probability of lensing is strongly suppressed 30,72 , which means that the galaxy (or galaxies) associated with a halo mass of more than about 10 15 M  (the M halo value of the exclusion curves for z = 1.5) would have to have a very high intrinsic (that is, unlensed) millimetre-wavelength flux density (more than about 20 mJy) to be included in the sample. Because of the effects of downsizing (that is, star formation is terminated at higher redshift in higher-mass galaxies than in lower-mass galaxies), it is unlikely that massive galaxies at z < 1.5 would have sufficiently high infrared luminosity to be detected by the SPT 73 . We furthermore assume that the survey is complete for M  > 10 11 M  . The assumption that the sample is complete to M halo > 10 11 M  is a conservative one because the galaxies hosted by such haloes (which would have M b  10 11 M  ) are unlikely to be sufficiently luminous to be detected without being very strongly lensed (µ > 10); erring on the side of overestimating the completeness yields a lower limit on the rareness. Substituting a minimum halo mass of, for example, 10 12 M  would make the value of the >ν rareness statistic less than that found for 10 11 M  ; that is, SPT0311-58 would be inferred to be even rarer.

The total area from which the SPT DSFG sample was selected is 2,500 deg 2 . However, the fact that most of the SPT DSFGs are strongly lensed implies that the effective survey area is potentially much less than 2,500 deg 2 because not only must a galaxy have a high intrinsic millimetre-wavelength flux density to be included in the sample but it also must be gravitationally lensed so that it exceeds the approximately 20-mJy threshold for inclusion in redshift follow-up observations. Properly accounting for the effects of lensing on the sample completeness would require defining an effective survey area as a function of halo mass and redshift: A eff (M halo , z) = 2,500 deg 2 × P(µ min | M halo , z), where P(µ | M halo , z) is the probability of a galaxy hosted by a halo of mass M halo at redshift z being lensed by a factor µ min , the minimum magnification necessary for a halo of mass M halo and redshift z to be detectable. However, given the large uncertainties in determining such a function, we opt for a simpler approach. Instead, in Fig. 3 we plot exclusion curves for the full sky (dotted line), for an area of 2,500 deg 2 (dashed line), which corresponds to the assumption that all haloes in the mass and redshift range specified above would be detected even if they were not lensed, and for an area of 25 deg 2 (solid line), which corresponds to the assumption that the survey area corresponds to only the approximately 1% of the SPT fields over which the magnification for sources at z > 1.5 will be at least 30,72