© ESO 2020

1. Introduction

Any structure on the ground that is exposed to the night sky experiences a radiation imbalance and therefore typically cools below the ambient air temperature (see Zhao et al. 2019 for an overview). Sky subcooling is strongest at sites with low water vapor column, hence under conditions commonly present at telescope sites (high and dry). As a consequence, the air surrounding a steel truss for example is cooled, causing a wake of cool air. Light rays passing the cooled air experience an optical path difference which induces a transient wavefront error in a telescope (Wilson 2013). This situation is most severe at small wind speeds and thus became known as the “low-wind effect”, and has recently been studied (Sauvage et al. 2015; Milli et al. 2018). A secondary detrimental consequence of subcooling is structural deflection due to thermal contraction, inducing mechanical stress and optical misalignment. Quantifying the effective sky temperature is therefore required when designing advanced ground-based telescopes and analyzing their image quality.

Previous studies (Berdahl & Fromberg 1982; Gliah et al. 2011; Sun et al. 2017; Li et al. 2017, 2018) used wide-angle integrating detectors subtending nearly the whole sky. Moreover, these detectors are tested and calibrated at sites with lower altitudes and higher humidity. To derive the effects of the sky temperature on a telescope structure and the optical elements such as the main mirror, a detailed understanding of sky irradiance as a function of the direction on the sky is required as well. To our knowledge, such a study has not yet been published.

In Sect. 2, we discuss the night-sky emissivity and define the effective bolometric sky temperature T_sky. Section 3 then sets up a thermal balance model of a structure on the ground exposed to sky radiation and forced convection. Furthermore, we discuss the effects of heat conduction and thermal inertia. Finally, we conclude in Sect. 4.

2. Sky radiance and sky temperature

The spectrum of the radiance L_sky for a given area on the sky irradiating a surface strongly depends on wavelength λ, zenith distance θ, and weather conditions. Its total bolometric radiance E_BOL, defined as the integral of the radiance spectrum over all wavelengths, is dominated by those fractions of the wavelength domain where the atmosphere is completely opaque. In these regions, the mean free path of the light is often a few hundred meters only; in optically very thick lines it may even be only a few meters. Thus the main emission is strongly bound to the local temperature near the ground. Those regions where the atmosphere is nearly transparent do not significantly contribute to the integrated radiance. We use here the common definition of a sky temperature T_sky as the equivalent temperature of a black-body radiation curve of the same bolometric radiance E_BOL (see Li et al. 2017, 2018, and references therein),

$$\begin{array}{c}\hfill E_{\mathrm{BOL}}={\displaystyle {\int }_0^{\infty }{L}_{\mathrm{sky}}(\lambda )\mathrm{d}\lambda :=\frac{\sigma }{\pi }{T}_{\mathrm{sky}}^4,}\end{array}$$(1)

where the unit is W m⁻² sr⁻¹, with the Stefan–Boltzmann constant σ  =  5.670367  ×  10⁻⁸ W m⁻² K⁻⁴.

2.1. General trend and climatology of the sky temperature

With our European Southern Observatory (ESO) Sky Model¹ (Noll et al. 2012; Jones et al. 2013) we are able to calculate detailed spectral dependencies and integrate them to derive the bolometric flux and the sky temperature T_sky. Although this sophisticated model contains a wide variety of components of illumination such as scattered moonlight, zodiacal light, scattered starlight, and atomic and molecular emission lines in the upper atmosphere (airglow), the vast majority of the bolometric flux originates from about 4 to 7 μm and longwards of 13.5 μm and is dominated by the emission of the lower atmosphere. Our model uses a spectral model of the Earth’s lower atmosphere in local thermal equilibrium calculated by means of the Line-By-Line Radiative Transfer Model (LBLRTM; Clough et al. 2005). This code is used with the spectral line parameter database aer, which is based on HITRAN (Rothman et al. 2009). Averaged atmospheric profiles for Cerro Paranal describing the chemical composition as a function of height (a combination of the ESO ambient Meteo Monitor, a standard atmosphere for low latitudes², and the 3D Global Data Assimilation System (GDAS) model) are also used as input. The GDAS model is provided by the National Oceanic and Atmospheric Administration (NOAA)³ and contains time-dependent profiles of temperature, pressure, and humidity. For the sky model, bimonthly averages over several years were built yielding reliable forecasts for typical conditions at Cerro Paranal. For the transmission and thermal emission of the atmosphere, the current model can be adapted to any site if appropriate GDAS profiles are used. Clear cloudless weather conditions are required with the LBLRTM standard setup.

We show the dependency of the thermal sky radiance as a function of precipitable water vapor (PWV) in Fig. 1. This quantity is equal to the integrated column density across the whole atmosphere, including all layers contributing to the opaqueness, and is given in units of equivalent thickness in millimeters of liquid water. The opaque fraction of the spectral energy distribution resembles a black body radiation curve. While the transparency in the window around λ = 10 μm is only marginally affected by the H₂O molecular bands, it varies strongly in the region from 15 to 28 μm (astronomical Q-band) with the water vapor content. Therefore, the weather-dependent variations and uncertainties mainly originate from this wavelength band. There are marginal transparency windows beyond these wavelengths (astronomical mid-infrared Z-band with 28 <  λ <  40 μm). As shown by Yoshii et al. (2010), this assumption still holds at sites of even higher altitude than Mauna Kea (e.g., as planned for the 6.5 m Tokyo Atacama Telescope project at 5640 m). Therefore, an approximation with a black body curve is normally sufficient. The ozone absorption in the center of the N-band does not follow the trend of the other opaque regions because of the high altitude of the absorbing molecules.

Fig. 1. Sky radiance Lsky(λ) for Cerro Paranal at different PWV conditions typical for observations in the direction of zenith. The astronomical observation bands and the main ozone absorption are marked. For further details on species that are not important for us here (e.g., CO2, N2O and CH4) see Fig. 1 and Table 1 in Smette et al. (2015).

The sky radiance (expressed as sky temperature as defined above) therefore depends significantly on the PWV. This circumstance was already found by Berdahl & Fromberg (1982), but these latter authors characterized the measured radiation as a function of the relative humidity near the ground, specifically the dew point temperature T_DP calculated from the humidity. Although radiation from low altitudes prevails, the variation is dominated by semi-transparent regions in the Q-band. Here we see radiation from higher layers having different temperatures and a water vapor content which is not necessarily linked to the relative humidity measured near the ground.

We find a phenomenological fit function of the sky temperature with the water vapor content for dry conditions (Fig. 2) which is precisely passing the points with deviations of less than 0.1 K. In more humid conditions (PWV > 15 mm) the relation still fits rather well, but has some larger uncertainties, on the order of 2 K. Although the fit is based on model data for the southern winter months June and July, it is also valid for the other months of the year. The small differences in the bimonthly mean profiles of pressure and temperature, as well as the relative vertical distribution of water vapor at Cerro Paranal, appear to have a negligible effect on the integrated radiance spectrum.

Fig. 2. Correlation of the sky temperature Tsky and PWV at zenith Only the range of conditions with PWV <  15 mm relevant for typical astronomical observations is used for the final fit function.

Moreover, previous studies (Berdahl & Fromberg 1982; Gliah et al. 2011; Sun et al. 2017; Li et al. 2017, 2018) used wide-angle integrating detectors subtending almost the whole sky. However, in order to calculate a radiation model for a structure, the irradiation of this structure as a function of its exposed direction to the sky has to be taken into account (see also Sect. 3.2). This effect is not as important for a structure such as a whole building or a solar cell, as studied by Berdahl & Fromberg (1982) or Li et al. (2017) for example. These latter authors therefore needed only integrated values for the entire sky. However, structures of a telescope inside a dome are exposed in a more complex way to different fractions of the sky. Thus, the exact spatial impacts are essential. When tilting the telescope to higher zenith distances, the radiation is increasingly dominated by low heights. Therefore, we derived the correlation of the irradiance with zenith distance for a wide variety of parameters (Fig. 3). A comparison of the summer and winter models shows a very good match of the data points. There are only marginal differences for θ = 60° and 70°. Therefore, again, a general season-independent fit can be derived. We find that this correlation does not depend on humidity after accounting for the zenith dependency mentioned above, and a direct empirical fit function with both parameters for good weather conditions (PWV < 15 mm_H2O and no clouds) is found:

$$\begin{array}{c}\hfill T_{\mathrm{sky}}(T_{\mathrm{amb}},\mathrm{PWV},\theta )=T_{\mathrm{amb}}-45.75+6.52ln(\mathrm{PWV})+0.00025{\theta }^{2.5},\end{array}$$(2)

Fig. 3. Correlation of the sky temperature Tsky and zenith distance at a constant level of PWV.

where T_amb is the ambient temperature near the ground in Kelvin, PWV is in mm_H2O, and the zenith distance θ is measured in degrees.

2.2. Individual solutions and local estimators

