© ESO, 2012

1. Introduction

One of the main characteristics of astronomical objects is their brightness in different spectral bands. The standard precision of ground-based photometry is adequate in most cases, but a number of astronomical problems require an even greater precision (Heasley et al. 1996; Everett & Howell 2001). One of the fundamental factors limiting the precision of ground-based photometry is the stellar scintillation occurring in the atmosphere as a result of its turbulent nature.

The fluctuations of the refractive index cause phase distortion in a plane light wave passing through the atmosphere to an entrance aperture of the telescope. As the wave propagates, the phase distortions lead to a redistribution of the amplitudes between different parts of the wavefront. Averaging within the aperture reduces the fluctuations caused by this mechanism, but does not eliminate them completely. In photometric practice, this effect is considered as an additional source of error. The scintillation noise expressed in stellar magnitudes does not depend on the object’s brightness, therefore it cannot be reduced by observing brighter stars.

Scintillation noise in high-precision and fast photometry has been studied for quite a long time (Young 1967, 1969; Dravins et al. 1997). The many methods proposed for reducing scintillation noise (Heasley et al. 1996; Dravins et al. 1998; Gilliland et al. 1993; Osborn et al. 2010) show that there is no perfect solution and that additional work is needed. Moreover, this same problem is often presented as an argument to do precise photometry from space.

Stellar scintillation is of interest not only to photometry. This phenomenon is a powerful tool for remote sensing of optical turbulence (OT) in the atmosphere. The mechanisms of scintillation generation and its characteristics are well studied theoretically and experimentally because they are closely related to the most important characteristic of the OT above astronomical sites, the seeing.

Unlike the seeing, scintillation received relatively little attention in the astro-climatic work; the error budget of high-precision photometry is still evaluated using the data from Young (1967) or other disparate estimates. Only recently such studies have been conducted in the general context of the characterization of optical turbulence above different astronomical observatories and prospective sites (Kenyon et al. 2006; Kornilov 2011). These data are useful for comparing sites in the context of photometry. Moreover, the scintillation noise of long-exposure photometric measurements depends not only on the intensity of high-altitude OT, but on the wind speed at altitudes above the tropopause (Kornilov 2011), which is important for the global dynamics of the atmosphere.

This paper presents scintillation noise measurements with the MASS instrument at observatories situated in different geographical areas. Section 2.1 recalls the theoretical description of the scintillation in two basic measurement regimes. In Sect. 2.2 the method used to estimate relevant parameters from the raw data is described. The next section describes the sites studied here, the original data, and the procedure for calculating the scintillation noise. Results and comparative analysis are presented in Sect. 4, Sect. 5 describes the temporal variability of scintillation noise, and the final section is a discussion and comparison with other available data.

2. Basic relations

2.1. Theoretical background

Scintillation is characterized by the variance of relative fluctuations of the flux I passing through a receiving aperture, the so-called scintillation index s²: (1)In the approximation of weak perturbations, the total scintillation index is the sum of scintillation produced by independent turbulent layers: (2)where W(z) is a weighting function that depends on the size and shape of the receiving aperture and does not depend on the distribution of the structural refractive-index coefficient . The weighting function has a simple physical meaning, it equals the scintillation index generated by a layer of unit intensity  located at a distance z.

In the general case of non-zero exposure, the weighting function also depends on the exposure time and on the wind speed at altitude h. The altitude and distance to the layer are trivially related as h = z   cosγ, where γ is the zenith angle of the observed object. Theoretical descriptions of the scintillation can be found, for example, in Young (1967); Roddier (1981); Dravins et al. (1997).

Using Taylor’s frozen-flow hypothesis (Taylor 1938), the wavefront evolution during exposure time τ for wind speed w(z) is reduced to a simple translation by w(z)τ. Invoking the Kolmogorov spectrum of refractive index perturbations, the final expression for the weighting function is (Tokovinin 2002; Kornilov 2011) (3)Here, the integration is performed over the modulus of the spatial frequency f, assuming that the functions in the integrand are either axisymmetric or are already averaged over the polar angle. The function A(f) is the aperture filter, which takes into account spatial averaging by the aperture. A_s(w,τ,f) is the spectral filter of the wind translation, which describes the temporal averaging, and S(z,f) is the Fresnel spectral filter, which describes the generation of amplitude distortions in the propagation of the wavefront.

For a circular aperture, the aperture D filter is A(f) =  [ 2   J₁(πDf)/πDf] ². The Fresnel filter S(z,f) depends on the wavelength λ of the detected radiation and can be calculated for any spectral energy distribution (Tokovinin 2003). For monochromatic radiation S(z,f) =  [ sin(πλzf²)/λ] ². The characteristic spatial scale of the filter is the Fresnel radius r_F = (λz)^1/2.

The wind-translation filter is not axisymmetric, but after averaging over the polar angle it is represented as , where the function  can be expressed in terms of special functions (Tokovinin 2002; Kornilov 2011). The important features of this function are its asymptotes for small ξ: and for large ξ: .

It is difficult to analyse the expression (3) in general form because it depends on many parameters. However, for certain relationships between the parameters the expression is considerably simplified. First of all, it is the definition of the small D ≪ r_F and large D ≫ r_F aperture regimes. In the first case we can set A(f) ≡ 1, in the second case we simplify the Fresnel filter by replacing the sine with its argument (Roddier 1981).

