Microwave background temperature at a redshift of 6.34 from H2O absorption

Distortions of the observed cosmic microwave background provide a direct measurement of the microwave background temperature at redshifts from 0 to 1 (refs. 1,2). Some additional background temperature estimates exist at redshifts from 1.8 to 3.3 based on molecular and atomic line-excitation temperatures in quasar absorption-line systems, but are model dependent3. No deviations from the expected (1 + z) scaling behaviour of the microwave background temperature have been seen4, but the measurements have not extended deeply into the matter-dominated era of the Universe at redshifts z > 3.3. Here we report observations of submillimetre line absorption from the water molecule against the cosmic microwave background at z = 6.34 in a massive starburst galaxy, corresponding to a lookback time of 12.8 billion years (ref. 5). Radiative pumping of the upper level of the ground-state ortho-H2O(110–101) line due to starburst activity in the dusty galaxy HFLS3 results in a cooling to below the redshifted microwave background temperature, after the transition is initially excited by the microwave background. This implies a microwave background temperature of 16.4–30.2 K (1σ range) at z = 6.34, which is consistent with a background temperature increase with redshift as expected from the standard ΛCDM cosmology4.

Distortions of the observed cosmic microwave background provide a direct measurement of the microwave background temperature at redshifts from 0 to 1 (refs. 1,2 ). Some additional background temperature estimates exist at redshifts from 1.8 to 3.3 based on molecular and atomic line-excitation temperatures in quasar absorption-line systems, but are model dependent 3 . No deviations from the expected (1 + z) scaling behaviour of the microwave background temperature have been seen 4 , but the measurements have not extended deeply into the matter-dominated era of the Universe at redshifts z > 3. 3. Here we report observations of submillimetre line absorption from the water molecule against the cosmic microwave background at z = 6.34 in a massive starburst galaxy, corresponding to a lookback time of 12.8 billion years (ref. 5 ). Radiative pumping of the upper level of the ground-state ortho-H 2 O(1 10 -1 01 ) line due to starburst activity in the dusty galaxy HFLS3 results in a cooling to below the redshifted microwave background temperature, after the transition is initially excited by the microwave background. This implies a microwave background temperature of 16.4-30.2 K (1σ range) at z = 6.34, which is consistent with a background temperature increase with redshift as expected from the standard ΛCDM cosmology 4 .
We used the Northern Extended Millimeter Array (NOEMA) to obtain a sensitive scan across the 3-mm atmospheric window towards the z = 6.34 massive dusty starburst galaxy HFLS3 (also known as 1HERMES S350 J170647.8+584623; see Methods) 5 . These observations reveal a broad range of emission features dominated by the CO, H 2 O and H 2 O + molecules and atomic carbon, on top of thermal dust continuum emission that is rising with frequency consistent with a dust temperature of T dust = 63.3 −5.8 +5.4 K (refs. 5,6 ) ( Fig. 1). The spectrum also shows a deep absorption feature due to the ortho-H 2 O(1 10 -1 01 ) ground-state transition at rest-frame 538 μm (observed at 3.95 mm, or 75.9 GHz). This absorption is about two times stronger than the continuum emission from the starburst at the same frequency (Fig. 2). For this effect to occur, a substantial population of the ortho-H 2 O 1 10 level (which lies 26.7 K above the 1 01 ground state) has to be excited by cosmic microwave background (CMB) photons as a basis for pumping of this level by the starburst infrared radiation field (see Extended Data Fig. 1). The effect becomes observable towards HFLS3 because of the warm CMB at this redshift, which is predicted to be T CMB = 20.0 K at z = 6.34 based on the standard ΛCDM cosmology (where T CMB (z) = T CMB (z = 0)*(1 + z) (1−β) , T CMB (z = 0) = 2.72548 ± 0.00057 K (ref. 7 ) and the power-law index β = 0). The absorption of photons from the CMB radiation field appreciably populates the H 2 O 1 10 level. The intense infrared radiation field from the starburst then preferentially de-populates the 1 10 level through radiative pumping, resulting in a deficit in the 1 10 level compared with 1 01 relative to a thermal distribution. In combination, these two processes result in an excitation temperature T ex of the H 2 O(1 10 -1 01 ) line that is lower than T CMB , such that the line becomes observable in absorption against the CMB. As the effect depends on the strength of the CMB radiation field, it can be used to measure T CMB for galaxies that have well-measured dust spectral energy distributions and dust continuum sizes, as is the case for HFLS3.
To understand the effect, we have calculated a series of spherically symmetric RADEX 8 models over a wide range of H 2 O column densities, assuming purely radiative excitation (Figs. 2 and 3; see Methods for additional details). Exposing a cold, H 2 O-bearing region associated with HFLS3 to the black-body CMB radiation field at T CMB (z = 6.34), the models suggest that 77.2% of the molecules will be in the 1 01 ground state and 20.3% will be in the upper 1 10 state, and all H 2 O transitions have an excitation temperature T ex equal to T CMB . As a result of this zero temperature contrast, no H 2 O emission or absorption would be observable, despite the fact that the H 2 O rotational ladder is excited by the CMB radiation. However, this picture changes when the same region is also exposed to the infrared radiation field emitted by the starburst nucleus of HFLS3, as the latter does not follow a single-black-body radiation pattern. Indeed, the infrared spectral energy distribution of HFLS3 reaches its peak intensity at 73.3 −1.3 +1.6 μm and can be approximated by a grey body with a Rayleigh-Jeans slope of β IR = 1.94 −0.09 +0.07 . This is due to the presence of dust at multiple temperatures and an increasing dust optical depth towards shorter wavelengths 5,6 . In this case, the level populations of the 1 01 and 1 10 states will deviate from the single-temperature thermal equilibrium population and change to 68.0% and 14.6%, respectively, for the solution shown in Fig. 2, resulting in an excitation temperature T ex of only 17.4 K for this transition. Owing to the ΔJ = 1 selection rule for photon emission/absorption, only three ortho-H 2 O transitions contribute to the modification of populations in the 1 01 and 1 10 levels, namely, the 538-μm 1 10 -1 01 and 180-μm 2 12 -1 01 transitions affecting the former, and the 108-μm 2 21 -1 10 transition affecting the latter (see Fig. 2; the 2 21 -1 01 transition is forbidden).  Fig. 1 to show that the line absorbs into the CMB (b, blue shading added for emphasis). The black curve is a fit to the spectrum. The red dashed curve is the best-fit radiative transfer model.