Our software molecfit (Smette et al. 2015; Kausch et al. 2015) was originally designed to correct astronomical spectra for telluric absorption. However, it can also be used to simulate individual thermal emission spectra for a certain date and time at a site. The code automatically derives the GDAS profile of that very day and time from the servers.

In order to test extreme weather conditions that are not well represented by the sky model, we calculated emission spectra for the warmest and the coldest nights with respect to local midnight at Cerro Paranal in July 2018 (Fig. 4). The reference temperatures were measured 30 meters above ground at the ambient monitoring station. They agree well with the temperatures that are derived from a fit of a black body to only those regions in the spectrum where the atmosphere is very opaque. For comparison, we also considered the warmest and coldest nights in 2018 at the high-altitude site at Mauna Kea observatory in Hawaii, as measured at the weather monitor of the Canada-France-Hawaii Telescope (CFHT). This observatory site has significantly lower air temperatures due to its high altitude of 4200m (Fig. 5). Therefore, the relative humidity at this altitude tends to be lower. However, the latter and the dew point temperature calculated from it (Buck 1981) are only loosely correlated with the integrated water vapor column, that is, the scatter is large (Otárola et al. 2010). Therefore, we have to use PWV (indicating the impact of the most important and highly variable trace gas H₂O) to optimally parameterize the atmospheric radiation (Otárola et al. 2010; Lakićević et al. 2016).

Fig. 4. Sky radiance at zenith on the coldest and the warmest days at Cerro Paranal (2600m) in July 2018. The fit of a black body is based on the marked points only, defining an upper boundary.

While we were able to take the PWV measured from the Low Humidity And Temperature PROfiler (L-HATPRO) microwave radiometer at Cerro Paranal (Kerber et al. 2012, see Sect. 2.3.1) for Fig. 4, the PWV values for Hawaii are only based on the (less accurate) GDAS profiles for the given nights. The PWV values and the resulting T_sky for the integrated sky spectrum and using Eq. (2) are provided in Table 1. The fit formula causes deviations in T_sky between −3.6 and +4.8 K for these extreme conditions, which suggests that its general uncertainty is not larger than a few Kelvin.

In cases where a sky model is not available at a given site, we tested the possibility of estimating the values of T_sky from the transmission curve. For the most part, the latter depends on the site altitude and the PWV (which can be obtained by means of GDAS data but might also be derived by various multi-wavelength observations without use of a model), and can therefore be easily tabulated for any site. Multiplying this transmission curve with a black body spectrum only depending on the ambient temperature T_amb results in the desired estimator $T_{\mathrm{sky}}^{\mathrm{est}}$ after integration over all wavelengths (Eq. (1)). As indicated in Figs. 4 and 5, such an estimate can almost perfectly match the accurate radiative transfer model result and confirms the parameterization in Eq. (2). Only the ozone band from 9.6 to 10.1 μm is strongly overestimated because most of the ozone emission originates from a higher layer in the atmosphere with a lower temperature. This kind of emission might be different for very low-altitude sites near populated areas, but the latter is beyond our study. Therefore, the ozone-related wavelength range was excluded from the integration. Other regions with special molecular features do not seem to be dominant enough to significantly contribute to our result. The resulting estimated sky temperature $T_{\mathrm{sky}}^{\mathrm{est}}$ in Table 1 does not deviate by more than 2.1 K from T_sky. For the more typical cases represented by the same parameters as used in Figs. 2 and 3, $T_{\mathrm{sky}}^{\mathrm{est}}$ lies about 1.5 K above the value of the full model for PWV <  15mm (Fig. 6).

Fig. 5. Same as Fig. 4, but for Mauna Kea, Hawaii (4200m)

Fig. 6. Comparison of the sky temperature $T_{\mathrm{sky}}^{\mathrm{full}}$ derived from the sky model for the full wavelength range, the temperature Tsky derived from sky models assuming the approximation of pure black body radiation above 30 μm, and the estimated $T_{\mathrm{sky}}^{\mathrm{est}}$ with the ambient temperature black body radiation and the atmospheric transmission curve only.

The ESO Sky Model yields radiation spectra only up to a wavelength of 30 μm. Above this region, pure black body emission is typically assumed in the standard online implementation (SkyCalc). The integration of that emission using Eq. (1) leads to T_sky used in the rest of the paper. This simplification leads to slight overestimation of the radiation for very dry conditions and high sites. However, a complete model by us using a custom patched version of the radiative transfer model ($T_{\mathrm{sky}}^{\mathrm{full}}$) shows that the correction to the formula given above is marginal for PWV > 1.0mm and even below. The discrepancies are of the order of only 0.2K (Fig. 6).

2.3. Ambient meteorological conditions at Cerro Paranal

2.3.1. Methods and data selection

At Cerro Paranal there are several astronomical site monitoring (ASM) systems installed which permanently record the atmospheric ambient conditions. We investigate data from a period of 4.5 years (2015-07-01T00:00:00 to 2019-12-31T23:59:59) based on the ESO Meteo Monitor, the L-HATPRO radiometer (Kerber et al. 2012), and the Differential Image Motion Monitor (DIMM). An overview of these data is given in Sarazin (2017). From the Meteo Monitor we extracted the wind direction, wind speed, and the ambient temperature T_amb measured at a height of 30 m (as in Sect. 2.2) as the telescope structure is more exposed to this part of the atmosphere than towards the ground (cf. Sect. 2.1).

The L-HATPRO radiometer database provides integrated PWV measurements at zenith (θ = 0°) and measurements of the infrared sky temperature T_IR achieved by a specific infrared radiometer mounted at the L-HATPRO, which operates in the 9.6 to 11.5 μm band. As mentioned in Sect. 2.1, the Earth’s atmosphere is widely transparent in this wavelength range. Therefore, T_IR measurements are very sensitive to (comparably warm) clouds at various heights. In particular, ice clouds (cirrus) can be detected, which is not possible with the microwave radiometer. The DIMM incorporates data on the seeing conditions and the sky transparency, the latter being measured by the relative root mean square (rms) of flux variations of reference stars along the line of sight during nighttime (normalized by the average flux).

We merged the data points of these three databases via their time stamps by finding the closest Meteo Monitor and DIMM neighbor to each individual radiometer measurement within Δt ≤ 300 s. By skipping those measurements with Δt >  300 s, we automatically restrict our merged data set to nighttime because the DIMM only provides data taken between dusk and dawn. We also omitted data points with invalid values and those with PWV  ≤ 0.02 mm.

The data set should represent typical observing conditions at Cerro Paranal, that is, omitting poor weather conditions, which is important because the subcooling of the telescope structure (Sect. 3) is strongest in the best (= photometric) conditions. In addition, the ESO Sky Model does not cover cloudy conditions, and therefore Eq. (2) is not applicable. In order to exclude cloudy periods, we therefore apply a DIMM rms cut of 2% (which is consistent with the criterion for possible clouds of the online ASM tool at ESO). This cut reduces our sample by only about 10%; good photometric conditions prevail at Cerro Paranal. In total we select 276 958 data points.

2.3.2. Statistics on the meteorological data

Investigations of the overall yearly variability reveal a strong trend towards higher PWV values in the southern summer months with an annual median value of PWV ≈2.47 mm (10% percentile: 1.06 mm; 25% percentile: 1.56 mm; 75% percentile: 4.08 mm; 90% percentile: 6.39 mm). In addition, weather situations with winds coming from a 90° cone centered on the north are the most frequent ones (56% of the overall time), with a tendency towards higher wind speeds, higher PWV values, and more clouds. As our data selection criteria exclude most of the poor conditions, the sensitivity of our investigations with respect to the wind direction is minimized. Our empirical equation for T_sky (Eq. (2)) depends on PWV and the ambient temperature T_amb. We assume these quantities to be independent. Our data set can be used to test this assumption. Figure 7 does not show a clear correlation. This finding can be tested by the Pearson coefficient r, which indicates whether or not a data set shows a linear correlation (r = ±1 means perfect linear correlation, r = 0 no linear correlation). As the absolute value is small in our case (r = 0.16), we conclude that our parameterization of T_sky is justified. The distribution of T_sky given in Fig. 8 reveals a median value of T_sky = 246.2 K and a standard deviation of σ = 5.7 K for the selected clear observing conditions.

Fig. 7. Two-dimensional histogram of the ambient temperature Tamb as a function of PWV for clear observing conditions (see text for more details). The bin widths of the histogram are 0.1 mm in PWV and 0.5 K in Tamb. The resulting Pearson coefficient r is indicated (see text for more details).

Fig. 8. Histogram of the calculated Tsky values; its mean value is 246.2 K, the standard deviation σ = 5.7 K.

Figure 9 shows the infrared temperature T_IR versus the sky temperature T_sky determined with Eq. (2), which was calculated by means of the T_amb values from the ESO Meteo Monitor at θ = 0°. As T_IR has a relatively weak dependency on PWV but is strongly sensitive to even thin cirrus clouds (due to the small wavelength range used for the measurements), the infrared radiometer temperature is related to the sky temperature T_sky only in the case of photometric observing conditions and T_sky ≳ 240 K. It is therefore not a good estimator for the sky temperature and its effect on the subcooling of telescope structures. In contrast, T_sky derived by Eq. (2) provides a robust estimate for that purpose only based on measurements of the PWV, the ambient temperature, and the zenith distance θ.