We can also identify two limiting cases by the value of the wind translation. At very short (zero) exposures wτ ≪ min(D,r_F), the filter A_s(wτ,f) ≡ 1. In the opposite case of long exposure wτ ≫ max(D,r_F), the A_s(wτ,f) = 1/πwτf.

These limiting situations are generally accepted, although in photometric practice, the small aperture approximation is not as interesting. In this paper, we will consider the evaluation of scintillation noise on a large telescope in the regimes of short (SE) and long (LE) exposures.

The scintillation index in a large telescope for short exposures is expressed by the well-known formula (see, e.g., Roddier 1981) (4)since in this case the weighting function is W(z) = 17.34   D ^− 7/3   z². The weighting function does not depend on the wavelength owing to the achromatism of the scintillation in a large telescope. This case is typical for fast photometry. The dependence on the zenith distance γ can be obtained from trivial geometric considerations: , where M_z = secγ. This formula ignores the central obscuration of astronomical telescopes, important in the short-exposure regime (Young 1967; Dravins et al. 1998).

Nevertheless, in conventional astronomical photometry typical exposures of tens of seconds and longer are used. In this case, the condition wτ ≫ D for long exposure is satisfied for typical wind speeds in the upper atmosphere, and (Kornilov 2011) (5)Naturally, the wavefront is translated by the wind’s component perpendicular to the line of sight. Given that, on average, the wind vector is directed horizontally, this effect also depends on the air mass M_z and on the azimuth of the wind with respect to the observed star (Young 1969). The dependence of the scintillation on the air mass therefore varies from  for transverse wind up to  in the case of longitudinal wind directed along the object azimuth. In the long-exposure regime the scintillation is little affected by the central obscuration (Young 1967), which has to be taken into account only in the extreme case of a ring-like aperture (Kornilov 2012).

Formulae (4) and (5) show that for a particular telescope and exposure time the scintillation index can be easily calculated if the integrals in these formulae are known. Using the notation introduced in Kenyon et al. (2006), (6)for the SE regime, and (7)for LE regime.

These scintillation noise parameters and  characterize the power of scintillation noise under given atmospheric conditions and represent the scintillation index in a 1-m telescope (with 1 s exposure for the LE regime). The quantities ℳ₂ and  are known as atmospheric moments of second degree (Tokovinin et al. 2003b).

2.2. Measurement of scintillation noise parameters with MASS

The Multi-Aperture Scintillation Sensor (MASS) is an instrument for measuring OT profile from scintillation (Kornilov et al. 2003, 2007). Photons from a single bright star collected by a small telescope are detected with four photo-multipliers that sample intensity in the pupil with four concentric annular apertures A, B, C, and D. The smallest circular aperture A at the centre has a typical diameter of 2 cm, the largest annular aperture D has an outer diameter of  ~10 cm. These diameters are determined by the optical magnification factor k, adjusted and measured for each instrument individually. Photon counts with a 1-ms exposure are processed statistically to derive four normal and six differential scintillation indices, which are related to the OT through weighting functions (Eq. (2)). Apart from the aperture geometry, these functions depend on the spectral energy distribution, i.e. on the spectral type of the star and on the instrument response. The functions W(h) are computed by numerical integration of the basic expression (3). In the calculations, it is assumed that the measurement occurs with “zero” exposure, neglecting the wind translation during 1 ms (Tokovinin et al. 2003b; Kornilov & Kornilov 2011).

The measured scintillation indices are used not only for the restoration of the vertical profile of OT, but also for the evaluation of several atmospheric moments required for computing such integral characteristics of the atmosphere as the free atmosphere seeing, isoplanatic angle, effective altitude of the turbulence, etc. (Tokovinin et al. 2003b; Kornilov et al. 2007). Among these moments, the second moment ℳ₂ is computed, so the evaluation of the parameter  is a trivial task.

The method of finding the parameter  has been proposed in Kornilov (2011). It relies on the fact that the scintillation index  for long exposure can be calculated from the mean fluxes measured with 1 s exposure that are stored in the MASS output files *.stm. It should be mentioned that for the small MASS apertures the LE regime begins at τ ≪ 1 s.

In this LE regime, the scintillation index (we call it LE-index to distinguish from the “fast” indices with 1 ms exposure) can be expressed as (8)where U′(z) are the MASS weighting functions for the LE regime. They are distinguished from the conventional weighting functions by the additional factor (πf)^-1 under the integral over the modulus of the spatial frequency. The typical behaviour of U′(z) for the whole set of MASS indices is presented in Kornilov (2011), so we only note that the asymptotic dependencies of these functions are U′(z) ∝ z^4/3 for infinitely small apertures, and U′(z) ∝ z² for large apertures.

We estimate  from the weighted sum (linear combination) of LE scintillation indices measured with MASS. The coefficients of this sum are adjusted to make the weighted sum of U′(z) as close as possible to 10.66   z². Then, as can be seen by comparing Eqs. (7) and (8), the weighted sum will approximate . Details of the procedure are given in Appendix A.

Figure 1 plots the relative error of these approximations for five typical MASS devices. It shows that for propagation distances between 4 and 32 km, the errors do not exceed  ± 5%. The contribution of closest layers (below 1 km) is underestimated, but it is negligible because of the factor z². The impact of the distant (beyond 32 km) OT is underestimated by  ≈ 10% as well. Here we exclude the measurements made far from the zenith (M_z > 1.3), so this under-estimation is relevant for altitudes above 25 km where the OT already has a low intensity.

This figure also shows the curve for the aperture D of the device with k = 16.3. It is evident that the result is greatly underestimated if we use this aperture only. The approximation is therefore better for devices with larger apertures (higher magnification coefficient) and a higher spectral response in the blue.

