© ESO, 2013

1. Introduction

The current trend in astronomy has been to build larger and larger telescopes. The operating costs for running these large telescopes are high and careful planning of observations is needed for telescope time is always in demand. This means that more accurate predictions and estimation of the sky background are needed to understand how long an exposure is necessary for a given observation with a certain signal to noise ratio. The brightest natural source of light in the night sky, and therefore the greatest contributor to the sky background noise, is scattered moonlight. Having a reliable model of the moonlight for sky background estimation is critical. Also, many observers are trying to characterize faint objects spectroscopically, and knowing accurately the spectrum of the background would allow astronomers to predict which spectral features are observable within a given exposure time. By improving the scattered moonlight model within a sky background model, we can increase telescope scheduling efficiency.

The long standing scattered moonlight model used by ESO (European Southern Observatory) was by Walker (1987) and provides a table of the magnitudes for five photometric bands of the night sky at five different moon phases. It does not depend on the positions of the Moon or target observation and was measured during solar minimum. This model is limited in producing a scattered moonlight spectrum which is accurate enough for the current and future telescope operations.

The current, widely used model was developed by Krisciunas & Schaefer (1991), with 52 citations (e.g. Davies et al. 2013; Knoetig et al. 2013; Trinh et al. 2013). It is an empirical fit to data in the V-band taken at the 2800 m level of Mauna Kea. Even though the fit was separated into various specific functions, such as initial intensity from the Moon, Rayleigh and Mie scattering, it was still a parametrization based on only 33 observations in one photometric band specifically for Mauna Kea. The accuracy is between 8 and 23% if not near full Moon. In a previous paper, Noll et al. (2012), we presented a spectroscopic extension of the Krisciunas & Schaefer (1991) model. It was optimized for Cerro Paranal and covered the optical regime. The lunar albedo was taken to be constant with respect to wavelength and scaling factors for the different functions were introduced to better fit data from Cerro Paranal.

Our advanced scattered moonlight model works from the UV to the near-IR, but has only been fully tested in the optical, due to a current lack of data. It has been calibrated and investigated with 141 optical spectra and has an overall accuracy of σ ≲ 0.2 mag. The model has been split into physically based modules which are given by either physical models or the best current fits. The present version is optimized for Cerro Paranal, but can be modified for any location with information about its atmospheric properties. Our model is fully 3D and can predict the amount of scattered moonlight for any configuration of the Moon and target. It includes higher order scattering and is therefore still reliable at high zenith distances (when either object is near the horizon). Because our scattered moonlight model produces a spectrum, it can be used for finding spectral features and trends as well as photometric magnitudes.

The original purpose of designing an advanced scattered moonlight model, as part of the full sky background model, was for the Austrian accession to ESO to improve the ETC (Exposure Time Calculator), and thus enhance telescope scheduling and efficiency. Our optical sky background model, but with an older version of the scattered moonlight model, is described in Noll et al. (2012). The full sky background model estimates the amount of background light from 0.3 to 30 micron for Cerro Paranal. It includes all relevant components, such as scattered moonlight and starlight, zodiacal light, airglow line emission and continuum, scattering and absorption within the Earth’s atmosphere, and thermal emission from the atmosphere and telescope. Each component was designed with the latest knowledge and results in the field and was thoroughly checked with archival ESO data. The new scattered moonlight model is the topic of this paper. The zodiacal light is found using the prescription from Leinert et al. (1998) and the airglow model is based on local observations and semi-empirical modeling (Patat 2008). The scattering and absorption are calculated with two codes, one that is fully 3D for the scattering and a line-by-line radiative transfer code for molecular absorption. The scattering and absorption will also be described in this paper, in the context of the scattered moonlight. The thermal emission is estimated as a gray body. The new sky background model has already been implemented by ESO.

In this paper we will discuss the optical portion of the advanced scattered moonlight model. In a subsequent paper, we will discuss the extension of this model into the near-IR using X-Shooter data (Vernet et al. 2011). For the remainder of the paper, the term “moon model” refers to the scattered moonlight model. In Sect. 2, we describe the moon model. In Sect. 3, we compare our scattered moonlight model with the observations and previous models, and finally in Sect. 4 we present our conclusions.

2. The model

To accurately calculate the scattered moonlight, there are several components that must be considered. The simplest way to account for the various pieces is to follow the path of light from the source to the instrument. The source of the scattered moonlight is the Sun. Then the light is reflected by the Moon which depends on the lunar albedo. This is mainly a function of the lunar phase. Next, the moonlight enters the Earth’s atmosphere, where it is scattered by molecules (Rayleigh) and aerosols (Mie), and absorbed. The light can be scattered multiple times, including off of the ground, before it reaches the telescope. This depends on the properties of the atmosphere and the positions of the Moon and target observation.

In the following subsections, we will discuss each of these steps in more detail. Section 2.1 introduces the calibration data set. In Sect. 2.2 we discuss the solar spectrum, in Sect. 2.3 the lunar albedo, and in Sect. 2.4 a more general discussion of the set-up of the radiative transfer equations. Then, we have a more detailed discussion about Rayleigh scattering in Sect. 2.5, Mie scattering in Sect. 2.6, and absorption in Sect. 2.7.

Fig. 1 Properties of the FORS1 data from Patat (2008). The data shown has been preselected to not have any clouds during observation and plots b), c), and d) also have ≤ 20% humidity. a) The distribution of the amount of humidity at the time of observation as given by ESO Ambient Conditions Database. Only observations with ≤ 20% humidity were analyzed and is shown here as the filled in part of the histogram. b) Shows the distribution of the relative moon distance. A moon distance of 1 corresponds to the average moon distance of 384 400 km. c) Here we have plotted the distribution of lunar phase angles g. Over-plotted is the number of spectra taken during a waxing and waning Moon, which effects the solar selenographic longitude. d) The distribution of the angular distance between the Moon and target observation ρ is displayed. See Sect. 2.1 for more details.