Fig. 9. Tsky vs. TIR at θ = 0° as two-dimensional histogram (bin width of 0.5 K) for clear observing conditions.

2.4. Ambient meteorological conditions at Cerro Armazones

In Sect. 3, we compute numerical quantities using typical environmental conditions at Cerro Armazones, the site of the Extremely Large Telescope (ELT; 3 060 m; Koehler 2015; Otárola et al. 2019). The median nighttime air temperature is T_amb = 9.1° C = 282.25 K and the median air pressure equals p = 712 hPa. Measurements of PWV such as those taken with L-HATPRO at Cerro Paranal (Sect. 2.3.1) are not available for Cerro Armazones. Nevertheless, estimates by Otárola et al. (2010) and Lakićević et al. (2016) suggest that the median value for clear sky conditions at Cerro Paranal, namely 2.47 mm (Sect. 2.3.2), should be relatively similar. Possible systematic deviations of the order of half a millimeter are not critical for our purposes. Applying Eq. (2) for the selected PWV value, we obtain a difference between T_sky and T_amb of −39.9 K at zenith and −37.8 K at the typical observing zenith angle θ = 37°. With the ELT median T_amb, this results in T_sky of 242.4 K and 244.5 K, respectively.

Otárola et al. (2019) presented a cumulative wind statistic for Cerro Armazones. However, that study was focused on radio astronomy and thus included daytime data. Fortunately, the complete weather data are publicly available⁴ and we selected the nighttime wind data between 22 : 30 and 9 : 00 h UTC for the years 2006–2009. We find a nocturnal median wind speed at 7 meters above the ground of 7.0 m s⁻¹, 33% of the time the wind speed is below 5.2 m s⁻¹, 25% of the time below 4.3 m s⁻¹ and 15% of the time below 3.0 m s⁻¹.

From a related data set covering the years 2004–2008, we find that the median nocturnal temporal gradient of T_amb is equal to −0.34 K h⁻¹ when the temperature is falling and +0.25 K h⁻¹ when it is rising (which is relevant for significantly shorter time periods). The corresponding quartiles are −0.17 and −0.60 K h⁻¹ (falling) versus +0.11 and +0.51 K h⁻¹ (rising), respectively.

3. Thermal modeling of telescope structures

3.1. Radiative heat flux

In this section, we derive the steady-state temperature of a structure that radiatively cools against the night sky and is heated by forced convection of the ambient air. We consider a small surface of area A that has a spectral emissivity of ϵ(λ, γ), where γ is the angle between the surface normal n and the observing direction. The radiative heat flux from the sky that is absorbed by the surface can be written in spherical coordinates (centered at zenith) according to Lambert’s Cosine Law (Pedrotti & Pedrotti 1993) as

$$\begin{array}{cc}\hfill {\dot{Q}}_{\mathrm{sky}}& =A{\displaystyle {\int }_0^{\pi /2}}\mathrm{d}\theta sin(\theta ){\int }_0^{2\pi }\mathrm{d}\varphi \kappa {\int }_0^{\infty }\mathrm{d}\lambda L_{\mathrm{sky}}(\lambda ,\theta ,\varphi )ϵ(\lambda ,\gamma ),\hfill \end{array}$$(3)

$$\kappa :=\max (cos(\gamma ),0),$$(4)

where L_sky is the spectral radiance from the sky (unit W m⁻² sr⁻¹m⁻¹), γ = γ(θ, ϕ) and ϕ is the azimuth. Both θ and ϕ are measured in radians, in contrast to the previous section. If the surface is oriented horizontally (n points towards zenith), we find γ = θ and κ = cos(γ).

Conversely, to calculate the emitted flux ${\dot{Q}}_{\mathrm{surf}}$ from the surface to the sky, one simply replaces L_sky(λ, θ, ϕ) in Eq. (3) by the black body spectral radiance

$$\begin{array}{c}\hfill L_{\mathrm{bb}}(\lambda ,T_{\mathrm{surf}})=\frac{2h_{\mathrm{Pl}}{c}^2}{{\lambda }^5}{[exp\left(\frac{h_{\mathrm{Pl}}c}{\lambda k_B{T}_{\mathrm{surf}}}\right)-1]}^{-1},\end{array}$$(5)

where T_surf is the surface temperature, h_Pl = 6.62618 × 10⁻³⁴ J s is Planck’s constant, k_B = 1.38065 × 10⁻²³ J/K is Boltzmann’s constant, and c = 2.99792 × 10⁸ m s⁻¹ is the speed of light in vacuum.

As shown in Figs. 4 and 5, the sky spectrum agrees with a black-body fit spectrum within several wavelength ranges (e.g., near 7 μm or near 15 μm). The temperature of this fit spectrum is close to the ambient air temperature at sites at higher altitudes. The radiation imbalance ${\dot{Q}}_{\mathrm{sky}}-{\dot{Q}}_{\mathrm{surf}}$ is therefore dominated by the wavelength bands where L_sky(λ)−L_bb(λ) has a large magnitude, specifically in the astronomical N-band (7 – 13.5 μm) and the Q-band (17 – 25 μm). The area between the spectra and their respective fits in Figs. 4 and 5 provides a graphical representation of this relationship.

The wavelength- and angle-resolved emissivity ϵ(λ, γ) is often unknown. Working with an average emissivity $\overline{ϵ}$ is possible, but it should be noted that its value is defined by the requirement to preserve the value ${\dot{Q}}_{\mathrm{sky}}$ of the wavelength integral in Eq. (3). Any emissivity weighted by the full black-body spectrum, as is sometimes provided by commercial measurement devices or listed in specification sheets, is not well suited to quantify subcooling. Instead, a better approximation is to only consider the emissivity around λ = 11 μm and, if available, around 18 μm.

3.2. Sky view factor

The relationship between spectral radiance and T_sky is (Zhao et al. 2019)

$$\begin{array}{cc}\hfill \sigma ⟨T_{\mathrm{sky}}^4⟩& ={\displaystyle {\int }_0^{\pi /2}\mathrm{d}\theta sin(\theta )cos(\theta ){\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\infty }\mathrm{d}\lambda L_{\mathrm{sky}}(\lambda ,\theta ,\varphi )}\hfill \\ \hfill & =\pi E_{\mathrm{BOL}},\hfill \end{array}$$(6)

in accordance with Eq. (1), where the angular brackets denote the geometrical averaging over the sky hemisphere.

Surfaces that are located inside buildings, such as a telescopes in a dome with open doors, can typically only see a solid angle Ω <  2π of the sky hemisphere. The sky view factor for a point q on the surface is defined by

$$\begin{array}{c}\hfill F_{\mathrm{sky}}=\frac{1}{\pi }{\displaystyle {\int }_0^{\pi /2}\mathrm{d}\theta sin(\theta )\kappa (\theta ,\varphi ){\int }_0^{2\pi }\mathrm{d}\varphi \mathcal{T}(\theta ,\varphi ),}\end{array}$$(7)

where the telescope transmission function 𝒯(θ, ϕ) is equal to 1 at directions θ,  ϕ for which an unobstructed path from q to the sky exists and zero if the path is blocked by any warm black body. However, we also allow for partly obscured ray paths like propagation through a windscreen mesh or indirect ray paths to the sky through reflection on a mirror (often the primary mirror M1), or any partially reflective surface like a bare metal for which 0 <  𝒯 <  1. The normalization factor π⁻¹ ensures that 0 ≤ F_sky ≤ 1.

The sky view factor for a horizontal surface inside a building with a circular aperture centered at n that is seen from q under the one-sided angle θ_a is equal to sin²(θ_a). This term assumes the value 1/2 at θ_a = 45°, although the solid angle of this aperture subtends only 29.3% of the hemisphere. If, on the other hand, the aperture is not centered at n, F_sky will be smaller. This reasoning shows that the orientation of surfaces inside a telescope dome is important, and not simply the distance from the dome slit. The difference is particularly large when comparing the sky view factors for the top side of the spider trusses (with rectangular cross section) or the secondary mirror housing (F_sky ≈ 0.8 – 1) with its side faces where n is perpendicular to the telescope axis. Under no circumstances can the side faces see more than half of the sky hemisphere (F_sky ≤ 1/2). Furthermore, the zenith angles θ at which the cosine weighting term κ is large are almost perpendicular to the telescope axis; more precisely, these rays may often point near the horizon (the side faces may even see part of the ground outside the telescope building). The average sky temperature is higher in these directions due to the θ^2.5 dependence; see Eq. (2).