Article
The over-proportional de-population of the 1 10 level occurs because the H 2 O(2 21 -1 10 ) transition at 108 μm dominates the modification of the level population 9, 10 . This transition lies near the peak of the dust spectral energy distribution, where the dust emission has a higher optical depth than at 538 μm, where the H 2 O(1 10 -1 01 ) transition occurs. The increase in dust optical depth with wavelength leads to an increased availability of 108-μm photons relative to 538-μm photons compared with the thermal equilibrium case of a single-black-body radiation field. This implies that the 2 21 -1 10 transition at 108 μm is exposed to a more intense infrared radiation field than the 1 10 -1 01 transition at 538 μm. For the infrared radiation field of HFLS3, the models therefore suggest that both H 2 O transitions should be found in absorption, but that the 2 21 -1 10 transition (which is not covered by our observations) will contribute a larger fraction to the 1 10 level de-population than that found for the thermal equilibrium case. The CMB photons, on the other hand, only provide the base population expected for the thermal black-body equilibrium case. This causes a reduced excitation temperature in the 1 10 -1 01 transition compared with the thermal equilibrium case, unless the 1 01 level is even more substantially de-populated due to the H 2 O(2 12 -1 01 ) transition at 180 μm. This, however, does not occur, as our models show that this transition is expected to be seen in emission (at a line strength approximately 3-5 times below a previously reported upper limit 5 ). This is because the upper-level population of the 2 12 -1 01 transition is also affected by the population of the 2 12 level through the 2 21 -2 12 line, which appears in emission due to the pumping of its upper 2 21 state by the 2 21 -1 10 absorption line. As such, the models suggest a net deficit in the upper-level population of the 1 10 -1 01 transition compared with the thermal equilibrium case, which causes the T ex of the line to end up below the CMB temperature. The RADEX models yield T ex = 17.4 K due to this level-population modification. To translate model-predicted temperature differences into an absorption-line flux that can be compared with the observations, the size of the emitting/absorbing region needs to be known. Based on NOEMA observations at rest-frame 122 μm (Extended Data Fig. 2), we estimate the dust continuum size of the emitting region at 108 μm (the wavelength of the pumping transition) of HFLS3 to be r 108μm = 1.62 ± 0.45 kiloparsecs (kpc). Within the uncertainties of the size estimate, the RADEX models suggest that the strength of the observed H 2 O absorption can be reproduced over about two orders of magnitude in H 2 O column density, with a lower limit of around 10 16 cm −2 . The minimum covering fraction of the dust continuum is about 60% when conservatively leaving T CMB as a free parameter (100% is assumed for the grid shown in Fig. 3a). The upper limit for the H 2 column density implied by the gas mass of HFLS3 of (1.04 ± 0.09) × 10 11 M sun (ref. 5 ) provides a lower limit to the gas phase [H 2 O]/[H 2 ] abundance of >2 × 10 −7 , which falls within the range of 10 −9 -10 −5 found for nearby starbursts 11 . The small difference ΔT = T ex − T CMB = −2.6 K is, therefore, sufficient to explain the observed strength of the H 2 O(1 10 -1 01 ) absorption line towards the CMB in HFLS3 when assuming a layer of cold, diffuse H 2 O-bearing gas with a high covering fraction in front of the warm dust continuum source associated with the H 2 O emission lines.
As the absorption line is observed in contrast to the CMB, we can use the strength of the absorption line to obtain a measurement of T CMB at the redshift of HFLS3. The RADEX models suggest that, to detect the H 2 O(1 10 -1 01 ) line in absorption against the CMB, T CMB (z = 6.34) must be >7-8 K, independent of the model assumptions (see Methods). The observed strength of the signal suggests 16.4 K < T CMB (z = 6.34) < 30.2 K (1σ, or 12.8 K < T CMB (z = 6.34) < 34.0 K 2σ) for HFLS3 when treating T dust , β IR and the wavelength where the dust optical depth reaches unity as free-fitting parameters for each dust continuum size sampled by the models. This explains why the effect has not been previously seen. T CMB must be sufficiently high to satisfy the requirement of a notable H 2 O 1 10 level population due to the CMB, such that a de-population by the infrared radiation field of the starburst will lead to a sufficiently important decrement to be observable in absorption against the CMB. This limits observability to z > 4.5 for dust spectral energy distribution shapes and dust continuum sizes of star-forming galaxies like HFLS3 (Fig. 3c), where only a few spectra at rest-frame 538 μm with sufficient signal-to-noise ratio to detect the effect exist. This differs from molecules like H 2 CO, for which absorption against the CMB has been predicted to occur at The overlapping region between the white boundary (that is, the minimum allowed absorption strength) and the size measurement (that is, the minimum required emitting area at 100% covering fraction) is the allowed parameter space for the absorption strength within 1σ r.m.s. The minimum required radius at N(H 2 O) ~ 10 17 cm −2 is due to a minimum in T ex in the models. b, Constraints on T CMB for the observed absorption strength (green line and shaded region) at the minimum size compatible with the observations (a), based on the same models (red/blue shaded regions are the allowed ranges within the source radius +1σ/+2σ r.m.s.). The source radius at face value (black line), as well as the −1σ and −2σ r.m.s. regions (not shown), are ruled out by the observations. The minimum filling factor of the dust emission region CF min is indicated for the +1σ and +2σ r.m.s. regions. The grey dashed line shows a model assuming a continuum radius of 5 kpc, which provides a conservative lower limit on T CMB . c, Observability of the H 2 O absorption as a function of redshift for three solutions allowed by the data without and with collisional excitation. The effect becomes observable at z ~ 4.5 and remains visible at similar strength to z > 12. The lower-redshift limit is higher in cases where collisional excitation is important, but the impact is minor below n(H 2 ) = 10 5 cm −3 .
any redshift up to the present day 12 , but for which no detections at high redshift currently exist 13 . For starbursts with dust as warm as HFLS3, its relative strength is expected to continue to increase with redshift all the way up to z ~ 7-8, and to remain observable back to the earliest epochs when such galaxies existed. The thermal Sunyaev-Zel'dovich (SZ) effect entails the modification of the thermal distribution of CMB photons by Thomson scattering off electrons at high temperatures in the intergalactic medium of galaxy (proto)clusters. The effect observed here requires the CMB photons to excite the H 2 O rotational ladder to create a thermal distribution of the lower energy levels, which is then modified through the absorption of far-infrared photons from the starburst radiation field permeating the interstellar medium. In both cases, the CMB is responsible for establishing unperturbed thermal distributions (of photons and H 2 O excitation, respectively), which are then modified by local conditions. The SZ effect is a broad-band modification of the thermal distribution of CMB photons via scattering, with an expected signal strength (relative to the CMB) that is independent of redshift. In contrast, the H 2 O absorption signal described herein is a narrow-band (spectral-line) absorption process of the CMB photons, catalysed (in part) by the CMB itself, with an absorption-line strength that increases with redshift owing to the increasing temperature of the CMB, relative to the fixed excitation temperature of the H 2 O(1 10 -1 01 ) transition.
Standard ΛCDM cosmologies predict a linear increase of T CMB with (1 + z). However, there are hypothetical physical mechanisms that could lead to departures from this linear behaviour, including the evolution of physical constants 14 , decaying dark-energy models 15 and axion-photon-like coupling processes 16,17 . Direct measurements of T CMB are, thus, a crucial test of cosmology, but they are currently limited to z < 1, due to the lack of sufficiently precise measurements of the thermal SZ effect in galaxy clusters at higher redshifts ( Fig. 4; see Methods for further details). A limited sample of additional constraints exists at z = 1.8-3.3 based on measurements of T ex for the ultraviolet transitions of CO, [CI] and [CII] in absorption-line systems along the lines of sight to quasars. These lines are not directly observed in contrast to the CMB, and they use the T ex of these lines as a proxy of T CMB , such that the resulting T CMB estimates are subject to model-dependent excitation corrections 3,17-28 . As an example, for the CO molecule, T ex typically already exceeds T CMB in the diffuse interstellar medium in the Milky Way due to collisional excitation, showing a rising excess with increasing CO optical depth due to photon trapping 18 . In contrast, our models suggest that collisional excitation of H 2 O becomes important only at very high densities, such that H 2 O-based measurements are probably only minimally affected by this effect. The H 2 O absorption against the CMB at z = 6.34 reported here thus provides the most direct constraint on T CMB currently available at z > 1. Indeed, the existence of this effect on its own directly implies that the CMB is warmer than at low redshift, because T CMB must be sufficiently high to notably excite the H 2 O 1 10 level, which lies 26.7 K above ground, as a basis for the observed decrement due to de-population of this level by the starburst radiation field.

