A&A 455, 385-393 (2006)
DOI: 10.1051/0004-6361:20064992

UBVRI twilight sky brightness at ESO-Paranal

F. Patat¹ - O. S. Ugolnikov² - O. V. Postylyakov³

1 - European Southern Observatory (ESO), K. Schwarzschildstr. 2, 85748, Garching b. München, Germany
2 - Space Research Institute, Russian Academy of Sciences, Profsoyuznaya ul., 84/32, Moscow, 117997 Russia
3 - A.M.Obukhov's Institute of Atmospheric Physics, Russian Academy of Sciences, Pyzhevsky per., 3, Moscow, 119017 Russia

Received 9 February 2006 / Accepted 10 April 2006

Abstract
Context. Twilight studies have proved to be important tools for analyzing the atmospheric structure with interesting consequences on the characterization of astronomical sites. Active discussions of this topic have started again recently in connection with the evaluation of Dome C, Antarctica, as a potential astronomical site and several site-testing experiments, including twilight brightness measurements, are being prepared.
Aims. The present work provides the first absolute photometric measurements of twilight sky brightness for ESO-Paranal (Chile). They are meant both as a contribution to the site monitoring and as reference values in the analysis of other sites, including Dome C.
Methods. The ${\it UBVRI}$ twilight sky brightness was estimated on more than 2000 FORS1 archival images, which include both flats and standard star observations taken in twilight and covering a Sun zenith distance range 94$^\circ$-112$^\circ$.
Results. The comparison with a low-altitude site shows that Paranal V twilight sky brightness is about 30% lower, implying that some fraction of multiple scattering has to take place at an altitude of a few km above sea level.

Key words: atmospheric effects - site testing - techniques: photometric

  
1 Introduction

The quality of an astronomical site is determined by several parameters, which may vary according to the wavelength range of interest. For the optical and near-IR domain, these include typical seeing, sky transparency, number of clear nights, humidity, night sky brightness, amount of precipitable water vapor, dust, and aerosols. While the seeing, extinction, sky brightness, and other quantities are commonly measured at most observatories, the twilight brightness is not. This is mainly because the relevant information on the typical atmospheric conditions can be derived from other measurements obtained during the night. Nevertheless, twilight observations provide an independent tool for probing the overhead atmosphere under much higher flux conditions, thus allowing more accurate results. The interested reader can find an extensive review of this topic in the classical textbook by Rozenberg (1966).

Very recently, the twilight has received particular attention due to the growing interest of the international astronomical community for what seems to be the new frontier of ground-based astronomy, i.e. Dome C - Antarctica. This site is exceptional in many respects. Besides the extremely good seeing conditions reported by Lawrence et al. (2004), several studies have shown very low amounts of precipitable water vapor that, coupled with a low sky emission, could imply that this is the best site for IR and sub-millimetric ground-based astronomy. A review of the characteristics of Dome C has been recently presented by Kenyon & Storey (2005).

One of the main concerns is related to the high latitude of this site. This, in fact, causes a significant reduction in the amount of dark time with respect to equatorial observatories, thus posing some doubts about the effective exploitation of the exceptional seeing in the optical. The possibility of opening spectral windows otherwise inaccessible from the ground is in itself a valid and sufficient scientific driver for Dome C. Nevertheless, arguments in favor of Dome C as a site for optical astronomy have been put forward. Among these, a smaller average contribution by the scattered moonlight to the global background and a cleaner atmosphere have been advocated as features that may possibly compensate for the reduced dark time (Kenyon & Storey 2005). In particular, since the last phases of twilight (also known as deep twilight) are dominated by multiple scattering (Rozenberg 1966), the amount of scattered sunlight is strongly dependent on the amount of aerosol in the lower atmospheric layers (see for example Ougolnikov et al. 2004). In a supposedly low aerosol content site like Dome C, this effect is expected to be very low, and it would, in turn, allow one to start the observations earlier than at normal sites. Even though the argument has a good physical ground, direct on-site measurements are still lacking. In this respect, it is worth mentioning that a couple of dedicated experiments for sky brightness measurements are currently being setup (A. Moore, J. Storey 2006, private communications).

In spite of the large number of investigations done in the past in this field, absolute twilight brightness measurements are rather rare, especially for large observatories placed in top-rated sites. To our knowledge, the only published work on twilight observations in the Johnson-Cousins standard system is the one by Tyson & Gal (1993) who, given their purposes, report only uncalibrated data for CTIO. In light of these facts, both with the purpose of also characterizing Paranal from this new point of view and of providing the community with absolute reference values obtained over a large time baseline, we present here ${\it UBVRI}$ twilight sky brightness measurements for the first time.

The paper is organized as follows. In Sect. 2 we introduce the basic concepts through a simplified model (which is discussed in more detail in Appendix A), while in Sect. 3 we describe the observations, data reduction, and calibration. The ${\it UBVRI}$ twilight sky brightness measured at ESO-Paranal is presented and discussed in Sect. 4, while Sect. 5 summarizes the results obtained in this work.

  
2 The twilight problem