The influence of the magnification coefficient on the results was verified a posteriori for the real data obtained with the device MD05 TMT (Armazones). When changing k from 15.9 to 15.4 (by  ≈ 3%), the median value of the parameter  decreased by 4%, so that the requirement for the accuracy of k is not stronger than in the typical processing of MASS data (Kornilov et al. 2007). The substitution of another spectral reaction curve with an effective wavelength shorter by  ≈ 14 nm changed the median by less than 0.2%. The dependence of the approximation on the star color is negligible.

3. Data processing

3.1. Data sets and site testing campaigns

In the past decade, a large amount of scintillation data from MASS instruments has been accumulated. The measurements were performed during various projects in different geographical locations. Some of these campaigns focused on the site selection for future large telescopes. Other studies were methodological in nature, or supported already operating telescopes. The list of astronomical sites studied in this paper is listed in Table 1.

They include all sites studied by the TMT site testing programme (Schöck et al. 2009), representing two sites in the northern hemisphere: Mauna Kea (Hawaii) and San Pedro Martir (Mexico) and three mountains in northern Chile: Cerro Tolar, Armazones and Tolonchar. The data from this campaign use a similar setup on all sites using MASS/DIMM units attached to custom-made telescopes capable of operating robotically within a wide range of wind conditions.

Measurements on Cerro Armazones in Chile were first carried out by TMT. After this site was chosen for the E-ELT, monitoring was taken over in 2010 by ESO. Shorter samples at the nearby sites Ventarrones and La Chira, studied in the frame of the E-ELT site-characterization campaign (Vernin et al. 2011), have been added for the purpose of a temporal correlation analysis.

The robotic MASS-DIMM site monitor at Cerro Pachón is a shared facility between the Gemini-S and SOAR telescopes, operational since 2004. Early results of turbulence monitoring at Pachón were reported in Tokovinin & Travouillon (2006). Here we use the much larger data set accumulated during eight years.

The longest series of observations were obtained at the Paranal and Pachón observatories. To control the stability of the results and for comparison with the two campaigns at Armazones, the Paranal data are also divided into two parts: 2004 − 2009 (Paranal A) and 2010 − 2011 (Paranal B). Because the measurements are carried out at Paranal in continuous mode, both samples are quite large. The time coverage of all campaigns is shown in Fig. 2. The overall volume of data amounts to more than 40 000 h.

The only site in the eastern hemisphere is represented by data from the Mt. Shatdzhatmaz in the northern Caucasus (Kornilov et al. 2010) where the Sternberg Astronomical Institute (SAI) observatory is building its 2.5-m telescope. Note that in operating observatories, the measurements with MASS device are continuing to provide operational information about OT in the atmosphere.

3.2. Computation of LE scintillation indices

For each of four MASS apertures, the LE scintillation index  is calculated as the mean square of the difference between adjacent 1-s mean fluxes, (9)where is the average of the series of N photon counts  { F_i }  with 1 s exposure. In the MASS output files, these counts are scaled to 1 ms exposure, which causes the factor of 10^-3 in the last term describing the contribution of the photon noise. Use of the difference effectively suppresses the contribution of the flux variations on time intervals greater than 1 s, unrelated to the scintillation. It is easy to show that the scintillation is uncorrelated at long exposures; it follows from the relation s²(τ) ∝ 1/τ in the LE regime.

Each estimate of the  uses all 1-s flux values during the MASS accumulation time. As a rule, this time is 1 min, hence N = 60. For an unbiased estimate of the variance of the difference, a factor of 1/(N − 1) is required; although the number of differences is also N − 1, the mathematical expectation is known and equal to zero.

For fluxes registered with 1-s exposure, the photon counting statistics can be considered as exactly Poissonian. Any difference of the PMTs from an ideal detector does not produce errors exceeding 1%, even for extremely small indices. Non-linearity of the photon counts is corrected with a fixed dead time of 20 ns, as typical for the MASS detectors.

In contrast to Kornilov (2011), the reduction to the zenith is performed after evaluation of the parameter  from the LE indices measured at a certain air mass M_z, by Eq. (A.1). Since we do not know the direction of high-altitude winds, two estimates of the scintillation parameter at zenith are calculated: the minimum  for the longitudinal wind, and the maximum  for the transverse wind (see Sect. 2.1). Data obtained at air masses M_z > 1.3 are not used to avoid the uncertainty in the reduction to the zenith caused by the unknown wind direction. In addition, as noted in Sect. 2.2, with a large air mass our method underestimates the contribution of high-altitude turbulence to the parameters  and .

Filtering used to remove invalid data and the estimates of the resulting errors are considered in Appendices B and C. The filtering procedure eliminates a small, though non-negligible, fraction of the data, on average about 10%. In some campaigns, a significant part of measurements was carried out far enough away from the zenith; in these cases, the proportion of the rejected measurements exceeds 20%.

Although we took great care in filtering out the bad data, the final results remain almost the same when all the data are used. This is so because the data were collected by automatic monitors for a long time and therefore contain only very few faulty measurements.

4. Results

4.1. The scintillation parameter S

As noted in Sect. 2.2, the calculation of the parameter  is made by rescaling the atmospheric moment ℳ₂ obtained in conventional MASS data processing and already reduced to the zenith (see Eq. (6)). Filtering of the original data consists only in the choice of the points corresponding to the correct estimates of the parameter . For some measurements, for various reasons, the ℳ₂ have not been obtained while the corresponding parameter  has been successfully evaluated. These situations are reflected by the fact that in Table 2 the total time T is somewhat shorter than in Table 3.