2.1. Calibration data set

We calibrated our full sky background model including the improved scattered moonlight model with a FORS1 (FOcal Reducer/low dispersion Spectrograph) data set from Patat (2008). FORS1 was a low resolution optical spectrograph at Cerro Paranal. The data were taken between 1999 and 2005 with four different observing modes and consist of long-slit spectra, where plain sky could be extracted. For the details about the different observing modes, see Table 1. There are 1186 spectra, but only 26% have moonlight. For the analysis, we considered only 141 spectra with good weather conditions according to ESO’s ambient weather conditions database¹. The condition of good weather depended on two criteria, humidity and clouds. First, we only considered spectra taken when the humidity was ≤20%. Second, there must not be any detection of possible clouds in the night sky, within 5 h of the observation. Since we are measuring scattered moonlight, the entire sky must be free of clouds to ensure stable atmospheric conditions. Figure 1 shows several properties of the data: the distribution of the humidity with the cutoff, moon distance, lunar phase angles with waxing and waning, and angular distance ρ.

2.2. Solar spectrum

The source of the scattered moonlight is the Sun. We use the spectrum from Colina et al. (1996), which covers wavelengths from 0.12 to 2.5 micron. The UV and optical portion are from measurements from satellites and ground based data. The near-IR spectrum is from data by the NASA CV-990 aircraft with a model spectrum. The optical and near-IR spectrum has an uncertainty of ≤5%, while the UV has an uncertainty of ~20% due to solar variability, with a resolution of 1 to 2 nm. Figure 2 shows the optical solar spectrum. It resembles the black body spectrum with typical absorption lines for a G2 star.

Fig. 2 Optical solar spectrum from Colina et al. (1996) with principal absorption lines and bands labeled (Wallace et al. 2011). This was used as the source for our scattered moonlight model (Sect. 2.2).

2.3. Lunar albedo

The lunar albedo (A) determines the amount of sunlight that is reflected off the lunar surface towards the Earth. It is directly related to the intensity of light that will enter the Earth’s atmosphere, called I^∗. They are related by (Kieffer & Stone 2005), (1)I_sol is the spectrum of the light intensity of the Sun (see Sect. 2.2), and Ω_M is the solid angle of the Moon, where we take Ω_M = 6.4177 × 10^-5 sr. For our set of sky observations (Patat 2008), we found the moon distances (M_dis) with the online tool called JPL Horizons². The distance of 384   400 km is the average moon distance.

The lunar albedo A depends on several factors, including lunar phase, solar selenographic longitude, and wavelength. We use an empirical fit based on the ROLO survey (Kieffer & Stone 2005). They used over 100 000 images of the Moon in 32 different photometric bands in the optical and near infrared for certain lunar phases. The fit given by Kieffer & Stone (2005) is, (2)The lunar phase is described by the parameter g, where g = 0° corresponds to a full Moon and g = 180° is a new Moon. Φ is the solar selenographic longitude which describes whether the Moon is waxing or waning. Because the maria are not uniformly distributed on the lunar surface, the overall reflectivity of the Moon varies depending on which portion of the Moon is illuminated. We used the Astronomical Ephemeris by NASA³ to find whether the Moon was waxing or waning. In Eq. (2) the parameters a_i,λ, b_j,λ, and d_x,λ for the 32 different wavelengths λ are provided in Table 4 and the constants p_n are in Eq. (11) of Kieffer & Stone (2005).

The last three terms in Eq. (2) are to correct for the opposition effect (exponential functions) and to minimize the residuals between the data and fit (cosine function). The opposition effect is where the Moon has an enhanced brightness near full Moon (Whitaker 1969). We have neglected the terms that depend on the observer selenographic latitude and longitude. They have a small impact, maximum of ~7% over a full Saros (libration) cycle, and are variables that are not easily available in astronomical data.

Fig. 3 Lunar albedo based on the Kieffer & Stone (2005) fit as a function of the phase angle g, with the extrapolation where g > 97°, for three different wavelengths. The variation due to the solar selenographic longitude is shown in yellow. See Sect. 2.3 for more details.

The ROLO fit only covers the phase angles between g = 1.55° and 97°. We extrapolated this fit to g = 180°. For the terms that depend on Φ, when g > 97° we took the values of Φ at g = 97°. This physically means that we are assuming that the ratio of maria to highlands is the same when g > 97° as g = 97°. The other terms in Eq. (2) are still well behaved when g > 97° and fit reasonably well with our data (Patat 2008). It is expected from previous studies that the general dependence of A on g is exponential for all values of g (outside the opposition effect) (e.g. Lane & Irvine 1973), which we see in our extrapolation as well. As g increases the overall flux from the Moon decreases, and becomes less important compared with the other components of the sky background. Thus, we believe the extrapolation to higher g is reasonable. In Fig. 3, we show the dependence of the lunar albedo on moon phase including the extrapolation when g > 97°. Also shown is the effect of Φ which is a bit larger at redder wavelengths and is over-plotted as a band of values in yellow. At g = 97°, Φ can cause deviations of ~6.8%.

The fit from the ROLO data was only done photometrically, however we require a full spectrum of I^∗ for the scattered moonlight model. To attain a fit for the entire optical spectrum, we performed a simple linear interpolation between the fit values provided. This interpolation is shown in Fig. 4 for several different moon phases. In general the wavelength dependence is fairly smooth and becomes less significant towards larger g. There are some small variations where the photometric values are clustered.