The calculation of scattered flux during twilight is a rather complicated problem that requires a detailed treatment of multiple scattering (see for example Blättner et al. 1974; Wu & Lu 1988) and an accurate description of the atmospheric composition and the physical phenomena taking place in the various layers (Divari & Plotnikova 1966; Rozenberg 1966). Notwithstanding the large amount of work done in the '60s and in the '70s, the problem is still a matter of investigations (see for example Anderson & Lloyd 1990; Ougolnikov 1999; Ougolnikov & Maslov 2002; Ekstrom 2002; Ougolnikov et al. 2004; Postylyakov 2004; Mateshvili et al. 2005). While it is well beyond the purposes of the present work to explore the problem from a theoretical point of view, we deemed it interesting to introduce a simple single-scattering model, which is useful both for understanding the basic principles of twilight and to providing a quick comparison to the observed data. The assumptions and the model itself are discussed in Appendix A, to which we refer the interested reader for the details, while here we concentrate on the model predictions only.

The calculated zenith ($\alpha=0$) ${\it UBVRI}$ sky brightness as a function of Sun zenith distance $\zeta$, computed using the average broad band extinction coefficients for Paranal (Patat 2003a), are plotted in Fig. 1. As one can immediately see, the single scattering component drops below the night sky brightness between $\zeta=99$$^\circ$ and $\zeta=100^\circ$, indicating that from this point on multiple scattering is the only contributor to the observed flux, as shown by Ugolnikov & Maslov (2002) on the basis of polarization measurements.

	Figure 1: Model twilight sky brightness at zenith. The thick curves include the night sky contribution, while the thin lines indicate only the scattered component. The vertical dashed-dotted line marks the Sun zenith distance when the lower boundary layer height is 120 km. The upper scale indicates the lower Earth's boundary layer height in km.
Open with DEXTER

Another aspect to be considered is that the transition to the flatter part of the atmospheric density profile (see Fig. A.2) definitely occurs during the multiple scattering-dominated phase. Since multiple scattering takes place with higher probability where the density is higher, i.e. in the lower atmospheric layers, the explanation given by Tyson & Gal (1993) for the observed rate of decline in brightness during twilight does not seem to be correct. In fact, these authors interpret the observed values as the pure consequence of the lower shadow boundary height change, neglecting extinction and multiple scattering. They conclude that, since their observations have been taken when $100\leq h_z \leq 400$ km (where h_z is the height of the lower Earth's shadow boundary along the zenith direction; see also Appendix A), the sky brightness rate is directly related to the slope of the density law in that region of the atmosphere. Nevertheless, calculation of the Sun's ephemeris for the site and epoch of Tyson & Gal's observations shows that, in the most extreme case (see their Table 1, R filter), it was 1$^\circ$ . $8 \leq\varphi \leq 7$$^\circ$ . 4 (where $\varphi=\zeta-90$). As the reader can figure out from Table A.1, this implies that h_z<80 km in all cases, i.e. well within the steep part of the density profile. Therefore, that the observed rate and the rate expected from pure single scattering in the higher atmospheric layers are consistent is just a coincidence.

  
3 Observations, data reduction and calibration

In order to measure the twilight sky brightness on Paranal, we used archival calibration data obtained with the FOcal Reducer/low-dispersion Spectrograph (hereafter FORS1), mounted at the Cassegrain focus of the ESO-Antu/Melipal 8.2 m telescopes (Szeifert 2002). The instrument is equipped with a $2048\times2048$ pixel (px) TK2048EB4-1 backside-thinned CCD and has two remotely exchangeable collimators, which give a projected scale of 0 $.\!\!^{\prime\prime}$2 and 0 $.\!\!^{\prime\prime}$1 per pixel (24 $\mu$m $\times$ 24 $\mu$m). According to this collimator, the sky area covered by the detector is 6$^\prime$ . 8 $\times$ 6$^\prime$ . 8 and 3$^\prime$ . 4 $\times$ 3$^\prime$ . 4. For this study we selected only the data obtained with the lower resolution collimator and the 4-port high-gain read-out mode, since this combination is used the most for imaging with FORS1. The read-out noise is 5.5 electrons (e^-) with this setup.

For our purposes, we selected two sets of data. The first is composed of broad-band ${\it UBVRI}$ twilight sky flats (hereafter TSF), which are regularly obtained as part of the calibration plan. In the current implementation, the observing software estimates with a simple algorithm the integration time for the first exposure in a series of 4 frames after taking a test exposure. Subsequent exposure times are adjusted on the basis of the previous exposure level and this allows one to obtain high signal-to-noise images with a rather constant count level, which is typically around 20 000 ADUs. This is achieved with exposure times that range from 0.25 s up to 5 min. Given these values, the sensitivity of FORS1 in the various passbands and the typical twilight sky brightness behavior (see for example Fig. 1), these observations are expected to approximately cover the Sun zenith distance range 94 $^\circ \leq \zeta \leq 101^\circ$, i.e. still within the nautical twilight. Since the most important part of this analysis concerns the deep twilight, it is clear that an additional set of data is required to complement the sky flats.

The calibration plan of FORS1 includes the observation of standard star fields (Landolt 1992) in ${\it UBVRI}$ passbands, which are regularly taken during twilight, typically just after the sky flat sequence is completed. For calibration purposes, a fraction of these exposures are obtained using relatively long integration times (typically 40 s for U and 20 s for ${\it BVRI}$) which, at an 8 m-class telescope, are sufficient to bring the sky background to exposure levels that are suitable for our purposes. In fact, the bulk of these observations covers the range 100 $^\circ \leq \zeta\leq 112^\circ$, i.e. well into astronomical twilight. For the sake of clarity we indicate them as long-exposure standards (LES).

In order to collect a statistically significant sample, we retrieved from the ESO Archive all suitable TSF obtained in ${\it UBVRI}$passbands from 1 Jan. 2005 to 30 Sep. 2005, for a total of 1083 frames (U: 148, B: 208, V: 226, R: 261, I: 240). Since a much higher night-to-night spread is expected in the deep twilight phase due to the natural fluctuations of the night sky emission and also because LES are less frequently obtained than TSF, a longer time interval must be considered. To this aim, we have put together the LES sample collecting all suitable images obtained from 1 Jan. 1999 (i.e. shortly after the beginning of FORS1 operations) to 30 Sep. 2005, thus covering almost 6 years, for a total of 3388 frames (U: 923, B: 611, V: 609, R: 635, I: 610). All images were processed within the xccdred package of IRAF. Due to the large amount of data and the purpose of this work, the bias subtraction was performed using only a pre-scan correction, while flat-fielding was achieved using a stack of all TSF in each given passband as a master flat, which was then adopted to correct each TSF and LES frame.

Due to the nature of the data, there is no need to take care of the possible presence of crowded stellar fields or bright extended objects, as is the case for night sky brightness estimates (see Patat 2003b). For this reason, the background in each image was measured using a simple and robust mode estimator. To avoid possible vignetting and flat-fielding problems, only the central $1024\times 1024$ pixels were considered. On these spatial scales and due to improper flat-fielding, FORS1 is known to show variations on the order of a few percent, while the gradients in the twilight sky are much smaller (see Chromey & Hasselbacher 1996) and can be safely neglected. Therefore, the mode <I> of pixel intensity distribution is assumed as the best estimate of the sky background. For each filter, this is converted into a surface brightness in the Johnson-Cousins system via the following relation

$$b=-2.5 \log <I> + 2.5 \log \left(p^2 t_{\rm exp}\right) + m_0$$

where p is the pixel scale (arcsec px^-1), $t_{\rm exp}$ is the exposure time (seconds), and m₀ the instrumental photometric zeropoint for the given passband. No color correction was applied for two reasons: a) color terms in FORS1 are very small (Patat 2003a); b) color correction depends on the intrinsic color that, for the twilight sky, changes according to time and position. Given this and the fact that night-to-night fluctuations are the dominating source of noise in the measurements, the color correction can be safely neglected. We have used the average values of the instrumental zeropoints computed for the period of interest using data taken only during photometric nights. Variations in the zeropoints are known to take place (Patat 2003a), mainly due to the aging of telescope reflective surfaces but, again, these are smaller than the inherent sky variations. Finally, no airmass correction was applied, which is expected to produce an additional spread in the data close to the end of astronomical twilight and a systematic increase in the average sky brightness in that region. Due to the large number of pixels used (N>10⁶), the typical internal photometric error is expected to be less than 1%.

	Figure 2: Distribution of twilight observations in Alt-Az coordinates for TSF (filled symbols) and LSE (empty symbols). The astronomical azimuth has been replaced with the difference in azimuth between the telescope pointing and the Sun, $\Delta a_\odot$.
