A multi-band map of the natural night sky brightness including Gaia and Hipparcos integrated starlight
Eduard Masana, Josep Manel Carrasco, Salvador Bará, Salvador J. Ribas
MMNRAS , 1–15 (2020) Preprint 6 January 2021 Compiled using MNRAS L A TEX style file v3.0
A multi-band map of the natural night sky brightness including
Gaia and
Hipparcos integrated starlight
Eduard Masana, ★ Josep Manel Carrasco, Salvador Bará, and Salvador J. Ribas , Departament Física Quàntica i Astrofìsica. Institut de Ciències del Cosmos (ICC-UB-IEEC), C Martí Franquès 1, Barcelona 08028, Spain Departmento Física Aplicada, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galicia, Spain Parc Astronòmic Montsec - Ferrocarrils de la Generalitat de Catalunya. Camí del Coll d’Ares, s/n, Àger (Lleida) 25691, Spain
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The natural night sky brightness is a relevant input for monitoring the light pollutionevolution at observatory sites, by subtracting it from the overall sky brightness determinedby direct measurements. It is also instrumental for assessing the expected darkness of thepristine night skies. The natural brightness of the night sky is determined by the sum ofthe spectral radiances coming from astrophysical sources, including zodiacal light, and theatmospheric airglow. The resulting radiance is modified by absorption and scattering before itreaches the observer. Therefore, the natural night sky brightness is a function of the location,time and atmospheric conditions. We present in this work GAMBONS (GAia Map of theBrightness Of the Natural Sky), a model to map the natural night brightness of the sky incloudless and moonless nights. Unlike previous maps, GAMBONS is based on the extra-atmospheric star radiance obtained from the
Gaia catalogue. The
Gaia -DR2 archive compilesastrometric and photometric information for more than 1.6 billion stars up to 𝐺 =
21 magnitude.For the brightest stars, not included in
Gaia -DR2, we have used the
Hipparcos catalogueinstead. After adding up to the star radiance the contributions of the diffuse galactic andextragalactic light, zodiacal light and airglow, and taking into account the effects of atmosphericattenuation and scattering, the radiance detected by ground-based observers can be estimated.This methodology can be applied to any photometric band, if appropriate transformationsfrom the
Gaia bands are available. In particular, we present the expected sky brightness for 𝑉 (Johnson), and visual photopic and scotopic passbands. Key words: light pollution - scattering - radiative transfer - atmospheric effects - instrumen-tation: photometers - site testing
The natural brightness of the night sky in different photometricpassbands is a relevant parameter for site characterization and lightpollution research, since it allows establishing a baseline againstwhich to evaluate the light pollution levels experienced in urban,periurban, rural and pristine dark sites.Brightness, in this context, is a shortcut name for the spectralradiance of the sky, weighted by the sensitivity function of thephotometric band, spectrally integrated over all wavelengths, andangularly integrated within the field-of-view of the detector.The measurement of the night sky brightness and the atmo-spheric conditions, in particular at sites of astronomical interest, hasa long history that continues to the present day (for example Turn-rose (1974) for Palomar and Mount Wilson Observatories; Benn & ★ E-mail: [email protected] (EM)
Ellison (1998) for La Palma; or Patat (2008) for Cerro Paranal). Thenatural sources of the night sky brightness were already outlinedin Roach & Gordon (1973). The reference paper by Leinert et al.(1998) studies and analyzes the different sources contributing to thesky brightness covering wavelengths from 𝜆 ≈
100 nm (ultravioletdomain) to 𝜆 ≈ 𝜇 m (far infrared). It includes airglow, zodiacallight, integrated star light, diffuse galactic light and extragalacticbackground light, and summarizes the best estimations at the timefor all of them. A similar work can be found in Noll et al. (2012)for the wavelength range from 𝜆 ≈
300 to 𝜆 ≈
920 nm, with anextended discussion on the airglow and atmospheric extinction atCerro Paranal. Duriscoe (2013) evaluates the artificial contributionto the sky brightness by subtracting the natural contributions fromthe overall sky radiance measured with all-sky imaging systems.Duriscoe (2013) used a model of integrated starlight (including dif-fuse galactic light) constructed from images made with the same © a r X i v : . [ a s t r o - ph . I M ] J a n E. Masana et al. equipment used for sky brightness observations, and modeling zo-diacal light and airglow from previously available data.In parallel, several models to compute the artificial light con-tribution to the overall night sky brightness have been developed,from the early works by Treanor (1973), Berry (1976) and Garstang(1986, 1989, 1991), to the more recent by Cinzano & Falchi (2012),Kocifaj (2007, 2018) and Aubé et al. (2005, 2018). The outcomeof these models are maps of the artificial brightness of the nightsky (for instance Cinzano et al. (2001) and Falchi et al. (2016)) towhich the natural brightness shall be added in order to obtain thetotal one, which is the physical observable available from directmeasurements.The publication of the
Gaia -DR2 archive (Gaia Collaborationet al. 2018), with high quality space-based global astrometric andphotometric information for more than 1.6 billion stars in all the skyup to 𝐺 =
21 mag (being 𝐺 the white light magnitude observed with Gaia in the astrometric instrument), offers a unique opportunity toimprove the computation of the integrated star light, by directlyintegrating the contributions of all the stars in the
Gaia catalogue.As a first result, we present in this work the first maps of the extra-atmospheric brightness of the night sky in the 𝑉 (Johnson) and thehuman vision photopic and scotopic bands, directly obtained fromthe Gaia and
Hipparcos (ESA 1997; van Leeuwen 2007) catalogues.
Hipparcos is used here to extend
Gaia catalogue to the bright end,not included in
Gaia -DR2. The inclusion of the remaining relevantsources, including zodiacal light and airglow, and the correctionsdue to attenuation and scattering into the field of view, allows tocompute the radiance at ground-level at any direction of the sky andany location and night time for the observer. This methodology canbe applied to any photometric band, if appropriate transformationsfrom the native
Gaia bands are computed. In particular, we havedetermined the sky brightness for the Johnson 𝑉 and human-visionphotopic and scotopic passbands. All together constitutes the firstversion of GAMBONS model of the natural night sky brightness.The novelty of this work is the determination of the radianceoutside the Earth’s atmosphere using Gaia and
Hipparcos cata-logues. However, a complete sky brightness model requires themodeling of the other components, as zodiacal light, airglow oratmospheric attenuation and scattering. We have implemented inGAMBONS the more updated and reliable models available in theliterature for those components, as described in the next sections,but it is not the scope of this work to discuss in detail each one ofthem. On the other hand, the final map of the night sky brightness ishighly dependent on the atmospheric conditions, including airglow.Thus, the accuracy of the results depends on the reliability of theatmospheric parameters used in the computation.The
Gaia and
Hipparcos catalogues are over-viewed in Sec-tion 2. Their photometric systems, as well as the other photometricbands used in this work are described in Section 3. Section 3 alsodeals with the problem of deriving the needed transformations be-tween different photometric bands. Section 4 describes the generalcharacteristics of our model. Sections 5 to 7 explain in detail thecomputation of the different sources of the natural night sky bright-ness and the models for atmospheric attenuation and scattering.Finally, we expose some applications of the model in Section 8.
Gaia (Gaia Collaboration et al. 2016) is an extremely ambitiousastrometric space mission of the scientific programme of the Euro-pean Space Agency (ESA).
Gaia measures, with very high accuracy, the positions, motions and parallaxes of a large number of stars andmany other kinds of objects as galaxies, asteroids and quasars.Consequently, a detailed three-dimensional map of more than1.6 billion stars of our Galaxy (approximately 1% of the stars popu-lating the Milky Way) is obtained. This map also includes, for almostall the objects, information about their brightness and colour, as wellas radial velocity, and several astrophysical parameters for a largefraction of them.
Gaia -DR2 photometry in 𝐺 band is almost complete up to 𝐺 ≈
20 mag, with a limiting magnitude slightly fainter than 𝐺 =
21 mag in some areas of the sky.
Gaia extends to the faintthe catalogue obtained by the
Hipparcos
ESA mission, launchedin 1989.
Hipparcos catalogue includes information about the po-sition, motion, brightness and colour of the sources. Although the
Gaia catalogue is several orders of magnitude more complete andprecise than
Hipparcos catalogue,
Gaia cannot observe very brightobjects (with
𝐺 <
Hipparcos catalogueaddresses the lack of
Gaia data. In particular, of the 118,218 starsin the
Hipparcos catalogue (ESA 1997), 83,034 are present in the
Gaia -DR2 -
Hipparcos best neighbour catalogue available on the
Gaia archive, and has therefore,
Gaia data available. The details ofthe
Gaia -DR2 -
Hipparcos cross match are given in Marrese et al.(2019). For the 35,000 remaining stars not having
Gaia informationwe have used photometry from the
Hipparcos catalogue. In terms offlux, the
Hipparcos stars account for around 20 per cent of the totalintegrated star light. After these considerations, our final cataloguecontains 1,692,953,899 stars.
Gaia photometric system
The
Gaia broad band photometric system is fully described inJordi et al. (2010). It is composed by three different passbands( 𝐺 , 𝐺 BP and 𝐺 RP ). The transmission curves represented in Fig. 1were derived by Maiz Apellaniz & Weiler (2018), as the most rep-resentative of the Gaia -DR2 data according to the
Gaia team . 𝐺 isa panchromatic band covering the wavelength range from about 350to 1000 nm. It is the result of the transmission of the mirrors andthe response of the instruments in Gaia . 𝐺 BP (blue photometer)and 𝐺 RP (red photometer) are two broad passbands in the wave-length ranges 330–680 nm and 640–1000 nm, respectively. Thecorresponding magnitudes are obtained from the integrated flux ofthe low-resolution spectra of the Gaia spectrophotometer. Thesephotometric data allow the classification of the sources by derivingthe astrophysical parameters, such as effective temperature, grav-ity, chemical composition and interstellar absorption. An additionaldescription of the photometric contents of
Gaia -DR2, includinguncertainty distribution, external comparisons and colour transfor-mations, can be found in Evans et al. (2018).A small systematic inconsistency in the BP photometric systemwas spotted by Maiz Apellaniz & Weiler (2018); Weiler (2018),likely caused by insufficient convergence of the 𝐺 BP calibrationin Gaia -DR2 for bright sources. Maiz Apellaniz & Weiler (2018)propose to mitigate this effect by using two different 𝐺 BP responsecurves ( 𝐺 𝑏 BP and 𝐺 𝑓 BP ) for the magnitude ranges brighter and fainterthan 𝐺 = .
87 mag, respectively. MNRAS , 1–15 (2020) ap of the natural sky brightness Figure 1.
Gaia passbands transmission used in this work (Maiz Apellaniz& Weiler 2018).
Figure 2. 𝑉 Johnson (Bessell & Murphy 2012), Hipparcos ( 𝐻 p ) (Weileret al. 2018), 𝑉 ’ scotopic (CIE 1951) and 𝑉 photopic (CIE 1990) passbandstransmission used in this work. Hipparcos catalogue provides photometric magnitudes in severalpassbands. Three of them, 𝐻 p , 𝐵 T and 𝑉 T , are the three nativeHipparcos–Tycho bands. Furthermore, the catalogue compiles ex-ternal values of the 𝑉 Johnson magnitude and ( 𝐵 − 𝑉 ) and ( 𝑉 − 𝐼 )colours in the Johnson-Cousins system (Johnson 1963).In this work we use 𝑉 and 𝐻 p magnitudes and ( 𝐵 − 𝑉 ) and( 𝑉 − 𝐼 ) colours to transform to other passbands, as described inSection 3.4. The passband transmissions used here for 𝑉 (Bessell& Murphy 2012) and 𝐻 p (Weiler et al. 2018) are shown in Fig. 2.In the case of the V band, the photonic response function definedin Bessell & Murphy (2012) was transformed back to its originalenergy response form, since the classical 𝑉 magnitude system isdefined in terms of in-band irradiances, not photon numbers. The scotopic 𝑉 (cid:48) ( 𝜆 ) and photopic 𝑉 ( 𝜆 ) spectral responses (Fig. 2)describe the response of the human eye fully adapted to low (<0.005 cd m − ) and high (> 5 cd m − ) environmental luminances,respectively. These photometric bands have been standardized by CIE (1951, 1990). The luminance of a radiance field seen withphotopic adaptation, in cd m − , is obtained by multiplying the in-tegrated radiance within the 𝑉 ( 𝜆 ) band (in W m − sr − ) by thephotopic luminous efficacy factor 𝐾 m =
683 lm W − . Similarly,the luminance of a radiance field seen under scotopic adaptation isgiven by the radiance in the 𝑉 (cid:48) ( 𝜆 ) band multiplied by the scotopicluminous efficacy factor 𝐾 (cid:48) m = − . Mesopic sensitivitybands can also be defined for intermediate adaptation states (0.005to 5 cd m − ) (CIE 2010). The luminance is the basic physical quan-tity at the root of the human perception of brightness. Note that the’visible’ 𝑉 Johnson band does not provide an accurate estimation ofthe visual brightness of the sky: skies with the same brightness in 𝑉 Johnson may correspond to very different luminances, dependingon the spectral composition of the sky radiance (Bará et al. 2020),and, of course, on the adaptation state of the observer.
The purpose of this work is to obtain the maps of the natural nightsky radiance in any photometric band, assuming that its spectral re-sponse is known. The input data available from
Gaia and
Hipparcos catalogues are not sky radiances (W m − sr − ), but the 𝐺 , 𝐺 BP , 𝐺 RP , 𝑉 magnitudes and the ( 𝐵 − 𝑉 ) and ( 𝑉 − 𝐼 ) colours of the individualstars. Therefore, we need to transform from these star magnitudesand colours in their native bands to the corresponding sky radiancein the desired target band. Let us recall that astronomical magni-tudes are a logarithmic measure of the in-band irradiances (W m − )at the entrance pupil of the observing instrument. The canonicaltransformation procedure consists then of three main steps: (i) inthe first one, the catalogue magnitudes are used to compute the irra-diances produced by each individual star in the original bands, (ii)the star irradiances in the original bands are used to estimate thestar irradiances within the target observation band, and finally (iii)the radiance 𝐿 of any patch of the sky in the target band is given bythe sum of the irradiances contributed by all stars contained withinthat patch, divided by the solid angle (in steradian) it subtends (Baráet al. 2020). This radiance is usually called surface brightness inthe astrophysical literature.In practice these steps can be performed in different ways, andin some cases they can be grouped into a single transformation. Thedetails depend on the available data. Note, for instance, that 𝑉 John-son is included in the
Hipparcos catalogue but not in
Gaia -DR2. So,we need to transform the
Gaia data from 𝐺 to 𝑉 . Furthermore, theinput provided by the Gaia -DR2 catalogue is given also in photonsper second within the instrument entrance pupil, which is a scaledversion of the actual in-band photon irradiance (photons s − m − )and hence of the irradiance (W m − ), being the details of this lasttransformation contingent on the stars’ spectra.In general, a transformation between two passbands ( 𝐴 and 𝐵 )is a function of, at least, one colour. For many cases, a polynomialexpression with only one colour 𝐶 , derived by subtracting the mag-nitude in two different passbands, is enough to obtain a good fitbetween 𝐴 and 𝐵 , in the form: 𝐴 = 𝐵 + ∑︁ 𝑖 𝑎 𝑖 𝐶 𝑖 (1)Taking into account the relation 𝐴 = 𝐴 − . ( 𝐸 A / 𝐸 ) ,between the magnitude 𝐴 and the corresponding irradiance 𝐸 A , MNRAS , 1–15 (2020)
E. Masana et al. with 𝐴 the zero point magnitude for the 𝐸 irradiance, we can geta relation for the ratio 𝐸 A / 𝐸 B from equation 1: 𝐸 A 𝐸 B = 𝑒 (cid:205) 𝑖 𝑎 𝑖 𝐶 𝑖 (2)The 𝐴 and 𝐵 constants, as well as the conversion from com-mon to natural logarithm are included in the 𝑒 𝑎 term.For instance, the transformation between 𝐺 and 𝑉 can be ex-pressed as function of the colour ( 𝐺 BP - 𝐺 RP ): 𝐸 V = 𝐸 G 𝑒 (cid:205) 𝑖 𝑎 𝑖 ( 𝐺 BP − 𝐺 RP ) 𝑖 (3)The methodology to get the coefficients 𝑎 𝑖 of the transforma-tions is the same as in Jordi et al. (2010). First, the set of stellarspectra in BaSeL-3.1 models (Westera et al. 2002) is used to com-pute the synthetic photometry in all the bands. BaSeL-3.1 coversa wide range of effective temperature, surface gravity and stellarmetallicity including all possible evolutionary stages in the stellarevolution. We also incorporate the effect of interstellar reddeningby using a wide range of possible absorption rates (from 0 to 11magnitudes at 550 nm) of the Gaia and
Hipparcos stars. FollowingJordi et al. (2010), a value equal to 0.03 mag has been assumed foreach synthetic magnitude of a Vega-like star. The synthetic pho-tometry was simulated in both photons per second within the
Gaia input pupil and in W m − , allowing the transformation of units tobe included in the fit (for instance from the radiant flux 𝐹 G in pho-tons per second within the Gaia input pupil to the correspondingirradiance 𝐸 V in W m − ; or from the irradiance 𝐸 V to the scotopicirradiance 𝐸 Sco , both in W m − ). Finally, a minimum least squaresfit is applied to get the 𝑎 𝑖 coefficients in equation 2.The conversion from the catalogue astronomical magnitude 𝑚 of each star to the irradiance it produces, 𝐸 (W m − ), is made usingthe standard definition of magnitudes, 𝐸 = 𝐸 r − . 𝑚 (4)where 𝐸 r is the reference irradiance associated with the ’zero-point’of the band ( 𝑉 or 𝐻 p for Hipparcos stars), that is, the irradiancecorresponding to 𝑚 = . 𝐸 r = ∫ ∞ 𝑆 ( 𝜆 ) 𝐸 Vega ( 𝜆 ) 𝑑𝜆 (5)where 𝐸 Vega ( 𝜆 ) is the STIS003 spectral irradiance of Vega fromBohlin & Gilliland (2004), normalized in such a way that 𝑚 Vega = .
03 mag for all the bands, and 𝑆 ( 𝜆 ) is the photometric passband.As pointed out above, the radiance 𝐿 of any region of the skyis computed as the sum of the irradiances produced by all sourcescontained in that region, divided by the region’s angular extent insteradians. This radiance can be expressed in the logarithmic scaleof magnitudes per square arcsecond within the 𝑆 ( 𝜆 ) band, 𝑚 S ,according to the definition: 𝐿 = 𝐿 r − . 𝑚 S (6)where 𝐿 r is the reference radiance given by 𝐿 r = 𝐾 ∫ ∞ 𝑆 ( 𝜆 ) 𝐸 Vega ( 𝜆 ) 𝑑𝜆 (7)being K=2.3504 10 − the steradian equivalent of 1 square arcsec-ond. The values of 𝐿 r for several bands used in this work are givenin Table 2. Figure 3.
Relationship between the irradiances at different bands used inthis work obtained using BaSeL-3.1 spectral library (Westera et al. 2002).The colour code represents the interstellar absorption (vertical scale). Theblack line is the fit of equation 2 with the coefficients in Table 1.
In principle, this transformation methodology can be appliedto any band, taking care to choose the colour that minimizes theresiduals of the fit. For
Gaia passbands this colour is ( 𝐺 BP - G RP ),while for 𝑉 , 𝐻 p and human-vision bands we choose ( 𝑉 − 𝐼 ) or( 𝐵 − 𝑉 ), depending on the transformation. The relationships betweenthe different bands used in this work are shown in Fig. 3. The 𝑎 𝑖 coefficients, the colours 𝐶 and the 𝜎 of the fits for all of them aresummarized in Table 1. We describe in this section the general expression for the in-bandintegrated radiance, 𝐿 obs ( 𝜶 , ℎ ) , detected by an observer in the di-rection 𝜶 = ( 𝑎, 𝑧 ) ( 𝑎 the azimuth and 𝑧 the zenith angle) relativeto its reference frame, at height ℎ above sea level. 𝐿 obs ( 𝜶 , ℎ ) musttake into account all the sources contributing to the natural night skybrightness: the integrated star light (ISL), the diffuse galactic light(DGL), the extragalactic background light (EBL) and the zodiacallight, that conform the radiance outside the Earth’s atmosphere; andthe atmospheric airglow. This extra-atmospheric radiance 𝐿 ( 𝜆, 𝜶 ) MNRAS , 1–15 (2020) ap of the natural sky brightness Table 1.
Coefficients and mean standard deviations of the transformations between irradiances in different bands, according to equation 2. 𝐸 (cid:48) denotes irradiancesin photon s − m − , whereas 𝐸 denotes irradiance in W m − . For the transformations from Gaia photometry, the last column shows the standard deviation forthe range ( 𝐺 BP - 𝐺 RP ) < Gaia -DR2 stars are.C 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝜎 all (%) 𝜎 Gaia (%) 𝐸 G / 𝐸 (cid:48) G 𝐺 𝑓 BP - 𝐺 RP -42.474 -0.1168 -0.006034 +0.005261 -6.496 10 − +2.521 10 − 𝐸 V / 𝐸 (cid:48) G 𝐺 𝑓 BP - 𝐺 RP -43.878 +0.028640 -0.20573 +0.017942 -5.444 10 − 𝐸 Sco / 𝐸 (cid:48) G 𝐺 𝑓 BP - 𝐺 RP -43.551 -0.201 -0.2284 +0.02308 -8.216 10 − 𝐸 Phot / 𝐸 (cid:48) G 𝐺 𝑓 BP - 𝐺 RP -43.749 +0.04745 -0.1972 +0.01678 -5.03 10 − 𝐸 G / 𝐸 V ( 𝑉 − 𝐼 ) − 𝐸 Sco / 𝐸 V ( 𝐵 − 𝑉 ) − 𝐸 Phot / 𝐸 H p ( 𝑉 − 𝐼 ) -1.0155 +0.39292 -0.12037 +0.0077708 -3.2423 10 − Table 2.
Reference radiances 𝐿 r based on STIS003 spectrum and 𝐺 = 𝑉 = 𝐻 p = .
03 mag for Vega.Passband 𝐿 r (W m − sr − ) 𝐺 𝑉 𝐻 p (formally including airglow) contributes to the overall map of thesky as seen by the observer as follows: • it produces a direct (attenuated) radiance 𝐿 d ( 𝜆, 𝜶 , ℎ ) = 𝐿 ( 𝜆, 𝜶 ) 𝑇 ( 𝜆, 𝑧, ℎ ) in the direction of observation 𝜶 . • it introduces a scattered radiance 𝐿 s ( 𝜆, 𝜶 s , 𝜶 , ℎ ) in the remain-ing pixels of the sky, including the pixel 𝜶 = 𝜶 s itself. 𝑇 ( 𝜆, 𝑧, ℎ ) is the effective atmospheric transmittance at wave-length 𝜆 , accounting for the radiance reduction due to the attenuationalong the beam path from the limits of the atmosphere to the locationof the observer.If we denote 𝐿 s ( 𝜆, 𝜶 , ℎ ) as the total scattered radiance reachingthe observer when looking at the direction 𝜶 : 𝐿 obs ( 𝜆, 𝜶 , ℎ ) = 𝐿 d ( 𝜆, 𝜶 , ℎ ) + 𝐿 s ( 𝜆, 𝜶 , ℎ ) (8) 𝐿 obs ( 𝜶 , ℎ ) = ∫ ∞ 𝑆 ( 𝜆 ) 𝐿 obs ( 𝜆, 𝜶 , ℎ ) 𝑑𝜆 (9)The effects of the atmospheric transmittance and scattering,i.e. the computation of 𝑇 ( 𝜆, 𝑧, ℎ ) and 𝐿 s ( 𝜆, 𝜶 , ℎ ) , are described indetail in Section 7. 𝐿 obs ( 𝜆, 𝜶 , ℎ ) can be expressed as the sum of the radi-ance coming from different contributors: 𝐿 obs , ISL ( 𝜆, 𝜶 , ℎ ) (in-tegrated starlight); 𝐿 obs , DGL ( 𝜆, 𝜶 , ℎ ) (diffuse galactic light); 𝐿 obs , EBL ( 𝜆, 𝜶 , ℎ ) (extragalactic background light); 𝐿 obs , zl ( 𝜆, 𝜶 , ℎ ) (zodiacal light); and 𝐿 obs , ag ( 𝜆, 𝜶 , ℎ ) (airglow): 𝐿 obs ( 𝜶 , ℎ ) = ∑︁ ∀ P (cid:20) ∫ ∞ 𝑆 ( 𝜆 ) 𝐿 obs , P ( 𝜆, 𝜶 , ℎ ) 𝑑𝜆 (cid:21) == 𝐿 obs , ISL ( 𝜶 , ℎ ) + 𝐿 obs , DGL ( 𝜶 , ℎ ) ++ 𝐿 obs , EBL ( 𝜶 , ℎ ) + 𝐿 obs , zl ( 𝜶 , ℎ ) + 𝐿 obs , ag ( 𝜶 , ℎ ) (10)In the next sections we will first describe in detail the compu-tation of the extra-atmospheric radiance for each component, andthereafter the effect of the atmospheric attenuation and scattering.Note that the different components of the radiance are usually given in different coordinate systems: equatorial ( 𝛼, 𝛿 ) or galactic ( 𝑙, 𝑏 ) coordinates for the astrophysical component; ecliptic ( Λ , 𝛽 ) for the zodiacal light; and horizontal ( 𝑎, ˆ ℎ ) coordinates for the air-glow and atmospheric attenuation, being ˆ ℎ the angular height abovethe horizon. Horizontal coordinates are linked to the local referenceframe of the observer and therefore the horizontal coordinates of anyextra-terrestrial source depend on the observer position and time.Moreover, both the airglow 𝐿 ag and the atmospheric transmittance 𝑇 depend only of the zenith angle 𝑧 = − ˆ ℎ . For the details aboutcoordinate transformation see for instance Green (1985). In the visual wavelength range, the integrated starlight (ISL) is oneof the most important contributors to the natural sky brightness.In some circumstances, as for lines of sight towards the GalacticCentre far from the Sun and for low solar activity, ISL can bethe brightest component. Early works on the night sky brightnessused simple models of the Galaxy to determine the ISL. Bahcall& Soneira (1980) constructed a model with the Galaxy consistingof an exponential disk and a power-law, spheroidal bulge, that waswidely used in the last years of the last century.Other early approach to the problem of the ISL was to usedata from imaging photopolarimeters (IPP ' s) on the Pioneer 10 and11 deep space probes (Weinberg et al. 1974). The full set of datacontains stars and background integrated light for almost all the skywith a spatial resolution around 2 ◦ in two bands, blue and red.From the early 2000s, the determination of the ISL takes ad-vantage of the emergence of large surveys, as Tycho-2, USNO-A2 or2MASS (see for instance Melchior et al. 2007). Following this line,the publication of Gaia -DR2 offers a new opportunity to determinethe ISL from the most complete photometric survey ever published.Before starting with the description of the use of the
Gaia data to establish the ISL, it is also worth mentioning the work ofDuriscoe (2013). He used more than 700 sky images taken frompristine mountaintop locations to get the contribution of the stel-lar contribution plus diffuse galactic and extragalactic light, afterremoving the contribution from airglow and zodiacal light.
Gaia -DR2 and
HipparcosWith the transformations described in Section 3.4, it is possibleto calculate the irradiance produced at the top of the atmospherein a given passband by each and every star in our catalogues, andthen to compute the sky map of the integrated starlight by first (a)adding these irradiances within small elementary patches of the sky,
MNRAS000
MNRAS000 , 1–15 (2020)
E. Masana et al. and (b) converting them into radiances by dividing the irradiancesby the solid angle (in sr) subtended by each patch. For that, thesky is first tessellated into 𝑁 pixels by using the HEALPix scheme(Górski et al. 2005). We used a resolution equal to 8 and therefore 𝑁 = − sr), but the methodology could be applied to anydesired HEALPix resolution.For each pixel we have collected all the stars in Gaia -DR2and
Hipparcos catalogues and added its irradiances for all the con-sidered passbands. However, there are more than 300 millions offaint sources without colour information in
Gaia -DR2. In orderto not lose these stars, we have assigned to them the mean 𝐺 BP - 𝐺 RP colour of the stars in its surrounded area, allowing in this waythe application of the transformation described in Section 3.4. Thein-band irradiances 𝐸 ∗ for each star are computed using equation 2with the appropriate coefficients and colour given in Table 3.4. 𝐸 p = 𝑁 p ∑︁ 𝑘 = 𝐸 ∗ 𝑘 (11)where 𝐸 p stands for irradiance in a given band in the pixel 𝑝 , and 𝑁 p stands for the number of stars inside the pixel. This total irradianceis divided by the solid angle Δ 𝜔 subtended by the pixel to obtainthe radiance outside the atmosphere from the ISL, 𝐿 ISL , p ( 𝑙, 𝑏 ) = 𝐸 p / Δ 𝜔 , with ( 𝑙, 𝑏 ) the galactic coordinates of the center of thepixel. Following the notation introduced in Section 4, with 𝐿 ∗ 𝑘 ( 𝜆 ) the contribution to the radiance at wavelength 𝜆 of the k-th star inthe pixel, 𝐿 ISL , p ( 𝑙, 𝑏 ) can be expressed as: 𝐿 ISL , p ( 𝑙, 𝑏 ) = 𝑁 p ∑︁ 𝑘 = ∫ ∞ 𝐿 ∗ 𝑘 ( 𝜆 ) 𝑆 ( 𝜆 ) 𝑑𝜆 (12)The integral in equation 12 can not be computed explicitly, aswe have not 𝐿 ∗ 𝑘 ( 𝜆 ) for each star. But it is equal to the flux in 𝐺 from Gaia -DR2 transformed to in-band irradiances expressed in W m − ,following equation 2 , and divided by the pixel solid angle, Δ 𝜔 .Stars fainter than 𝐺 = . Gaia − 𝐷𝑅 lost flux . The modelprovides a realistic description of the stellar content of the MilkyWay, including its kinematic and dynamics, the mass distribution,the star-formation rate and evolution of different stellar populations.It has been used to simulate the stars in the Milky Way up to 𝐺 = . ( 𝑙, 𝑏 ) over the sky.These stars have been classified in two groups: bright stars (10 . <𝐺 < . faint stars (21 . < 𝐺 < . 𝐺 between 20.0 and 21.0 mag have been assigned toone group or another following a probability linear distribution thatapproximately reproduce the selection function of the Gaia -DR2catalogue. Then, the ratio between the total flux of both groupswas computed (Fig. 4). Only a few points on the galactic plane andaround the galactic center have a contribution from faint stars > 𝐿 ISL , p ( 𝑙, 𝑏 ) by adding the contribution of the very faint stars notpresent in Gaia -DR2, computed from the ratios obtained from the
Figure 4.
Ratio between flux in photons sec − in the 𝐺 band of faint stars(20 < 𝐺 < . . < 𝐺 <
21 mag) obtainedusing the Besançon Galaxy Model.
Besançon Galaxy Model model in order to improve the accuracy ofour model.
The diffuse galactic light (DGL) is the diffuse background radiationproduced by the scattering of the starlight in the dust grains presentin the interstellar space. It contributes typically between 20 and 30per cent of the total integrated light from the Milky Way (Leinertet al. 1998).DGL is very difficult to map because it is masked by the lightcoming from the unresolved stars and, for ground based observa-tions, airglow and zodiacal light. Despite all this, several works tocharacterize and even model the DGL in the optical range have beendone. A simple estimation of the DGL can be done using the rela-tion between the ISL and DGL intensities in a given line of sight.It is supported by the fact that DGL is mainly originated from theforward scattering on interstellar grains, tracking in this way thestarlight in a given sky direction. Leinert et al. (1998) give the ratiosDGL to ISL for different galactic latitudes based on Toller (1981).A second and more accurate approach is to use the data from Pio-neer probes (Arai et al. 2015). The Imaging Photopolarimeter (IPP)instruments on Pioneer 10/11 collected data in blue (395 nm – 495nm) and red (590 nm – 690 nm) bands. The measured radiances atheliocentric distances greater than 3 AU are assumed as not affectedby the zodiacal light. After removing the contribution of the ISLusing stellar counts or synthetic models, it is possible to get the DGLplus the Extragalactic Background Light (the integrated radiationfrom all light sources outside the Galaxy). A summary of the useand limitations of this methodology can be found in Toller (1990).In this work, we have adopted a different approach. From theearly work of Laureijs et al. (1987), several authors have pointed outthe relation between the diffuse emission of the dust at 100 𝜇 m andits emission at visible wavelength. The relation can be written, usingthe notation in Matsuoka et al. (2011) and Kawara et al. (2017), as: 𝐼 𝜈,𝑖 ( DGL ) = 𝑏 𝑖 𝐼 𝜈, − 𝑐 𝑖 𝐼 𝜈, (13) 𝐼 𝜈, = 𝐼 𝜈, SFD − . − (14)where 𝐼 𝜈, is the spectral radiance at 100 𝜇 m from the In-terstellar Medium (ISM), 𝑏 𝑖 and 𝑐 𝑖 are free parameters, and 𝐼 𝜈, SFD is the 100 𝜇 m spectral radiance from the diffuse emission map ofSchlegel et al. (1998) (SFD hereafter). Equation 14 accounts for MNRAS , 1–15 (2020) ap of the natural sky brightness the Extragalactic Background Light (EBL) that must be subtractedfrom the SFD map. The EBL emission at 100 𝜇 m is ≈ − (Matsuoka et al. 2011). Their optical and 100 𝜇 m emissions are notcorrelated.The negative quadratic term in equation 13 reflects the ob-served saturation in the DGL radiance when regions with high100 𝜇 m emission become optically thick (Ienaka et al. 2013). Theuse of equation 13 is therefore limited to 𝐼 𝜈, SFD <50 MJy sr − . Thisrestriction applies mainly to low galactic latitudes | 𝑏 | (cid:47) ◦ , wherehigh 100 𝜇 m emission at very optical thick regions are found. In thiscases, equation 13 could give preposterous DGL radiances severaltimes grater than ISL ones. To avoid this, we have imposed an upperlimit in the DGL to ISL radiance ratio equal to 0.35, according withthe highest values reported by Toller (1981).Our computation of DGL is based on the 100 𝜇 m spectralradiance SFD map. It is a combination of IRAS and COBE/DIRBEdata, with a resolution of few arcminute. It is the main source toderive the dust temperature, opacity and extinction. Values are givenin MJy sr − (1MJy=10 − W m − Hz − ). The zodiacal foregroundemission and bright stars have been removed from the map, but notLMC, SMC and M31 extragalactic sources.The coefficients 𝑏 𝑖 and 𝑐 𝑖 are taken from Kawara et al. (2017),who gives their values for several optical wavelengths from 0 . 𝜇 mto 0 . 𝜇 m. This range of wavelengths almost fully covers the opticalbands considered in this work, allowing the interpolation of the 𝑏 𝑖 and 𝑐 𝑖 given in Kawara et al. (2017) for the whole range and thereforethe computation of the radiance of the DGL, 𝐿 DGL ( 𝜆 ) , followingthe same steps as for the ISL. The resultant map of the DGL in the 𝑉 band is shown in Fig. 5. The Extragalactic Background Light (EBL) is a minor contributorto the sky brightness. EBL is the sum of all extragalactic sources,mainly resolved and unresolved galaxies and intergalactic matter. Itis assumed as isotropic. Due to its low intensity, its observation andmeasurement is strongly disturbed by the zodiacal light and airglow.In spite of its low contribution, we have included ELB in our model,based on the data from Driver et al. (2016). In that paper, the EBLis derived for a wide range of wavelengths, including the
𝑈𝐵𝑉 𝐼 bands, from a combination of wide and deep galaxy number-countdata from the Galaxy And Mass Assembly, COSMOS/G10, HubbleSpace Telescope (HST) Early Release Science, HST UVUDF, andvarious near-, mid-, and far-IR data sets from ESO, Spitzer, andHerschel. We used the data of table 2 of the paper to interpolate atany wavelengths in the optical range. With this procedure, we obtaina radiance of the EBL at the 𝑉 band equal to 1.1 nW m − sr − .The final map of the radiance in 𝑉 Johnson passband outsidethe atmosphere including the integrated starlight, the diffuse galacticlight and the extragalactic background light is shown in Fig. 6. Thisdata, together with the radiances for the 𝐺 , scotopic and photopicpassbands, is available online. A sample of the data is shown inTable 3. Zodiacal light is originated by the scattering of the Sun light in thedust particles near the ecliptic plane. It represents a significant termof the natural night sky brightness. The zodiacal light decreases
Figure 5.
Diffuse Galactic Light radiance in the 𝑉 Johnson passband. with the angular distance to the Sun, with the exception of the anti-solar point, where we find the gegenschein contribution caused bythe backward scattering of the solar light. It also decreases withthe heliocentric distance. The values for the optical emission (at 𝜆 =
500 nm) of the zodiacal light at 1 A.U. were compiled byLevasseur-Regourd & Dumont (1980) and updated by Leinert et al.(1998). Kwon et al. (2004) presents a new set of data, with a smallersky coverage near the Sun than the previous work. Both sets of datahave a high degree of agreement and offer confident values for thezodiacal light brightness. In this work we have used the data givenin Leinert et al. (1998).The spectrum of the zodiacal light is slightly reddened withrespect to the solar spectrum, and therefore a color correction mustbe taken in to account when computing its radiance in other bandsdifferent to 𝑉 . Leinert et al. (1998) accounts for the effect of thecolor in the zodiacal light through the factor 𝑓 co , a measure of thequotient of the zodiacal light and the solar radiance, normalized at 𝜆 =
500 nm, function of the wavelength and elongation ( 𝜖 ): 𝑓 co ( 𝜆, 𝜖 ) = 𝐿 zl ( 𝜆 )/ 𝐿 (cid:12) ( 𝜆 ) 𝐿 zl ( )/ 𝐿 (cid:12) ( ) (15)From the values of 𝑓 co ( 𝜆, 𝜖 ) in Leinert et al. (1998), the zodia-cal light spectrum at 1 A.U. can be obtained from the solar spectrum: 𝐿 zl ( 𝜆, Λ , 𝛽 ) = 𝑓 co ( 𝜆, 𝜖 ) 𝐿 zl ( ) 𝐸 (cid:12) ( 𝜆 ) 𝐸 (cid:12) ( ) (16)where 𝐸 (cid:12) is the STIS002 spectral irradiance of the Sun from theCALSPEC library (Bohlin et al. 2014), and the elongation can becomputed from the ecliptic coordinates ( Λ , 𝛽 ) .As in the previous sections, the in-band radiance 𝐿 zl ( Λ , 𝛽 ) isobtained by integrating 𝐿 zl ( 𝜆, Λ , 𝛽 ) multiplied by the photometricpassband transmission 𝑆 ( 𝜆 ) .The effect of the variation of the visual brightness of the zodi-acal light with the heliocentric distance of the Earth 𝑟 (in A.U.) ismodelled (see Leinert et al. (1980)) by the factor 𝑓 R = 𝑟 − . .Finally, the influence of the Earth position relative to the planeof the interplanetary dust cloud is also taken into account for highecliptic latitudes ( | 𝛽 | ≥ ◦ ). It introduces a sinusoidal variationof ±
10 per cent in the zodiacal light brightness, having the mostextreme values when Earth is above or below this plane, and meanvalues at the nodes, placed at
Ω = ◦ (Leinert et al. 1998). Thisfactor 𝑓 S takes the form: 𝑓 S ( 𝛽 ) = (cid:40) + . ( Λ E − Ω ) , if | 𝛽 | ≥ ◦ , otherwise (17) MNRAS000
Ω = ◦ (Leinert et al. 1998). Thisfactor 𝑓 S takes the form: 𝑓 S ( 𝛽 ) = (cid:40) + . ( Λ E − Ω ) , if | 𝛽 | ≥ ◦ , otherwise (17) MNRAS000 , 1–15 (2020)
E. Masana et al.
Figure 6.
Sky map of the radiance outside the Earth atmosphere (integrated star light, diffuse galactic light and extragalactic background light) in the 𝑉 Johnsonpassband.
Table 3.
Radiances (W m − sr − ) outside the Earth atmosphere including the integrated starlight, the diffuse galactic light and the extragalactic backgroundlight for the 𝑉 Johnson, 𝐺 , scotopic and photopic passbands. Only the first rows are showed. The full table is available online.HEALPix Id 𝑙 gal ( ◦ ) 𝑏 gal ( ◦ ) 𝐿 V 𝐿 G 𝐿 Sco 𝐿 Phot − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − where Λ E is the ecliptic longitude of the Earth.The in-band zodiacal radiance outside the Earth atmosphere ata point of ecliptic coordinates ( Λ , 𝛽 ) can then be expressed by: 𝐿 zl ( Λ , 𝛽 ) = 𝑓 R 𝑓 S ( 𝛽 ) 𝐿 zl ( Λ , 𝛽 ) (18) Aiglow is a faint light emission originated in the upper atmosphere.It is caused by chemiluminescence, i.e. the emission of light (lu-minescence) from the decay of excited states of the products of achemical reaction. In the case of the airglow, the chemiluminescenceis triggered by the high-energy solar radiation. Airglow is emittedfrom several high altitude atmosphere layers, starting around 90 km(mesopause), where the bright OI and OH emissions, together withfainter O and NaD emissions, concentrate. Between ≈
250 km and ≈
300 km we found the emission of several OI transitions. Theairglow spectrum in the visible wavelengths is dominated by the OIgreen line at 558 nm (see Hart (2019a)), the OI red lines at 630nm and 636 nm (both produced in the layer between 200 km and300) and the FeO pseudo-continuum around 590 nm (Saran et al.(2011), Unterguggenberger et al. (2017)). Finally, the hydrogen atgeocorona contributes to the airglow in the 𝐻 𝛼 line. Other lines, like 𝐿 𝛼 and 𝐿 𝛽 are also generated in the geocorona by fluorescence, butthey fall out our spectral range of interest. For a detailed analysis ofthe complex processes involved in the airglow generation see Hart(2019a) or Noll et al. (2012).Airglow is a highly variable source of the natural night skybrightness. The variability of the emissions due to the different con-stituents of the airglow (atomic and molecular O , Na, OH, . . . ) atdifferent time scales (nightly, yearly, long term, . . . ) is thoroughlyanalyzed by Hart (2019a,b). For some of the main components itcan reach 100 per cent (maximum minus minimum divided by themedian). Also, airglow depends on the solar activity cycle (Patat2008). Finally, the airglow emission could change with the geo-graphic latitude, specially for lines originating in the mesopause(OH, O , Na I D, FeO or most of the green OI lines), and also withgeomagnetic latitude, in particular for ionospheric lines such as theOI lines at 630 and 636 nm. This makes the prediction of the actualvalue of the airglow emission difficult to model, and it could beconsidered as a free parameter if comparisons with real measure-ments are done. The basic model described here is only intended toprovide a reference value and it shall be judiciously applied to anyparticular observation.The most common model for the airglow brightness, i.e. the in-band integrated airglow radiance 𝐿 ag ( 𝑧 ) at different zenith angles 𝑧 , MNRAS , 1–15 (2020) ap of the natural sky brightness is the one based on the factorization of the zenith brightness 𝐿 ag ( ) and the zenith angle dependence given by the van Rhijn function(Leinert et al. 1998): 𝐿 ag ( 𝑧 ) = 𝐿 ag ( ) (cid:16) − (cid:2) 𝑅𝑅 + 𝐻 (cid:3) 𝑠𝑖𝑛 𝑧 (cid:17) (19)where 𝑅 = 𝐻 is the height above theEarth’s surface of an assumed thin homogeneously emitting layerresponsible for airglow. Different values have been quoted in theliterature for 𝐻 . Here we adopt 𝐻 =
87 km (Hart 2019a).In general,the choice of 𝐻 is only crucial at high zenith angles, but we needto bear in mind that, as mentioned above, different emissions areoriginating at different altitudes, and consequently, our choice of 𝐻 could be inaccurate if strong emissions are originating at highaltitude. For more accurate results, in particular in case of strongline emissions at different heights, a more complex airglow model,beyond the van Rhijn approximation with fixed H, should be used.Some remarks should be made before applying this equationto our calculations. First, note that equation 19 is expressed in termsof in-band integrated radiances, 𝐿 ag ( 𝑧 ) , not of spectral radiances 𝐿 ag ( 𝜆, 𝑧 ) , although it could be deemed approximately valid (withthe same values for 𝑅 and 𝐻 ) for 𝐿 ag ( 𝜆, 𝑧 ) . Second, the van Rhijnformulation is meant to describe the at-the-layer radiance, 𝐿 ag ( 𝑧 ) ,as seen at zenith angles 𝑧 from the reference frame of an observerlocated at sea level ( ℎ = ℎ > 𝑧 from sea level will appear to be at somewhat larger angles 𝑧 (cid:48) = 𝑧 (cid:48) ( 𝑧, ℎ ) , if seen by an observer at ℎ > ℎ is much smaller than 𝑅 and 𝐻 , thissmall correction can be ignored, and equation 19 can be used to de-scribe the at-the-layer radiance in their particular reference frames,irrespective from ℎ .The nadir-oriented spectral radiance at the layer, 𝐿 ag ( 𝜆, ) ,can be determined by observational measurements or from syn-thetic models. In the case of observational determinations, as inHart (2019a), some other sky radiance components can contributeto some extent to the continuum value. For calculating the resultspresented in Section 8 of this work we have used the Cerro ParanalAdvanced Sky Model (Noll et al. 2012 and Jones et al. 2013) toget a synthetic airglow line and continuum emission spectrum. Thespectrum, computed from ESO’s SkyCalc web interface, is calcu-lated for the Cerro Paranal altitude (2640 m above the sea level)in the wavelength range from 350 nm to 1050 nm. For a referencespectrum, we set the value of the Monthly Averaged Solar RadioFlux equal to 100 sfu (1 sfu = 10 − W m − s − ), the approximateaverage value in the period 2009-2020 (one solar cycle) according tothe data of the Canadian Space Weather Forecast Centre (CSWFC).This value can be set to other value for any specific application ofthe model. The out-coming spectrum is shown in Fig. 7. After cor-recting for the vertical atmospheric transmittance at Cerro ParanalObservatory, 𝑇 CPO ( 𝜆, ℎ ) , also provided by the ESO’s SkyCalctool, we can determine the airglow radiance at the emission layer: 𝐿 layer ( 𝜆, ) = 𝐿 CPOag ( 𝜆, ℎ )/ 𝑇 CPO ( 𝜆, ℎ ) (20)As stated above, the final spectral airglow radiance is highlyspatially and temporally variable and usually it should be consideredas a free parameter. Its value could be set if more information about Figure 7.
The airglow spectrum computed with the ESO’s SkyCalc webinterface for an altitude of 2640 m. Logarithm binning 𝜆 / Δ 𝜆 = the airglow emission is available. This issue will be discussed inSection 8.1. In the previous sections we have described the contributions of theastrophysical sources to the night sky brightness. Although slightlydifferent models could be adopted for the background light or thezodiacal light, it is expected that the actual radiance outside theEarth’s atmosphere does not differ significantly from the valuesgiven in this work. We have also discussed the contribution of theairglow, a highly variable source that should be considered as a freeparameter of the model.However, in order to obtain the local map of the night skybrightness, we need to deal with the effects of the Earth atmosphereon the outside radiance. In particular, the effect of the atmosphericattenuation and the scattering must to be added to the radiance com-ing from the astrophysical sources. Both effects are highly variablein time, since they strongly depend on the particular atmosphericconditions at the place and time of the observation. For the numer-ical examples presented in Section 8 of this work we adopted somestandard values for the main atmospheric parameters. Different at-mospheric states will give rise to different results.The spectral radiance is modified by the terrestrial atmosphereby means of two interrelated processes. On the one hand, the radi-ance propagating toward the observer is attenuated by absorptionand scattering by the atmospheric constituents. On the other hand,some amount of light initially propagating in different directionsgets scattered into the observer’s line of sight, being added to therecorded brightness. These two opposite effects stem from the samebasic interaction processes at the atomic and molecular levels.The attenuation of a beam along an atmospheric path of length 𝑑 can be expressed in terms of the optical thickness 𝜏 ( 𝑑 ) or, equiv-alently, of the transmittance 𝑇 ( 𝑑 ) , as: 𝐿 ( 𝑑 ) = 𝐿 ( ) 𝑒 − 𝜏 ( 𝑑 ) = 𝐿 ( ) 𝑇 ( 𝑑 ) (21)where 𝐿 ( ) is the initial radiance and 𝐿 ( 𝑑 ) the radiance aftertravelling the distance 𝑑 .Henceforth we closely follow the formulation by Kocifaj(2007). Let us denote by 𝑘 ext ( ℎ (cid:48) ) the volume extinction coefficientat height ℎ (cid:48) above sea level, due to aerosols (A) and molecules (M),such that: 𝑘 ext ( ℎ (cid:48) ) = 𝑘 𝑀 ext ( ℎ (cid:48) ) + 𝑘 𝐴 ext ( ℎ (cid:48) ) (22) MNRAS000
The airglow spectrum computed with the ESO’s SkyCalc webinterface for an altitude of 2640 m. Logarithm binning 𝜆 / Δ 𝜆 = the airglow emission is available. This issue will be discussed inSection 8.1. In the previous sections we have described the contributions of theastrophysical sources to the night sky brightness. Although slightlydifferent models could be adopted for the background light or thezodiacal light, it is expected that the actual radiance outside theEarth’s atmosphere does not differ significantly from the valuesgiven in this work. We have also discussed the contribution of theairglow, a highly variable source that should be considered as a freeparameter of the model.However, in order to obtain the local map of the night skybrightness, we need to deal with the effects of the Earth atmosphereon the outside radiance. In particular, the effect of the atmosphericattenuation and the scattering must to be added to the radiance com-ing from the astrophysical sources. Both effects are highly variablein time, since they strongly depend on the particular atmosphericconditions at the place and time of the observation. For the numer-ical examples presented in Section 8 of this work we adopted somestandard values for the main atmospheric parameters. Different at-mospheric states will give rise to different results.The spectral radiance is modified by the terrestrial atmosphereby means of two interrelated processes. On the one hand, the radi-ance propagating toward the observer is attenuated by absorptionand scattering by the atmospheric constituents. On the other hand,some amount of light initially propagating in different directionsgets scattered into the observer’s line of sight, being added to therecorded brightness. These two opposite effects stem from the samebasic interaction processes at the atomic and molecular levels.The attenuation of a beam along an atmospheric path of length 𝑑 can be expressed in terms of the optical thickness 𝜏 ( 𝑑 ) or, equiv-alently, of the transmittance 𝑇 ( 𝑑 ) , as: 𝐿 ( 𝑑 ) = 𝐿 ( ) 𝑒 − 𝜏 ( 𝑑 ) = 𝐿 ( ) 𝑇 ( 𝑑 ) (21)where 𝐿 ( ) is the initial radiance and 𝐿 ( 𝑑 ) the radiance aftertravelling the distance 𝑑 .Henceforth we closely follow the formulation by Kocifaj(2007). Let us denote by 𝑘 ext ( ℎ (cid:48) ) the volume extinction coefficientat height ℎ (cid:48) above sea level, due to aerosols (A) and molecules (M),such that: 𝑘 ext ( ℎ (cid:48) ) = 𝑘 𝑀 ext ( ℎ (cid:48) ) + 𝑘 𝐴 ext ( ℎ (cid:48) ) (22) MNRAS000 , 1–15 (2020) E. Masana et al. (which of course are wavelength-dependent). For a layered at-mosphere as implicitly assumed in equation 21, the optical thicknessappearing in equation 22 can be computed as: 𝜏 ( 𝑑 ) = 𝜏 ( ℎ , ℎ ) = ∫ ℎ ℎ 𝑘 ext ( ℎ (cid:48) ) 𝑑ℎ (cid:48) (23)where 𝑑ℎ (cid:48) is a function of the trajectory followed by the lightrays. Under the assumption of nearly-rectilinear propagation be-tween atmospheric layers, at an angle 𝑧 with respect to the localzenith, we can write: 𝑑ℎ (cid:48) = 𝑀 𝑀 ( 𝑧 ) 𝑑ℎ , and 𝑑ℎ (cid:48) = 𝑀 𝐴 ( 𝑧 ) 𝑑ℎ ,where 𝑀 𝑀 ( 𝑧 ) and 𝑀 𝐴 ( 𝑧 ) are the molecular and aerosol airmasses,respectively, corresponding to 𝑧 , and 𝑑ℎ is measured along the ver-tical. Then: 𝜏 ( ℎ , ℎ ) = 𝑀 𝑀 ( 𝑧 ) ∫ ℎ ℎ 𝑘 𝑀 ext ( ℎ ) 𝑑ℎ + 𝑀 𝐴 ( 𝑧 ) ∫ ℎ ℎ 𝑘 𝐴 ext ( ℎ ) 𝑑ℎ (24)Under the assumption of exponential number density profilesfor molecules and aerosols, the extinction coefficients take the form: 𝑘 𝑀 ext ( ℎ ) = 𝜏 𝑀 ℎ m 𝑒 − ℎ / ℎ m ; 𝑘 𝐴 ext ( ℎ ) = 𝜏 𝐴 ℎ a 𝑒 − ℎ / ℎ a (25)where 𝜏 𝑀 and 𝜏 𝐴 are the (wavelength-dependent) molecular andaerosol optical thicknesses of the whole atmosphere, measuredalong a vertical path ( 𝑧 =0), for which 𝑀 𝑀 ( ) = 𝑀 𝐴 ( ) =
1, and ℎ m and ℎ a are the molecular and aerosol scale heights, respectively.For our present calculations we have taken ℎ m = 8 km and ℎ a =1.54km. The molecular - Rayleigh component 𝜏 𝑀 ( 𝜆 ) can be approxi-mated analytically (Teillet 1990), with 𝜆 in 𝜇 m, by: 𝜏 𝑀 ( 𝜆 ) = . 𝜆 − . (26)And the aerosol component by: 𝜏 𝐴 ( 𝜆 ) = 𝜏 𝐴 ( 𝜆 ) (cid:16) 𝜆𝜆 (cid:17) − 𝛼 (27)being 𝛼 the Ångstrom exponent.Data on atmospheric aerosols (aerosol optical depth andÅngstrom exponent) can be obtained from the AERONET network https://aeronet.gsfc.nasa.gov . For the results presented be-low, we chose 𝜏 𝐴 (
550 nm ) = . 𝛼 = 1. The choice of other val-ues will affect the resulting night sky brightness, especially at highzenith angles, where the effect of atmospheric attenuation is moresignificant. For a highly accurate modelling, it is recommended touse the available local values.The airmasses, at a first approximation, do not depend onwavelength. Then, the spectral attenuation across a path at zenithangle 𝑧 , between the layers at ℎ and ℎ is: 𝑇 ( 𝑧, 𝜆 ; ℎ , ℎ ) = exp (cid:110) − 𝑀 𝑀 ( 𝑧 ) 𝜏 𝑀 ( 𝜆 ) [ 𝑒 − ℎ / ℎ m − 𝑒 − ℎ / ℎ m ]− 𝑀 𝐴 ( 𝑧 ) 𝜏 𝐴 ( 𝜆 ) [ 𝑒 − ℎ / ℎ a − 𝑒 − ℎ / ℎ a ] (cid:111) (28) And the overall atmospheric attenuation, from ℎ = ℎ to ℎ = ∞ is: 𝑇 ( 𝑧, 𝜆 ; ℎ ) = exp (cid:110) − 𝑀 𝑀 ( 𝑧 ) 𝜏 𝑀 ( 𝜆 ) 𝑒 − ℎ / ℎ m − 𝑀 𝐴 ( 𝑧 ) 𝜏 𝐴 ( 𝜆 ) 𝑒 − ℎ / ℎ a (cid:111) (29)For the airmass we use the expression of Kasten & Young(1989), defined for all zenithal distances as: 𝑀 ( 𝑧 ) = ( 𝑧 ) + . ( . − 𝑧 ) − . (30)with 𝑧 in degrees.Light from all other sky directions is also absorbed and scat-tered. When interacting with the atmospheric constituents locatedalong the line of sight a fraction of this light gets scattered into theobserver’s field of view, adding to the detected radiance. Denot-ing by 𝐿 ( 𝜆, 𝜶 s ) the extra-atmospheric radiance from a differentialpatch of the sky of solid angle 𝑑𝜔 around the generic direction 𝜶 s ,the total radiance scattered into the observer line of sight, 𝐿 s ( 𝜆, 𝜶 , ℎ ) can be calculated as 𝐿 s ( 𝜆, 𝜶 , ℎ ) = ∫ Ω Ψ ( 𝜆, 𝜶 , 𝜶 s , ℎ ) 𝐿 ( 𝜆, 𝜶 s ) 𝑑𝜔 (31)where Ψ ( 𝜆, 𝜶 , 𝜶 s , ℎ ) is the function describing the spectral radiancescattered toward the observer along the direction 𝜶 , per unit 𝑑𝜔 ,due to a unit amplitude radiant source located at 𝜶 s (all directionsmeasured in the observer reference frame). The integral is extendedto the whole hemisphere above the observer, Ω . For this work wehave used the Kocifaj-Kránicz Ψ ( 𝜆, 𝜶 , 𝜶 s , ℎ ) described in equa-tion (18) of Kocifaj & Kránicz (2011) with an effective scatteringphase function composed of aerosol and molecular (Rayleigh) com-ponents weighted by the corresponding optical depths. The aerosolscattering phase function is described by a Henyey-Greenstein func-tion with asymmetry parameter 𝑔 = Gaia -Hipparcos map,requires performing an integration over the whole celestial hemi-sphere for every direction of observation, 𝜶 . Several simpler butless accurate approaches have been proposed in the literature toaccount for the radiance scattered into the line of sight without cal-culating this integral. One of them is based on replacing 𝜏 𝑀 / 𝐴 ( 𝜆 ) in equation 29 by an effective optical depth 𝜏 𝑀 / 𝐴 , eff ( 𝜆 ) = 𝛾 𝜏 𝑀 / 𝐴 ( 𝜆 ) with 𝛾 <
1. The value of 𝛾 depends on the aerosol albedo andasymmetry parameter, and typical values are in the range 0.5-0.9(Hong et al. 1998). This is the approach used for diffuse sources inDuriscoe (2013), with 𝛾 = .
75, based on the empirical results ofKwon (1989). It can be a practical option for obtaining reasonablyaccurate results when computing time is a constraint, especially ifthe night sky brightness is evaluated in sky pixels of sufficient sizeto spatially average the contributions of the brightest stars.
In this section we outline some applications of the GAMBONSmodel. The examples presented here have been calculated for par-
MNRAS , 1–15 (2020) ap of the natural sky brightness Figure 8.
Comparison of GAMBONS maps (left) and SQC images (right) in 𝑉 Johnson passband for Senyús ( 𝜆 = ◦ (cid:48) ” E, 𝜙 = ◦ (cid:48) ” N, h=1142 m,Catalonia, Spain). The colour scale represents the sky brightness in Johnson V mag arcsec − . Bright patches near the horizon in the SQC images are due tolight pollution. Top: 2018 August 10, 21:46 UT. Bottom: 2019 March 23, 20:28 UT. ticular choices of the atmospheric parameters (e.g. aerosol opticaldepth) and airglow zenith radiance. Quantitative comparisons withfield measurements must take into account the actual state of theatmosphere and the airglow at the precise time of taking them.All-sky maps like the ones described below can be computedfor any location and time and freely downloaded from the GAM-BONS website ( http://gambons.fqa.ub.edu ). As a first application, the GAMBONS model has been used to ob-tain all-sky maps of the brightness of the natural night sky. Thesemaps provide a realistic description of what one can expect to ob-serve in pristine locations free from light pollution sources, underdifferent atmospheric conditions. The multi-band capability of theGAMBONS approach allows to generate all-sky images in severalbands.The maps presented in this section were calculated for differentphotometric bands and are displayed in different units. In lightpollution studies it is still frequent to work in the classical Johnson 𝑉 band, reporting the brightness (i.e., the in-band radiances) inthe logarithmic scale of magnitudes per square arcsecond. Let usremind that a region of the sky is said to have a brightness of 𝑚 V magnitudes per square arcsecond if each square arcsecond ofits visual field gives rise -at the entrance pupil of the observinginstrument- to the same irradiance a star of magnitude 𝑚 = 𝑚 V would produce Bará et al. (2020). In this section we use the 𝑉 bandas defined in Section 3.2. The magnitudes per square arcsecond arethen calculated from the in-band radiances 𝐿 (in W m − sr − ) as 𝑚 V = − . 𝑙𝑜𝑔 ( 𝐿 / 𝐿 r ) (32)with 𝐿 r = . − sr − as indicated in Table 2.For human visual observations, the quantity of choice is theluminance, measured in SI units cd m − , corresponding to the pho-topic, mesopic, or scotopic adaptation states of the eye.The all-sky images captured with some specific devices canbe compared with the GAMBONS maps. One of these devices isthe Sky Quality Camera (SQC) from Euromix Ltd. (Slovenia), acommercial DSLR (Digital Single-Lens Reflex) camera with fish-eye lens to evaluate the night sky brightness in the whole sky. Thiskind of devices are a versatile solution to measure the night skybrightness (see Hänel et al. (2018)), providing directional informa-tion. The SQC software allows processing raw all-sky images andprovides different products, including calibrated data in the 𝑉 band(Jechow et al. 2018) or (Vandersteen et al. 2020)). In Fig. 8 we showthe comparison of two of those images with our model. The images MNRAS000
Comparison of GAMBONS maps (left) and SQC images (right) in 𝑉 Johnson passband for Senyús ( 𝜆 = ◦ (cid:48) ” E, 𝜙 = ◦ (cid:48) ” N, h=1142 m,Catalonia, Spain). The colour scale represents the sky brightness in Johnson V mag arcsec − . Bright patches near the horizon in the SQC images are due tolight pollution. Top: 2018 August 10, 21:46 UT. Bottom: 2019 March 23, 20:28 UT. ticular choices of the atmospheric parameters (e.g. aerosol opticaldepth) and airglow zenith radiance. Quantitative comparisons withfield measurements must take into account the actual state of theatmosphere and the airglow at the precise time of taking them.All-sky maps like the ones described below can be computedfor any location and time and freely downloaded from the GAM-BONS website ( http://gambons.fqa.ub.edu ). As a first application, the GAMBONS model has been used to ob-tain all-sky maps of the brightness of the natural night sky. Thesemaps provide a realistic description of what one can expect to ob-serve in pristine locations free from light pollution sources, underdifferent atmospheric conditions. The multi-band capability of theGAMBONS approach allows to generate all-sky images in severalbands.The maps presented in this section were calculated for differentphotometric bands and are displayed in different units. In lightpollution studies it is still frequent to work in the classical Johnson 𝑉 band, reporting the brightness (i.e., the in-band radiances) inthe logarithmic scale of magnitudes per square arcsecond. Let usremind that a region of the sky is said to have a brightness of 𝑚 V magnitudes per square arcsecond if each square arcsecond ofits visual field gives rise -at the entrance pupil of the observinginstrument- to the same irradiance a star of magnitude 𝑚 = 𝑚 V would produce Bará et al. (2020). In this section we use the 𝑉 bandas defined in Section 3.2. The magnitudes per square arcsecond arethen calculated from the in-band radiances 𝐿 (in W m − sr − ) as 𝑚 V = − . 𝑙𝑜𝑔 ( 𝐿 / 𝐿 r ) (32)with 𝐿 r = . − sr − as indicated in Table 2.For human visual observations, the quantity of choice is theluminance, measured in SI units cd m − , corresponding to the pho-topic, mesopic, or scotopic adaptation states of the eye.The all-sky images captured with some specific devices canbe compared with the GAMBONS maps. One of these devices isthe Sky Quality Camera (SQC) from Euromix Ltd. (Slovenia), acommercial DSLR (Digital Single-Lens Reflex) camera with fish-eye lens to evaluate the night sky brightness in the whole sky. Thiskind of devices are a versatile solution to measure the night skybrightness (see Hänel et al. (2018)), providing directional informa-tion. The SQC software allows processing raw all-sky images andprovides different products, including calibrated data in the 𝑉 band(Jechow et al. 2018) or (Vandersteen et al. 2020)). In Fig. 8 we showthe comparison of two of those images with our model. The images MNRAS000 , 1–15 (2020) E. Masana et al.
Figure 9.
All sky luminance map ( 𝜇 cd m − ) for observers with scotopic(top) and photopic (bottom) luminance adaptation. 2019 March 23, 20:28TU, Senyús ( 𝜆 = ◦ (cid:48) ” E, 𝜙 = ◦ (cid:48) ” N, h=1142 m, Catalonia,Spain). were taken by one of the authors (SJR) from Senyús ( 𝜆 = ◦ (cid:48) ” E, 𝜙 = ◦ (cid:48) ” N, h=1142 m), a pristine dark site with very lowlight pollution near the Montsec Protected Area (Catalonia, Spain)at two different epochs. The agreement between the SQC data andthe model is fairly good. The main features of the sky (zodiacallight, clearly visible in the image of the bottom panel in the Westdirection, or Milky way) in these SQC images are coincident withthose predicted by GAMBONS.Note that this is a qualitative comparison. As mentioned inSection 6, the radiance of the airglow is very variable both in shortand large time scales, and for many practical applications it may beconsidered as a free parameter. We have used the value of the SolarRadio Flux from the Canadian Space Weather Forecast Centre toget the airglow spectra for the both epochs, but it could not accountfor local or very short variations of the airglow radiance. The sameapplies to the atmospheric parameters related with the attenuation.After some tests, the best match is obtained with 𝜏 𝐴 = .
15 and 𝛼 =
1. Furthermore, SQC images record light pollution, and its 𝑉 band could not coincide with the 𝑉 Johnson band used in thiswork, as it is derived from the RGB image taken with a DSLRcamera. This could introduce some differences between SQC andGAMBONS values, function of the colour of the sky.The GAMBONS model can also be used to determine the
Table 4.
Mean contributions to zenith natural night sky brightness for alatitude equal to 40 ° N.Component Mean radiances Percentage(nW m − sr − ) (%)Airglow 165 47.7Zodiacal Light 95 27.5Integrated Star Light 75 21.7Diffuse Galactic Light 10 2.9Extragalactic Background Light 0.8 ≈ amount of artificial light polluting the night sky. As in the previouscase, it requires the fitting of the parameters not available from directmeasurements, as e.g. the airglow, and in some cases the aerosolproperties, if they are unknown. The fit could be laborious due tothe complex spatial features of the actual airglow and attenuation.The model can be easily extended to other photometric bands usedin light pollution studies, as e.g. the SQM and TESS-W (Bará et al.2019) or RGB (Sánchez de Miguel et al. 2019; Kolláth et al. 2020).Both topics will be addressed in forthcoming papers. By way ofexample, in Fig. 9 we show the visual scotopic and photopic allsky maps for the same location and date. They are expressed inluminance units 𝜇 cd m − . The model allows to compute the expected contribution of the dif-ferent sources (ISL, DGL+ELB, zodiacal light and airglow) to thetotal night sky brightness (moonless and cloudless) under differentconditions (i.e. different locations, times or airglow intensities). Thecontribution also depends on the extent of the field of view that weconsider. We refer here as zenith values to the average radiance ina circular region of radius 10 degrees around the zenith with noweights applied to the different points.We have run the model for a mid latitude location ( 𝜙 = ◦ )and the atmospheric conditions given in Section 7 and the reference(constant) airglow spectrum described in Section 6, and averagedthe contribution of the different sources over a year. The modelwas run in one hour intervals during the astronomical night (i.e.Sun more than 18 degrees below the horizon). The result for the 𝑉 band is shown in Table 4. Note that these values are sensitive tothe airglow intensity, which is highly variable.The contributions also depend on the observer latitude, asshown in Fig. 10. Zodiacal light has a larger impact in locationsclose to the Equator when we observe near the zenith, as the planeof the ecliptic reaches greater heights above the horizon at theselatitudes.Finally, the radiances are also a function of the observationtime. The presence or absence of the Milky Way and the height ofthe ecliptic above the horizon are the two main factors that modulatethe relative contributions to the sky brightness throughout the year.As shown in Fig. 11, for the midnight radiance around the zenithfor an observer at mid latitudes, the Milky Way could become thebrightness component. As it can be observed in both figures, weassume a constant airglow radiance, regardless of the location andtime. For a given location, the natural night sky brightness is highlydependent on the observing time. As we have already mentioned, the
MNRAS , 1–15 (2020) ap of the natural sky brightness Figure 10.
Annual average radiance at midnight in the Johnson 𝑉 band atzenith for the different contributors to the natural sky brightness as functionof the observer latitude. Figure 11.
Midnight radiance at zenith in the Johnson 𝑉 band for thedifferent contributors to the natural sky brightness as function of the day ofthe year for an observer at latitude equal to 40 ° N. presence or absence of the Milky Way or the altitude of the eclipticplane above the horizon strongly determine the brightness of the sky.We must bear this in mind when trying to characterize the darknessof a given place. According to our model, for zenith measurements,the natural variation of the night sky brightness along the year dueto the variation of the astrophysical contributions can reach morethan 0.6 magnitudes. As an example, Fig. 12 shows the variationof the 𝑉 magnitude around the zenith (i.e. a 10 degrees circularregion around it) for an observer at 𝜙 = ° latitude and h=1000m.a.s.l. The figure plots only the times when the height of the Sunabove the horizon is < − ° . The presence of the Milky Way intwo different epochs of the year is clearly visible. The zodiacal lightis also modulated by the position of the Earth above the eclipticplane and the Sun–Earth distance. The variation is smoother if awider area of the sky is considered. But even for relatively largefields of view of tens of degrees, the variation of the natural nightsky brightness is up to several tenths of mag V arcsec − . Again, weconsider constant airglow. Figure 12.
Variation of the natural night sky brightness (mag V arcsec − atzenith) along the year for an observer at 𝜙 = ◦ and h=1000 m.a.s.l. Notethat the colour scale is zoomed with regard to Fig. 8. Figure 13.
Darkest value of the natural night sky brightness (in mag V arcsec − ) as a function of the latitude for three different fields of view: 10degrees around the zenith (zenith angle 𝑧 = (0, 10)); 30 degrees around thezenith ( 𝑧 = (0, 30)); and all-sky ( 𝑧 = (0, 90)). The atmospheric and airglowconditions are the ones assumed in this section. There might be a question about the darkest value the natural nightsky can attain, in absence of any source of artificial light pollution.This value is significant for light pollution studies, as it establishesan upper limit to the actual measurements of the sky darkness.Darker measurements could indicate higher atmospheric attenu-ation or sporadic lower airglow radiances, not necessarily beingindicative of less polluted (’darker’) skies. The darkest value of thezenith natural night sky depends on the observer latitude and thetime of observation, as well as on the atmospheric conditions, al-titude of the observer above sea level, and airglow intensity. Dueto attenuation, less stars are expected to be perceived from sites atlower altitudes above sea level, as well as at higher zenith angles,as recently analyzed in Cinzano & Falchi (2020). For a given lati-tude, using our model (with the particular airglow and atmosphericconditions assumed in this section), we are able to determine theminimum value of the night sky brightness (i.e. the darkest sky ata given latitude). It is shown in Fig. 13. The darkest skies (aver-aged over 10 degrees around the zenith) are found at mid-latitudesin both hemispheres. The presence of the ecliptic plane near thezenith for observers near the Equator prevents them to reach thedarkest values. This effect almost disappears if we considered thefull sky dome, but it is well visible even for a wide field of view(0 ° < 𝑧 < ° ). MNRAS000
Darkest value of the natural night sky brightness (in mag V arcsec − ) as a function of the latitude for three different fields of view: 10degrees around the zenith (zenith angle 𝑧 = (0, 10)); 30 degrees around thezenith ( 𝑧 = (0, 30)); and all-sky ( 𝑧 = (0, 90)). The atmospheric and airglowconditions are the ones assumed in this section. There might be a question about the darkest value the natural nightsky can attain, in absence of any source of artificial light pollution.This value is significant for light pollution studies, as it establishesan upper limit to the actual measurements of the sky darkness.Darker measurements could indicate higher atmospheric attenu-ation or sporadic lower airglow radiances, not necessarily beingindicative of less polluted (’darker’) skies. The darkest value of thezenith natural night sky depends on the observer latitude and thetime of observation, as well as on the atmospheric conditions, al-titude of the observer above sea level, and airglow intensity. Dueto attenuation, less stars are expected to be perceived from sites atlower altitudes above sea level, as well as at higher zenith angles,as recently analyzed in Cinzano & Falchi (2020). For a given lati-tude, using our model (with the particular airglow and atmosphericconditions assumed in this section), we are able to determine theminimum value of the night sky brightness (i.e. the darkest sky ata given latitude). It is shown in Fig. 13. The darkest skies (aver-aged over 10 degrees around the zenith) are found at mid-latitudesin both hemispheres. The presence of the ecliptic plane near thezenith for observers near the Equator prevents them to reach thedarkest values. This effect almost disappears if we considered thefull sky dome, but it is well visible even for a wide field of view(0 ° < 𝑧 < ° ). MNRAS000 , 1–15 (2020) E. Masana et al.
We present GAMBONS (the GAia Map of the Brightness Of theNatural Sky) model of the natural night sky brightness. The modelconsiders each of the contributors to the sky brightness (i.e. the in-tegrated star light, the diffuse galactic light, the extragalactic back-ground light, the zodiacal light and the airglow), plus the atmo-spheric attenuation and scattering into the field-of-view. The mainnovelty with respect to previous models is the use of the
Gaia -DR2 catalogue to evaluate the integrated star light. Furthermore,the photometric data in
Gaia -DR2 ( 𝐺 , 𝐺 BP and 𝐺 RP bands) allowto transform from Gaia photon fluxes to radiances in any otherphotometric band. In this way, GAMBONS becomes a multi-bandmodel to study the natural night sky brightness.In particular, although it amounts less than the 3 per cent of thesky brightness, we have put an effort to evaluate the contribution ofthe astrophysical diffuse light (diffuse galactic light plus extragalac-tic background light), by using the relation between the emission at100 𝜇 m and the optical emission. The spectrum template of the air-glow used by GAMBONS is based on the Cerro Paranal AdvancedSky Model, through the ESO SkyCalc tool.Two highly variable inputs to the model are the absolute zenithvalue of the airglow radiance and the aerosol properties at the timeof observation. They vary from place to place and in different timescales. Information on aerosol properties at a large network of mon-itoring sites worldwide can be obtained from AERONET.GAMBONS is expected to enable multiple applications in thestudy and characterization of the natural night sky brightness. Wehave pointed out several of them. In particular we have shownthe variation of the sky darkness as function of the latitude ofthe observer and time, as well as the variation of the contributionof the different sources of natural light with these two variables.Furthermore, GAMBONS could help, if reliable determinationsof airglow intensity and atmospheric conditions are available, toestablish a reference value of the sky darkness, for any given locationin a moonless and cloudless night, against which to evaluate themeasurements reported by the observers. In a forthcoming paperwe will deal with the use of GAMBONS to remove the natural nightsky brightness from the raw all-sky images in order to quantitativelyestimate the levels of light pollution, and to put to the test thepredictions obtained from models of atmospheric propagation ofartificial light at night.GAMBONS and the data related are accessible via web at http://gambons.fqa.ub.edu . ACKNOWLEDGEMENTS
We thank Dr. M. Hart for kindly providing the airglow spectrataken at Apache Point Observatory. We thank the reviewer for theirthorough and useful comments.This work was supported by the MINECO (SpanishMinistry of Economy) through grant RTI2018-095076-B-C21(MINECO/FEDER, UE). EM and JMC acknowledges financial sup-port from the State Agency for Research of the Spanish Ministryof Science and Innovation through the “Unit of Excellence Maríade Maeztu 2020-2023” award to the Institute of Cosmos Sciences(CEX2019-000918-M). SB acknowledges Xunta de Galicia, grantED431B 2020/29. SJR thanks Cal Joanet family for allowing him totake SQC measurements from their home in the village of Senyús.
DATA AVAILABILITY
The data underlying this article are available in the article and in itsonline supplementary material.
REFERENCES
Arai T., et al., 2015, ApJ, 806, 69Aubé M., Franchomme-Fossé L., Robert-Staehler P., Houle V., 2005, in At-mospheric and Environmental Remote Sensing Data Processing and Uti-lization: Numerical Atmospheric Prediction and Environmental Moni-toring. p. 589012Aubé M., Simoneau A., Wainscoat R., Nelson L., 2018, Monthly Notices ofthe Royal Astronomical Society, 478, 1776Bahcall J. N., Soneira R. M., 1980, ApJS, 44, 73Bará S., Tapia C., Zamorano J., 2019, Sensors, 19, 1336Bará S., Aubé M., Barentine J., Zamorano J., 2020, MNRAS, 493, 2429Benn C. R., Ellison S. L., 1998, New Astron. Rev., 42, 503Berry R. L., 1976, Journal of the Royal Astronomical Society of Canada,70, 97Bessell M., Murphy S., 2012, PASP, 124, 140Bohlin R. C., Gilliland R. L., 2004, AJ, 127, 3508Bohlin R. C., Gordon K. D., Tremblay P. E., 2014, PASP, 126, 711CIE 1951, Proceedings of the Commission Internationale de l’Éclairage,Vol. 1, Sec 4; Vol 3, p. 37, Bureau Central de la CIE, ParisCIE 1990, Commission Internationale de l’Éclairage 1988 2 · Spectral Lu-minous Efficiency Function for Photopic VisionCIE 2010, Commission Internationale de l’Éclairage. Recommended systemfor mesopic photometry based on visual performance.Cinzano P., Falchi F., 2012, MNRAS, 427, 3337Cinzano P., Falchi F., 2020, J. Quant. Spectrosc. Radiative Transfer, 253,107059Cinzano P., Falchi F., Elvidge C. D., 2001, MNRAS, 328, 689Driver S. P., et al., 2016, ApJ, 827, 108Duriscoe D. M., 2013, PASP, 125, 1370ESA 1997, The HIPPARCOS and TYCHO catalogues. Astrometric andphotometric star catalogues derived from the ESA HIPPARCOS SpaceAstrometry Mission. ESA Special Publication Vol. 1200Evans D. W., et al., 2018, A&A, 616, A4Falchi F., et al., 2016, Science Advances, 2, e1600377Gaia Collaboration et al., 2016, A&A, 595, A1Gaia Collaboration et al., 2018, A&A, 616, A1Garstang R. H., 1986, PASP, 98, 364Garstang R. H., 1989, PASP, 101, 306Garstang R. H., 1991, PASP, 103, 1109Górski K. M., Hivon E., Banday A. J., Wand elt B. D., Hansen F. K.,Reinecke M., Bartelmann M., 2005, ApJ, 622, 759Green R. M., 1985, Spherical AstronomyHänel A., et al., 2018, J. Quant. Spectrosc. Radiative Transfer, 205, 278Hart M., 2019a, PASP, 131, 015003Hart M., 2019b, AJ, 157, 221Hong S. S., Kwon S. M., Park Y. S., Park C., 1998, Earth, Planets, andSpace, 50, 487Ienaka N., Kawara K., Matsuoka Y., Sameshima H., Oyabu S., TsujimotoT., Peterson B. A., 2013, ApJ, 767, 80Jechow A., Ribas S. J., Domingo R. C., Hà ¶ lker F., Kolláth Z., Kyba C. C.,2018, Journal of Quantitative Spectroscopy and Radiative Transfer, 209,212Johnson H. L., 1963, Photometric Systems. p. 204Jones A., Noll S., Kausch W., Szyszka C., Kimeswenger S., 2013, A&A,560, A91Jordi C., et al., 2010, A&A, 523, A48Kasten F., Young A. T., 1989, Appl. Opt., 28, 4735Kawara K., Matsuoka Y., Sano K., Brand t T. D., Sameshima H., TsumuraK., Oyabu S., Ienaka N., 2017, PASJ, 69, 31Kocifaj M., 2007, Applied optics, 46, 3013Kocifaj M., 2018, J. Quant. Spectrosc. Radiative Transfer, 206, 260MNRAS , 1–15 (2020) ap of the natural sky brightness Kocifaj M., Kránicz B., 2011, Lighting Research & Technology, 43, 497Kolláth Z., Cool A., Jechow A., Kolláth K., Száz D., Tong K. P., 2020,J. Quant. Spectrosc. Radiative Transfer, 253, 107162Kwon S. M., 1989, Journal of Korean Astronomical Society, 22, 141Kwon S. M., Hong S. S., Weinberg J. L., 2004, New Astron., 10, 91Laureijs R. J., Mattila K., Schnur G., 1987, A&A, 184, 269Leinert C., Richter I., Pitz E., Hanner M., 1980, in Halliday I., McIntoshB. A., eds, Vol. 90, Solid Particles in the Solar System. pp 15–18Leinert C., et al., 1998, A&AS, 127, 1Levasseur-Regourd A. C., Dumont R., 1980, A&A, 84, 277Maiz Apellaniz J., Weiler M., 2018, VizieR Online Data Catalog, ppJ/A+A/619/A180Marrese P. M., Marinoni S., Fabrizio M., Altavilla G., 2019, A&A, 621,A144Matsuoka Y., Ienaka N., Kawara K., Oyabu S., 2011, ApJ, 736, 119Melchior A. L., Combes F., Gould A., 2007, A&A, 462, 965Noll S., Kausch W., Barden M., Jones A. M., Szyszka C., Kimeswenger S.,Vinther J., 2012, A&A, 543, A92Patat F., 2008, A&A, 481.2, 575Roach F. E., Gordon J. L., 1973, The Light of the Night SkyRobin A. C., Reylé C., Derrière S., Picaud S., 2003, A&A, 409, 523Robin A. C., Marshall D. J., Schultheis M., Reylé C., 2012, A&A, 538, A106Sánchez de Miguel A., Kyba C., Aubé M., Zamorano J., Cardiel N., TapiaC., Bennie J., Gaston K. J., 2019, Remote Sensing of Environment, 224,92–103Saran D. V., Slanger T. G., Feng W., Plane J. M. C., 2011, Journal ofGeophysical Research (Atmospheres), 116, D12303Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525Teillet P. M., 1990, Appl. Opt., 29, 1897Toller G. N., 1981, PhD thesis, State University of New York, Stony Brook.Toller G. N., 1990, in Bowyer S., Leinert C., eds, IAU Symposium Vol. 139,The Galactic and Extragalactic Background Radiation. p. 21Treanor P. J., 1973, The Observatory, 93, 117Turnrose B. E., 1974, PASP, 86, 545Unterguggenberger S., Noll S., Feng W., Plane J. M. C., Kausch W.,Kimeswenger S., Jones A., Moehler S., 2017, Atmospheric Chemistry& Physics, 17, 4177Vandersteen J., Kark S., Sorrell K., Levin N., 2020, Remote Sensing, 12,1785Weiler M., 2018, A&A, 617, A138Weiler M., Jordi C., Fabricius C., Carrasco J. M., 2018, A&A, 615, A24Weinberg J. L., Hanner M. S., Beeson D. E., DeShields L. M. I., GreenB. A., 1974, J. Geophys. Res., 79, 3665Westera P., Samland M., Bruzual G., Buser R., 2002, The BaSeL 3.1 models:metallicity calibration and application. p. 166van Leeuwen F., 2007, A&A, 474, 653This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000