According to Velikodsky et al. (2011), there are some discrepancies in the overall flux calibration amongst the various Moon observations. They find that the ROLO data at 603 nm is 13% too faint. They used an interpolation of the ROLO fit to find the albedo at 603 nm to compare it with their observations. Since there is no information about the wavelength dependence for this correction, we divide the albedo at all wavelengths by 0.87.

We currently neglect the effects of polarization. We assume that the change in the resulting scattered moonlight due to polarization is smaller than our current errors. In general the polarization of the lunar surface is small, less than 10% at maximum. However, the maria can have a decent amount of polarization (~30%) (Dollfus & Bowell 1971).

Fig. 4 Values of the lunar albedo given by ROLO (Kieffer & Stone 2005) presented as tick marks and the interpolation for the lunar albedo as a function of wavelength for several lunar phases g (Sect. 2.3).

2.4. Radiative transfer

The light from the Moon can be scattered many times or absorbed in the Earth’s atmosphere before reaching the telescope. We have developed a fully 3D single scattering code with estimates for double and higher order scattering. At the various scattering points, an effective airmass is calculated for absorption.

The transmission of light t can be directly related to the optical depth at zenith τ₀ or the extinction coefficient k and airmass X by the following, (3)We will use this formalism in the following sections to describe the effects of scattering.

Fig. 5 Geometry of the scattering in the Earth’s atmosphere (cf. Wolstencroft & van Breda 1967; Noll et al. 2012). Point M, where the moonlight enters the top of the atmosphere, is not in the same plane as the other points. The azimuth of M as seen from S is AM (not shown). S and T are at an azimuth A0 (also not shown) for an observer at O. Single scattering is shown at point S along paths s1 and s2. Double scattering occurs at point D or off the ground G (Sect. 2.4).

2.4.1. Single scattering

The single scattering can be done fully in 3D (Wolstencroft & van Breda 1967; Staude 1975; Bernstein et al. 2002). We wish to obtain the integrated scattered light towards the azimuth A₀ and zenith distance z₀. We take the Moon to be a point source. We consider the scattering path elements S of density n(σ), with σ being the radius vector from the center of Earth C to S. S is located between the top of atmosphere T to the observer O, at height H₀ above the surface (see Fig. 5). The distance between O and C is σ₀ = H₀ + R, where R is the radius of Earth (6371 km for the mean radius). For each path element S at distance s₂ from O, the contributions of radiation from M, where the moonlight enters the atmosphere, to the intensity at S along s₁ are considered. This intensity includes the scattering of light out of the path (and possible absorption), which depends on the effective column density of the scattering/absorbing particles. This is given by, (4)The scattering intensity also depends on the wavelength-dependent extinction cross section C_ext(λ) of the various processes, Rayleigh, Mie and absorption. They are related to the optical depth τ (see Eq. (3)) and Eq. (4) at the zenith by, (5)After scattering at S the light travels along path s₂ and is extincted by B₂(z₀,σ). The total scattering intensity per solid angle I_scat at (A₀,z₀) for the flux of the Moon F^∗, related to I^∗ (Eq. (1)), at (A_M,z_M) as seen by O, for Rayleigh and Mie scattering separately, is given by, (6)where (7)Equation (7) is for Rayleigh (R), Mie (M), and absorption (A). In Eq. (6), C_scat is the wavelength-dependent scattering cross section, which will deviate from C_ext if absorption occurs. The scattering phase function P depends on the scattering angle θ, the angle between the paths s₁ and s₂, and is related to the zenith distances and azimuths at S by, (8)Here, , , , and are the zenith distances and azimuths of the target and Moon at S.

In Eq. (7), the various C_ext are calculated by Eq. (4) for their respective densities. C_ext,A is calculated assuming the same particle distribution as Rayleigh scattering. This is optimal for molecular oxygen absorption.

We neglect the effect of polarization on the scattering phase function. Scattered moonlight does have some degree of polarization, similar to scattered sunlight, and has large areas of the sky with very little polarization (Horváth et al. 1998; Gál et al. 2001). The distribution of polarized light is not fully understood and would be difficult to implement.

For the vertical distribution of the scattering molecules, we use the standard barometric formula, (9)Here, h = σ − R, the sea level density n₀ = 2.67 × 10¹⁹ cm^-3, and the scale height h₀ = 7.99 km above the Earth’s surface (Staude 1975; Bernstein et al. 2002). For the troposphere and the lower stratosphere, where most of the scattering occurs, this is a good approximation. Cerro Paranal is at an altitude of H₀ = 2.64 km. For the thickness of the atmosphere, we take H_max = 200 km. For Rayleigh scattering, C_scat = C_ext, i.e. no absorption is involved.

For the height distribution of aerosols, we also use Eq. (9) with n₀ = 1.11 × 10⁴ cm^-3 and h₀ = 1.2 km. This is the tropospherical distribution of Elterman (1966). Dust, in particular soot, also absorbs radiation, thus C_scat is lower than C_ext. We used the OMI (Ozone Monitoring Instrument) satellite data to find the median ratio of 0.97 (range between 0.90 and 0.99) for the area around Cerro Paranal (Levelt et al. 2006). The model is not very sensitive to changes in the molecular or aerosol distribution, because we scale them with the extinction curve.

2.4.2. Double scattering