Open with DEXTER

  
4 Twilight sky brightness at ESO-Paranal

Since the data were not obtained specifically for twilight brightness measurements, they are inhomogeneously distributed on the sky. In fact, in an Alt-Az plot where the ordinary azimuth is replaced by the difference in azimuth between the sky patch and the Sun ( $\Delta a_\odot\equiv a-a_\odot$), the data points tend to cluster in two regions, which correspond to evening and morning observations (Fig. 2). Besides target azimuth and altitude we computed a series of other quantities for each data point, which are relevant for the subsequent analysis. These include Sun azimuth and altitude, Sun-target angular separation, Moon phase, Moon altitude, and Moon-target angular separation. To avoid contamination from scattered moonlight in the LSE, we selected only those data points for which the Moon was below the horizon.

Since the twilight sky brightness for a given Sun zenith distance changes with the position on the sky, it is necessary to make a selection on the Alt-Az coordinates to study its behavior as a function of $\zeta$. Given the nature of the available data, which appear to be rather concentrated (Fig. 2), it seems reasonable to restrict the analysis only to zenith region. In order to have a sufficient amount of measurements, we used all data points with zenith distance $\vert\alpha \vert\leq 40^\circ$, which is of course expected to cause some additional spread in the observed relation. The results are presented in Figs. 3-7 for the ${\it UBVRI}$ passbands.

Figure 3: Zenith twilight sky surface brightness in the U passband from TSF (empty symbols) and LSE (filled symbols) with $\vert\alpha \vert\leq 40^\circ$. The vertical dashed lines mark the end of nautical ( left) and astronomical ( right) twilight. The solid curve traces the simple model described in the text. The horizontal lines indicate the U dark-time night-sky brightness average value (dashed) for Paranal (Patat 2003a) and the $\pm$1$\sigma$ levels (dotted).

Open with DEXTER

	Figure 4: Same as Fig. 3 for the B passband.
Open with DEXTER

	Figure 5: Same as Fig. 4 for the V passband.
Open with DEXTER

	Figure 6: Same as Fig. 3 for the R pass band.
Open with DEXTER

	Figure 7: Same as Fig. 3 for the I passband.
Open with DEXTER

As expected, the single scatter model drops much faster than the actual observations. While for ${\it VRI}$ the model deviates from the data around $\zeta\sim 97^\circ$, for B and especially for U the model already underestimates the surface brightness at $\zeta\leq96^\circ$, indicating that multiple scattering is more efficient at shorter wavelengths. This agrees with the findings of Ougolnikov & Maslov (2002), who show that the contribution of single scattering in the phases immediately following sunset is about 40%, 60%, 70%, and 80% in U, B, V, and R, respectively. These fractions remain roughly constant until $\zeta\leq 95^\circ$, after which the role of single scattering becomes weaker and weaker and multiple scattering rapidly takes over. In all passbands, the night sky brightness level is reached at around $\zeta=105^\circ{-}106^\circ$.