Note that the algorithm of calculating ℳ₂ implemented in the early versions of the MASS software had some instability, which caused artefacts in the distribution and a systematic overestimation of the second moment. Therefore, in almost all cases we used the results of data reprocessing carried out later with the latest version of the software atmos-2.93 (Kornilov & Kornilov 2011). For measurements on Tololo and Pachón we used the ℳ₂ calculated from the restored vertical profiles of the OT.

The cumulative distributions of the index  (in m^7/3) are shown in Fig. 3, and their characteristic points are given in Table 2. Of all the distributions, the Mauna Kea site stands out, as expected because of its high altitude and low turbulence in the upper atmosphere. The rest of the curves differ only in details, especially if we take into account the sampling effects for those campaigns with non-uniform seasonal coverage (we discuss this in the following section) and some differences in the MASS devices. That the Tololo and Pachón sites show a stronger scintillation noise is likely a result of evaluating  at these sites from the OT profiles. Nevertheless, this systematic difference does not exceed 10%.

The comparison of the Paranal A and Paranal B sub-samples shows that the period 2010–2011 does not differ in scintillation power from the previous years. The same should be expected for the measurements at the Armazones summit due to its geographical proximity to the Paranal (30 km). However, the parameter  for Armazones B is higher by 15% than for Armazones A. Apparently, this discrepancy is caused by a significant difference in the geometry of the entrance apertures of the instruments used in both cases. A similar systematic difference is observed for the parameter of , see Table 3.

The differential distributions are clearly asymmetric. For all studied sites, the modes of the distributions coincide with the 23 − 26% quantiles, so one can use the first quartile from Table 2 as the most probable value of . The typical skewness γ₁ amounts to 2 − 3 and the excess kurtosis is γ₂ ~ 10. Here we do not give more precise values, because adequate calculation of the higher moments of a statistical distribution requires careful removal of outliers and long-term trends.

The distributions of  are close to the normal distributions with slightly different means and widths. It is also noticeable that the asymmetry of the distributions are somewhat different for different observatories. Typical skewness and excess kurtosis of  are small, γ₁ ~ 0.2−0.4 and γ₂ ~ 0−0.2.

4.2. The S scintillation parameter

The relative errors of estimates are quite large,  ~0.3−0.4, widening the distributions somewhat. The significance of this effect was estimated by comparing the distributions of the initial results obtained with 60 s and with 240 s accumulation time (the latter are averages of four consecutive 60-s measures). As expected, the distribution of 240-s estimates is somewhat narrower than that of the 60-s estimates, but its median is higher by about 5% because of the asymmetry.

In Fig. 4 the cumulative distributions of the 240 s estimates of  are shown for seven sites. For clarity, only the curves corresponding to the assumption of the longitudinal wind (i.e. the lower limit of ) are plotted. In Table 3 the characteristic points of the distributions are listed for both wind directions; they can be considered as lower and upper limits of the  parameter.

Both the curves in Fig. 4 and the data in Table 3 show that the distributions of the  at various observatories differ even less than the distributions of . The curve for Mauna Kea again shows a lower scintillation. Other observatories, except the SPM, have approximately the same characteristics. The scintillation parameters measured at SPM are greater than the typical values by about 20%. However, we emphasize again that these results reflect each particular set of data without taking into account the seasonal distribution of the observations. It will be demonstrated in Sect. 5.2 that the  parameter is subject to seasonal variations.

The differential distributions of the  parameter are very similar to the differential distributions for the ; they have roughly the same asymmetry and a somewhat higher excess. Just as for , the distributions of  are close to normal; typical values of their skewness and excess kurtosis are γ₁ ~ 0.2−0.4 and γ₂ ~ 0−0.5.

5. Temporal variations of the S₂ and S₃ parameters

5.1. Nightly variations

The study of temporal variations of the astroclimatic parameters on time scales of  ~1 h (in course of a night) is important in terms of short-term predictions for operational scheduling of observations (Racine 1996; Skidmore et al. 2009). Of course, the variability of the parameters S₂ and S₃ is not as critical as the variability of the seeing. However, for a long time-series of photometric measurements, control of the scintillation noise is also highly desirable.

To quantify the nature of the short-term variability, different approaches are used. These are the calculation of the auto-correlation function (ACF) in the usual sense (Tokovinin et al. 2003a) and the average relative or absolute differences of a parameter for a certain time delay. The fractional difference (FD) proposed in Racine (1996) is not very amenable to a rigorous mathematical analysis (being a ratio, it has a non-trivial distribution function). The absolute difference (AD) used in Skidmore et al. (2009); García-Lorenzo et al. (2010) is restrictedly suitable for non-stationary processes where the variance depends strongly on the average.

To achieve the objective outlined by Racine (1996) in his introduction of the FD, one could analyse logarithms of the parameters instead of the parameters themselves. This method is recommended by Jenkins & Watts (1969) for non-stationary processes of multiplicative nature, where the dispersion is nearly a linear function of the mean. We have already used such a transformation when evaluating random errors of  and . In Sect. 4 we pointed out that  and  are distributed almost normally.

Thus, to describe the temporal behaviour of the scintillation parameters we use the structure functions of their logarithms. For the  parameter this is (10)where the averaging is performed over all available pairs for the given delay (lag) T. Similarly, we define the function  for the parameter . To keep the analysis free from normalization, we consider non-normalized functions and define the characteristic correlation time of the scintillation power using an absolute criterion.