The above scattering equations only deal with single scattering. We also considered double scattering. For this, we divided the path s₁ into two different possible paths s_1,1 and s_1,2 where scattering and absorption can occur at the path elements D (for double scattering), in between the two paths (see Fig. 5). At D the scattering is treated in a similar way as in Eq. (6), but over the additional path elements. To simplify this integral and speed up the computing time, we made an approximation that allowed the integration to become an analytical expression with percent level accuracy (Eq. (4)). We assumed that the Moon coordinates where absolute coordinates with respect to the observer, such that the path s₂ ≪ R and the Moon is sufficiently far enough away to disregard slight angular shifts to the Moon coordinates for the different scattering points S and D. Near the horizon and when scattering occurs at a far distance from the telescope, this assumption breaks down. However, most of the scattering occurs in the lower troposphere and a few km from the observer (with mean distances of ~10 km). The calculation was performed along a grid, which exponentially grows with distance. Where θ is small and therefore the contribution from the forward scattering peak of aerosols is high, we increased the resolution of the grid. For the higher order scattering, including double scattering, absorption by molecules via an effective airmass was no longer considered.

In the case of double scattering, scattering can also occur at the ground, labeled as point G in Fig. 5. Since Cerro Paranal is located at 2635 m in a mountainous region, we took G to be at a height of 2 km, although the exact height has little influence. To estimate the ground reflection, we used the A1 profile from Sutter et al. (2007) based on soil samples taken in the Chilean desert. We then scaled it to the values from OMI (Levelt et al. 2006) at the three provided wavelengths for the region around Cerro Paranal. We followed the procedure of Kleipool et al. (2008) to obtain a proper average ground reflection, since there are many errors associated with this calculation. This gave us a wavelength dependent ground reflection at our desired location, as shown in Fig. 6.

Double scattering is the sum of two different single scattering events, which can be off of molecules, aerosols, or the ground. We sum the contributions from the three different types of scattering to arrive at a final double scattering intensity, called I_DS. Because the moon model was intended for the ETC, it needed to be computationally quick and the exact position of the Moon is irrelevant. For this purpose, we made additional approximations for I_DS. We computed weighted averages over both the zenith distance to the target z₀ and to the Moon z_M. The weighting for z₀ was based on all the sky spectra from Patat (2008) and for z_M from theoretical modeling of the positions of the Moon. I_DS is then only a function of the angular separation between the Moon and target ρ and the optical depth τ. The average uncertainty for this simplification is on the order of 5%, and is higher at large zenith distances but quickly decreases towards smaller zenith distances.

Fig. 6 Wavelength-dependent ground reflectance used for double scattering (see Sect. 2.4 for more details).

2.4.3. Multiple scattering

We now consider the intensity for multiple scattering, anything higher than two, called I_MS. The intensity for single scattering will be labeled as I_SS.

For I_MS, we compared the contributions for I_DS with I_SS. We assumed that the ratio between I_DS and I_SS would be the same for each consecutive order of scattering. We then summed over the geometric series and multiplied it by I_DS, (10)We imposed an upper limit for I_DS/I_SS of 0.9, so I_MS would not diverge. Similar to the weighted averages for I_DS, we calculated a weighted average over the difference between the azimuths for the target and Moon A₀ − A_M, which were weighted evenly. With Eq. (8), this gave a weighted average over the angular separation ρ. In the end, I_MS is only a function of optical depth τ, which can be related to wavelength.

The total scattering intensity for the Moon I_tot is simply the sum of the contributions from I_SS, I_DS, and I_MS. With the weighting over the various parameters, we created a table of correction factors f for the added amount of scattering due to I_DS(τ,ρ) and I_MS(τ). (11)The higher order scattering, I_DS and I_MS, typically contribute only a few percent and is strongest at high optical depth. Some typical scattering curves are shown in Fig. 7 for three different wavelengths. At the redder wavelengths, Mie scattering is dominating and the scattering function approaches an exponential. At bluer wavelengths, the characteristic Rayleigh scattering phase function is the most influential.

We have compared our scattering intensities I_tot with a radiative transfer code called libRadtran (Mayer & Kylling 2005) for a grid of zenith distances of the Moon and target, and the angular separation ρ in 10° steps for many different optical depths (corresponding to the full optical wavelength range). We have made this comparison to check the validity and accuracy of our scattering code. libRadtran is a widely used third party code for calculating radiative transfer in the atmosphere. It has several different numerical solvers, including the default which uses the plane parallel approximation. With this approximation, libRadtran is not accurate when either the source (Moon) or target are near the horizon. For all of the cases, 85% of I_tot from our code and libRadtran agreed within a relative error of 20%. For I_tot with zenith angles less than 70°, 75% agreed to within 10% relative error. The comparison between our code and libRadtran is shown in Fig. 8. With zenith angles ≤50°, there is very good agreement. As expected, at higher zenith angles where libRadtran is no longer accurate, the two codes diverge.

In the next three subsections, we will describe how Rayleigh, Mie and absorption are treated within the context of the scattering and absorption equations.

Fig. 7 Relative intensity to the light entering the atmosphere for scattering as a function of the scattering angle for a given geometry at three different wavelengths. The bluer wavelength is dominated by Rayleigh scattering and the redder wavelength by Mie scattering (Sect. 2.4).

Fig. 8 Comparison of our scattering code used in the advanced scattered moonlight model with libRadtran (commonly used radiative transfer code), including single, double and multiple scattering. In different colors we have plotted the relative intensities for each based on the maximum zenith distance of the Moon (source) or target (observer). When both zenith distances are below 50° the agreement is excellent. See Sect. 2.4 for more details.

2.5. Rayleigh scattering

The scattering off of molecules in the Earth’s atmosphere can be well described by Rayleigh scattering, which assumes the particles are much smaller than the wavelength. Rayleigh scattering is characterized with having a steep dependence on wavelength ~λ^-4 and at Cerro Paranal it is quite stable (Noll et al. 2012). For the extinction due to Rayleigh scattering we use the following parametrization from Liou (2002), (12)Here, λ is in μm, and P is the pressure at a height H of the observer, which we take to be 744 ± 1.5 hPa and 2.64 km, respectively (Noll et al. 2012).