In order to give a more quantitative description of the observations, we fitted the surface-brightness data in the range 95 $^\circ \leq \zeta \leq 105^\circ$ using second order polynomials of the form $a_0 + a_1 (\zeta-95) + a_2 (\zeta-95)^2$, with $\zeta$ expressed in degrees and the surface brightness in mag arcsec^-2. The results are presented in Table 1, where we also report the RMS deviation from the fitted function $\sigma$ and the slope $\gamma$deduced from a linear fit to the data in the range 95 $^\circ \leq \zeta \leq 100^\circ$, i.e. during the interval typically used to obtain TSF exposures, when the contribution by the night sky is still moderate. A first aspect to be noticed is the spread shown by the data points around the mean laws, which is due to the night-to-night variations in the atmospheric conditions. The dispersion becomes particularly large in the I passband, where the fluctuations appear to be quite pronounced. We notice that the decay rate during nautical twilight tends to decrease for increasing values of wavelength.

To convert the values of $\gamma$ reported in Table 1 into surface brightness variation per unit time, one has to multiply them by ${\rm d} \varphi/{\rm d}t$, which is given by:

$$\frac{\rm d\varphi}{{\rm d}t} = \cos \varphi \sin H_\odot \cos \delta_\odot \cos \phi \;\frac{{\rm d}H_\odot}{{\rm d}t}$$

where $H_\odot$ and $\delta_\odot$ are the hour angle and declination of the Sun, respectively, while $\phi$ is the site latitude (see for example Smart 1977). Since d $H_\odot/{\rm d}t\simeq0.25$ deg min^-1, for $\varphi \sim0$ (i.e. at the time of sunrise and sunset) for Paranal ($\phi=-24$$^\circ$ . 6) one obtains ${\rm d} \varphi/{\rm d}t=0.23$ deg min^-1 at the equinoxes ( $\delta_\odot=0$). Applying this factor to the values reported in the last column of Table 1, one obtains brightness decline rates that range from 0.28 (U) to 0.22 (I) mag arcsec^-2min^-1. These values can be directly compared with those reported by Tyson & Gal (1993). The data available to these authors did not allow them to quantify the differences between the various filters, so they report a mean value of $\gamma=0.23\pm 0.02$ mag arcsec^-2 min^-1 that is within the range defined by our data.

Table 1: Twilight sky brightness fitted parameters in the range 95 $^\circ \leq \zeta \leq$ 105$^\circ$. All values are expressed in mag arcsec^-2.

	Figure 8: Broad-band zenith twilight-sky colors. The curves have been computed using second oder polynomials fitted to the observed data. For comparison, the colors of the Sun are U-B=0.13, B-V=0.65, V-R=0.52 and V-I=0.81, while those of the night sky at Paranal are U-B=-0.36, B-V=1.03, V-R=0.74, and V-I=1.90(Patat 2003a).
Open with DEXTER

With the aid of the second-order best-fit relations, we computed the color curves presented in Fig. 8. Due to the dispersion of the observed data, these colors can only be regarded as indicative, especially in the region $\zeta>102^\circ$, where the inherent night-sky brightness fluctuations start to be significant, particularly in the red passbands. It is interesting to note that while U-B and V-R colors do not change very much as the Sun sinks below the horizon, significant changes take place in B-V and especially in V-I. One expects that since multiple scattering boosts the light at shorter wavelengths with respect to the pure single scattering component, the overall color gets bluer and bluer as the Sun deepens below the horizon. Then, at some point, the night sky glow, which has completely different colors, starts to contribute and the colors progressively turn to those typical of the night sky. The observed U-B, B-V, and V-R curves indeed show this behavior (see Fig. 8), while V-I turns steadily redwards. This is due to the interplay between input Sun spectrum, scattering efficiency, extinction, multiple scattering and the emergence of the night glow, which combine in quite a complicated way. This is illustrated clearly in Fig. 9, where we present twilight spectra obtained on Paranal with FORS1 and his twin instrument FORS2. The sky spectra were extracted from spectrophotometric standard stars observations taken during twilight (see Table 2) and were wavelength and flux-calibrated with standard procedures in IRAF. The exposure times ranged from 10 s to 120 s and the signal-to-noise ratio was increased, thereby meshing all pixels in the direction perpendicular to the dispersion, after removing the region of the detector affected by the well-exposed standard star spectrum. For comparison, Fig. 9 also shows the typical dark time spectrum of Paranal (Patat 2003a) and the solar spectrum.

Figure 9: Twilight spectra obtained at Paranal at different Sun depression angles. For convenience, the spectra have been corrected to the corresponding V flux shown by the photometric data. Spectral resolution and sky patch differ from spectra to spectra (see Table 2). For presentation the spectra have been vertically shifted by the following amounts: +0.50 ( $\varphi =9.4$), +0.35 ( $\varphi =11.9$), -0.10 ( $\varphi =14.8$), -1.0 ( $\varphi >=18$). The dotted line traces the solar spectrum, normalized to the continuum of the $\varphi =9.4$ spectrum at 4500 Å.

Open with DEXTER

Table 2: Basic data for twilight sky spectra shown in Fig. 9.