When dealing with noncircular apertures to the sky (e.g., the open door slit of a telescope dome) as sketched in Fig. 10, it can be more convenient to work with Cartesian coordinates x,  y,  z. The effect of κ in Eq. (7) is simply to project the area of the sky hemisphere that can be seen from the point q onto the unit disk on the tangential surface (perpendicular to n). If we let ρ = (x, y, z) be a vector of unit length that points from q to the sky, and x, y, z are local coordinates, meaning that any vector (x, y, 0) lies on the tangential surface and n = (0, 0, 1), we can write

$$\begin{array}{c}\hfill F_{\mathrm{sky}}=\frac{1}{\pi }{\displaystyle \int }{{\displaystyle \int }}_{x^2+y^2<1}\mathrm{d}x\mathrm{d}y\mathcal{T}(x,y).\end{array}$$(8)

Fig. 10. Evaluation of sky view factors by projection. Red arrow: normal vector n above the point q on a surface. Blue quadrilaterals: two windows of solid angle Ω ≈ 0.16 sr through which the sky can be seen from q (i.e., 𝒯 = 1); the right window is tilted against n by γ = 60°. Green rectangles: projection of the windows onto the tangential plane with areas π Fsky ≈ Ω and Ω × cos(60° ) = Ω/2, respectively.

In buildings, the sky windows often have much more complicated shapes than in Fig. 10. One can numerically evaluate F_sky in this case by describing the edges of the sky windows by a closed set of straight lines. The projected window then forms an (irregular) polygon on the unit disk whose area can be computed easily by the “shoelace formula” (Braden 1986). Dividing this area by π yields F_sky. Any window edges whose rays ρ would point below the tangential plane, that is, if γ ≥ 90° and thus κ = 0, generate projected points that lie on the unit circle in this scheme.

A collection of formulas to calculate F_sky in special geometries is found in Chapter 10 of Lienhard (2011).

If a structure is exposed to some part of the sky through a mirror, or in general through a series of reflective optics, F_sky is defined by the ray bundle emanating from q that ultimately reaches the sky, hence with 𝒯 >  0. In other words, the effective aperture is the image of the real aperture (say, the rim of the dome slit) as seen from q. One can compute F_sky by evaluating the angles γ of the ray segments between q and the first reflective surface. This method is correct even if the mirrors are powered (i.e., have nonzero curvature) by the principle of étendue conservation (Chaves 2017) and assuming that the sky is a Lambertian radiator (setting aside the dependence of T_sky on θ).

3.3. Thermal balance equation

With the help of Eqs. (3)–(7), we can write the radiative heat loss due to sky subcooling as Lienhard (2011)

$$\begin{array}{cc}\hfill {\dot{Q}}_{\mathrm{sky}}-{\dot{Q}}_{\mathrm{surf}}& =A\overline{ϵ}{F}_{\mathrm{sky}}\sigma ({⟨T_{\mathrm{sky}}^4⟩}_{\mathcal{T}}-T_{\mathrm{surf}}^4)<0,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill T_{\mathrm{surf}}& =T_{\mathrm{amb}}+\mathrm{\Delta }{T}_{\mathrm{surf}},\hfill \end{array}$$(10)

where T_amb is the ambient air temperature inside the dome on the upwind side of the telescope structure and ΔT_surf quantifies the surface subcooling. The expression ${⟨T_{\mathrm{sky}}^4⟩}_{\mathcal{T}}$ denotes the average sky temperature with the same transmission 𝒯 as in Eqs. (7) and (8). As discussed in the previous section, this average geometrically depends on the patch of sky seen from q and the telescope pointing.

On the other hand, convection heats the surface at a rate of

$$\begin{array}{c}\hfill {\dot{Q}}_{\mathrm{conv}}=Ah(T_{\mathrm{amb}}-T_{\mathrm{surf}})=-Ah\mathrm{\Delta }{T}_{\mathrm{surf}}>0,\end{array}$$(11)

where h is the interface heat transfer coefficient (in units of W K⁻¹ m⁻²). In thermal equilibrium and neglecting thermal conduction, we require that

$$\begin{array}{c}\hfill {\dot{Q}}_{\mathrm{sky}}-{\dot{Q}}_{\mathrm{surf}}+{\dot{Q}}_{\mathrm{conv}}=0,\end{array}$$(12)

and therefore the heat fluxes are balanced. Our goal is to quantify the steady-state plate subcooling ΔT_surf <  0 that satisfies

$$\begin{array}{cc}\hfill \mathrm{\Delta }{T}_{\mathrm{surf}}& =\sigma \omega ({⟨T_{\mathrm{sky}}^4⟩}_{\mathcal{T}}-(T_{\mathrm{amb}}+\mathrm{\Delta }{T}_{\mathrm{surf}}{)}^4),\hfill \\ \hfill \omega & :=\frac{\overline{ϵ}{F}_{\mathrm{sky}}}{h}.\hfill \end{array}$$(13)

Expression (13) is a quartic equation with two real and two complex roots. One of the two real roots lies below −1300 K, and so there is only one meaningful solution. We calculated the latter in closed form, but the expression is rather complicated. Instead, we carry out a Taylor expansion to first order in ΔT_surf about ΔT_surf = δT. Solving for ΔT_surf, we obtain

$$\begin{array}{cc}& \mathrm{\Delta }{T}_{\mathrm{surf}}=\frac{\sigma \omega }{1+4T_0^3\sigma \omega }({⟨T_{\mathrm{sky}}^4⟩}_{\mathcal{T}}-T_0^3(T_{\mathrm{amb}}-3\delta T)),\hfill \\ \hfill & T_0:=T_{\mathrm{amb}}+\delta T.\hfill \end{array}$$(14)

In the remainder of this paper, we choose δT = −3 K so that we cover the typical range of subcooling in telescopes of about −6 K ≤ ΔT_surf ≤ 0 with good accuracy.

Following Lienhard (2011), we define the radiation heat transfer coefficient

$$\begin{array}{c}\hfill h_{\mathrm{rad}}:=4T_0^3\sigma \overline{ϵ}{F}_{\mathrm{sky}}\approx 4.9\mathrm{W}{\mathrm{K}}^{-1}{\mathrm{m}}^{-2}\times \overline{ϵ}{F}_{\mathrm{sky}},\end{array}$$(15)

and, rewriting Eq. (14), finally express ΔT_surf as

$$\begin{array}{cc}& \mathrm{\Delta }{T}_{\mathrm{surf}}=\frac{\eta }{1+\eta }\mathrm{\Delta }{T}_D,\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill & \mathrm{\Delta }{T}_D:=\frac{{⟨T_{\mathrm{sky}}^4⟩}_{\mathcal{T}}}{4T_0^3}-\frac{T_0-4\delta T}{4}\approx 0.89\mathrm{\Delta }{T}_{\mathrm{sky}},\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill & \eta :=\frac{h_{\mathrm{rad}}}{h},\hfill \end{array}$$(18)

where the number η quantifies the subcooling efficiency. We thus find that ΔT_surf scales linearly with η as long as η ≪ 1. We highlight that ΔT_surf ∝ 1/h in this unsaturated regime, which applies to most cases of interest in telescope subcooling. At ELT median conditions (Sect. 2.4), η is of the order of 0.1 and the subcooling is driven by the temperature offset ΔT_D which equals −33.9 K for θ = 37°. In the opposite limit of η​ → ​∞, we obtain ΔT_surf​ → ​ΔT_sky := T_sky − T_amb. In practice, this extreme can probably only be reached on a black surface facing the sky inside a vacuum vessel with a window that transmits infrared radiation.

3.4. Convective heat transfer

The thermal coupling of a structure to the ambient air through forced convection gives rise to a boundary flow problem (Martinez 2020a). The convective heat transfer coefficient h depends on the air speed v near the surface (which is often much lower than the wind speed outside the dome), the air pressure p, and, to a great extent, on the shape of the structure. The structures closest to or even inside the telescope beam are usually steel trusses supporting the secondary and possibly the tertiary or quaternary mirrors. The simplest truss shape is cylindrical with diameter d. An approximate expression for h_cyl, modeling the net heat exchange of airflow perpendicular to a cylinder, was derived by Churchill & Bernstein (1977). Their original formula is expressed in terms of the Reynolds, Nusselt, and Prandtl numbers (Martinez 2020b) and we checked that the values of these numbers for typical telescope structures are within the validity range of the formula. We do not reproduce the general formula here but instead simplify it to

$$\begin{array}{cc}& h_{\mathrm{cyl}}=c_1\frac{\sqrt{\beta }}{s}{(1+c_2{\beta }^{5/8})}^{4/5},\hfill \\ \hfill & \beta :=mspv,\hfill \end{array}$$(19)

where c₁ = 0.0179 W K⁻¹ m⁻¹, c₂ = 1.5 × 10⁻⁴ for airflow at T_amb = 282 K, and s = π d is the cylinder circumference. We introduced the Reynolds number scaling quantity m, which equals 1.0 m W⁻¹ if the incoming airflow is perfectly laminar and the cylinder has a smooth surface, and m ≈ 1.5 m W⁻¹ for partly turbulent inflow. However, if a precision of more than about 20% in h_cyl is demanded, one has to resort to computational fluid dynamics (CFD) simulations.