The phase function for Rayleigh scattering is also well defined and is given by, (13)where θ is the scattering angle.

The Rayleigh scattering is taken at each path element S and D, and is a well described and stable component of the model.

2.6. Mie scattering

For the aerosol scattering, we have gathered information from several different sources and tried to build the most physical description. We take the Ångström Law fit derived by Patat et al. (2011) for Cerro Paranal in the optical wavelengths. We then decomposed this fit into several aerosol types given by Warneck & Williams (2012), and produced the phase function using a Mie scattering code for log normal distributions based on Bohren & Huffman (1983).

A typical way to describe the optical depths τ and wavelength dependence λ of aerosols in the optical is with the Ångström Law, τ = βλ^−α, where α and β are the fitting parameters (Ångström 1929). The fit by Patat et al. (2011) is based on 600 observations of 8 different spectrophotometric standard stars taken with FORS1. The data were taken within a short time span of 6 months in 2009, so there could be some unknown, long-time variations. They compared the leftover aerosol signature with what was predicted by the radiative transfer code LBLRTM (Line-By-Line Radiative Transfer Model; Clough et al. 2005). The tropospheric aerosols (height <10 km) needed to be scaled down by 25% in LBLRTM to match the observations. The optical depth at Cerro Paranal is very low and the aerosols are almost purely background aerosols, with no urban component (Patat et al. 2011). In the optical, their fit is given by, (14)where k₀ = 0.013 ± 0.002 and a = −1.38 ± 0.06 mag airmass^-1, valid from 0.4 to 0.8 μm. We use the best fit k₀ value quoted in Patat et al. (2011) of 0.014 mag airmass^-1. The aerosol component can fluctuate significantly, as can be seen in data taken at Mauna Kea (Buton et al. 2013). They have the Ångström exponent of a = −1.3 ± 1.4 and state that at 330 nm it can vary as much as 0.4 mag airmass^-1 during photometric nights.

We are interested in extending the scattered moonlight model into the infrared, as well as wanting a more physical basis for the aerosol extinction to more accurately calculate the phase function. Therefore, we decomposed the Ångström Law. We used the size distributions given in Warneck & Williams (2012) for the remote continental tropospheric and stratospheric aerosols. The tropospheric aerosols are split into three types: nucleation, accumulation, and coarse modes. Nucleation consists of newly produced particles from either direct emissions of combustion products or gas-phase condensation reactions. Accumulation is from coagulation of nucleation particles, condensation of products from the gas-phase chemical reactions onto particles, and chemical reactions in the aqueous phase of clouds. The mass of the coarse mode comes from mineral dust, sea salt, and biogenic material. The stratospheric aerosols are in a layer around 25 km above sea level and consist mostly of sulfuric acid particles with an admixture of nitrosyl sulfates and solid granules containing silicates (Warneck & Williams 2012). We use a log normal distribution given by, (15)N is the cumulative number density distribution in particles cm^-3, and r is the particle radius in μm. Table 2 provides the values for the parameters n, R, and log s of the various aerosol modes, which represent the number of particles, average radius, and a distribution parameter, respectively. To calculate the extinction curve for each aerosol type we used a Mie scattering IDL code for log normal distributions written by G. Thomas⁴ based on Bohren & Huffman (1983) for single particles (Grainger et al. 2004). With the extinction curves for the various aerosol types along with the fit for Cerro Paranal, we determined how much each curve needed to be scaled to match the overall fit. However, there are degeneracies. The refractive index N is not known for these background aerosols, so we varied it between 1.3 and 1.5. Varying N changes the scaling of the different aerosol components needed to match the fit. There is also a degeneracy between two of the aerosol types, remote continental tropospheric coarse mode and stratospheric, because their optical extinction curves are similar. In Fig. 9a we have plotted the various Mie extinction curves. The extinction from the different aerosol types with a range in N between 1.3 and 1.5 are shown, along with the total curve that is consistent with the observed one. It can easily be seen that the tropospheric nucleation mode has little to no effect, whereas the tropospheric accumulation has the greatest influence.

Fig. 9 Different properties of the various types of aerosols. In both plots the width of the line is due to the uncertainty of the refractive index (N ranges between 1.3 and 1.5). a) Shows the extinction coefficient of Mie Scattering at 500 nm. The range for the total Mie extinction curve are given by the following values for tropospheric nucleation, accumulation, and coarse, and stratospheric modes, respectively: 100, 100, 60, and 60% for N = 1.3 and 100, 45, 5, and 100% for N = 1.5. The latter is what we have used for the analysis. Also plotted is the Patat et al. (2011) fit with the reported errors. The total Mie extinction curve (after scaling the various components) is in good agreement with the fit in the optical regime. b) This shows the Mie phase functions for the different modes and the total Mie phase function. See Sect. 2.6 for more details.

With the current FORS1 data set (Patat 2008), there is little to no difference in the overall results between N = 1.3 or 1.5 and the various amounts of the different aerosol components. We chose N = 1.5, with 45% tropospheric accumulation, 5% tropospheric coarse modes, and 100% stratospheric aerosols, where the total matches the Patat et al. (2011) fit. The tropospheric nucleation contributed extremely little to the overall aerosol scattering, and so was left at 100%. Since Cerro Paranal is located at 2635 m, we lowered the amount of tropospheric aerosols and kept all the stratospheric particles.