With the exception of the first spectrum, which was obtained with a very low resolution ($\sim$130 Å FWHM), the remaining data allow one to detect quite a number of details. For $\varphi<12^\circ$, i.e. during the nautical twilight, the spectrum is rather different from that of the Sun, even though it shows clear solar features, like the CaII H&K lines and the G-band at about 4300 Å. The Rayleigh-scattered Sun flux clearly contributes in the region bluewards of 5000 Å down to $\varphi=15^\circ$, after which the pseudo-continuum of the night sky emission takes over.

For $\varphi<9^\circ$, the contribution by night-sky emission lines is very weak. Characteristic lines like [OI]5577 Å and [OI]6300, 6364 Å (the latter overimposed on a O₂ absorption band) are barely visible, while the OH Meinel bands start to appear above 8000 Å. One remarkable exception are the NaI D lines, which are known to be present during the so-called sodium flash (Rozenberg 1966). A similar phenomenon is present for the [OI]6300, 6364 Å doublet (see for example Roach & Gordon 1973), which is indeed already visible in Fig. 9 at $\varphi=9$$^\circ$ . 4. With the onset of astronomical twilight, the spectrum in the red is more and more dominated by the OH bands. Another remarkable fact is the behavior of the H₂O (6800, 7600 Å) and O₂ (7200 Å) molecular absorption bands. During the bright twilight, when single scattering is still relevant ($\varphi=6$$^\circ$ . 5), they already appear to be significant, but they become even deeper in the multiple scattering-dominated phase ($\varphi=9$$^\circ$ . 4), due to the longer optical path traveled by the multiply scattered photons. For higher values of $\varphi$, they progressively disappear due to the weakening of the scattered Sun's continuum.

	Figure 10: Comparison between the zenith twilight-sky brightness measured at Paranal (small symbols) and the Crimean Astrophysical Observatory (large symbols) for the V passband.
Open with DEXTER

In order to estimate the effects of site altitude on the twilight sky brightness, we compared the results presented here with data obtained at the Southern Laboratory of the Sternberg Astronomical Institute (Moscow, Russia) during three morning twilights on December 9-11, 2002 (see Fig. 10). This facility is located within the Crimean Astrophysical Observatory (CrAO), at a latitude of 44$^\circ$ . 7 North and 600 m above sea level (hereafter asl). The observations were performed using a wide-field CCD-camera with a field of view of 8 $^\circ\times6^\circ$ and exposure times that ranged from 0.01 to 18 s. Photometric calibration was achieved using field stars included in the Tycho-2 Catalog (Hog et al. 2000). The photometric passband of this instrument is fairly similar to the Johnson-Cousins V, with a color correction in the order of 0.01 mag for the (B-V) color range shown by twilight data.

The two data sets clearly show that the twilight background at CrAO is systematically brighter; the difference is constant during the dark twilight period and vanishes at nightfall. More precisely, the comparison between V band data in moonless conditions at CrAO (zenith) and ESO-Paranal ( $\vert\alpha\vert\leq 20^\circ$) shows that the mean difference in the Sun depression range 5$^\circ$ . 5 $\leq \varphi \leq 11$$^\circ$ . 0 is $\Delta V = 0.27\pm0.03$. On the other hand, the typical atmospheric pressure value for ESO-Paranal is P₁=743hPa, and P₂=961 hPa for CrAO during the observations. Interestingly, the magnitude difference implied by the pressure ratio at the two sites is $-2.5\log~ (P_1/P_2)=0.28$, which is indeed very similar to the measured difference $\Delta V$. Therefore, we can conclude that the deep-twilight sky brightness is proportional to the atmospheric pressure or, equivalently, to the atmospheric column density above the observer. In turn, this implies that the light undergoes multiple scattering throughout the whole atmosphere and not only in the upper layers. Given that the difference in altitude between Paranal and CrAO is only 2 km, the observations we present here suggest that some fraction of multiple scattering has to take place in the first few km above sea level.

  
5 Discussion and conclusions

In this paper we have presented for the first time absolute ${\it UBVRI}$twilight brightness measurements for the ESO-Paranal Observatory (Chile) spanning almost 6 years. These measurements will serve as reference values for the similar studies that will be soon conducted at Dome C, Antarctica, as part of the site-testing campaign. The planned in situ spectrophotometric measurements will finally clarify whether this exceptional location shows a lower aerosol content, as expected from both the icy soil and its large distance from the sea coast (Kenyon & Storey 2005).

The twilight sky brightness measurements presented here were obtained from VLT-FORS1 archival data not specifically taken for this kind of study. Also, the initial twilight phases (0 $^\circ \leq \varphi \leq6^\circ$) are not covered, due to the large telescope diameter. In a sense this is quite unfortunate, since for these low Sun depression angles the lower shadow's boundary passes through the atmospheric layers below $\sim$30 km (see Table A.1), where the ozone and aerosol stratospheric concentration is at a maximum. These phases are in fact used to retrieve ozone and aerosol profiles, using both intensity and polarization measurements (see for example Wu & Lu 1988; Ugolnikov et al. 2004; Mateshvili et al. 2005). Nevertheless, during the deep twilight, when the direct Sun radiation illuminates only the upper atmospheric layers and single scattering on air molecules becomes progressively less important, the amount of aerosols and ozone plays a relevant role through multiple scattering. Therefore, even though they have a much more complicated interpretation, deep twilight observations may still give some insight into the conditions in the lower atmosphere. In the discussion of the supposedly shorter twilight duration at Dome C (see Kenyon & Storey 2005), what really matters is the behavior during the deep twilight. An example of this kind of analysis is shown in Fig. 11, where the data obtained at Paranal are compared to the MCC++ model calculations for a site at 2.6 km asl (see Postylyakov 2004 for a detailed description). This code treats the radiative transfer in a spherical atmosphere including Ozone, aerosol, and molecular scattering, also taking the backscattering by the Earth's surface into account. As one can see, the overall behavior is reproduced fairly well. The deviations are possibly due to the differences in real and model aerosol, since multiple scattering is very sensitive to it. The model adopts a urban microphysical model for the first 10 km of the atmosphere, which is certainly different from what is expected for a desertic area close to the sea, as is the case of Paranal. Dedicated instruments for twilight sky-brightness monitoring coupled with detailed modeling may indeed give usefully contribute to the already existing site-testing tools by providing independent indications of the overhead aerosol profile.