At v = 7.6 m s⁻¹ and s = 0.6 π = 1.88 m, one reaches c₂ β^5/8 = 1. For much lower air speeds (more precisely, if β ≪ 1.3 × 10⁶), the scaling law h_cyl ∝ (m p v/s)^1/2 applies. We note that, although v = 7.6 m s⁻¹ is close to the median wind speed on Cerro Armazones (cf. Sect. 2.4), the air speed inside the dome is often much lower. This holds in particular for structures below the M2 level, and therefore this regime is normally valid.

In the opposite limit of large β, for example for high air speed and/or large structures, h_cyl ∝ β/s = m p v. The Reynolds number then becomes so large that the boundary layer is turbulent throughout, and thus the structure size is no longer relevant for h_cyl. An example for this limit is subcooling of the outer dome cladding, where v approaches the free air wind speed and the interface length s can amount to tens of meters for streamlines that follow the dome surface.

Figure 11 shows a plot of h_cyl for three different values of s representative of structures in the VLT and ELT telescopes and the ELT dome door (all for the ELT environmental conditions and m = 1.25 m W⁻¹). Sundén et al. (2010) contains an overview of h for cylinders with noncircular cross sections.

Fig. 11. Convective heat transfer coefficient hcyl vs. local air speed v for three different interface lengths s representative of structures in the VLT and ELT telescopes, plus the ELT dome door.

Figure 12 shows a plot of ΔT_surf as a function of $\overline{ϵ}{F}_{\mathrm{sky}}$ according to Eq. (16), where the colors indicate four different local air speeds v. For these calculations and the remainder of this paper, we assume ELT median conditions as defined in Sect. 2.4 for the zenith angle θ = 37°. The colored fans in Fig. 12 are meant to give an impression of the influence of the interface length s : As s grows from s = 0.4 π m = 1.26 m to s = 4.3 m, h_cyl diminishes and therefore η grows, and thus the subcooling intensifies. Consequently, subcooling becomes more severe as telescopes grow.

Fig. 12. Steady-state surface subcooling ΔTsurf as a function of $\overline{ϵ}{F}_{\mathrm{sky}}$ for four different local air speeds v. The upper edges of the colored fans correspond to s = 1.26 m, the dashed lines to s = 1.88 m, and the lower edges to s = 4.3 m, respectively.

Figure 13 is a two-dimensional plot of T_surf(F_sky, v) for s = 4.3 m (ELT top ring truss circumference) and $\overline{ϵ}=0.2$, which is typical for dedicated low-emissivity paints; see Yarbrough (2010) for an overview of some respective commercial products. The tube structures of several large telescopes such as Gran Telescopio Canarias (GTC, La Palma), Keck I+II (Mauna Kea, Hawaii), and the two Gemini telescopes (Cerro Pachón, Chile, and Mauna Kea) are covered with LO/MIT⁵ paint by the company Solec (a silvery paint with small aluminum particles). On surfaces that scatter visible stray light, one may apply selective foil such as Acktar NanoBlack^{® 6} or Maxorb (Mason et al. 1982). These foils are black at visible wavelengths but highly reflective and thus low in emissivity ($\overline{ϵ}<10\%$) in the thermal infrared. Combining Eq. (19) with Eq. (16), we find a simplified fit function near $\overline{ϵ}=0.2$ and F_sky = 0.5 given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{T}_{\mathrm{surf}}\approx 0.636\mathrm{\Delta }{T}_{\mathrm{sky}}\overline{ϵ}{F}_{\mathrm{sky}}{s}^{0.355}{v}^{-0.551},\end{array}$$(20)

Fig. 13. Steady-state surface subcooling ΔTsurf as a function of Fsky and v for s = 4.3 m and $\overline{ϵ}=0.2$. Contours: increments of −1 K; thick contour: commonly accepted surface subcooling limit in telescopes of ΔTsurf = −1.5 K.

where s is entered in meters and v is in meters per second, and again ΔT_sky = T_sky − T_amb. We found the power laws of s and v using the rule b = x f′/f for a function f(x) = a x^b with derivative f′, applied at s = 0.6 π m = 1.88 m and v = 1 m s⁻¹.

3.5. Optical path difference

Temperature differences ΔT in an air volume cause the optical path length difference (OPD) of

$$\begin{array}{c}\hfill \mathrm{OPD}=\frac{\mathrm{d}n}{\mathrm{d}T}{\displaystyle \int }\mathrm{\Delta }T(\mathit{q}(l))\mathrm{d}l,\end{array}$$(21)

where q is a point in the air volume and l measures geometrical distance along the ray path. The refractive index of air n(λ, T, p) as a function of temperature T and pressure p is given by Edlen’s modified formula (Bönsch & Potulski 1998). Its derivative by T can be written in simplified form as

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}T}=\frac{1}{{\rho }_0}\frac{\mathrm{d}\rho }{\mathrm{d}T}=-\frac{p}{{\rho }_0{R}_{\mathrm{air}}{T}^2},\end{array}$$(22)

where ρ = p/(R_air T) denotes the air density (R_air = 287.06 J kg⁻¹ K⁻¹ is the gas constant of dry air, yielding ρ = 0.879 kg m⁻³ at ELT median conditions) and the constant ρ₀ = 4450 kg m⁻³ in the wavelength range 0.4 μm ≤ λ ≤ 3 μm.

At ELT median conditions, we obtain dn/dT = −9.828 × 10⁻¹² K⁻¹ Pa⁻¹ × p which is equal to −7.0 × 10⁻⁷ K⁻¹. For the ray path from the top of the ELT secondary mirror (M2) crown level to M1, then back to M2 and finally to the prime focus of total length 36 m + 31 m + 13 m = 80 m, an air cooling as small as ΔT = 25 mK causes an accumulated OPD of 1.4 μm, equivalent to about 0.64 λ in K-band (2 – 2.4 μm). The thermal OPD is thus significant and may actually become the leading wavefront error contribution in 30-meter-class telescopes and greater.

In the following, we establish the link between the subcooling temperature offset ΔT_surf and the OPD. We rewrite the surface area in Eq. (11) as A = s L, where L is the length of a cylindrical truss that we assume to be oriented perpendicular to the local airflow vector v and perpendicular to the telescope optical axis, as would be the case for example for part of the top ring, as shown in Fig. 14.

Fig. 14. Slice through a cylindrical truss of diameter d = s/π running perpendicular to v and to the telescope axis. Red arrow: ray of light traversing the cool wake of thickness δ and temperature difference ΔTwake.

On the downwind side of the truss, a wake of cooled air is generated. We consider now a ray of light in a thin air slice perpendicular to v and assume for simplicity that the wake has a thickness along the ray of δ with constant temperature offset ΔT_wake <  0. The total volume of cooled air flowing away from the truss per time Δt is then equal to V/Δt = δ L v and transports the heat flow

$$\begin{array}{c}\hfill -{\dot{Q}}_{\mathrm{conv}}=\frac{V}{\mathrm{\Delta }t}{c}_{\mathrm{p}}\rho \mathrm{\Delta }{T}_{\mathrm{wake}}=c_{\mathrm{p}}\rho vL\delta \mathrm{\Delta }{T}_{\mathrm{wake}},\end{array}$$(23)

where c_p = 1005 J kg⁻¹ K⁻¹ is the specific heat capacity of air at constant pressure.

We can then derive the OPD imposed on the ray by the cooled air generated by a single truss of circumference s in a single ray pass by combining Eq. (23) with Eq. (11)

$$\begin{array}{c}\hfill \mathrm{OPD}=\frac{\mathrm{d}n}{\mathrm{d}T}\delta \mathrm{\Delta }{T}_{\mathrm{wake}}=\frac{\mathrm{d}n}{\mathrm{d}T}\frac{sh\mathrm{\Delta }{T}_{\mathrm{surf}}}{c_{\mathrm{p}}\rho v}.\end{array}$$(24)

In reality, the wake is of course not isothermal, but instead the temperature distribution depends on the aerodynamics and is nonstationary. However, the actual temperature distribution along the ray is immaterial for the problem at hand since the path integral in Eq. (21) is a conserved quantity equal to δ ΔT_wake when averaged over time and along the truss length L. In the absence of this spatial and temporal averaging, transient vortexes can distribute air pockets of different temperature in the plane orthogonal to the rays (Martinez 2020b) generating a transient wavefront error; however, this aspect is beyond the scope of this work.

When h_rad ≪ h (e.g., with low-emissivity coatings and under most airflow conditions inside the dome) meaning that ΔT_surf is proportional to ΔT_D/h, and as dn/dT is proportional to p and thus to ρ, we can simplify Eq. (24) to

$$\begin{array}{c}\hfill \mathrm{OPD}\approx v_0\frac{\overline{ϵ}{F}_{\mathrm{sky}}s}{v}.\end{array}$$(25)

For ELT median conditions, v₀ ≈ 10⁻⁷ m/s. We point out that the OPD does not depend on h : In steady state, the heat power lost by the radiation imbalance is exactly compensated by convective heat transfer from the inflowing air so that h cancels.