The structure functions for seven sites are presented in Figs. 5 and 6, where the contribution of uncorrelated noise was subtracted. For this calculation we used 1-min measurements of both  and . Note that all structure functions have a similar behaviour for delays of less than 2 h. At longer delays, the time coverage of the individual data sets influences the computed structure functions and their behaviour becomes irregular. The initial part of the functions  and  is well approximated by the formula (11)The main features of this dependence are infinite derivative at zero and saturation at the level of σ² when T → ∞. The point T = τ_c corresponds to . The time constant τ_c therefore characterizes relative changes in the power of the scintillation. The structure function (11) corresponds to the temporal spectral density proportional to 1/ω (flicker noise), which saturates at low frequencies ω ≪ 1/πτ_c.

The approximation of the measured structure functions by Eq. (11) on the initial segment T < 2 h shows that the parameters σ² and τ_c are similar at all sites, with typical values τ_c ~ 1 h. The values of τ_c are given in Tables 4 and 5. The parameter σ² lies in the range of 0.25−0.4. The small additional parameter δ ~ 0.05 compensates for the imperfect subtraction of uncorrelated noise; it is not considered here.

For practical purpose, another way to evaluate the stability of the scintillation noise may be more convenient. To do this, we transform the difference of the logarithms in (10) to a logarithm of ratio and define the statistics (12)and the similar function R₃(T) for the  parameter. The functions R₃(T) are plotted in Fig. 7. The time lags τ_1.5 corresponding to changes of  and  by 1.5 times are given in Tables 4 and 5.

At almost all sites, the  parameter is changing faster than the . Apparently, in this case the variations of the wind speed are important. The fastest changes of the scintillation power are observed at Mauna Kea: 0.2 and 0.3 h. The atmosphere above Mt. Shatdzhatmaz is changing most slowly: 0.7 and 0.8 h. However, note that most measurements at this site have been performed from October to December, which is the most stable season.

Variation of the scintillation power over longer time scales can be estimated by comparing its variability within a night to the variability during the measurement season. To do this for both scintillation parameters, we calculated the average  ⟨ lnS² ⟩ _j and the corresponding variances  within each night j. Then, the variance of the average nightly quantities  was evaluated. The intra-night variability is quantified by its mean over the whole campaign, . Nights containing less than ten measurements of the  and  were discarded in this analysis.

The , corrected for the contribution of random measurement errors, are presented in Tables 4 and 5. On average, the variance of the random noise is 0.01 for the  parameter, and about 0.1 for the  (see Tables 2 and 3). Evidently, the relative variability within a night is almost identical at all observatories. As expected, the variability of the  is slightly higher than the variability of the . There is one statistically insignificant exception for the measurements carried out at Tololo. The average value of the variability within nights for  is 0.13    ±    0.01, for  is 0.16    ±    0.03.

The variability of the nightly averages is a combination of regular seasonal changes (see the next section) and the sporadic variability from night to night. We did not isolate the seasonal behaviour, so the full variances  are listed in Tables 4 and 5. Not surprisingly, the maximum values are observed for Shatdzhatmaz and S. Pedro Martir, because the seasonal variability at these sites is as high as 0.1 − 0.14 (about 2.5 times higher than the intra-night variability).

Measurements at Tolonchar also show a significant long-term variability. The data from this site are a poor fit to Eq. (11) and are therefore not represented in Fig. 5. However, it is possible to compare the temporal variability of Tolonchar to the other sites using the τ_1.5 metric as shown in Tables 4 and 5. Many parameters describing the scintillation noise at this site are substantially different from the other sites. The median variability from night to night at Tolonchar is 0.31    ±    0.07 for  and 0.24    ±    0.08 for  (less than for , what is atypical).

5.2. Seasonal variability

For all observatories, the long-term evolution of the parameters  and  is seasonal. Seasonal variability is evident in the raw data, but for a more detailed and quantitative study, we calculated the statistical characteristics by months of the year. These medians and quartiles (in the form of bars) for the parameter  are shown in Fig. 8 and for  in Fig. 9.

Note that almost all samples on which these statistics were determined are sufficiently large and their volumes vary by no more than two times from month to month. The exception are the Mt. Shatdzhatmaz data, where the vast majority of clear weather happens during from October to December. The seasonal dependence of  for Armazones B (not plotted) is similar to that of Armazones A, with slightly higher values.

In these plots, the noticeable difference between the northern and southern hemispheres is clearly visible, especially for . There is also some dependence of the seasonal variability on the latitude of the observatory. The lowest seasonal variability is found at Mauna Kea (for both  and ). The most pronounced seasonal variability of  can be seen on the curve for Shatdzhatmaz from May to August; the scintillation noise increased by almost two times, although this has little effect on the overall median (see Table 3). In contrast, the data for the SPM observatory contain many measurements from May to August, which overestimates the median for the entire data set.

Conversely, in the southern hemisphere the median of  is minimal in the period from May to September, and the median  is maximal from July to September, while in the northern hemisphere there is a minimum in this season. However, in general the seasonal variability does not exceed the amplitude of distributions defined by the quartiles. This means that at any time of the year there is a reasonable probability of both good and poor conditions for high-precision photometry.

Seasonal variations in the scintillation noise are caused by the redistribution of seasonal winds in the upper atmosphere. Moreover, an increase in the wind speed leads to an increase in the  parameter, but also to a decrease of the . Data on the wind speed above the tropopause may also be obtained from our measurements.