Some interesting results are obtained by comparing the estimates obtained for Paranal (2600 m asl) with those of a significantly lower site like CrAO (600 m asl). Even though the bright twilight and night sky brightnesses are very close at the two sites, Paranal is about 30% darker than CrAO in the V passband during the deep twilight (see Fig. 10). Due to the higher altitude, Paranal suffers from a lower extinction that would then turn into a brighter twilight sky if all other atmospheric properties are identical and multiple scattering takes place mostly in the troposphere (5-10 km, Ougolnikov & Maslov 2002). The observations actually show the opposite behavior, and the brightness ratio is fairly consistent with the atmospheric pressure ratio (see previous section). A natural interpretation is that a fraction of the multiple scattering events takes place at heights that are lower than was originally thought, say below 3 km asl.

Whether this is due to the lower density of air molecules, to a smaller amount of ground-level aerosols, or to a combination of the two needs further investigation, as does a comparison with other astronomical sites at even higher altitudes, such as Mauna Kea.

	Figure 11: Comparison between Paranal zenith sky brightness and the MCC++ model (Postylyakov 2004). The dashed curve traces the pure model solution without the contribution of night sky emission.
Open with DEXTER

  
Appendix A: A simple semi-analytical model for twilight brightness

A.1 Basic assumptions

The model is based on the following simplifying assumptions: (i) Earth is a sphere with radius R₀=6380 km; (ii) the atmosphere extends up to $\Delta R=400$ km and the numerical density n(h) of the scattering particle density is given by the MSIS-E-90 model profile (Hedin 1991); (iii) the effect of atmospheric refraction can be neglected; (iv) the Sun is a point source and all incoming sun rays are parallel; (v) only single scattering, plus attenuation, is considered; (vi) Rayleigh scattering by air molecules is the only source of sunlight diffusion. While some of these assumptions are reasonable, (v) and (vi) are a bit crude and will certainly lead to discrepancies between the model result and the actual observations. For a more sophisticated single scattering model taking refraction and the presence of ozone and aerosols into account, the reader is referred to Adams et al. (1974).

A.2 Model geometry

As said, we assume that the Earth is a sphere of radius R₀ and that the observer is placed at an elevation $h_{\rm s}$ asl. Since we consider only low values of $h_{\rm s}$ (<3 km), we make the further simplifying assumption that the horizon is a plane tangent to the sphere of radius $R_0+h_{\rm s}$ in O (see Fig. A.1), i.e. neglecting the horizon depression. We indicate with $\varphi$ the Sun depression ($\varphi>0$) and the zenith distance of the generic sky patch under consideration with $\alpha$. For the sake of simplicity we derive the sky brightness only along the great circle passing through the zenith (Z) and the Sun. This angular distance is counted positively in the direction of the Sun, so that negative angles indicate sky patches in the anti-Sun direction. Under these simplifying assumptions, the lower boundary of the Earth's shadow is described by a straight line, which is tangent to the sphere in $\rm H_0$, and which is indicated by a dotted-dashed line in Fig. A.1. When the observer is looking into the generic direction $\alpha$, the corresponding line of sight crosses the lower shadow boundary in $\rm P_0$, and the contribution to the observed flux comes from all the scattering elements along the segment $\rm P_0 P_1$, where $\rm P_1$ indicates the intersection between the line of sight and the upper atmospheric boundary, which is placed at an altitude $\Delta R$ asl.

Figure A.1: Geometry of the problem. For the sake of clarity, the observer's elevation $h_{\rm s}$ has been exaggerated.

In this geometry, a volume element placed in P receives the sun light, which is attenuated along its path OP, and it scatters the photons in the observer's direction, with a scattering angle $\theta=\pi/2-\alpha+\varphi$. according to the scattering phase function $\Phi(\theta)$, which obeys the usual normalization condition

$$\int_{4\pi} \Phi(\theta) \; {\rm d}\Omega =1 .$$ (A.1)

Before reaching the observer, it undergoes extinction along OP. The total flux will finally result from the integration along the illuminated segment $\rm P_0 P_1$. In order to compute the required quantities, we must first derive a series of useful relations. Since one has to evaluate the particle density along the generic light path, it is fundamental to know the height $h=\rm GP$ asl for any given point along the trajectory. This is particularly simple for the travel between O and $\rm P_1$, for which one can easily derive the following relation:

$$h_{\rm OP_1}=\sqrt{l^2 +2l(R_0+h_{\rm s})\cos\alpha + (R_0+h_{\rm s})^2} - R_0$$ (A.2)