Because the average air cooling is inversely proportional to the air mass passing the telescope structure per time, OPD ∝ v⁻¹. This fundamental relationship holds independently of the actual values of the structural subcooling ΔT_surf that various parts of the telescope converge to in steady state and explains the name “Low Wind Effect” (Sauvage et al. 2015; Milli et al. 2018); at very low wind speeds, forced convection gives way to natural convection so that the limit v → 0 is never reached.

By the same argument, the OPD does not depend on air pressure either, and therefore telescopes at high altitude are in principle afflicted by the same magnitude of thermal OPD as at the ground level (more specifically, at high altitude the PWV column tends to be smaller, and therefore ΔT_sky grows in magnitude and v₀ increases).

The structures in a telescope causing large OPD are those with high-emissivity coatings (e.g., any regular industrial paint, irrespective of color) and those with large dimensions and/or high ratios F_sky/v. The sky view factor F_sky tends to be larger for structures closer to the dome slit than those near the primary mirror, but the same holds for v (Cullum et al. 2000). Moreover, we remind the reader that structures facing M1 may see a significant fraction of the sky in reflection. Therefore, it is not necessarily clear which telescope structures are the most critical in terms of causing subcooling OPD.

Figure 15 shows a two-dimensional plot of the OPD as a function of F_sky and v for a single ray pass through the wake of a single truss of circumference s = 4.3 m (ELT top ring) and $\overline{ϵ}=0.2$, again at ELT median conditions. The OPD from a single truss of typical size therefore amounts to only a few tens of nanometers. Moreover, the OPD caused by the top ring for example may tend to spread rather uniformly across the pupil, in which case it would not strongly distort the wavefront of starlight. However, the situation is different for vertically oriented structures such as vertical trusses in the telescope tube or a tertiary/quaternary mirror tower. Assuming that the cold wake thickness is similar to the truss diameter δ ≈ d, the OPD for structures running parallel to the optical axis must be scaled by a factor of up to the truss aspect ratio, for instance L/d = 22/0.6 ≈ 37 for one of the vertical ELT telescope tube trusses.

Fig. 15. Steady-state OPD from subcooling as a function of Fsky and v for a single ray pass through the wake of a single truss of circumference s = 4.3 m and $\overline{ϵ}=0.2$ (with ELT median conditions). Contours: OPD in nanometers in increments of factors of 101/5 ≈ 1.585, starting at 10 nm.

To estimate how the OPD from subcooling scales with telescope size in general, we assume that telescopes of increasing sizes are supported by trusses whose circumference s and length L in the direction parallel to the telescope axis both scale linearly with the primary mirror diameter D, causing the OPD from the wake of each truss to scale like D². This relationship holds for the structure above the M1 level only, as long as we assume a roughly identical focal ratio of the primary mirrors. However, the tube structure of the upcoming generation of 30-meter-class telescopes and greater is becoming more delicate than that of the preceding 8- to 10-meter-class telescopes: The number of trusses increases, while each truss circumference grows only modestly. In the somewhat unrealistic limit that s does not grow at all and L rises proportionally to D, the OPD scales only linearly with D (but spreads more widely across the pupil). In practice, the scaling power law will lie somewhere in between these extremes.

The rays from a star at center field propagating through a rotationally symmetric telescope travel in a plane defined by the entry point of the ray in the pupil and the optical axis. The closer the entry point is to the central obscuration (CO), the closer the ray segments between reflections are to each other. The air volume that in vertical projection is close to the rim of the CO is therefore traversed at least three times, boosting the OPD on the downwind side of the CO. Unfortunately, the M2 crown and the tertiary mirror tower introduce a lot of vertical surface area in this region with high F_sky. Moreover, any vertically extended tertiary mirror tower and the M2 crown are optically conjugated to a large range of heights above the telescope, which can induce field-dependent errors in the adaptive optics correction. A quantitative analysis requires extensive CFD simulations of the entire dome volume with temperatures (Cho et al. 2011; Vogiatzis et al. 2014, 2018; Ladd et al. 2016; Oser et al. 2018).

3.6. Heat conduction

So far, we have neglected heat conduction within the structures, which will tend to equalize the structure temperature and thus spread the subcooling effect. In Appendix A, we derive the characteristic length scale due to heat conduction ξ with typical values of a few tens of centimeters in telescope trusses (see Table A.1).

The OPD is approximately proportional to $\overline{ϵ}{F}_{\mathrm{sky}}s$ as Eq. (25) shows. Therefore, it is fortunately equivalent to average the temperature around the perimeter of trusses when estimating the OPD in its wake, or for rectangular cross sections for example, to sum the OPD contributions of the four side faces. However, one should not assume that the subcooling spreads evenly along the length of a truss, which typically far exceeds ξ.

3.7. Transient effects

Telescope structures have a significant thermal inertia that delays temperature changes. In Appendix B, we derive the characteristic thermal relaxation timescale τ with typical values of a few hours (see Table A.1).

Sky subcooling is therefore a slow process that generally increases during the night. We note that any incorrect temperature conditioning of the telescope at sunset, including overheating (ΔT_surf >  0), will generate additional OPD, starting already at the beginning of the observations (Wilson 2013).

We point out that in contrast to the characteristic length ξ, the relaxation time τ is linearly dependent on w/h, where w is the truss wall thickness, and thus varies more strongly among different structures in a telescope, even if they are made of the same material. Additionally, the wind veers, the sky temperature drifts, there is variable cloud coverage, the telescope pointing zenith angle varies and so on, and therefore sky subcooling is an inherently transient process and any steady-state computation is only an approximation.

A serious source of wavefront errors in some telescopes is “mirror seeing” (Wilson 2013) caused by ΔT_surf on the optical mirror surface, which can sometimes become even more detrimental than “dome seeing” (Racine et al. 1991). However, mirror seeing is not normally linked to radiative cooling because typical mirrors have broadband coatings, and therefore very low emissivity such as $\overline{ϵ}<2$%. However, subcooling can become a problem if the uncoated backside of a secondary mirror is (partly) exposed to the night sky (the emissivity of typical mirror substrates made of Zerodur^® is ≈ 0.9).

To compound the complexity, the ambient air temperature T_amb also fluctuates. During a typical night, T_amb drops quickly after sunset and decreases later more slowly, meaning that the radiation imbalance from the cold sky initially accelerates the thermalization of the telescope. Conversely, a rise in T_amb exacerbates the subcooling. The median temporal gradients of T_amb given in Sect. 2.4 are smaller than typical subcooling gradients of the order of Ṫ_amb ≈ −1.5 K/τ ≈ −1 K h⁻¹, but are not negligible.

One can shorten the subcooling relaxation time whilst lowering the emissivity and slowing down the structural deflection of the telescope structure by covering the truss surfaces with a thermal insulation layer topped by aluminum foil or thin sheet metal; a strategy chosen at the Subaru Telescope on Mauna Kea, Hawaii (Bely 2006). However, we note that applying an insulation layer on top of a large thermal mass does not diminish the subcooling (i.e., does not mitigate ΔT_surf in the steady state) because the radiation imbalance on its surface is not reduced.

4. Conclusions

Subcooling of structures on the ground is caused by a thermal-radiation imbalance between the night sky and a structure on the ground at ambient temperature which can reach 150 W m⁻² and above. Structural elements in telescopes therefore cool by a few Kelvin (see Eqs. (16) and (20)) and in turn cool the surrounding air. Wakes of cooled air flowing from elements in or near the telescope beam, such as spider trusses, the top ring, the tertiary mirror tower, and so on, cause an optical path length difference (OPD), introducing a transient wavefront error. While the size of modern telescopes continues to grow, the budgets of acceptable wavefront error are diminishing for diffraction-limited science cases. Therefore, subcooling poses a serious challenge to the generation of 30-meter telescopes and beyond, and second-generation instruments with adaptive-optics correction in existing 10-meter-class telescopes.

Additionally, we find the following:

• The integrated thermal sky radiation can be characterized by the effective bolometric sky temperature T_sky. This temperature is dominated by the lower atmosphere and T_sky therefore floats with the ambient air temperature T_amb near the ground. Based on spectra from the Cerro Paranal Sky Model, we derived a simple fit formula for T_sky as a function of T_amb, the PWV column height, and the zenith angle θ; see Eq. (2). The formula is valid for clear sky conditions, which are present at Cerro Paranal for about 90% of the time, and can be applied to any relatively dry site. For the investigated range of PWV values from 0.5 to 20 mm, the difference between T_sky and T_amb varies from −50 K to −25 K for θ = 0°. The discrepancy can decrease by several Kelvin for higher zenith angles.

• For Cerro Armazones, the site of the ELT, we estimate a median T_sky of about 244 K at a zenith angle of 37°, which lies about 38 K below T_amb.