5.3. Correlation between close observatories

Strong correlation of scintillation power on time scales of the order of 0.5−1 h (Sect. 5.1) implies that the size of the spatial region in which the power of turbulence, wind, and their distribution over altitude can be considered constant, is  ~100 km. In this case, a strong correlation of the scintillation power should be observed for closely spaced observatories (<50 km).

Some observatories on our list are geographically located very close to each other. The distance between Tololo and Pachón is about 10 km. The distance between Paranal and Armazones is no more than 25 km. Additionally, in an area around Paranal, the ESO team investigated two sites: La Chira (15 km) and Ventarrones (32 km). Two other summits in northern Chile were studied by the TMT team: Tolar and Tolonchar are spaced from Paranal by about 300 km, first to the north and the other to the west.

The cross-correlation of the scintillation parameter  (or ) measured at two observatories, was calculated as follows. For the required time delay T, the pairs of measurements for which the difference between their acquisition time is less than T/2 (30 s) were selected from the two data sets. Then, using the classical method, the correlation coefficient ρ(T) of the first sample on the second for the logarithms of the scintillation parameters lnS² was calculated. Finally, we corrected the correlation coefficient for the contribution of random measurement errors by its renormalization with corrected variances.

Naturally, the correlation reflects not only the short-term variations of the scintillation power, but also its variability on the time scale of days, including seasonal variability. This can be clearly seen in Fig. 10 which shows the cross-correlation between Tololo and Pachón with time lags  −5 < T <  + 5 h. In addition to the correlation peak of  ~1 h width caused by the intra-night variability, there is a wide pedestal at the level of ρ(T) ≈ 0.5. The cross-correlation between measurements at Tololo and Paranal, also depicted in the figure, shows no correlation peak, while its constant level of ρ(T) ≈ 0.2 reflects the common nature of seasonal variability at observatories in Chile. In contrast, a typical cross-correlation between Paranal and Shatdzhatmaz is  ≈−0.2 (the seasonal variability is anti-correlated).

The width of the correlation peak is comparable to the width of the peak of the auto-correlation function (ACF) for each individual site. For illustration, the ACF for Tololo is plotted in Fig. 10. A similar situation is observed for other close observatories, see Fig. 11, except for the weak cross-correlation peak between Paranal and Ventarrones. The maximum cross-correlation is larger for the parameter  than for , but not by a significant amount. The height of the peak above the extended pedestal is approximately 0.2. Generally, when there is a significant correlation, the ratio of the cross-correlation peak to the pedestal agrees well with the relative variability during the night and the variability from night to night (Tables 4 and 5).

The lack of a significant delay in the variations of the scintillation power indicates that at large spatial scales (10 − 100 km), the wind transportation of the OT is not dominant compared to the intrinsic evolution of the turbulence.

6. Discussion and conclusion

The resulting characteristics of the scintillation noise in the regime of short exposures are shown in Table 6, for long exposures – in Table 7. Instead of the squared parameters, those tables list S₂ and S₃. The S₃ was calculated as the average between the values obtained under the assumptions of longitudinal and transverse wind. For the Armazones site, the average values computed with a weight proportional to the number of measurements for campaigns Armazones A and Armazones B are presented.

We recall that the scintillation noise at the zenith σ_S and σ_L for a telescope with diameter D is calculated by the formulae (13)for short exposures, and (14)for a long exposure τ.

These expressions follow directly from the Eqs. (4) and (5). To convert the values of Tables 6 and 7 into stellar magnitudes, they should be multiplied by the constant 1.086, i.e., the amplitude of the scintillation noise in magnitudes is σ_mag = 1.086σ.

As we pointed out in Sect. 2.1, for the SE regime the central obscuration effect should be considered. A good approximation is presented in Kornilov (2012). For extra-large telescopes (D ≳ 8 m), the effect of the turbulence outer scale is substantial in both the SE and LE regimes, and should also be taken into account (Kornilov 2012).

The results shown in these tables agree well with the estimates of the scintillation parameters S₂ and S₃ for Tololo and Pachón obtained in Kenyon et al. (2006) by numerical calculation of the moments ℳ₂ and  by Eqs. (6) and (7) on the basis of the measured OT vertical profiles and the modelled wind profile. Our method gives a similar estimate of S₃ = 0.0030 for the Maidanak observatory (Kornilov 2011). Unfortunately, to our knowledge other statistically reliable estimates of the scintillation noise are not available.

The results indicate that scintillation noise is mostly defined by the global turbulence at altitudes of 10–15 km, at the tropopause and above. The only site in our list where the local effects are noticeable and make a difference with other sites is the Tolonchar. In this case, the proximity to the main ridge of the Andes, which are quite high and perpendicular to the global air circulation, could create a quasi-stationary vertical vortex with a scale of tens of kilometres.

The stronger scintillation noise (about 10% excess) at the SPM observatory is likely a result of a biased seasonal distribution of the observations together with a significant variability from night to night (see Table 5). On the other hand, the most probable value of S₃ at this observatory exceeds the value at Paranal by only 3%.

6.1. Comparison of the S₃ parameter with Young’s equation