where l is the coordinate along $\rm OP_1$ (l=0 in O). Another useful relation is the one that gives the lower limit for the integral along the line of sight, i.e. the optical path l₀ between O and $\rm P_0$. In order to get the required expression, one first needs to find the length of $\rm H_0P_0$ and h₀, i.e. the height asl of $\rm P_0$. These can be obtained after some trigonometry and are respectively,

$${\rm H_0P_0}=R_0 \left [\tan \varphi - \frac{\sin \alpha}{\co... ...\varphi}{\cos \varphi} -\frac{h_{\rm s}}{R_0}\right) \right ]$$

and

$$h_0=\sqrt{{\rm H_0P_0}^2 + R_0^2} - R_0 .$$

With these two equations at hand, one can easily find l₀:

$$\begin{array}{lll} l_0 & = & \sqrt{(R_0+h_{\rm s})^2\cos^2 \a... ...0-h_{\rm s})} \\ & & -(R_0+h_{\rm s})\cos \alpha. \end{array}$$

The upper limit l₁ along the line of sight turns out to be:

$$\begin{array}{lll} l_1 &=&\sqrt{(R_0+h_{\rm s})^2\cos^2 \alph... ...^2-h_{\rm s}^2}\\ & &-(R_0+h_{\rm s})\cos \alpha. \end{array}$$

The next step is the calculation of $h_{\rm QP}$, the height along the sun ray. In order to do so, we introduce the perigee height $\delta$, i.e. the minimum distance between the Earth's surface and the sun ray passing through P, which can be expressed as a function of the lcoordinate along $\rm P_0 P_1$ as:

$$\delta=(l-l_0) \cos (\alpha-\varphi).$$

Having this in mind, one can derive the length of QH as follows:

$$QH=\sqrt{2R_0(\Delta R-\delta) + \Delta R^2 - \delta^2}.$$

If we then introduce a coordinate q along QP (with q=0 in Q), considering the fact that $({\rm QH}-q)^2 + (R_0+\delta)^2=(R_0+h)^2$, we finally obtain

$$h_{QP}=\sqrt{(QH-q)^2 + (R_0+\delta)^2} - R_0$$ (A.3)

with $0\leq q \leq \rm QP$. The optical path of the unscattered sun ray QP is given by the sum of QH and HP, the latter being ${\rm HP}={\rm H_0P_0}-(l-l_0) \sin~(\alpha-\varphi)$. By means of Eqs. (A.2) and (A.3) one can now compute the height asl (and hence the particle number density) along the light path. The height along the zenith direction, h_z=TZ₀, can be readily derived and is given by

$$h_z=R_0 \frac{1-\cos \varphi}{\cos \varphi} \cdot$$

From this one can figure out how fast h_z grows when the sun depression $\varphi$ increases (see also Table A.1). This, coupled to the rapid decrease of the atmospheric density as a function of height, is the cause for the very quick drop in the twilight sky surface brightness.

Table A.1: Height of lower boundary of Earth's shadow at zenith (h_z) and at 60$^\circ$ from zenith along the solar meridian in the Sun's direction.

Figure A.2: Normalized density (solid line) and temperature (dashed line) profiles according to the MSIS-E-90 model (Hedin 1991). For comparison, the dots trace the values of the US Standard Atmosphere (McCartney 1976, Table 2.6).

A.3 Density profile and optical depth

As anticipated, for the density profile we adopted the MSIS-E-90 model density profile (Hedin 1991), which is presented in Fig. A.2. As the plot shows, the global profile can be roughly described by two laws: one exponential (for h<120 km, the so-called homosphere) and a power law (for 120<h<400 km, the so-called thermosphere). Clearly, as the Sun depression increases, the lower Earth's shadow boundary will pass through more and more tenuous atmospheric layers and, therefore, the change in density slope should turn into a change in the twilight sky surface brightness decline. Under the assumptions made for this simplified model, this should happen when $\varphi\simeq101^\circ$, i.e. close to the end of nautical twilight.

Since in the following section we are interested in the product between the number density and the extinction cross section, we can derive this for a given passband by assuming the measured extinction coefficient $\kappa(\lambda)$ (in mag airmass^-1) and integrating the previous density profile along the vertical (i.e. at airmass 1). In fact, assuming that all the extinction comes from Rayleigh scattering, one can write:

$$\tau_z(\lambda) = n_0 \; C_{\rm ext}(\lambda) \int_{h_{\rm s}} ^{\Delta R} \frac{n(h)}{n_0} \; {\rm d}h.$$

Then, considering that $\tau_z(\lambda)=1.086\kappa(\lambda)$, the product $n_0 \; C_{\rm ext}(\lambda)$ can be readily derived from the previous relation and used in all further optical depth calculations.

A.4 Scattering cross section and phase function

The $\lambda^{-4}$ wavelength dependency of the Rayleigh scattering cross section is implicitly taken into account by the extinction coefficients, which are to be considered as input data to the model and not as free parameters. As for the scattering phase function we used the canonical expression for air molecules (McCartney 1976):

$$\Phi(\theta) = \frac{3}{16\pi} \left[1+0.9324\cos^2~(\theta)\right]$$

where the multiplicative constant comes from the normalization condition expressed by Eq. (A.1).

A.5 Scattered flux