• The OPD due to subcooling of a single truss of circumference s is proportional to $\overline{ϵ}{F}_{\mathrm{sky}}s/\mathit{v}$, where $\overline{ϵ}$ denotes surface emissivity, F_sky is the sky view factor averaged around the truss cross section, and v is the local air speed; see Eqs. (24) and (25). As large telescopes are generally supported by trusses with both larger circumference and higher (vertical) length, the OPD scales between linearly and quadratically with the primary mirror diameter.

• It is crucial to apply low-emissivity coatings to all telescope structures above the primary mirror level that are exposed to a large fraction of the sky. Good choices are either bare aluminum or paints with embedded aluminum particles. On surfaces that may scatter visible stray light, one can apply selective foil that is black in the visible but reflective in the thermal infrared such as Acktar NanoBlack^® or Maxorb. Conversely, regular industrial paint usually has high emissivity ($\overline{ϵ}>0.9$), even if it is white.

• Sufficient airflow should be maintained in the telescope dome from sunset until the end of the night, for example by opening louvers. Because $\mathrm{OPD}\propto \overline{ϵ}{F}_{\mathrm{sky}}/\mathit{v}$, the wavefront error from subcooling is not just caused by surfaces near the level of the dome slit, but also by any structures closer to the primary mirror if they are poorly ventilated and/or have high emissivity.

• Special attention should be devoted to the outer dome cladding because of the combination of high sky view factor (F_sky ≈ 1) and large surface area. It is advisable to use a low-emissivity coating, additionally with low solar absorptivity to limit daytime heating. Again, aluminum alloys or paints with aluminum pigments are good choices, while any colored paints should be avoided.

• The most critical areas in a telescope structure in terms of wavefront distortion due to subcooling tend to be the vertical ones (i.e., running parallel to the optical axis), in particular surfaces near the projected rim of the central obscuration and the lateral spider truss faces.

• The ratio of heat capacity to convective efficiency of any telescope structure sets the timescale of thermal relaxation, but these quantities do not influence the steady-state subcooling temperature, as discussed in Sect. 3.7. The time to reach thermal equilibrium can be hours, and therefore subcooling often becomes more serious in the second half of the night.

• Handheld commercial thermal cameras are not suitable to measure the bolometric sky temperature. These cameras normally use microbolometer detectors (Rogalski 2019), which are only sensitive in the wavelength range of about 8–14 μm which largely overlaps with the astronomical N-band, a region with particularly high atmospheric transmission. On the other hand, if such cameras are employed to estimate subcooling of telescope structures, it is best to take the measurement at the end of the night directly after closing the dome to avoid measuring spuriously low temperatures due to partial sky reflection.

Acknowledgments

W. Kausch is funded by the Hochschulraumstrukturmittel provided by the Austrian Federal Ministry of Education, Science and Research (BMBWF). S. Noll was financed by the project NO 1328/1-1 of the German Research Foundation (DFG). R. Holzlöhner wishes to thank M. Brinkmann (ESO) for stimulating discussions on fluid dynamics and A. Otárola (TMT) for help with the Cerro Armazones environmental data statistics. The authors further thank J. Vinther (ESO) for his continued support of SkyCalc.

References

Appendix A: Heat conduction length scale

We consider a plate and a pipe-shaped truss, both with wall thickness w, letting y be the local surface tangent coordinate along the truss length L and x run along the orthogonal tangential direction. We can write the two-dimensional heat equation for the temperature field u(x, y):=ΔT_surf(x, y) as

$$\begin{array}{cc}& -k_m(u_{\mathit{xx}}+u_{\mathit{yy}})=q_{\mathrm{sky}}-q_{\mathrm{surf}}+q_{\mathrm{conv}}\hfill \\ \hfill & =-\frac{h}{w}[(1+\eta )u-\eta \mathrm{\Delta }{T}_D],\hfill \end{array}$$(A.1)

where u_xx, u_yy denote the respective second partial derivatives of u by x and y, $q=\dot{Q}/(A\mathit{w})$ is the specific lateral heat flow per bulk wall volume (unit W m⁻³) and we have inserted the heat flow terms on the left-hand side of Eq. (12). In Eq. (A.1), we use Eq. (15) where now h_rad = h_rad(x, y), k_m is the thermal conductivity of the truss wall material and we assume that the heat conduction normal to the surface is so efficient that the wall becomes isothermal along z, thus u(z) = const and k_m/w ≫ h.

We define the characteristic length ξ which quantifies the ratio of lateral heat conduction in the truss wall to convection and rewrite Eq. (A.1) as

$$\begin{array}{cc}& {\xi }^2(u_{\mathit{xx}}+u_{\mathit{yy}})-u(x,y)+U_{\mathrm{rad}}(x,y)=0,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill & \xi :=\sqrt{\frac{w{k}_m}{h}},\hfill \end{array}$$(A.3)

where we assume again h_rad ≪ h as discussed after Eq. (15) and the radiative forcing term U_rad(x, y) is equal to η ΔT_D.

We now consider two different simple structures: The first is a rectangular plate of dimensions L × s × w with L ≫ s ≫ w. We require that there be no heat flux across the plate edges, hence the temperature becomes stationary there. Further, the radiative forcing term contains a step function centered at x = 0, for example because the plate is bent into an L-shaped form with a 90° corner running parallel to the long axis y where the sky view factor changes abruptly:

$$\begin{array}{cc}\hfill u_x(\pm s/2)& =0,\hfill \\ \hfill U_{\mathrm{rad}}(x,y)& :=\mathrm{\Delta }U\mathrm{sgn}(x)+U_0,\hfill \end{array}$$(A.4)

where u_x is the partial derivative of u by x and sgn(x) = 1 for x >  0, sgn(x) = − 1 for x <  0 and 0 otherwise. The solution of Eq. (A.2) with the boundary conditions in Eq. (A.4) is

$$\begin{array}{c}\hfill u(x^{\prime })=\mathrm{\Delta }U[tanh(s^{\prime })sinh(x^{\prime })-\mathrm{sgn}(x)(cosh(x^{\prime })-1)]+U_0,\end{array}$$(A.5)

where tanh, sinh, and cosh denote the hyperbolic equivalents of the tan, sin, and cos functions, respectively, and we have introduced the reduced variables x′: = x/ξ and s′: = s/(2ξ).

The blue and red curves in Fig. A.1 show u(x), where s = π 0.6 m ≈ 1.88 m, which has the shape of an inverted letter S with its inflection point at (x​ = ​0,  u​ = ​U₀). For this plot we used the material properties of construction steel S355 as listed in Table A.1 and chose w = 10 mm, ΔU = −1 K, and U₀ = −2 K. We assumed the two different conductive heat transfer coefficients h₁ = 6 W K⁻¹ m⁻² and h₂ = 3 W K⁻¹ m⁻², corresponding to the scaling lengths ξ₁ = 0.27 m and ξ₂ = 0.39 m (blue and red curves, respectively).

Fig. A.1. Example of steady-state temperature distribution u(x) with heat conduction. Blue and red curves: plate of size π 0.6 m with step-like radiative forcing; ξ = 0.27 m and 0.39 m, respectively. Turquoise and orange: same ξ values, but sinusoidal forcing around pipe-shaped truss of diameter d = s/π = 0.6 m. The transition region of the blue curve has the width Δx = 0.55 m.

We can define a transition region width Δx satisfying u_x(0) Δx/2 = ΔU, as indicated by the blue dotted line in Fig. A.1. Using u_x(0) = ΔU tanh(s′)/ξ, we find

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }x}{s}=\frac{1}{s^{\prime }tanh(s^{\prime })}\approx 0.33-0.11(s^{\prime }-1)+O(s^{\prime 2}),\end{array}$$(A.6)

where the last term is the Taylor expansion about the value s′ = 3 relevant for the presented examples. We conclude that the transition region can cover a significant fraction of a typical steel plate.

The second structure is a pipe-shaped truss of diameter d = s/π with sinusoidal radiative forcing around its circumference instead of the step-like forcing, which could be a model of a vertical truss in the telescope tube, where F_sky peaks on the side facing the telescope axis and becomes minimal on the opposite side oriented towards the inside of the dome. We modify the plate boundary conditions so that the problem becomes periodic in the circumferential coordinate x, as already analyzed by Joseph Fourier (1822)

$$\begin{array}{cc}& u(s/2)=u(-s/2),\hfill \\ \hfill & u_x(s/2)=0,\hfill \\ \hfill & U_{\mathrm{rad}}(x,y)=\mathrm{\Delta }Ucos(2\pi jx/s)+U_0,\hfill \end{array}$$(A.7)

where j is any integer. As in the plate problem, U_rad(x, y) lies between the extremes U₀ + ΔU and U₀ − ΔU, which are assumed now at x = 0 and x = ±s/2, respectively. The solution is

$$\begin{array}{cc}& u(x^{\prime })={\mathrm{\Psi }}_j\mathrm{\Delta }Ucos(\pi j{x}^{\prime }/s^{\prime })+U_0,\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill & {\mathrm{\Psi }}_j:=\frac{s^{\prime 2}}{s^{\prime 2}+{\pi }^2{j}^2}=\frac{1}{1+(2\pi j\xi /s{)}^2},\hfill \end{array}$$(A.9)