Article
A combined fit to the available data (Fig. 4) is consistent with the redshift scaling expected from ΛCDM. Fitting for the adiabatic index γ in the equation of state between pressure P and energy density ρ for the sum of baryonic and dark matter and radiation-that is, P rm = (γ − 1)ρ rm with a standard formalism (see Methods)-we find γ = 1.328 −0.007 +0.008 , which agrees with the standard value of γ = 4/3 expected in ΛCDM. At the same time, we find an effective dark energy equation of state parameter w eff = P de /ρ de = −1.011 −0.017 +0.018 , which is consistent with the w = −1 expectation for a dark energy density that does not evolve with time.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-04294-5.

NOEMA observations
The target was observed in the 3-mm wavelength band 1 (rest-frame 400 μm) with NOEMA as part of project S20DA (Principal Investigators: D. A. Riechers, F. Walter). Three partially overlapping spectral setups were observed under good weather conditions between 26 July 2020 and 25 August 2020 with ten antennas in the most compact D configuration, using a bandwidth of 7.7 GHz (dual polarization) at 2-MHz spectral resolution per sideband. We also included previously  Fig. 2. We also make use of previously published 5 rest-frame 158-μm NOEMA data, which were adopted without further modification.

Line and continuum parameters
The flux of the H 2 O(1 10 -1 01 ) line was extracted by simultaneous Gaussian fitting of the line and continuum emission (including a linear term for the continuum) in the one-dimensional spectrum shown in Fig. 1, which was extracted from the image cube. The source is unresolved at the frequency of the H 2 O(1 10 -1 01 ) line, such that the main uncertainties are due to the slope of the continuum emission and the appropriate fitting of other nearby lines, in particular, CO(5-4). The uncertainties in these parameters are part of the quoted uncertainties. We find a line peak flux of −818 ± 145 μJy at a line full width half maximum (FWHM) of 507 ± 111 km s −1 , centred at a frequency of 75.8948 GHz (±46 km s −1 ; the calibration uncertainties on the line FWHM and centre frequency are negligible and that on the line peak flux is <10%-that is, minor compared with the measurement uncertainty). Given the rest frequency of the line of 556.9359877 GHz, this corresponds to a redshift of 6.3383, which is consistent with the systemic redshift of HFLS3 (z = 6.3335 and 6.3427 with uncertainties of ±14 and ±54 km s −1 at Gaussian FWHM of 243 ± 39 and 760 ± 152 km s −1 , respectively, for the two velocity components detected in the 158-μm [CII] line) 5 . For comparison, the H 2 O(2 02 -1 11 ) and H 2 O(2 11 -2 02 ) emission lines in HFLS3 have FWHM of 805 ± 129 and 927 ± 330 km s −1 , respectively 5 -that is, only marginally broader than the 1 10 -1 01 line at the current measurement uncertainties. The continuum flux at the line frequency is 396 ± 15 μJy, corresponding to 48% ± 9% of the absorption-line flux (the relative flux calibration uncertainty between the line and continuum emission is negligible). We also measured the 335.5-GHz continuum flux by two-dimensional fitting to the continuum emission in the visibility plane. We find a flux of 33.9 ± 1.1 mJy, which agrees with previous lower-resolution observations at the same wavelength 5 . The major (minor) axis FWHM diameter of the source is 0.617 ± 0.074 arcsec (0.37 ± 0.20 arcsec). This yields the physical source size quoted in the main text at the redshift of HFLS3.