The scattered flux can be computed in the same way as done, for example, by Krisciunas & Schaefer (1991) for the Moon light. If $L_\odot$ is the luminosity of the sun at a given wavelength (expressed in photons per unit time and per unit wavelength), the flux received by the Earth at the top of the atmosphere is $F_\odot^0 = L_\odot/4\pi d^2$, where d=1 AU. If we now consider an infinitesimal volume element dV placed along the line of sight at a distance lfrom O, the number of scattered photons received by the observer per unit time, unit area, and unit wavelength is given by

$${\rm d}f = F_\odot^0 {\rm e}^{-\tau({\rm QP})}\; n[h(l)] \; C... ...ac{\Phi(\theta)}{l^2}\; {\rm e}^{-\tau ({\rm OP})}\;{\rm d}V .$$ (A.4)

Given the geometry of the problem, the infinitesimal volume element can be written as d $V={\rm d}S \; {\rm d}l\equiv \pi l^2 \phi^2 \; {\rm d}l$, where $\phi$ is the semi-amplitude of the angle subtended by dS at the distance of the observer. Since the solid angle subtended by dS is $\rm d\Omega=\pi \phi^2$, the surface brightness produced by the volume element is simply

$${\rm d}b=\frac{{\rm d}f}{\rm d\Omega}= F_{\odot}^0 {\rm e}^{-... ...lambda) \; \Phi(\theta) \; {\rm e}^{-\tau\rm (OP)} \; {\rm d}l$$

and the total surface brightness is finally obtained integrating along the line of sight within the illuminated region:

$$b=\int_{l_0}^{l_1} {\rm d}b\; .$$

Finally, to take the contribution by the night sky emission into account, we added to the computed flux the one implied by the average values measured for Paranal ( $h_{\rm s}=2.6$ km) and reported by Patat (2003a, Table 4). With this, the model is completely constrained and there are no free parameters.

Acknowledgements

We wish to thank S.L. Kenyon and J.W.V. Storey for inspiring this work and the ESO Archive Group for the support received during the data retrieval. O. Ugolnikov is supported by a Russian Science Support Foundation grant. Our gratitude also goes to the referee, Dr. A. Tokovinin, for his useful suggestions and comments.
This paper is based on observations made with ESO Telescopes at Paranal Observatory.

References

• Adams, C. N., Plass, G. N., & Kattawar, G. W. 1974, J. Atm. Sci., 31, 1662 [NASA ADS] (In the text)
• Anderson, D. E., & Lloyd, S. A. 1990, JGR, 95, 7429 [NASA ADS] (In the text)
• Blättner, W. G., Horak, H. G., Collins, D. G., & Wells, M. B. 1974, Applied Optics, 13, 534 [NASA ADS] (In the text)
• Chromey, F. R., & Hasselbacher, D. A. 1996, PASP, 108, 944 [NASA ADS] [CrossRef] (In the text)
• Divari, N. B., & Plotnikova, L. I. 1966, Sov. Astron., 9, 840 [NASA ADS] (In the text)
• Ekstrom, P. 2002, SPIE, 4815-14 (In the text)
• Hedin, A. E. 1991, J. Geophys. Res., 96, 1159 [NASA ADS] (In the text)
• Hog, E., et al. 2000, The Tycho-2 Catalogue on CD-ROM, Copenhagen University Observatory (In the text)
• Kenyon, S. L., & Storey, J. W. V. 2005, PASP, in press [arXiv:astro-ph/0511510] (In the text)
• Krisciunas, K., & Schaefer, B. E. 1991, PASP, 103, 1033 [NASA ADS] [CrossRef] (In the text)
• Landolt, A. U. 1992, AJ, 104, 340 [NASA ADS] [CrossRef] (In the text)
• Lawrence, J. S., Ashley, M. C. B., Tokovinin, A., & Travouillon, T. 2004, Nature, 431, 278 [NASA ADS] [CrossRef] (In the text)
• Mateshvili, N., Fussen, D., Vanhellemont, F., et al. 2005, J. Geophys. Res., 110, D09209 [CrossRef] (In the text)
• McCartney, E. J. 1976, Optics of the Atmosphere (New York: John Wiley & Sons) (In the text)
• Ougolnikov, O. S. 1999, Cosmic Research, 37, 159 [NASA ADS] (In the text)
• Ougolnikov, O. S., & Maslov, I. A. 2002, Cosmic Research, 40, 224 [NASA ADS] [CrossRef] (In the text)
• Ugolnikov, O. S., Postylyakov, O. V., & Maslov, I. A. 2004, J. Quant. Spec. Radiat. Transf., 88, 233 [NASA ADS] (In the text)
• Patat, F. 2003a, A&A, 400, 1183 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Patat, F. 2003b, A&A, 401, 797 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Postylyakov, O. V. 2004, J. Quant. Spec. Radiat. Transf., 88, 297 [NASA ADS] (the same issue with paper by Ugolnikov, Postylyakov, Maslov) (In the text)
• Roach, F. E., & Gordon, J. L. 1973, The light of the night sky (Boston, Dordrecht: Reidel) (In the text)
• Rozenberg, G. V. 1966, Twilight (New York: Plenum Press) (In the text)
• Smart, W. M. 1977, Textbook on Spherical Astronomy, Sixth Edition (Cambridge: Cambridge University Press) (In the text)
• Szeifert, T. 2002, FORS1+2 User's Manual, VLT-MAN-ESO-13100-1543, Issue 2.3
• Tyson, N. D., & Gal, R. 1993, AJ, 105, 1206 [NASA ADS] [CrossRef] (In the text)
• Wu, B., & Lu, D. 1988, Appl. Opt., 27, 4899 [NASA ADS] (In the text)