With the log normal Mie scattering code we also calculated the Mie phase function of the different aerosol components and the overall Mie phase function. This method of determining the phase function gives a realistic forward scattering peak, but can poorly approximate the back scatter. This is caused by the fact that the aerosol particles are not spherically symmetric, and therefore the scattering phase function tends to flatten out at higher angles, rather than dip to a minimum as in the Mie approximation (e.g. Horvath et al. 2006). Since the forward scattering dominantes the back scatter by a few orders of magnitude, we prefer the Mie scattering phase function over the Henyey-Greenstein (Henyey & Greenstein 1941) approximation, which depends on the average asymmetry parameter of the aerosols. Additionally, the majority of the stratospheric aerosols in this region are sulfates, which are in droplets, and so are fairly round. The amount of coarse mode aerosol is low, which is the least round type of aerosol. This will also provide another method to break the degeneracies amongst the various aerosol modes. Finally, in astronomical observations it is rare to have the target and Moon at angular distances >120°, where the back scatter peak becomes important. In Fig. 9b, the phase functions for the different aerosol types are plotted along with the total Mie phase functions. Once again each mode is shown with a range of N between 1.3 and 1.5.

2.7. Absorption

The absorption is calculated using the radiative transfer code LBLRTM (Clough et al. 2005) with a merged atmospheric profile for Cerro Paranal and the effective airmass X_eff at point O calculated from the single scattering at S. LBLRTM is a widely used radiative transfer code in atmospheric sciences and it uses the atomic line database HITRAN (High-Resolution TRANsmission; Rothman et al. 2009). This gives the amount of absorption from the various atmospheric molecules. It requires an atmospheric profile, and we use a merged one, which combines the temperature, pressure, and chemical composition as a function of height from three different sources. We combined MIPAS (Michaelson Interferometer for Passive Atmospheric Sounding) equatorial profile (see Seifahrt et al. 2010) with one from GDAS (Global Data Assimilation System)⁵, and finally scale the quantities of the lower atmosphere to the local meteorological data taken at Cerro Paranal. The MIPAS profile contains height information on 30 molecular species. The GDAS database, which is based on measurements, provides modeled profiles of temperature, pressure, and the relative humidity up to a height of about 26 km on a 3 h basis with a spatial grid of 1° × 1°. The local meteorological data from the ESO MeteoMonitor⁶ gives the local temperature, relative humidity, and pressure. We currently use bimonthly mean, merged profiles for all the nights.

For each single scattering point S, an effective airmass X_eff relative to the direct path from the zenith to observer is calculated to scale the amount of absorption due to the given geometry of the Moon and target. We do not calculate X_eff from the higher order scattering, since this contribution is quite small. X_eff then modifies the transmission from the moonlight as calculated by LBLRTM.

For determining X_eff, we compare the intensities with and without extinction, I⁽⁺⁾ and I⁽⁻⁾, respectively. We first look at a specific path element for a given S. We find the column density of that path element for both Rayleigh B_n,R and Mie scattering B_n,M (similar to Eq. (4)). Then with C_scat and phase function P(θ), we calculate I⁽⁻⁾ by, (16)where, (17)To calculate I⁽⁺⁾, we use I⁽⁻⁾ and τ from Eq. (7), and is given by, (18)Here, τ_n is given at a reference wavelength, where the absorption is strong and the scattering is weak. Within this constraint, we arbitrarily choose a wavelength of 1 μm for the scattering at the zenith, which gives τ_R + M = 0.02.

Then with I⁽⁺⁾ and I⁽⁻⁾ along with the optical depth at zenith τ₀, X_eff is, (19)Here, τ₀ is the reference τ_R + M = 0.02 plus τ_A.

Fig. 10 Two examples of the quality of the model. The observed sky spectrum is shown in black with the full sky background model in red. The scattered moonlight portion of the model is given in blue. In both cases the model fits the observed spectrum fairly well. a) This observation was done with moderate amount of moonlight. The moonlight comprises roughly half of the overall sky background flux. b) This one was with a significant amount of moonlight, as can be seen by the increase in overall flux. The scattered moonlight completely dominates the amount of background light. See Sect. 3 for more details.

Finally, the scaled scattering intensity I_scat,A for absorption is given by, (20)The three main absorbing molecules in the optical are O₂, H₂O, and O₃. LBLRTM (Clough et al. 2005) with the merged atmospheric profile produces a transmission curve for the absorption. We use the particle distribution from Rayleigh scattering as an approximation for the absorbing molecules, which is optimal for O₂. The water vapor, however, is distributed differently from the main air molecules. The absorption from H₂O mostly occurs at low altitudes. Since the distribution of particles is scaled with the extinction curve, the amount of H₂O is fairly accurate. Changing the profile did not significantly change the resulting modeled spectrum.

The ozone is treated separately using the van Rhijn formula with a height of 25 km, which is given by, (21)Here X_oz is the effective airmass for ozone and z_M is the zenith distance to the Moon. There are two assumptions made. The first one is that we treat the ozone only at a height of 25 km and not at a distribution of heights. The thickness can vary with season and geographic location. This assumption allows for a much simpler calculation of X_oz. The second is that we are only considering the light that is coming along the direct path from the Moon, by using the Moon’s zenith distance. We assume that the majority of the light has not been scattered prior to passing through the ozone layer. Since most of the scattering takes place in the lower troposphere, which is distinctly below 25 km, this appears to be a safe assumption. Following Patat et al. (2011), we multiply the amount of ozone given in the standard atmospheric profile by 1.08.

Fig. 11 a) Mean and b) σ of the difference between the observed and modeled spectra for the full sky background of all the spectra, those with moonlight, and only those without moonlight. See Sect. 3 for more details.