Brightness temperature contrast
The H 2 O(1 10 -1 01 ) line leads to a decrement in continuum photons from the starburst and, as such, is observed as a lack of continuum emission at its frequency at the position of the starburst. It therefore appears as negative flux in an image where starburst continuum emission has been subtracted. In addition, (sub)millimetre-wavelength interferometric images reveal structure against a flat sky background defined by the large-scale CMB surface brightness, which the interferometer does not detect itself due to its limited spatial sampling. Therefore the fraction of the signal due to the decrement in CMB photons at the position of the starburst not only appears as negative flux without subtracting any further signal but it also corresponds to a lack of continuum emission at the line frequency in practice. As the mere presence of an absorption-line signal stronger than the measured continuum emission implies absorption against the CMB, this interpretation is not limited by uncertainties in the galaxy continuum flux or uncertainties in the absolute flux calibration.

Line-excitation modelling
RADEX is a radiative transfer program to analyse interstellar line spectra by calculating the intensities of atomic and molecular lines, assuming statistical equilibrium and considering collisional and radiative processes, as well as radiation from background sources. Optical depth effects are treated with an escape probability method 8 . Studies of nearby star-forming galaxies show that the observed absorption strengths of the ground-state H 2 O and H 2 O + transitions are due to cooler gas that is located in front of, and irradiated by, a warmer background source that is emitting the infrared continuum light that also excites the higher-level H 2 O emission lines 11,29 . We therefore adopt the same geometry for the modelling in this work, which is adequately treated within RADEX (that is, treating the dust continuum plus the CMB as background fields for the absorbing material) 8 . The dust continuum emission is modelled as a grey body with treating T dust , β IR and the wavelength where the dust optical depth reaches unity as free-fitting parameters for each dust continuum size and T CMB sampled by the models. The observed spectral energy distribution of HFLS3, including all literature 5 photometry and the measurements presented in this work, is then treated as the contrast between the dust continuum and CMB background fields, such that the resulting fit parameters for the dust continuum source change with T CMB in a self-consistent manner. In the RADEX models, we derive the H 2 O peak absorption depth into the CMB. We then multiply the best matching peak absorption depth found by RADEX with a Gaussian matched to the fitted line centroid and line width obtained from the observed line profile in Fig. 2 to determine the model line profile. In this approach, the shallower absorption in the line wings either corresponds to a lower filling factor of the H 2 O layer at the corresponding velocities or to lower H 2 O column densities. Although collisions of H 2 O molecules with H 2 is another mechanism that can modify the level populations especially at very high gas densities (which is an important mechanism for the cooling of low-excitation-temperature transitions of molecules like H 2 CO to below T CMB ) 12,30 , the RADEX models show that they do not affect our findings (see Fig. 3c). We therefore adopt models with essentially no collisions by assuming a very low gas density of n(H 2 ) = 10 cm −3 . We then compare our findings to those obtained when adopting conditions that are similar to those found in local starburst galaxies 11 and to those found for high-density environments with n(H 2 ) > 10 5 cm −3 . The cross sections for collisions out of the 1 01 level are always larger than those out of the 1 10 level, independent of the collision partner and the temperature at which the collisions take place [31][32][33] . Therefore collisions cannot be responsible for an over-proportional de-population of the 1 10 level relative to the 1 01 ground state, and the net effect of including collisions is a decrease in the absorption depth into the CMB by reducing the T CMB − T ex temperature difference at very high gas densities compared with cases without collisions. For reference, the effect of collisions on the determination of T CMB is negligible for the typical conditions found in local starbursts (that is, n(H 2 ) ~ 10 4 cm −3 ; T kin = 20-180 K) 11 and only starts to have an impact for very high densities n(H 2 ) > 10 5 cm −3 . For a given continuum source size, the constraints on T CMB would therefore be tighter (that is, would more quickly become inconsistent with the observations) for the high-density case than for the case without collisions, such that the latter approach is more conservative (see Fig. 3b).
The overall impact of collisional excitation would therefore be more stringent requirements on the source size, covering fraction and water column, such that their inclusion would only further strengthen our conclusions. We note that this is the opposite effect to the case of the studies of ultraviolet lines 3,17,[19][20][21][22][23][24][25][26][27][28] , where neglecting collisional excitation results in less conservative constraints on T CMB . If we were to assume that the H 2 O absorption were to emerge from within the infrared continuum-emitting region, a larger source size would probably be required to obtain the same absorption-line strength due to a reduced effective radiation field strength from the starburst. Previous modelling attempts of nearby galaxies assuming such geometries have not been able to produce H 2 O(1 10 -1 01 ) line absorption on the scales necessary to explain the observations of HFLS3, which may indicate that even more complex assumptions would be required 11 . Thus, the resulting constraints would, once again, be less conservative, perhaps acting in a similar manner as the high-density case. Excluding both of these effects from the models leads to a maximally conservative estimate of T CMB and its uncertainties. Assuming a plane-parallel or similar geometry instead of a spherical geometry would only have a minor impact on our findings 8 . The models shown in Fig. 3 assume a filling factor of unity, which is the most conservative possible assumption. A more clumpy geometry with a lower covering fraction remains possible for all T CMB values for which the predicted absorption strength exceeds the observed value (see shaded regions in Fig. 3b). For reference, the minimum covering fractions consistent with the continuum size at the observed signal strength are shown for the different cases considered in Fig. 3. The line absorption is also found to be optically thick, with an optical depth of τ H2O = 21.1 for the solution shown in Fig. 2b. To determine the redshift above which the effect becomes observable (Fig. 3c), we fixed r 108μm , T dust , β IR and M dust to the observed values and the H 2 O column density to the value corresponding to the model spectrum. H 2 O line absorption into the dust continuum of HFLS3 would already become visible at z > 2.9, but absorption into the CMB only becomes observable at z > 4.5 (or higher for H 2 densities of >10 5 cm −3 ). These values account for changes in the shape of the dust grey-body spectrum (that is, changes in the relative availability of 538-μm and 108-μm photons) due to changes in T CMB with redshift. To better quantify the impact of different modelling parameters, we have varied T dust and β IR beyond their previously estimated uncertainties (nominal reference values without considering variations in T CMB from the literature are T dust = 63.3 −5.8 +5.4 K and β IR = 1.94 −0.09 +0.07 ) 5,6 . This is necessary because both parameters are dependent on the varying T CMB in our models (and therefore are changing parameters in Fig. 3b, c), such that their true uncertainties need to be re-evaluated. We independently varied β IR in the 1.6-2.4 range and T dust in the ±20 K range as functions of T CMB around the best-fit values. This shows that β IR > 2.0 and T dust lower by more than 10 K from the best fits yield very poor fits to the spectral energy distribution data, whereas Δβ IR > −0.1 below the best-fit value would require a larger continuum size than the measured r 108μm + 1σ and therefore are disfavoured by the size constraint. Excluding these ranges, the extrema across this entire range would extend the uncertainty range in the predicted T CMB by only −1.7 and +5.4 K and −0.8 and +4.4 K for the r 108μm + 1σ and r 108μm + 2σ cases, respectively. For comparison, the difference between the +1σ and +2σ uncertainty ranges is −3.6 and +3.8 K). This shows that the impact of the uncertainties in the dust spectral energy distribution fitting parameters on those in T CMB are subdominant to those in the continuum size measurement. Conversely, we have studied the impact of changes in T CMB on the best-fit T dust and β IR . For the values corresponding to r 108μm + 1σ and r 108μm + 2σ ranges, T dust typically changes by <0.5 K and β IR typically changes by <0.1-0.2 when varying the parameters independently. These changes are larger than the actual uncertainties, because the fit to the dust spectral energy distribution becomes increasingly poorer with these single-parameter variations. At the same time, these changes are subdominant to those induced by changes in dust continuum size within the +1σ and +2σ uncertainty ranges, which is consistent with our other findings.