which is plotted in Fig. A.1 in the turquoise and orange curves for j = 1 and s = 0.6 π m = 1.88 m and again setting ξ₁ = 0.27 m and ξ₂ = 0.39 m, respectively (as u(x′) is periodic in the pipe case, we have chosen to shift it along x so that its inflection point lies at x = 0 in the plot, as for the plate solution). The factor Ψ_j compresses the amplitude of the sinusoidal temperature variation around the pipe, in particular if s <  2π j ξ. We note that we can approximate any periodic boundary condition by a Fourier series, hence the linear combination of solutions in Eq. (A.9) with different spatial frequencies j.

If a plate is thermally insulated on one side, such as can be the case for the outer dome cladding, the convective heat exchange on the insulated side is almost turned off. Conversely, for plates whose two side faces are exposed to the ambient air, the interface surface area A is twice that of the insulated case and F_sky should be averaged over the two sides. However, in such a case the average F_sky cannot exceed a value of one-half because of geometric constraints, and therefore the product F_sky A may in some cases not differ strongly from the insulated plate case.

Inside a sealed air-filled structure such as a pipe truss, the enclosed air exchanges heat energy with the internal walls and there is also internal thermal radiation. However, as the internal airflow is usually not forced and because the wall temperature variation is normally too small to induce a significant radiation imbalance, we can neglect these factors.

Summarizing the plate and the pipe examples, we find that u(x) becomes compressed if the characteristic length ξ = (w k_m/h)^1/2 is large compared to s, that is, if lateral conduction dominates over convection. Irrespective of this condition, the average of u(x) around the truss circumference is always equal to the average of U_rad, thus U₀ = −2 K in our two examples. The limit s ≪ 2π ξ may be referred to as the isothermal regime, while the opposite limit s ≫ 2π ξ may be thought of as thermally disconnected. For a typical telescope pipe steel truss with a wall thickness of w = 10 mm, neither limit clearly applies around its circumference. However, the thermally disconnected regime applies along a typical length L ≫ 2π ξ.

Appendix B: Thermal inertia

Until now, we have discussed only the steady-state solution. We now finally consider the full time-dependent heat equation for a three-dimensional body in air with temperature offset u = u(x, y, z, t):=T_body(x, y, z, t)−T_amb(t) with time t,

$$\begin{array}{c}\hfill c_m{\rho }_m(u_t+{\dot{T}}_{\mathrm{amb}})-k_m(u_{\mathit{xx}}+u_{\mathit{yy}}+u_{\mathit{zz}})=q_{\mathrm{sky}}-q_{\mathrm{surf}}+q_{\mathrm{conv}},\end{array}$$(B.1)

where c_m and ρ_m are the specific heat capacity and the density of the body material, respectively.

We can consider the pipe problem of the previous section with j = 1 and the initial condition u(x, y, z, 0) = 0 and, setting u_zz = Ṫ_amb = 0, find the solution

$$\begin{array}{c}\hfill u(x^{\prime },t\prime )={\mathrm{\Psi }}_1\mathrm{\Delta }Ucos(\pi x^{\prime }/s^{\prime })(1-e^{-t^{\prime }/{\mathrm{\Psi }}_1})+U_0(1-e^{-t^{\prime }})\end{array}$$(B.2)

as a function of the reduced time t′ = t/τ with the lateral thermal relaxation time

$$\begin{array}{c}\hfill \tau :=\frac{w{c}_m{\rho }_m}{h}.\end{array}$$(B.3)

We now evaluate the time t = τ after which the term (1 − exp(−t/τ)) grows from zero to 1 − 1/e ≈ 0.63. We assume again the material properties of construction steel S355 listed in Table A.1 and w = 10 mm to arrive at τ₁ = 1.64 h and τ₂ = 3.27 h for h₁ = 6 W K⁻¹ m⁻² and h₂ = 3 W K⁻¹ m⁻², respectively.

Table A.1 shows an overview of some material properties applicable to telescopes. We do not provide an emissivity value for S355 steel since it depends strongly on the surface preparation and regular steel is normally coated anyway to protect from rust. The aluminum alloy 5754 is a candidate for the ELT outer dome cladding ($\overline{ϵ}=6$% corresponds to a clean and 20% to a contaminated surface, respectively).

Zerodur is a low-thermal-expansion glass ceramic used for various telescope mirrors. The thermal conductivity of Zerodur and other glasses and ceramics is often so low that the relation k_m/w ≫ h no longer holds for mirrors with typical thickness values of w = 50 mm (ELT M1 segments) to hundreds of millimeters. As a consequence, transient thermal phenomena become three-dimensional problems (Banyal & Ravindra 2011) with the relaxation time along z of τ_z = w² c_m ρ_m/k_m.

All Tables

All Figures

	Fig. 1. Sky radiance Lsky(λ) for Cerro Paranal at different PWV conditions typical for observations in the direction of zenith. The astronomical observation bands and the main ozone absorption are marked. For further details on species that are not important for us here (e.g., CO2, N2O and CH4) see Fig. 1 and Table 1 in Smette et al. (2015).
In the text

	Fig. 2. Correlation of the sky temperature Tsky and PWV at zenith Only the range of conditions with PWV <  15 mm relevant for typical astronomical observations is used for the final fit function.
In the text

	Fig. 3. Correlation of the sky temperature Tsky and zenith distance at a constant level of PWV.
In the text

	Fig. 4. Sky radiance at zenith on the coldest and the warmest days at Cerro Paranal (2600m) in July 2018. The fit of a black body is based on the marked points only, defining an upper boundary.
In the text

	Fig. 5. Same as Fig. 4, but for Mauna Kea, Hawaii (4200m)
In the text

	Fig. 6. Comparison of the sky temperature $T_{\mathrm{sky}}^{\mathrm{full}}$ derived from the sky model for the full wavelength range, the temperature Tsky derived from sky models assuming the approximation of pure black body radiation above 30 μm, and the estimated $T_{\mathrm{sky}}^{\mathrm{est}}$ with the ambient temperature black body radiation and the atmospheric transmission curve only.
In the text

	Fig. 7. Two-dimensional histogram of the ambient temperature Tamb as a function of PWV for clear observing conditions (see text for more details). The bin widths of the histogram are 0.1 mm in PWV and 0.5 K in Tamb. The resulting Pearson coefficient r is indicated (see text for more details).
In the text

	Fig. 8. Histogram of the calculated Tsky values; its mean value is 246.2 K, the standard deviation σ = 5.7 K.
In the text

	Fig. 9. Tsky vs. TIR at θ = 0° as two-dimensional histogram (bin width of 0.5 K) for clear observing conditions.
In the text

Fig. 10. Evaluation of sky view factors by projection. Red arrow: normal vector n above the point q on a surface. Blue quadrilaterals: two windows of solid angle Ω ≈ 0.16 sr through which the sky can be seen from q (i.e., 𝒯 = 1); the right window is tilted against n by γ = 60°. Green rectangles: projection of the windows onto the tangential plane with areas π Fsky ≈ Ω and Ω × cos(60° ) = Ω/2, respectively.

In the text

	Fig. 11. Convective heat transfer coefficient hcyl vs. local air speed v for three different interface lengths s representative of structures in the VLT and ELT telescopes, plus the ELT dome door.
In the text

	Fig. 12. Steady-state surface subcooling ΔTsurf as a function of $\overline{ϵ}{F}_{\mathrm{sky}}$ for four different local air speeds v. The upper edges of the colored fans correspond to s = 1.26 m, the dashed lines to s = 1.88 m, and the lower edges to s = 4.3 m, respectively.
In the text

	Fig. 13. Steady-state surface subcooling ΔTsurf as a function of Fsky and v for s = 4.3 m and $\overline{ϵ}=0.2$. Contours: increments of −1 K; thick contour: commonly accepted surface subcooling limit in telescopes of ΔTsurf = −1.5 K.
In the text

	Fig. 14. Slice through a cylindrical truss of diameter d = s/π running perpendicular to v and to the telescope axis. Red arrow: ray of light traversing the cool wake of thickness δ and temperature difference ΔTwake.
In the text

	Fig. 15. Steady-state OPD from subcooling as a function of Fsky and v for a single ray pass through the wake of a single truss of circumference s = 4.3 m and $\overline{ϵ}=0.2$ (with ELT median conditions). Contours: OPD in nanometers in increments of factors of 101/5 ≈ 1.585, starting at 10 nm.
In the text

	Fig. A.1. Example of steady-state temperature distribution u(x) with heat conduction. Blue and red curves: plate of size π 0.6 m with step-like radiative forcing; ξ = 0.27 m and 0.39 m, respectively. Turquoise and orange: same ξ values, but sinusoidal forcing around pipe-shaped truss of diameter d = s/π = 0.6 m. The transition region of the blue curve has the width Δx = 0.55 m.
In the text