3. Comparison with FORS1 data and previous models

We have compared our model with the FORS1 data (see Sect. 2.1 for details) (Patat 2008) and with the previous moon model from Noll et al. (2012) based on Krisciunas & Schaefer (1991). We found that overall, for the complete sky background model, an uncertainty of σ ~ 0.2 mag with the FORS1 data set. For the moon model, the uncertainty is ~0.15 mag.

Figure 10 shows two examples of how well the model fits the observed data. The spectrum in Fig. 10a has moderate moonlight and is around 1st quarter and (b) has significant moonlight and is near full Moon. The full sky model is consistent with the observed spectrum in both scenarios. Also shown is the scattered moonlight model which is a large portion of the overall sky background flux. For (b), the model tends to slightly overestimate the amount of background light.

To better understand the errors of the moon model, we made a mean and σ spectrum of the difference between the FORS1 observations and the full sky model for all the spectra, those with moonlight, and those without moonlight (Fig. 11). The mean spectrum is centered around zero, with increased fluctuations at the redder wavelengths due to airglow emission. The σ spectrum without moonlight is slightly lower, but all three groups tend to lie around 0.2 mag, and increase slightly at the red end. Loosely comparing the spectrum with moonlight to those without, gives a rough estimate of the quality of our advanced moon model.

Fig. 12 Means and uncertainties for the scattered moonlight versus observed data at several 4 nm wide continuum bands. The y-axis is the average of the observed minus modeled fluxes for data with good weather conditions and a significant amount of moonlight. Over-plotted is the same analysis with the previous model from Noll et al. (2012) based on Krisciunas & Schaefer (1991), labeled as KS91. The numbers below each point are the number of spectra considered. See Sect. 3 for more details.

To better evaluate only the scattered moonlight model, we subtracted the other components of the sky background, using the sky background model, from the observations. This, in principle, should provide only an observed scattered moonlight spectrum. We then compare this observed moon-only spectrum with the scattered moonlight model. Such a comparison contains all the errors associated with the full sky background model as well as the errors from the moon model. For an estimate of the uncertainty in the other sky background components, we looked at the full sky model without moonlight. We then computed a weighted σ to represent the uncertainty from the other sky background components, by comparing the relative flux from the moonlight with full sky background for the Moon observations used in this analysis. If there is less moonlight, then σ from the other components should be larger, versus an observation dominated by moonlight where σ should mostly come from the moon model. We did this comparison using all the data with a significant amount of moonlight (>100 phot s^-1   m^-2   μm^-1 arcsec^-2) and decent weather conditions. Also, we ignored one of the observing modes, Grism 600B, which had only 8 data points and the flux calibrations appeared to greatly deviate from the other observing modes. This left us with a total of 82 spectra. The results can be seen in Fig. 12. There is a jump in the mean around 500 nm where the majority of the observations changes from one observing mode to another. We believe this comes from errors in the flux calibration of the Patat (2008) data. For the mean of these spectra to be centered at zero, we needed to multiply the moon model by a factor of 1.2 to correct for the flux calibration. The error bars shown here are the uncertainties from the moon model squared minus the square of the weighted σ from the other components. At blue wavelengths σ is small, and it gets larger towards the red. The estimated weighted σ of the other sky background components is inflated in the UV and red wavelengths. The airglow/residual continuum model probably overestimates the flux, and so the residuals, especially in the red, are large. This causes the Moon contribution to be underestimated and σ to be inflated, so the last three continuum nodes are not shown for this reason. In Fig. 11a, this effect due to the airglow model can be seen from the downturn in the mean difference spectra. Also plotted in Fig. 12 with the same analysis is the previous moon model discussed in Noll et al. (2012), which is an extension of Krisciunas & Schaefer (1991) moon model and scaled to Patat (2008) data, labeled as KS91. The mean and σ for this previous model are worse than our new advanced moon model. KS91 has a larger σ at blue wavelengths and the mean is more off centered at redder wavelengths compared with this work. The mean and σ for our new advanced moon model and the previous one, before being corrected for the errors from the other sky background components are shown in Table 3, along with the weighted σ from the rest of the sky background.

To further show the improvement of our new advanced scattered moonlight model, we have plotted the full sky background model using the extrapolated Krisciunas & Schaefer (1991) moon model (labeled again as KS91) and our new moon model with an observed spectrum in Fig. 13a. It can be clearly seen that the new moon model fits the observed spectrum much better than the previous model. The relative residuals are plotted on panel (b).

Fig. 13 Comparison of the full sky background model with the previous moon model based on Krisciunas & Schaefer (1991), labeled as KS91, and then one presented in this paper. a) The observed spectrum along with the full sky background model with the old (KS91) approach and this work. b) This shows the relative residuals for the previous model and the work presented in this paper. See Sect. 3 for more details.

4. Conclusion

We have developed an advanced scattered moonlight model. It is based on physical processes and data, and has fewer empirical parametrizations than previous models. It is spectroscopic and traces spectral trends seen in observations. Currently, it has been evaluated for the optical range from 0.36 to 0.89 micron. The uncertainty in the moon model is around 0.15 mag, and the full sky background model has σ ~ 0.2 mag. There are several advantages of our new moon model:

• More physical than previous models (Krisciunas & Schaefer 1991; Noll et al. 2012).

   ⋄ Uses a fit for the lunar albedo based on more than 100 000 observations of the Moon (Kieffer & Stone 2005).

   ⋄ Uses fully 3D single scattering calculations.

   ⋄ Calculates fully 3D double scattering calculations, including ground reflection, then simplifies for faster computing.

   ⋄ Approximates multiple scattering by comparing double and single scattering.

   ⋄ Calculates the absorption and single scattering simultaneously.

   ⋄ Decomposes the Mie scattering Ångström Law into typical aerosol distributions for Cerro Paranal.

   ⋄ Calculates the Mie phase function based on the aerosols found by decomposing the Ångström Law.