Other H 2 O transitions in HFLS3
Five H 2 O lines were previously detected towards HFLS3 (2 02 -1 11 , 2 11 -2 02 , 3 12 -2 21 , 3 12 -3 03 and 3 21 -3 12 ) and two additional lines were tentatively detected (4 13 -4 04 and 4 22 -4 13 ) 5 . The J up = 3 transitions are due to ortho-H 2 O and all other transitions are due to para-H 2 O. All of these transitions appear in emission. Given the high critical densities of these transitions, our RADEX models cannot reproduce the strength of these lines as the same time as the observed ortho-H 2 O(1 10 -1 01 ) absorption strength, which suggests that they emerge from different gas components. For reference, to reproduce the strength of the H 2 O(2 11 -2 02 ) in Fig. 1 alone with collisional excitation, n(H 2 ) = 2 × 10 7 cm −3 and T kin = 200 K would be required, but the H 2 O(1 10 -1 01 ) would no longer appear in absorption against the CMB if it were to emerge from the same gas component. This is consistent with the picture that the H 2 O absorption is due to a cold gas component along the line of sight to the warm gas that gives rise to the emission lines 11 . Observations of the para-H 2 O(1 11 -0 00 ) ground state do not currently exist for HFLS3, but our models do not show this line in absorption towards the CMB.