Until now, astronomers (see, e.g., Everett & Howell 2001; Mann et al. 2011) used Young’s formula (Young 1967) to estimate the contribution of scintillation noise to the accuracy of photometric measurements, (15)where D is telescope diameter in meters, and h_obs is the observatory altitude above sea level. Since the original expression contains a bandwidth rather than exposure, there was a misunderstanding in the translation of one to the other. Equation (15) is taken from Gilliland et al. (1993), where it was corrected after the intervention of Young (1993). The dependence on the observatory altitude was originally proposed by Reiger (1963), who assumed an exponential dependence  from altitude, but the scale height of h₀ = 8 km was established by Young. The numerical coefficient was determined from observations mainly with the 0.9 m telescope.

For comparison with our data, we rewrite the previous formula as S₃ = 0.0030   e ^{− h_obs/h₀}. Using this formula and the site altitudes from Table 1, we obtain S₃ estimates from 0.0018 for the highest summit of Mauna Kea to 0.0023 for the Shatdzhatmaz. A comparison with Table 7 shows that Young’s formula underestimates the median amplitude of the scintillation noise by a factor of 1.5 (scintillation power by two times). The values of S₃ inferred from Young’s formula are similar to the first quartiles of its actual distributions derived both here and in Kenyon et al. (2006); Kornilov (2011).

The dependence of the scintillation on the observatory altitude is ambiguous because of many factors affecting it. Our results do not show this dependence for most of the studied sites located between 2000 and 3000 m a.s.l. because this effect does not exceed the accuracy of the method. Only for the highest site, Mauna Kea, we can possibly relate lower scintillation noise to higher altitude. On the other hand, the Mt. Ventarrones data show a domination of local effects.

6.2. Conclusions

We presented the results of the evaluation of scintillation noise in observations on telescopes with large diameter D ≫ r_F or D ≳ 1 m for the optical and near-infrared in the Earth’s atmosphere. These estimates were obtained by an indirect method based on the data of the measurement with the MASS instruments, without involving measurements on large telescopes. However, 1) this method has a reliable theoretical basis; 2) identical instruments were involved and measurements were obtained using the same technique; 3) data for each studied site were obtained over a long time period and have a large statistics.

The scintillation noise at short exposures is quantified by the parameter S₂ proposed in Kenyon et al. (2006), while the parameter S₃ characterizes observations with long exposures. Using these parameters, one can calculate the scintillation noise for a telescope of any reasonable diameter for any long exposure. For example, 8-m telescopes such as Gemini or LSST at Cerro Pachón will have scintillation noise of 50 μmag for τ = 5 min on a typical night.

The data from seven of the 11 sites studied here are of particular interest to astronomers because they describe the conditions at the existing observatories Tololo, Pachón, Paranal, S. Pedro Martir, Mauna Kea, or the observatories soon to come into operation (Armazones, Shatdzhatmaz). Scintillation noise measurements at four other summits in northern Chile are more interesting from the point of view of general turbulence behaviour in the upper atmosphere.

In addition to the general statistical characteristics of the scintillation noise parameters, their temporal variability was investigated, from short-time scales of minutes to seasonal variations. Characteristics of the short-term (during the night) changes are important for optimising observational strategy in the high-precision photometry. Usually, we can assume that the scintillation noise power is sufficiently stable for 0.5 − 1 h, but sometimes it may change significantly on timescales of 10 − 15 min.

The seasonal variations are significant and in some observatories they reach a factor of two in power. They are directly related to the latitude of the observatory. The minimum amplitude is observed for Mauna Kea, which is close to the equator, and the maximum for the mid-latitude observatory at Mt. Shatdzhatmaz.

The main conclusion from the comparison of scintillation noise at different observatories is that there are no major differences. For the purpose of choosing the best site for high-precision photometry, all sites are essentially equal. A much larger effect can be achieved by choosing the best season for observations and by using the real-time information about the power of the scintillation in operational planning of the photometric observations.

Apart from the high-precision photometry, the parameters measured in the paper can be interesting for error-budget evaluation of high-precision differential astrometry (Kenyon et al. 2006; Cameron et al. 2009).

Acknowledgments

The authors thank the Sternberg Astromomical Institute MASS group: N. Shatsky, B. Safonov, S. Potanin, M. Kornilov, who provided data from Mt. Shatdzhatmaz. The data for Cerro Tololo and Pachon were obtained from the robotic site monitor operated by the CTIO site testing team. Data from Armazones, Tolar, Tolonchar and Mauna Kea were taken in the framework of the TMT site testing program.

References

Appendix A: The weighting function approximation

The best approximation is found in the same way as for other atmospheric moments (Tokovinin et al. 2003b). We solve a linear system of equations U′d = q, where U′ is the weighting function matrix of dimension m × n, d is the vector of unknown coefficients, q is the desired function 10.66z². The number of nodes n = 50 of the distance grid  { z_i }  is substantially greater than the maximum number of indices m = 10. As always, a log-uniform distance grid is used, with higher density at low altitudes.

The system is solved by singular value decomposition, discarding singular values lower than 5    ×    10^-4. The quality of the approximation is controlled by the noise amplification factor . This system is weighted  ∝ z^-1 to ensure a best fit in the range of 5 − 25 km. Here we use only the four normal indices because it became clear that an approximation that includes the differential indices is better only at low altitudes, the least interesting zone from the standpoint of the scintillation in large telescopes.

Using the linear relationship between the weighting functions and the scintillation indices, we can finally write (A.1)where  are the LE indices (with 1-s exposure) measured in the MASS apertures A, B, C, and D.