• Provides a scattered moonlight spectrum, not just photometric magnitudes (Walker 1987; Krisciunas & Schaefer 1991).

• The uncertainties across the entire optical range are much lower than previous (Krisciunas & Schaefer 1991; Noll et al. 2012).

Acknowledgments

This study is carried out in the framework of the Austrian ESO In-kind project funded by BM:wf under contracts BMWF-10.490/0009-II/10/2009 and BMWF-10.490/0008-II/3/2011.

References

All Tables

All Figures

Fig. 1 Properties of the FORS1 data from Patat (2008). The data shown has been preselected to not have any clouds during observation and plots b), c), and d) also have ≤ 20% humidity. a) The distribution of the amount of humidity at the time of observation as given by ESO Ambient Conditions Database. Only observations with ≤ 20% humidity were analyzed and is shown here as the filled in part of the histogram. b) Shows the distribution of the relative moon distance. A moon distance of 1 corresponds to the average moon distance of 384 400 km. c) Here we have plotted the distribution of lunar phase angles g. Over-plotted is the number of spectra taken during a waxing and waning Moon, which effects the solar selenographic longitude. d) The distribution of the angular distance between the Moon and target observation ρ is displayed. See Sect. 2.1 for more details.

In the text

	Fig. 2 Optical solar spectrum from Colina et al. (1996) with principal absorption lines and bands labeled (Wallace et al. 2011). This was used as the source for our scattered moonlight model (Sect. 2.2).
In the text

	Fig. 3 Lunar albedo based on the Kieffer & Stone (2005) fit as a function of the phase angle g, with the extrapolation where g > 97°, for three different wavelengths. The variation due to the solar selenographic longitude is shown in yellow. See Sect. 2.3 for more details.
In the text

	Fig. 4 Values of the lunar albedo given by ROLO (Kieffer & Stone 2005) presented as tick marks and the interpolation for the lunar albedo as a function of wavelength for several lunar phases g (Sect. 2.3).
In the text

Fig. 5 Geometry of the scattering in the Earth’s atmosphere (cf. Wolstencroft & van Breda 1967; Noll et al. 2012). Point M, where the moonlight enters the top of the atmosphere, is not in the same plane as the other points. The azimuth of M as seen from S is AM (not shown). S and T are at an azimuth A0 (also not shown) for an observer at O. Single scattering is shown at point S along paths s1 and s2. Double scattering occurs at point D or off the ground G (Sect. 2.4).

In the text

	Fig. 6 Wavelength-dependent ground reflectance used for double scattering (see Sect. 2.4 for more details).
In the text

	Fig. 7 Relative intensity to the light entering the atmosphere for scattering as a function of the scattering angle for a given geometry at three different wavelengths. The bluer wavelength is dominated by Rayleigh scattering and the redder wavelength by Mie scattering (Sect. 2.4).
In the text

Fig. 8 Comparison of our scattering code used in the advanced scattered moonlight model with libRadtran (commonly used radiative transfer code), including single, double and multiple scattering. In different colors we have plotted the relative intensities for each based on the maximum zenith distance of the Moon (source) or target (observer). When both zenith distances are below 50° the agreement is excellent. See Sect. 2.4 for more details.

In the text

Fig. 9 Different properties of the various types of aerosols. In both plots the width of the line is due to the uncertainty of the refractive index (N ranges between 1.3 and 1.5). a) Shows the extinction coefficient of Mie Scattering at 500 nm. The range for the total Mie extinction curve are given by the following values for tropospheric nucleation, accumulation, and coarse, and stratospheric modes, respectively: 100, 100, 60, and 60% for N = 1.3 and 100, 45, 5, and 100% for N = 1.5. The latter is what we have used for the analysis. Also plotted is the Patat et al. (2011) fit with the reported errors. The total Mie extinction curve (after scaling the various components) is in good agreement with the fit in the optical regime. b) This shows the Mie phase functions for the different modes and the total Mie phase function. See Sect. 2.6 for more details.

In the text

Fig. 10 Two examples of the quality of the model. The observed sky spectrum is shown in black with the full sky background model in red. The scattered moonlight portion of the model is given in blue. In both cases the model fits the observed spectrum fairly well. a) This observation was done with moderate amount of moonlight. The moonlight comprises roughly half of the overall sky background flux. b) This one was with a significant amount of moonlight, as can be seen by the increase in overall flux. The scattered moonlight completely dominates the amount of background light. See Sect. 3 for more details.

In the text

	Fig. 11 a) Mean and b) σ of the difference between the observed and modeled spectra for the full sky background of all the spectra, those with moonlight, and only those without moonlight. See Sect. 3 for more details.
In the text

Fig. 12 Means and uncertainties for the scattered moonlight versus observed data at several 4 nm wide continuum bands. The y-axis is the average of the observed minus modeled fluxes for data with good weather conditions and a significant amount of moonlight. Over-plotted is the same analysis with the previous model from Noll et al. (2012) based on Krisciunas & Schaefer (1991), labeled as KS91. The numbers below each point are the number of spectra considered. See Sect. 3 for more details.

In the text

Fig. 13 Comparison of the full sky background model with the previous moon model based on Krisciunas & Schaefer (1991), labeled as KS91, and then one presented in this paper. a) The observed spectrum along with the full sky background model with the old (KS91) approach and this work. b) This shows the relative residuals for the previous model and the work presented in this paper. See Sect. 3 for more details.

In the text