Origin of the lower and upper limits on T CMB
Our models show that the lower limit on T CMB at a given redshift based on the observed H 2 O absorption is due to the minimum 'seed' level population due to the CMB black-body radiation field. To determine a conservative lower limit, we have calculated models with continuum sizes up to r 108μm = 5 kpc (see Fig. 3b), corresponding to a +7.5σ deviation from the observed continuum size, and recorded the temperatures at which such weakly constrained models turn into absorption. We find that this results in a lower limit of T CMB > 7-8 K, independent of the model assumptions. This finding alone does not explain the existence of an upper limit in Fig. 3b. For a given size of the dust continuum emission, an increase in T CMB also requires an increase in M dust to still reproduce the observed dust spectral energy distribution, which leads to an effective increase in the dust optical depth at a given wavelength. The result of a rising optical depth is that the grey-body spectrum between 538 and 108 μm increasingly resembles a black-body spectrum and, hence, a decrease in the H 2 O absorption against the CMB. This effect is responsible for the upper limit in allowed T CMB for a given dust continuum size and absorption strength.

Uncertainties of T CMB measurements
The uncertainties shown for the literature data in Fig. 4 are adopted from the literature sources without modification, and they typically represent the statistical uncertainties from the individual measurements or sample averages. Individual cluster measurements of the thermal SZ effect may be affected by dust associated with foreground galaxies or the Milky Way, the galaxy clusters or background galaxies that may be amplified by gravitational lensing, uncertainties in the reconstruction of the Compton-y parameter maps due to flux uncertainties, radio emission due to active galactic nuclei and/or relics, the kinetic and relativistic SZ effects, and general bandpass and calibration uncertainties 17 . Furthermore, uncertainties on the cluster geometry-and therefore line-of-sight travel distance of the CMB photons through the cluster-and on the temperature of the intra-cluster gas limit the precision of individual SZ measurements. Sample averages may also be affected by systematics in the stacking procedures. Individual data points deviate by up to at least two standard deviations from the trend, which may indicate residual uncertainties beyond the statistical error bars provided, such that the error bars shown in Fig. 4 are underestimated. The main source of uncertainty for the ultraviolet absorption-line-based measurements are due to the assumption of no collisional excitation, which is not taken into account in the statistical uncertainties shown in Fig. 4. Attempts to take this effect into account appear to suggest substantially larger uncertainties than indicated by individual error bars 27 (Fig. 4). To expand on earlier estimates 21 , we have calculated RADEX models for typical T kin , n(H) and column densities found from [CI] measurements in the diffuse interstellar medium 34 , which suggests that collisional excitation contributes to the predicted T ex of the lower fine-structure transition. Although we show the original unmodified data, the ultraviolet-based measurements are therefore subject to uncertainties due to model-dependent excitation corrections in addition to the statistical uncertainties. Furthermore, the fine-structure levels of tracers like the [CI] lines can be excited by ultraviolet excitation and following cascades. To constrain T CMB based on these measurements, the kinetic temperature, particle density and local ultraviolet radiation field must be known, and are typically determined based on tracers other than the species used to constrain T CMB . Also, some measurements are based on spectrally unresolved lines, which limits the precision of kinetic temperature measurements based on thermal broadening 21 . Owing to these uncertainties, the ultraviolet absorption-line-based measurements are probably consistent with the standard ΛCDM value, but they do not constitute a direct measurement of T CMB without notable further assumptions. For reference, the median T CMB /(1 + z) estimate based on the [CI] measurements alone (excluding upper limits) is 3.07 K, with a median absolute deviation of 0.09 K and a standard deviation of 0.31 K. Therefore the current sample median deviates from the ΛCDM value by about one standard deviation. A combination of the (uncorrected) measurements based on CO, [CI] and [CII] provides a median value of 2.84 K, with a median absolute deviation of 0.15 K and a standard deviation of 0.25 K. This highlights the importance of the corrections discussed above and in the literature and the value of measurements with systematic uncertainties that differ from this method to obtain a more complete picture. The main source of uncertainty of the H 2 O-based measurements, beyond the caveats stated in the line-excitation-modelling section, are the statistical uncertainties on the source size, the lack of a direct measurement of the absorbing H 2 O column density, variations in the dust mass absorption coefficient and the filling factor. Given the high metallicity suggested by other molecular line detections, the limitation to high filling factors due to the source size and the constraint on the gas mass from dynamical mass measurements, the main source of uncertainty resides in the source size due to limited spatial resolution in the current data. As such, major improvements should be possible by obtaining higher, (sub-)kpc resolution (that is, <0.2") imaging with the Atacama Large Millimeter/submillimeter Array (ALMA; for other targets) and planned upgrades to NOEMA, and, in the future, with the Next Generation Very Large Array (ngVLA). Statistical uncertainties will also be greatly reduced by observing larger samples of massive star-forming galaxies over the entire redshift range where measurements are possible, closing the gap to SZ-based studies, which are currently limited to z < 1. The resulting improvement in precision will provide the constraints that are necessary to confirm or challenge the evolution of the CMB temperature with redshift predicted by standard cosmological models.

Accessibility of the line signal
The frequency range currently covered by NOEMA is 70.4-119.3, 127.0-182.9 and 196.1-276.0 GHz (with greatly reduced sensitivity above about 115 and 180 GHz in the first two frequency ranges). ALMA covers the 84-500-GHz range with gaps at 116-125 and 373-385 GHz, with a future extension down to 65 GHz (with greatly reduced sensitivity below about 67 GHz). The ngVLA is envisioned to cover the 70-116-GHz range. Excluding regions of poor atmospheric transparency, the H 2 O(1 10 -1 01 ) line is therefore observable in these frequency ranges at redshifts of z = 0.1-0.4, 0.5-2.0, 2.1-3.4 and 3.8-6.9 in principle, but the detectability of the line in absorption against the CMB is probably limited to the z ~ 4.5-6.9 range if the spectral energy distribution shape of HFLS3 is representative. At lower frequencies, the Karl G. Jansky Very Large Array and, in the future, ALMA and the ngVLA also provide access to the <52-GHz range, such that the signal also becomes observable at z > 9.7 in principle. In conclusion, the absorption of the ground-state H 2 O transition against the CMB identified here could be traced from the ground towards star-forming galaxies across most of the first approximately 1.5 billion years of cosmic history.

Detectability of the line signal for different spectral energy distribution shapes
To investigate whether the effect is expected to be detectable for different galaxy populations, we have applied our modelling to the z = 3.9 quasar APM 08279+5255, for which the dust spectral energy distribution is composed of a dominant 220-K dust component and a weaker 65-K dust component, contributing only 10-15% to the far-infrared luminosity [35][36][37][38][39][40][41][42][43][44][45][46] . The models suggest that the line is expected to occur in emission and that it would not be expected to be detectable in absorption at any redshift out to at least z = 12 in galaxies with similar dust spectral energy distributions. Other far-infrared-luminous, high-redshift, active galactic nucleus host galaxies typically show a stronger relative contribution of their lower-temperature dust components, such that the effect may remain detectable in less extreme cases. For galaxies with lower dust temperatures than HFLS3, the effect may be present even at lower redshifts, but is typically expected to be weaker in general and to disappear at redshifts where T CMB approaches their T dust . For a dust spectral energy distribution shape resembling the central region of the Milky Way but otherwise similar properties, the effect is expected to be reduced by more than two orders of magnitude at its redshift peak, and to become virtually unobservable at the redshift of HFLS3. Thus, dusty starburst galaxies appear to be some of the best environments to detect the effect.

Derivation of equation of state parameters
To determine the adiabatic index, we assume a standard Friedmann-Lemaitre-Robertson-Walker cosmology with zero curvature and a matter-radiation fluid that follows the standard adiabatic equation of state quoted in the main text. This would correspond to a redshift scaling T CMB (z) = T CMB (z = 0)*(1 + z) 3(γ − 1) in the presence of a dark energy density that does not scale with redshift. The dark energy density is parameterized to scale with a power law (1 + z) m , where m = 0 corresponds to a cosmological constant. With standard assumptions, this yields a redshift scaling of T CMB (ref. 15 ):