The weighting functions U′(z) and, therefore, the coefficients d_j, depend on the aperture dimensions and spectral sensitivity of the MASS detectors. Although all MASS/DIMM devices are almost identical, they have been used with different feeding optics and somewhat differ in the spectral response of the detectors and in the magnification coefficient k for the overall optical system “MASS+telescope” that determines the aperture diameters in plane of the entrance pupil. Therefore, the coefficients d_j were computed individually for each instrument considered in this work. These values are given in Table A.1.

Appendix B: Filtering of the LE scintillation indices

During the computation of the indices using Eq. (9) and the data stored in the file *.stm, a set of the flux values  { F_i }  is omitted if 1) there is no information about the sky background or the background is too high; 2) the flux ratio between MASS apertures B and C does not satisfy the condition F_B < F_C; 3) the number of 1-s flux values in a 1-min. accumulation period N < 10; and 4) the flux in aperture D is low, F_D < 100. These criteria eliminate obviously wrong or inaccurate data. However, they do not exclude all situations with incorrect data and we developed additional data-filtering criteria:

• g₁ is the ratio of the maximum term in Eq. (9) to the total sum for aperture D. This parameter detects strong outliers;

• g₂ the ratio of the flux range in the aperture D between its 95% and 5% quantiles (i.e. ignoring the two lowest and two highest values in series of 60 points) to its rms fluctuations . This indicator is sensitive to a trend;

• g₃ is calculated similarly g₂ using the full flux range. It monitors isolated overshoots or, in their absence, the general trend.

We determined the filtering thresholds for these indicators from their empirical distributions.

As an example, the differential and complementary cumulative distributions of g₂ for the data obtained at the Shatdzhatmaz are shown in Fig. B.1. Similar distributions are found at other sites. The inflection point on the right (descending) branch of the differential curve is chosen as the threshold for all sites (at higher g₂ the distributions clearly change their character). In this example, the cut-off is set at g₂ > 6.7, which discards  ≈ 8% of the data. Visual inspection confirms that the discarded data are affected by either clouds or large tracking errors.

Using differences for calculating the LE scintillation indices does not completely suppress the low-frequency flux variation, which increases the estimated variance. Most of these variations have the form of a trend during the accumulation time. We denote the relative change in the flux due to the trend by Δ. It is easy to show that a trend adds Δ²/2N² to the true value of the LE index . If the relative error of the index estimate must not exceed p, the following condition must be satisfied: (B.1)If p = 0.01 and N = 60, the allowable trend must be . Since the LE indices in the apertures A, B and C enter into the result (the parameter ) with lower coefficients, and are themselves higher than the index in D, it is sufficient to check the condition (B.1) for the aperture D only. In addition, monitoring of the absolute value of the trend is also needed, since strong systematic changes in flux induce such as an increase in the variance that the condition (B.1) again begins to be fulfilled.

We used yet another powerful criterion based on the theoretical relation between the measured values. The behaviour of the weighting functions U′(z) restricts the ratios of the LE indices. Clearly, the ratios U′(z) for apertures A, B, and C to the weighting function for the aperture D do not equal the ratio of the indices themselves, but for any profiles of OT and wind the ratios of  must lie between the theoretical minimum and maximum.

The ratio of the LE indices  is weakly dependent on the distance, which is reflected in the form of the observed distribution of the ratio of  in apertures C and D shown in Fig. B.2. The cumulative distribution is very steep and the differential distribution is very narrow. The ratio of indices in apertures A and D typically has a broader distribution, which is explained by the behaviour of the functions U′(z). This ratio is most sensitive to the specific characteristics of the MASS device and to the minimum altitude of the turbulence, which produces noticeable scintillation. The minimum ratio of LE indices is reached when the whole turbulence is at  ~ 32 km.

Figure B.2 indicates a good match between the measured ratios and their theoretical estimates, so we boldly used the experimental distributions to determine the lower and upper limits (approximately at the level of 0.1%) for rejecting the outliers. Typical instrumental causes of such outliers are a wrong sky background value or an accidental partial vignetting of the entrance pupil. Of course, random errors in the LE indices also widen the observed distribution.

Appendix C: Error estimation of LE indices

The errors of LE indices  were estimated in two ways. First of all, assuming stationarity over one minute, normal distribution, and uncorrelated 1-s fluxes, we estimated the probable relative error of the sample variance ε = (2/N)^1/2 as 0.18 for all apertures.

The second method consists of computing the mean squared difference between the logarithms of successive estimates: . Assuming quasi-stationarity, the relative error is ε = (0.5 ⟨ δ ⟩ )^1/2. The probability distribution of δ is similar to the χ distribution with one degree of freedom. The medians of these estimated errors are thus 0.175, 0.174, 0.172 and 0.172 for apertures A, B, C and D (for the Shatdzhatmaz data).

Assuming the same relative errors of  for all apertures, ε = 0.18, and assuming uncorrelated errors in all apertures, we calculated the relative error E₁ of the parameter . Because noise amplification by the linear combination (A.1) depends on the coefficients, which differ between individual MASS instruments, the errors of  are also different. The median values of these errors are 0.28−0.36.

The errors in the different apertures are partially correlated (see Kornilov & Kornilov 2011), so the estimate E₁ of the errors was additionally checked by calculating the mean difference of the logarithms of the adjacent  values, as for the LE indices. These estimates E₂ are somewhat smaller than E₁, but they agree well.

The average errors of the 4-min  estimates obtained in the this way are listed in Table 3. Clearly, they are practically identical for all sites. The same method was used to estimate the errors of the parameter , which are typically 0.06 − 0.09 (see Table 2). These errors are smaller because the estimates are calculated from  ~1000 samples.

All Tables

All Figures