text:
# The artificial night sky brightness mapped from DMSP satellite Operational Linescan System measurements

## 1 Introduction

The night sky is a world heritage. In recent decades there has been a rapid increase in the brightness of the night sky in nearly all countries. This has degraded astronomical viewing conditions. The increase in night sky brightness is one of the most noticeable effects of light pollution, which can be defined as the alteration of natural light levels in the outdoor environment due to artificial light sources. The widespread use of artificial lighting with little regard to fixture shielding or energy conservation is the primary source of anthropogenic light contributing to the brightness of the night sky. The astronomical community has expressed its concern over the growth of the sky brightness in a number of official documents and positions (e.g. the Resolutions of the General Assemblies of the International Astronomical Union (IAU) XVI/9 1976, XIX/B6 1985, XX/A2 1988, XXIII/A1 1997 and the Positions of the American Astronomical Society). The UNESCO, the United Nations and the Commission Internationale de l’Eclairage (CIE) give consideration to this concern. The Commission 50 of IAU (The protection of existing and potential astronomical sites) is working actively in order to preserve the astronomical sky and many studies and meetings have been dedicated to this topic (e.g. Crawford 1991; Kovalevsky 1992; McNally 1994; Isobe & Hirayama 1998; Cinzano 2000a; see also Cinzano 1994 for a large reference list). Laws, bills, standard rules, ordinances, regulations limit in many countries the direct wasting of light in the sky from lighting fixtures and, in some cases, also the quantity of light reflected by lighted surfaces. The International Dark-Sky Association is active worldwide in this battle of culture and intelligence with the aim of building awareness and saving for mankind the possibility of feeling part of the universe.

An effective battle against light pollution requires the knowledge of the condition of the night sky across large territories, recognition of vulnerable areas, the determination of the growth trends and the identification of most polluting cities. Therefore a method to map the artificial sky brightness across large territories is required. This is also useful in order to recognize areas with low level of light pollution and potential astronomical sites. In the past, mappings of sky brightness for extended areas has been performed based on population density data with some simple modelling. Examples include the works of Walker (1970, 1973) in California and Arizona, Albers  in USA, Bertiau et al.  in Italy, Berry  and Pike & Berry  in Ontario (Canada). These authors used population data of cities to estimate their upward light emission and a variety of empirical or semi-empirical propagation laws for light pollution in order to compute the sky brightness produced by them. Recently advances in the availability and gain control of the DMSP satellite now provides direct measurements of the upward light emission from nearly the entire surface of the Earth (60°South to 72°North). The first global DMSP image of the Earth at night was produced at 10 km resolution (Sullivan 1989, 1991) using DMSP film strips, the only data available at that time. Beginning in 1992, the availability of digital DMSP data has enabled a new set of higher spatial resolution products (Elvidge et al. 1997a, 1997b, 1997c). These data have been used to model population distribution (Sutton et al. 1997) and urban sprawl impacts on food production (Imhoff et al. 1997a, 1997b). The data have also been used to document light pollution as expressed in the increase in the upward flux of light over time for Japan (Isobe & Hamamura 1998).

In this paper we present a method to map the artificial sky brightness across large territories. In order to bypass errors arising when using population data to estimate upward flux, we construct the maps based on direct measuraments of the upward flux as observed from space using DMSP satellite night-time images and compute the downward flux to the Earth’s surface with detailed modelling of light pollution propagation in the atmosphere. We also present, as an application, detailed maps of artificial sky brightness in Europe in V and B astronomical photometrical bands with a resolution of approximately 1 km. In section 2 we describe the OLS-DMSP satellite observations, their reduction and analysis. We also discuss the relation upward flux versus city population. In section 3 we describe the mapping technique and in section 4 we discuss the application to the maps of Europe. The maps are presented in section 5 together with our comments. Section 6 contains our conclusions.

## 2 Observations and data analysis

### 2.1 Satellite data

U.S. Air Force Defense Meteorological Satellite Program (DMSP) satellites are in low altitude (830 km) sun-synchronous polar orbits with an orbital period of 101 minutes. With 14 orbits per day they generate a global night time and day time coverage of the Earth every 24 hours with their main purpose to monitor the distribution of clouds and to assess navigation conditions. The U.S. Department of Commerce, NOAA National Geophysical Data Center (NGDC) serves as the archive for the DMSP program. The digital archive was initiated in 1992. Starting in 1994, NGDC has embarked on the development of night time lights processing algorithms and products from the DMSP-OLS (Elvidge et al. 1997a, 1997b, 1997c).

The Operational Linescan System (OLS) is an oscillating scan radiometer with low-light visible and thermal infrared (TIR) imaging capabilities which first flew on DMSP satellites in 1976. At night the OLS uses a Photo Multiplier Tube (PMT) to intensify the visible band signal. The purpose of this intensification is to observe clouds illuminated by moonlight. The PMT data have a broad spectral response from 440 to 940 nm (485 - 765 nm FWHM) with highest sensitivity in the 500 to 650 nm region (see Fig. 1). This covers the range for primary emissions from the most widely used lamps for external lighting: Mercury Vapour (545 nm and 575 nm), High Pressure Sodium (from 540 nm to 630 nm), Low Pressure Sodium (589 nm). The sensitivity of the PMT, combined with the OLS-VDGA (variable digital gain amplifier) and fixed gain OLS pre- and post- amplifiers allows measurement of radiances down to $`10^{10}`$ W cm-2 sr-1 $`\mu `$m-1 (Elvidge et al. 1999). This implies that the OLS-PMT could detect radiation with an effective wavelength near 550 nm down to a luminance of approximately 0.2 mcd m-2. The TIR detector is sensitive to radiation from 10.0 to 13.4 $`\mu `$m (10.3 - 12.9 FWHM).

The OLS acquires swaths of data that are 3000 km wide. The OLS sinusoidal scan maintains a nearly constant along-track pixel-to-pixel ground sample distance (GSD) of 0.56 km. Likewise, the electronic sampling of the signal from the individual scan lines maintains a GSD of 0.56 km. The effective instantaneous field of view (EIFOV) is larger than the 0.56 km GSD and varies with scan angle. At nadir the EIFOV is 2.2 km and expands to 5.4 km at the edge of scan 1500 km out from nadir. Most of the data received by NGDC has been ”smoothed” by on-board averaging of 5 $`\times `$ 5 pixel blocks, yielding data with a GSD of 2.8 km. Other details of the OLS are described by Lieske .

DMSP platforms are stabilized using four gyroscopes (three axis stabilization), and platform orientation is adjusted using a star mapper, an Earth limb sensor, and a solar detector (Elvidge et al. 1997a). Daily radar bevel vector sightings of the satellites provided by Naval Space Command allow Air Force orbital mechanics models to compute geodetic subtrack of each orbit giving the position of satellites every 0.4208 seconds. This position together with OLS scan angle equations, an oblate ellipsoid Earth sea level model and digital terrain elevations from U.S. Geological Survey, EROS Data center, allow geolocation algorithms to estimate latitude and longitude for each pixel center.

Normally the OLS is operated at high gain settings for cloud detection. City lights are saturated in these data and radiances can not be extracted. In 1996 and 1997, NGDC made special requests to the Air Force for collection of OLS PMT data at reduced gain settings. This request was granted for the darkest nights of lunar cycles in March of 1996 (8 nights) and January and February (10 nights each) of 1997. During these experimental data collections on board algorithms which adjust the visible band gain were disabled. The major on-board algorithm which effects the night time visible band data is the along scan gain control (ASGC), which adjust the gain dynamically in response to scene brightness. There is also a on-board bi-directional reflectance (BRDF) algorithm designed to reduce the brightness of the image ”hot-spot” which occurs where the observation angles matches the illumination angle. The BRDF algorithm has minimal effect when lunar illumination is low, as was the case during the 28 nights in which the gain controlled data were acquired. The two on-board gain control algorithms were turned off to simplify the retrieval of radiances from the special data collections.

During the special data acquisition the OLS VDGA gain, which is normally operated at 60 dB, was reduced to avoid saturation in urban centres. On one set of nights the gain was operated at a setting of 24 dB. This produced data that avoided saturation on major on major urban centres, but did not permit detection of city edges and lighting in smaller towns. To overcome this dynamic range limitation, data was also acquired at gain settings of 40 and 50 dB.

### 2.2 Data reduction

These data were used to assemble a cloud-free composite image calibrated to top-of-atmosphere (TOA) radiances. The composition provides additional advantages in the removal of ephemeral lights sources, such as fire and lightning, plus the retrieval of lights from small towns that are near the detection limits of the sensor and processing algorithms. The primary processing steps include: 1) establishment of a reference grid with finer spatial resolution than the input imagery using the one kilometre equal area Interrupted Homolosine Projection (Goode 1925; Steinwand 1993) developed for the NASA-USGS Global 1km AVHRR project; 2) identification of the cloud free section of each orbit based on OLS-TIR data; 3) identification of lights, removal of noise and solar glare, cleaning of defective scan lines and cosmic rays; 4) projection of the lights from cloud-free areas from each orbit into the reference grid; 5) calibration to radiance units using prior to launch calibration of digital number for given input telescope illuminance and VDGA gain settings in simulated space conditions; 6) tallying of the total number of light detections in each grid cell and calculation of the average radiance value; 7) filtering images based on frequency of detection to remove ephemeral events.

The final image was transformed into latitude/longitude projection with 30”$`\times `$30” pixel size with data in 8-bit byte format and power law scaling (see Elvidge et al. 1999 for details). The maps of Europe described below were obtained from an image extracted from the global image. The image of Europe is 4800$`\times `$4800 pixels in size, starting at longitude 10°30’ west and latitude 72°north. Orbits involved are listed in appendix B.

Radiance range of final composite image goes from a minimum of $`1.54\times 10^9`$ to a maximum of $`3.17\times 10^7`$ W cm-2 sr-1 $`\mu `$m-1. The minimum luminance detectable for light with effective wavelength of 550 nm (which power is $`1.47\times 10^3`$ W lm-1) is approximately 3 mcd m-2. Assuming an average vertical extinction of $``$0.3 mag in visual band, the minimum detectable luminance on the ground would be of the order of 4 kcd/km2. Two un-shielded fixtures of the ”globe” kind, with clean transparent glass (fixture efficiency $``$80 per cent), equipped with a 250 W high pressure sodium lamp with 125 lm W-1 efficiency, and placed every square kilometer, could be sufficient to produce this luminance.

### 2.3 Data analysis

We first analysed the composite image in order to evaluate the emission versus population relationship. We chose a number of European cities of various populations111Population data refer to 1991 for Italy, Spain and Greece, 1990 for France, 1996 for Germany and have been provided by their national bureaux of statistic. and measured their relative upward flux per unit solid angle summing the counts of all pixels pertinent to each city and multiplying for pixel size at that latitude. Fig. 2 show the measured upward flux of a sample of European cities normalized, for display purposes, to the average flux of a city of $`10^5`$ inhabitants in the same country.

The upward emission increases linearly with the population. One possible source of error is the on-board averaging of the 5$`\times `$5 pixel blocks during the smoothing process. During smoothing it is possible for saturated pixels in the cores of urban centres to be averaged in with unsaturated pixels to produce an unsaturated smoothed pixel value. This phenomenon will be addressed through the inclusion of OLS data acquired at even lower gain settings in the updated nighttime lights map NGDC is preparing for the 1999-2000 time period. Another uncertainty in the city analysis is that we did not precisely match the outline of the lights to match the population data reporting area.

Even if the upward emission vs. city population relationship depends on the local lighting habits, our results are in agreement with relations successfully applied or measured in other countries. Elvidge et al.  found a linear relation between composite radiance and population for a sample of cities of North America and obtained a smaller scatter when correlating radiance with electric consumption. Walker  found linear proportionality between population and street light emission for a number of California cities, with a few departures above or below the mean depending on the industrial or residential character of the city. He also measured the sky illumination produced by three cities, obtaining a dependence on the power of 0.8 of their population. Bertiau et al.  used successfully a linear model with brightness proportional to population to predict light pollution, and so did Walker (1970, 1973) and Garstang (1987, 1988, 1989a, 1989b, 1989c, 1991a, 1991b, 1991c, 1992, 1993, 2000) who obtained a good fit with many observations including Walker  population - distance data. Berry (1976) fitted well observations of sky brightness in city centers in Ontario with a propagation law for light pollution based on the approach of Treanor  but assuming a dependence on the power of 0.5 of the population. However Garstang’s linear models fit Berry’s  observations well, suggesting that the power of 0.5 found by him was produced by the extinction of light emitted by outskirts of large cities in propagating to the center and does not depend on the upward flux versus population relationship (Garstang 1989a).

Bertiau et al.  in the early 1970’s found that the city upward emission in Italy depended on its economic and commercial development, so they were forced to include in their model a development factor. Twentyfive years later we are unable to identify an affluence effect on the brightness of italian cities. Cities in southern Italy have the same light output as comparable size cities in northern Italy, even if the former have a per capita income that is nearly half than the latter.

The proportionality between satellite data and population would allow us to replace the population distribution from census data, usually adopted in computations of night sky brightness, with satellite data, independently by the fact that we would be observing the light coming from the external night-time lighting or whatever else. This replacement constitutes an improvement because census data are based on city lists which, though they can be associate with the city geographical position, do not provide spatially explicit detail of the population geographical distribution.

## 3 Mapping technique

Scattering from atmospheric particles and molecules spreads the light emitted upward by the sources. If $`e(x,y)`$ is the upward emission per unit area in $`(x,y)`$, the total artificial sky brightness in a given direction of the sky in a site in $`(x^{},y^{})`$ is:

$$b(x^{},y^{})=e(x,y)f((x,y),(x^{},y^{}))𝑑x𝑑y,$$
(1)
where $`f((x,y),(x^{},y^{}))`$ give the artificial sky brightness per unit of upward light emission produced by the unitary area in $`(x,y)`$ in the site in $`(x^{},y^{})`$. The light pollution propagation function $`f`$ depends in general on the geometrical disposition (altitude of the site and the area, and their mutual distance), on the atmospheric distribution of molecules and aerosols and their optical characteristics in the chosen photometric band, on the shape of the emission function of the source and on the direction of the sky observed. In some works this function has been approximated with a variety of semi-empirical propagation law like Treanor Law (Treanor 1973; Falchi & Cinzano 2000; Cinzano et al. 2000), Walker Law (Walker 1973; Joseph et al. 1991), Berry Law (Berry 1976; Pike 1976), Garstang Law (Garstang 1991b).

We obtained the propagation function $`f((x,y),(x^{},y^{}))`$ for each couple of points $`(x,y)`$ and $`(x^{},y^{})`$ with detailed models for the light propagation in the atmosphere based on the modelling technique introduced and developed by Garstang (1986, 1987, 1988, 1989a, 1989b, 1989c, 1991a, 1991b, 1991c, 1992, 1993, 2000) and also applied by Cinzano (2000b, 2000c, 2000d). The models assume Rayleigh scattering by molecules and Mie scattering by aerosols and take into account extinction along light paths and Earth curvature. These models allow associating the predictions with well-defined parameters related to the aerosol content, so the atmospheric conditions at which predictions refer can be well known. Here we will describe only the main lines of the models and our specific implementation, leaving the reader to the cited papers for details.

A telescope of area $`\pi d^2/4`$ situated in the observing site O collects from within an infinitesimal section $`dV=\left(\pi ϵ^2u^2du\right)`$ of the cone of angle $`2ϵ`$ around the line-of-sight at a distance $`u`$ and with thickness $`du`$, a luminous flux $`d\mathrm{\Phi }`$ given by:

$$d\mathrm{\Phi }=\frac{\pi d^2}{4}\frac{1}{u^2}M_s(u)\xi _1(u)\left(\pi ϵ^2u^2du\right),$$
(2)
where $`M_s(u)`$ is the luminous flux scattered in unit solid angle toward the observer from particles of aerosol and molecules inside unit volume of atmosphere at the distance $`u`$ along the line of sight, and $`\xi _1(u)`$ is the extinction of the light along its path to the telescope.

Calling $`e_s`$ the upward flux of the source and $`S=M_s/e_s`$ the scattered flux per unit solid angle per unit upward flux, the propagation function $`f`$, expressed as total flux per unit area of the telescope per unit solid angle per unit total upward light emission, is found integrating eq. (2) from the site to infinity:

$$f=_{u_0}^{\mathrm{}}S(u)\xi _1(u)𝑑u,$$
(3)
with:

$$\xi _1=\mathrm{exp}\left[_0^u\left(N_\mathrm{m}(h)\sigma _\mathrm{m}+N_\mathrm{a}(h)\sigma _\mathrm{a}\right)𝑑x\right],$$
(4)
where $`N_\mathrm{m}(h)`$ and $`N_\mathrm{a}(h)`$ are respectively the vertical number densities of molecules and aerosols, $`\sigma _\mathrm{m}`$ and $`\sigma _\mathrm{a}`$ are their scattering cross sections. The altitude $`h`$ depends on the integration variable $`x`$, on the zenith distance and azimuth of the line of sight, on the altitudes of site and source, and on their distance.

The luminous flux per unit solid angle per unit upward flux coming directly from the source and scattered toward the observer from unit volume along the line of sight is:

$$S_d(u)=\left(N_\mathrm{m}(h)\sigma _\mathrm{m}f_\mathrm{m}(\omega )+N_\mathrm{a}(h)\sigma _\mathrm{a}f_\mathrm{a}(\omega )\right)i(\psi ,s),$$
(5)
where $`i(\psi ,s)`$ is the direct illuminance per unit flux produced by each source on each infinitesimal volume of atmosphere along the line-of-sight of an observer and $`f_\mathrm{m}(\omega )`$, $`f_\mathrm{a}(\omega )`$ are the normalized angular scattering functions of molecular and aerosol scattering respectively. The scattering angle $`\omega `$, the emission angle $`\psi `$, the distance $`s`$ of the section from the source and the altitude $`h`$ of it, depend on the altitudes of the site and the source, their distance, the zenith distance and the azimuth of the line-of-sight and the distance $`u`$ along the line of sight, through some geometry.

If $`I(\psi )`$ is the normalized light flux per unit solid angle emitted by the considered source at the zenith distance $`\psi `$ and $`s`$ is the distance between the source and the considered infinitesimal volume of atmosphere, the illuminance per unit flux is:

$$i(\psi ,s)=I(\psi )\xi _2/s^2$$
(6)
in the range where there is no shielding by Earth curvature and zero elsewhere. The extinction $`\xi _2`$ along the path is:

$$\xi _2=\mathrm{exp}\left[_0^s\left(N_\mathrm{m}(h)\sigma _\mathrm{m}+N_\mathrm{a}(h)\sigma _\mathrm{a}\right)𝑑x\right].$$
(7)
A single scattering model is not sufficient to describe the artificial sky brightness produced by a source. In a real atmosphere several scatterings may occur during the travel of a photon from the source to the telescope. The optical thickness $`\tau =k𝑑r`$, where $`k`$ is an attenuation coefficient, determines how important secondary and higher scattering is. If $`\tau >>1`$ (thick layer) multiple scattering is dominant. The fraction of incident radiation which has been scattered once is $`(1\mathrm{e}^\tau )`$ and the fraction which it is scattered again is of order $`(1\mathrm{e}^\tau )^2`$. If $`(1\mathrm{e}^\tau )`$ is sufficiently small, which happens when $`\tau `$ is small, secondary and higher order scattering can be neglected. In absence of aerosol the optical thickness of the atmosphere at wavelength of 550 nm is about 0.1 (Twomey 1977). The aerosol optical thickness can be 0.05 in cleaner regions of the globe, but it can grow to higher values, even depending on seasonal changes (Garstang 1988). Then single scattering is the major contributor to scattered radiation but secondary scattering is not negligible. The error in neglecting third and higher order scattering can be significant for optical thickness higher than about 0.5. Therefore in the computation of the total light flux $`S`$ that molecules and aerosols in the infinitesimal volume scatter toward the observer we take into account light already scattered once. We assumed, as Garstang (1984, 1986), that the light coming to the considered infinitesimal volume along the line of sight after a scatter has approximately a direction near that of the direct light, so that the scattering angle $`\omega `$ in first approximation is always the same. In this case the total illuminance $`S`$ can be written:

$$S=S_dD_S,$$
(8)
where $`D_S`$ is a correction factor which take into account the illuminance due at light already scattered once from molecules and aerosols which can be evaluated with the approach of Treanor  as extended by Garstang (1984, 1986). Details on assumptions can be found in the quoted papers.

## 4 Application

In practice, we divided the surface of Europe into land areas with the same positions and dimensions as projections on the Earth of the pixels of the satellite image. We assumed each land area be source of light pollution with an upward emission $`e_{x,y}`$ proportional to the radiance measured in the corresponding pixel multiplied by the surface area. The total artificial sky brightness at the centre of each area, given by the expression (1), becomes:

$$b_{i,j}=\underset{h}{}\underset{k}{}e_{h,k}f((x_i,y_j),(x_h,y_k)).$$
(9)
We obtained the propagation function $`f((x_i,y_j),(x_h,y_k))`$ for each couple of points $`(x_i,y_j)`$ and $`(x_h,y_k)`$, the positions of the observing site and the polluting area, from eq. (3) after inserting eq. (5) and eq. (6) multiplied for the double scattering factor $`D_S`$ computed as below. We considered every land area as a point source located in its centre except when $`i=h`$, $`j=k`$ in which case we used a four points approximation (Abramowitz & Stegun 1964). The resolution of the maps, depending on results from an integration over a large zone, is greater than resolution of the original images and is generally of the order of the distance between two pixel centres. However where sky brightness is dominated by contributions of nearest land areas, effects of the resolution of the original image could became relevant.

We obtained the maps for B and V photometric astronomical bands (Johnson 1955).

### 4.1 Atmospheric model

#### 4.1.1 Molecular atmosphere vertical distribution

We assumed the atmosphere in hydrostatic equilibrium under the gravitational force. Neglecting the Earth curvature, the force per unit surface supporting a molecular layer of thickness $`dh`$ is $`dp=g\rho dh`$ where $`g`$ is the gravitational acceleration and $`\rho `$ the density of the layer. Replacing $`dp`$ with the equation of state of perfect gas, working for dry air, $`dp=d\rho R\overline{T}/\overline{M}`$, where $`R`$ is the gas constant, $`\overline{M}`$ the mass of a mole of dry air and $`\overline{T}`$ the average temperature, and integrating, we find that in first approximation the number density $`N_\mathrm{m}`$ of the gaseous component of atmosphere decreases exponentially with altitude $`h`$. If $`\overline{n}`$ is the average number of molecules in a mole of dry air:

$$N_\mathrm{m}(h)=\frac{\rho _0\overline{n}}{\overline{M}}\mathrm{exp}\left(\frac{\overline{M}\overline{g}}{R\overline{T}}h\right)=N_{\mathrm{m},0}\mathrm{e}^{ch}.$$
(10)
Measurements show that this is a good approximation for the first 10 km (see e.g. Fig. 3). We adopted this simple model as Garstang , neglecting improvements done by Garstang \[1991a\] for higher altitudes.

As Garstang  we assumed an inverse scale altitude $`c=0.104\mathrm{km}^1`$ and a molecular density at sea level $`N_{\mathrm{m},0}=2.55\times 10^{19}\mathrm{cm}^3`$ .

#### 4.1.2 Haze aerosol vertical distribution

We are interested in average atmospheric conditions, better if typical and not for the particular conditions of a given night, so a detailed modelling of the local aerosol distribution for a given night is beyond the scope of this paper. As Garstang  and Joseph et al.  we assumed an exponential decrease of number density for the atmospheric haze aerosols:

$$N_\mathrm{a}(h)=N_{\mathrm{a},0}\mathrm{e}^{ah}.$$
(11)
Measurements show (see e.g. Fig. 3) that for the first 10 km this is a reasonable approximation. To account for presence of sporadic denser aerosol layers at various heights or at ground level as Garstang \[1991b\] is beyond the scope of this work. We also neglected the effects of the Ozone layer and the presence of volcanic dust studied by Garstang (1991a, 1991c).

We take into account changes in aerosol content as Garstang  introducing a parameter $`K`$ which measures the relative importance of aerosol and molecules for scattering light in V band:

$$K=\frac{N_{\mathrm{a},0}\sigma _\mathrm{a}}{N_{\mathrm{m},0}\sigma _\mathrm{m}11.11\mathrm{e}^{cH}},$$
(12)
where $`H`$ is the altitude over sea level of the ground level. Changing the parameter $`K`$ we are able to compute the map for different aerosol contents, i.e. for different products $`N_\mathrm{a}\sigma _\mathrm{a}`$. As Garstang , the inverse scale altitude of aerosols was assumed to be $`a=0.657+0.059K`$. Effects of changes of aerosol scale altitude has been checked in sec. 4.7.1. More detailed atmospheric models could be used whenever available.

### 4.2 Angular scattering functions

We take into account both Rayleigh scattering by molecules and Mie scattering by aerosols. For molecular Rayleigh scattering the angular scattering function is :

$$f_\mathrm{m}(\omega )=3(1+\mathrm{cos}^2(\omega ))/16\pi .$$
(13)
The integrated Rayleigh scattering cross section in V band was assumed to be $`\sigma _\mathrm{m}=1.136\times 10^{26}`$ cm-2 sr-1 and in B band $`\sigma _\mathrm{m}=4.6\times 10^{27}`$ cm-2 sr-1.

The normalized angular scattering function for atmospheric haze aerosols can be measured easily with a number of well known remote-sensing techniques like classical searchlight probing (see e.g. Elterman 1966), modern bistatic lidar probing, measurements of the day-light or moonlight sky scattering function (see e.g. Hulburt 1951; Volz 1987; Krisciunas & Schaefer 1991). Nevertheless in this paper we are not interested in a specific function for a given site at a given time but in the typical average function so we adopted the function tabulated by Mc Clatchey et al.  as interpolated by Garstang \[1991a\] and we neglected geographical gradients:

$$\begin{array}{c}\mathrm{For}\mathrm{\hspace{0.33em}0}\omega 10^{}:\hfill \\ f_\mathrm{a}(\omega )=7.5\mathrm{exp}\left(0.1249\omega ^2/(1+0.04996\omega ^2)\right),\hfill \\ \\ \mathrm{For}\mathrm{\hspace{0.33em}10}^{}<\omega 124^{}:\hfill \\ f_\mathrm{a}(\omega )=1.88\mathrm{exp}\left(0.07226\omega +0.0002406\omega ^2\right),\hfill \\ \\ \mathrm{For}\mathrm{\hspace{0.33em}124}^{}<\omega 180^{}:\hfill \\ f_\mathrm{a}(\omega )=0.025+0.015\mathrm{sin}(2.25\omega 369.0).\hfill \end{array}$$
(14)
where $`\omega `$ is the scattering angle in degrees. The total integrated scattering cross section $`N_\mathrm{a}\sigma _\mathrm{a}`$ is given by eq. (12) for a given $`K`$. In B band $`\sigma _\mathrm{a}`$ is $`1.216`$ times the $`\sigma _\mathrm{a}`$ in V band.

### 4.3 Upward emission function

The normalized emission function of each area gives the relative upward flux per unit solid angle in every direction. It is the sum of the direct emission from fixtures and the reflected emission from lighted surfaces, normalized to its integral and is not known. The high number of sources contributing to the sky brightness of a site with casual distribution and orientation, smooth the shape of the average normalized emission function which can be considered in first approximation axisymmetric. It is possible to measure the average upward emission of a chosen area at a number of different elevation angles when a large number of satellite measurements from very different orbits is available. It will be possible to obtain it directly by integrating upward emission from all lighting fixtures and all lighted surfaces on the basis of lighting engineering data and models as soon they will be available.

In this paper we assumed that all land areas have the same average normalized emission function. This is equivalent to assuming that lighting habits are similar on average in each land area and that differences from the average are casually distributed in the territory. The average normalized emission function can be constrained from radiance measurements of cities at various distance from satellite nadir. It can also be obtained from comparison of Earth-based observations and models predictions (Cinzano, in prep.). We chose to assume this function and check its consistency with satellite measurements, rather than directly determine it from satellite measurements because at very low elevation angles the spread is too much large to constrain adequately the function shape.

We adopted for the average normalized emission function the normalized city emission function from Garstang  :

$$I(\psi )=\frac{1}{2\pi }\left[2a_1\mathrm{cos}\left(\psi \right)+0.554a_2\psi ^4\right],$$
(15)
where $`a_1`$ and $`a_2`$ are shape parameters. Here we assumed $`a_1=0.46`$, $`a_2=0.54`$ for the typical average function, corresponding to Garstang parameters $`G=0.15`$ and $`F=0.15`$. This function was tested by Garstang with many comparisons between model predictions and measurements. This author assumed this function be produced by the sum of direct emission from fixtures at high zenith distances and Lambertian emission from lighted horizontal surfaces at higher zenith angles. Nevertheless, upward flux can be emitted at all zenith angles both from fixtures and vertical or horizontal surfaces, so we preferred to consider Garstang’s function like a parametric representation with $`a_1`$ and $`a_2`$ as shape parameters without any meaning of fraction of direct and reflected light. We tested also the normalized city emission function of Cinzano (2000b, 2000c) which assumes a slightly higher emission at intermediate angles in respect to the function (15):

$$I(\psi )=\frac{1}{2\pi }\left[a_1+0.554a_2\psi ^4\right],$$
(16)
assuming $`a_1=0.79`$ and $`a_2=0.21`$. Comparison between these functions are presented by Cinzano \[2000b\].

We checked these functions by studying the relation between the upward flux per unit solid angle per inhabitant of a large number of cities and their distance from the satellite nadir in a single orbit satellite image taken the 13th January 1997 h20:27 from satellite F12 which was chosen for its large cloudfree area. Taking into account the average orbital distance $`R_S`$ and the Earth curvature radius $`R_T`$, it is possible, with some geometry, to relate the distance $`D`$ from satellite nadir with the emission angle $`\psi `$:

$$\psi =\frac{D}{R_T}\mathrm{arcsin}\frac{R_T\mathrm{sin}(D/R_T)}{s},$$
(17)
with:

$$s=\sqrt{R_T^2+(R_T+R_S)^22R_T(R_T+R_S)\mathrm{cos}(\frac{D}{R_T})}.$$
(18)
The flux per unit solid angle per inhabitant in relative units is obtained dividing the measurements for the correspondant extinction coefficient $`\xi _3(\psi )`$ computed for curved-Earth from eq. (20) of Garstang \[1989a\]:

$$\begin{array}{c}\xi _3=\mathrm{exp}\left[N_\mathrm{m}\sigma _\mathrm{m}\left(A1+11.778KA2\right)\right]\hfill \\ A1=\frac{1}{c}\mathrm{sec}\psi \left(1\mathrm{e}^{cs\mathrm{cos}\psi }+\frac{16}{9\pi }\frac{\mathrm{tan}^2\psi }{2cR_T}B1\right)\hfill \\ B1=\left(c^2s^2\mathrm{cos}^2\psi +2cs\mathrm{cos}\psi +2\right)\mathrm{e}^{cs\mathrm{cos}\psi }2\hfill \\ A2=\frac{1}{a}\mathrm{sec}\psi \left(1\mathrm{e}^{as\mathrm{cos}\psi }+\frac{16}{9\pi }\frac{\mathrm{tan}^2\psi }{2aR_T}B2\right)\hfill \\ B2=\left(a^2s^2\mathrm{cos}^2\psi +2as\mathrm{cos}\psi +2\right)\mathrm{e}^{as\mathrm{cos}\psi }2.\hfill \end{array}$$
(19)
Fig. 4 shows the relative flux per unit solid angle per inhabitant averaged over ranges of 100 km.

In Fig. 5 we show the $`I(\psi )`$ obtained with eq. (17) and (19) compared with the Garstang  function (solid line) and with the Cinzano \[2000b\] function (dashed line) with the assumed shape parameters. The fits are good with both. Errorbars are not necessarily related to fluctuations in function shape but rather with fluctuations in the total flux per inhabitant. Given that the extinction effects seems to balance the changes in measured flux at angles off from nadir (see Fig.4), we did not need to correct the input images from single orbits for these off-nadir effects when computing the composite image.

Snow reflects approximately 60 per cent of downlight changing the shape of the upward emission function. Lighted roads usually are cleaned in a few days so reflection of street lighting by snow on road surfaces is unlikely to be important, but reflection of the artificial sky light by snow on the rest of the land could enhance noticeably the upward flux in more polluted areas. Assuming roughly that 10 per cent of the upward flux be scattered downward by atmospheric particles and molecules and that 60 per cent of it be reflected upward again by the snow covered terrain, the increase of artificial sky brightness by snow reflexion of sky light could reach 6 per cent. The upward flux due to snow reflection in some zones of Europe (our images are taken in winter) is likely detected by the OLS-PMT, but the different shape of the upward emission function could produce small errors. We plan to use specific satellite surveys to detect snowed areas where the upward emission function must be corrected. Even the offshore lights in the North Sea, where oil and gas production sites are active, could have a different upward emission function.

### 4.4 Geometric relations

In this paper we are more interested in understanding and comparing light pollution distributions in the European territory rather than in predicting the effective sky brightness for observational purposes, so we computed the artificial sky brightness at sea level, in order to avoid introduction of altitude effects in our maps. We plan to take account of altitudes in a forecoming paper devoted to mapping the naked-eye star visibility which requires the computation of star-light extinction and natural sky brightness for the altitude of each land area.

With the hypothesis of sea level, geometrical relations from Garstang \[1989a\] between quantities shown in Fig. 6 taking into account Earth curvature, became simpler. They are listed in the appendix A. In eq. (12) now is $`H=0`$.

Equations (4), (7) have been integrated by Garstang (1989a: eq. (18), (19), (22) and eq. (20), (21)). They are listed in appendix A too.

Correction for double scattering of Garstang \[1989a\] became:

$$D_S=1+N_{\mathrm{m},0}\sigma _\mathrm{m}(11.11Kf_2+\frac{f_1}{3}),$$
(20)
where $`f_1`$ and $`f_2`$ are given in eq. (33). Integration of eq. (3) must be done only where the Earth curvature does not shield the line of sight from the source. We need to start integration from $`u_0`$ as given by (from eq. (12) and (13) of Garstang 1989a):

$$u_0=\frac{2R_T\mathrm{sin}^2(D/2R_T)}{\mathrm{sin}z\mathrm{cos}\beta \mathrm{sin}(D/R_T)+\mathrm{cos}z\mathrm{cos}(D/R_T)}.$$
(21)
We neglected the presence of mountains which might shield the light emitted from the sources from a fraction of the atmospheric particles along the line-of-sight of the observer. Assuming flat Earth, mountains between source and observatory shield the light emitted from the source with an angle less than $`\theta =arctg\frac{H}{p}`$ where $`H`$ is the height of the mountain and $`p`$ its distance from the source. The ratio between the luminance in the shielded and not shielded cases is given, in first approximation, by the ratio between the number of particles illuminated in the two cases:

$$\frac{b_s}{b_{ns}}\frac{\sigma _\mathrm{a}f_\mathrm{a}(\omega )_{\frac{Hq}{p}}^{\mathrm{}}N_\mathrm{a}(h)𝑑h+\sigma _\mathrm{m}f_\mathrm{m}(\omega )_{\frac{Hq}{p}}^{\mathrm{}}N_\mathrm{m}(h)𝑑h}{\sigma _\mathrm{a}f_a(\omega )_0^{\mathrm{}}N_\mathrm{a}(h)𝑑h+\sigma _\mathrm{m}f_m(\omega )_0^{\mathrm{}}N_\mathrm{m}(h)𝑑h},$$
(22)
where $`q`$ is the distance of the site from the source, $`N_\mathrm{m}(h)`$ and $`N_\mathrm{a}(h)`$ are respectively the vertical number densities of molecules and aerosols at the altitude $`h`$, $`\sigma _\mathrm{m}`$ and $`\sigma _\mathrm{a}`$ are their scattering sections and roughly $`\omega \pi \theta `$. This expression shows that, given the vertical extent of the atmosphere in respect to the highness of the mountains, the shielding is not negligible only when the source is very near the mountain and both are quite far from the site (Garstang 1989a; see also Cinzano 2000c): $`\frac{q}{p}>>\frac{h_\mathrm{a}}{H}`$ where $`h_\mathrm{a}`$ is the vertical scale height of the atmospheric aerosols. This condition in the considered territories usually applies to poorly lighted areas only, which produce little pollution. Earth curvature emphasizes this behaviour.

### 4.5 Relation with atmospheric conditions

The adopted modelling technique allows us to assess the atmospheric conditions for which a map is computed giving observable quantities like the vertical extinction at sea level in magnitudes (Garstang 1989a):

$$\mathrm{\Delta }m=1.0857N_{\mathrm{m},0}\sigma _\mathrm{m}\left(\left(\frac{550}{\lambda }\right)^4\frac{1}{c}+\left(\frac{550}{\lambda }\right)\frac{11.778K}{a}\right),$$
(23)
where for V band $`\lambda =550`$ and for B band $`\lambda =440`$. Neglecting Earth curvature, the horizontal visibility defined as the distance at which a black object would show a brightness of 0.98 of the brightness of background horizon due to scattered light is (Garstang 1986):

$$\mathrm{\Delta }x=\frac{3.91}{N_{\mathrm{m},0}\sigma _\mathrm{m}}\left(\left(\frac{550}{\lambda }\right)^4+11.778K\left(\frac{550}{\lambda }\right)\right)^1.$$
(24)
Effects of curvature have been discussed by the cited author. Other relations exist with measurables like the Linke turbidity factor for total solar radiation received on Earth (Garstang 1988). The optical thickness $`\tau `$ is (from eq. (22) of Garstang 1986):

$$\tau =\mathrm{\Delta }m/1.0857.$$
(25)
With $`K`$=1 we obtain $`\tau =0.3`$ so double scattering approximation is adequate in our map computations. In order that $`\tau 0.5`$ we need $`K`$2.2.

### 4.6 Spectral emission

DMSP satellites do not carry a spectrograph able to obtain spectra of the upward light emitted by each land area. However it is possible to recover information about the spectral emission from the differential effects of extinction on light of different wavelengths. Calling $`I(\lambda )`$ the specific radiance of each land area at wavelength $`\lambda `$, $`T(\lambda )`$ the PMT sensitivity curve and $`\xi _4(\lambda ,K)`$ the vertical extinction for a given clarity parameter $`K`$, the radiance $`r(\lambda )`$ measured by the OLS-PMT is:

$$r(K)=_0^{\mathrm{}}T(\lambda )I(\lambda )\xi _4(\lambda ,K)𝑑\lambda .$$
(26)
The vertical extinction can be obtained from eq. (23):

$$\begin{array}{c}\xi _4(\lambda ,K)=\mathrm{exp}\left(p_1\lambda ^4\right)\mathrm{exp}\left(p_2K\lambda ^1\right)\hfill \\ p_1=N_{\mathrm{m},0}\sigma _\mathrm{m}\lambda _0^4/c\hfill \\ p_2=11.778N_{\mathrm{m},0}\sigma _\mathrm{m}\lambda _0/a,\hfill \end{array}$$
(27)
where the molecular scattering cross section $`\sigma _\mathrm{m}`$ at $`\lambda _0=550`$ nm, the molecular density at sea level $`N_{\mathrm{m},0}`$ and the inverse scale heights of molecules and aerosols $`c`$ and $`a`$ are given in section 4.1. Equation (26) could be inverted with a Richardson-Lucy algorithm in order to recover $`I(\lambda )`$ from $`r(K)`$, a function which can be obtained comparing many radiance calibrated images obtained in different atmospheric conditions. The average clarity parameter $`K`$ for each image can be obtained comparing the variations of measured radiance with increasing distance from satellite nadir.

Fig. 7 shows in the lower panel the $`r(K)`$ curves produced by the example spectra $`I(\lambda )`$ in the upper panel. Curves have been scaled so that $`r(K=1)=1`$. Spectra with a similar effective wavelength after convolution with the PMT sensitivity curve, like e.g. a constant (dashed curve), a blackbody at 4000°K (dotted curve) and a narrow gaussian centered at 650 nm (long dashed curve), give almost the same $`r(K)`$ curve. However spectra with lower or higher effective wavelength, like a gaussian centered at 440 nm (dot-long dashed curve), 550 nm (solid curve) or at 800 nm (dot-dashed curve), give different $`r(K)`$ curves.

We obtained the maps for B and V photometric astronomical bands. The relative B-V color index of each land area could be inferred from the $`r(K)`$ curves but requires some assumptions on the shape of the spectra. In order to be simple in this paper we assumed it to be constant everywhere and we only took into account the different propagation of the light in the atmosphere in the two bands. We plan to study differences in city spectra in future works.

### 4.7 Calibration

We calibrated the maps on the basis of both (i) accurate measurements of sky brightness together with extinction from the Earth-surface and (ii) analysis of before-fly radiance calibration of OLS-PMT.

#### 4.7.1  Calibration with Earth-based measurements

A detailed calibration requires sky brightness measurements at a large number of sites on sea level, taken in nights with the same vertical extinction and horizontal visibility of the map under calibration, averaged over many nights in order to smooth atmospheric fluctuations. Observations need to be under the atmosphere, i.e. as actually observed from the ground without any extinction correction applied. To obtain the artificial sky brightness it would be necessary to measure the natural sky brightness in some sites where the maps suggest that the artificial one is negligible, and subtract the mean from all measurements. Moreover given the fast growth rate of artificial sky brightness, which e.g. in Italy reaches 10 per cent per year (Cinzano 2000d), measurements have to be taken in the same period and the calibration will refer to that time.

Measurements of sky brightness in Europe in V and B bands in the period 1996-1999 are scarce, so we calibrated our maps with all available measurements in the chosen bands taken in 1998 and in 1999 in clean or photometric nights even if extinction was not available or not exactly the required one (Catanzaro & Catalano 2000; Cinzano 2000d; Della Prugna 1999; Favero et al. 2000; Piersimoni, Di Paolantonio & Brocato 2000; Poretti & Scardia 2000; Zitelli 2000). Some data has been taken by one of us for this purpose with a small portable telescope and a CCD device (Falchi 1999). Most of the sites are at sea level but we included also a few sites at altitudes under 1300 m o. s. l.. In lack of measurements of natural sky brightness in Europe at sea level, we assumed it at minimum solar activity $`B=22.7`$ mag arcsec-1 in B band, and $`V=21.6`$ mag arcsec-1 in V band, estimating an incertitude of at least $`\pm 0.1`$ mag arcsec-1. Natural sky brightness increases when solar activity increases (Walker 1988) and the solar activity in 1998 was close to minimum but not at the minimum, so it could be underestimated and consequently the artificial brightness in darker sites considerably overestimated. The sky brightness has been transformed into photon radiance with formulae of Garstang (1989a: eq. (28) and eq. (39)).

A least square fit of a straight line $`y=a+x`$ over the logarithmic measured radiances versus the logarithmic predicted radiances gives the logarithmic calibration coefficients $`a_B=0.63\pm 0.04`$ and $`a_V=0.00\pm 0.04`$. We assumed unavailable the incertitudes of measurements as given by atmospheric conditions and emission function fluctuations. The uncertainty of the calibration coefficients produces an incertitude of about 10 per cent in the calibrated predicted radiances. However single data points show differences even of about 60 per cent, so that we consider safer to assume this last value as estimate of the uncertainty of our calibrated predictions for a given site. More precise calibrations will be possible when a large number of measurements of sky brightness at sea level together with extinction become available. A large CCD measurement campaign is being set up.

In Figs 8 and 9 we compared our calibrated predictions with the available measurements of artificial sky brightness respectively in V and in B bands. Photon radiance in V band is expressed in units of 0.3419 1010 ph s-1 m-2 sr-1, corresponding approximately to a luminance of one $`\mu `$cd m-2, and in B band in units of 1010 ph s-1 m-2 sr-1. Errorbars indicate measurement errors which are much smaller than the effects of fluctuations in atmospheric conditions.

We checked the effects on predictions of Figs 8 and 9 of changes of (i) atmospheric conditions and (ii) the shape of $`I(\psi )`$. An increase of the aerosol content parameter $`K`$ increases exponentially the extinction and linearly the scattering along the line of sight, so that, excluding a small area very near the source where the extinction is negligible, the predicted artificial sky brightness decreases everywhere: a well known result (Garstang 1986, 1989a, 2000; Cinzano 2000c). The decrease is larger for areas farther from sources. These usually are darker, so in the log-log diagram, increasing haze darker sites move toward low predicted brightness more than less dark sites. A change of aerosol scale length has a similar effect. A change of the shape of the upward emission function $`I(\psi )`$ has a different effect. Decreasing the relative emission at low elevations (large $`\psi `$) the sky brightness is decreased approximately proportionally in almost all sites, excepts very near the source (see also Cinzano & Diaz Castro 2000), so in a log-log diagram all sites move toward low predicted brightness of nearly the same quantity. This opens a way to detect cities with more light-wasting installations. In fact if a series of measurements of sky brightness at increasing distances from a city lies on a straight line parallel to the calibration line but displaced toward lower predicted values, likely the relative emission of the city at low elevation is higher than assumed in making the map. This is a typical result of poorly shielded or too much inclined street light.

Shifts in measurements obtained with different instrumental setups could arise because, as will be shown in Fig. 10, the average emission spectrum has typically its maximum at a side of the V band sensitivity curve, so that the resulting measurements are quite sensitive to little differences between the instrumental curve and the standard V band curve.

#### 4.7.2 Calibration with pre-fly irradiance measurements

Map calibration based on pre-fly irradiance calibration of OLS PMT requires the knowledge, for each land area $`(i,j)`$, of (i) the average vertical extinction $`\mathrm{\Delta }m`$ during satellite observations and (ii) the relation between the radiance in the chosen photometrical band and the radiance measured in the PMT spectral sensitivity range, which depends on the emission spectra. Both of them are unknown.

If $`\overline{r}`$ is the energetic radiance in $`\left[10^{10}\mathrm{Wcm}^2\mathrm{sr}^1\right]`$, the upward energy flux in $`\left[W\right]`$ in the PMT photometric band is:

$$e(i,j)=\overline{r}(i,j)\frac{A(i,j)}{I(\psi =0)}\mathrm{\hspace{0.17em}10}^{0.4\mathrm{\Delta }m},$$
(28)
where $`I(\psi )`$ is the upward emission function and A is the surface area in $`km^2`$ of the land area $`(i,j)`$:

$$A(i,j)=\left(\frac{2\pi \mathrm{\Delta }x}{3606060}R_T\right)^2\mathrm{cos}(l),$$
(29)
with $`l`$ latitude of the land area, $`R_T`$ average Earth radius in km and $`\mathrm{\Delta }x`$ pixel size in arcsec.

The photon radiance in the V photometric band corresponding to the energetic radiance measured by the PMT is:

$$\begin{array}{c}\overline{f_r}(i,j)=\frac{_0^{\mathrm{}}T_\mathrm{V}(\lambda )I(\lambda )\frac{\lambda }{hc}𝑑\lambda }{_0^{\mathrm{}}T_{\mathrm{P}\mathrm{M}\mathrm{T}}(\lambda )I(\lambda )\frac{\lambda }{hc}𝑑\lambda }\frac{_0^{\mathrm{}}T_{\mathrm{P}\mathrm{M}\mathrm{T}}(\lambda )S(\lambda )\frac{\lambda }{hc}𝑑\lambda }{_0^{\mathrm{}}T_{\mathrm{P}\mathrm{M}\mathrm{T}}(\lambda )S(\lambda )𝑑\lambda }\overline{r}(i,j)\hfill \\ =\frac{_0^{\mathrm{}}T_\mathrm{V}(\lambda )I(\lambda )\frac{\lambda }{hc}𝑑\lambda }{_0^{\mathrm{}}T_{\mathrm{P}\mathrm{M}\mathrm{T}}(\lambda )I(\lambda )\frac{\lambda }{hc}𝑑\lambda }\frac{<\lambda >}{hc}\overline{r}(i,j),\hfill \end{array}$$
(30)
where $`T_V`$ and $`T_{\mathrm{P}\mathrm{M}\mathrm{T}}`$ are the sensitivity curves respectively of V band and PMT detector, $`I(\lambda )`$ is the energy spectrum of the emission from the chosen land area, $`S(\lambda )`$ is the energy spectrum of the pre-fly calibration source, $`h`$ is the Plank constant and $`c`$ is the velocity of light. The second ratio is the effective wavelength $`<\lambda >`$ of the combination of the sensitivity curves of the PMT and the calibration source, divided by $`hc`$.

In order to check the consistency of our Earth-based V-band calibration with pre-fly radiance calibration of PMT images, we obtained a tentative map calibration assuming for all land areas an average vertical extinction $`\mathrm{\Delta }m=0.3`$ mag V at sea level and constructing a synthetic spectrum for a typical night-time lighting. We very roughly assumed that 50 per cent of the total power emitted by each land area be produced by HPS lamps (SON standard) and 50 per cent by Hg vapour lamps (HQI). The composite spectrum is visible in Fig. 10 (solid line) together with the V band (dashed line).

The result of integration of eq. (30) is $`\frac{\overline{f_r}}{\overline{r}}1.48\times 10^{18}`$ ph s-1 W-1, much less than the photon flux per unit power at 550 nm, $`2.79\times 10^{18}`$ ph s-1 W-1. Eq. (28) and eq. (30), taking in account the internal constants of the program and omitting the $`\mathrm{cos}(l)`$ already included in it, give the V-band logarithmic calibration coefficient $`a_V=0.01`$. In spite of the uncertainties both in the extinction and in the average emission spectrum, this calibration agree very well with the Earth-based calibration.

## 5 Results

Figs 11, 12, 13 and 14 show the artificial sky brightness in Europe at sea level in V and B bands. Colours correspond to ratios between the artificial sky brightness and the natural sky brightness of: $`<`$0.11 (black), 0.11-0.33 (blue), 0.33-1 (green), 1-3 (yellow), 3-9 (orange), $`>`$9 (red). Original maps are 4800$`\times `$4800 pixel images saved in 16-bit standard fits format with fitsio Fortran-77 routines developed by HEASARC at the NASA/GSFC. Images have been analysed with ftools 4.2 analysis package by HEASARC and with quantum image 3.6 by Quantum Image Systems. Maps have been computed for clean atmosphere with aerosol clarity $`K=1`$, corresponding to a vertical extinction of $`\mathrm{\Delta }m=0.33`$ mag in V band, $`\mathrm{\Delta }m=0.56`$ mag in B band, horizontal visibility $`\mathrm{\Delta }x=26`$ km, optical depth $`\tau =0.3`$. We limited our computations to zenith sky brightness even

though our method allows determinations of brightness in other directions. This would be useful to predict visibility in large territories of particular astronomical phenomena, like e.g. comets. A complete mapping of the artificial brightness of the sky of a site, like Cinzano \[2000c\], using satellite data instead of population data is possible (Cinzano 2000, in prep.). Falchi & Cinzano  and Cinzano et al.  obtained in 1998 the first maps of the artificial sky brightness from satellite data using a DMSP single-orbit image and replacing $`f`$ in eq. (1) with the Treanor Law, a semi-empirical law which assumes a very simplified model with homogeneous atmosphere, vertical heights small in relation to the horizontal distances, scattering limited to a cone of small angle co-axial with the direct beam and flat Earth. Differences with our maps mainly arise where Earth curvature play a role in limiting the propagation of light pollution. Our study constitute the natural improvement of their seminal work.

Recommendation 1 of the IAU Commission 50 (Smith 1979) states that the increase in sky brightness at 45° elevation due to artificial light scattered from clear sky should not exceed 10 per cent of the lowest natural level in any part of the spectrum between wavelengths 3000Å and 10000Å. So this is the level over which the sky must be considered ”polluted”. The maps shows that only a few areas in Europe are under the limit of 10 per cent at zenith and some of them could still result in being quite polluted at higher zenith distances.

## 6 Conclusions

We presented a method to map the artificial sky brightness in large territories in astronomical photometric bands with a resolution of approximately 1 km. We computed the maps with detailed models for the propagation in the atmosphere of the upward light flux measured with DMSP satellites Operational Linescan System. The use of this modelling technique allows us to (i) assess the atmospherical conditions for which the maps are computed giving observable quantities and (ii) take into account Earth curvature. This cannot be done properly when using semiempirical propagation laws. The use of satellite data constitutes an improvement over the use of population data to estimate upward flux because (i) it allows spatially explicit detail of the geographic distribution of emissions and (ii) some polluting sources, like e.g. industrial areas, ports and airports, are not well represented in population data.

We presented, as an application of the described method, the maps of artificial sky brightness in Europe at sea level in V and B bands. We are extending the maps to the rest of the World in a fore coming World Atlas of Artificial Sky Brightness and preparing predictions for the state of the night sky in future years.

## Acknowledgments

We are indebted to Roy Garstang of JILA-University of Colorado for his friendly kindness in reading and commenting on this paper, for his helpful suggestions and for interesting discussions.

## Appendix A Geometrical relations

With the hypotesis of sea level, geometrical relations from Garstang \[1989a\] between quantities in Fig. 6 became:

$$\begin{array}{c}s=\sqrt{(ul)^2+4ul\mathrm{sin}^2(D/2R_T)}\hfill \\ \mathrm{with}l=\sqrt{4R_T^2\mathrm{sin}^2(D/2R_T)}\hfill \\ h=R_T\left(\sqrt{1+\frac{u^2+2uR_T\mathrm{cos}(z)}{R_T^2}}1\right)\hfill \\ \omega =\theta +\varphi \hfill \\ \mathrm{with}\theta =\mathrm{arccos}(\frac{q_1}{\sqrt{4R_T^2\mathrm{sin}^2(D/2R_T)}})\hfill \\ \varphi =\mathrm{arccos}(\frac{l^2+s^2u^2}{2ls})\hfill \\ \psi =\mathrm{arccos}(\frac{q_2}{s})\hfill \end{array}$$
(31)
$$\begin{array}{c}q_1=R_T(\mathrm{sin}(D/R_T)\mathrm{sin}(z)\mathrm{cos}(\beta )+\mathrm{cos}(D/R_T)\mathrm{cos}(z))\hfill \\ q_2=u\mathrm{sin}(z)\mathrm{cos}(\beta )\mathrm{sin}(D/R_T)+q_3\hfill \\ q_3=u\mathrm{cos}(z)\mathrm{cos}(D/R_T)2R_T\mathrm{sin}^2(D/2R_T)\hfill \end{array}$$

where $`R_T`$ is the Earth curvature radius.

Equations (4), (7) have been integrated by Garstang (1989a: eq. (18), (19), (22) and eq. (20), (21)). For sea level and $`z<90`$° and $`\psi <90`$°, they are:

$$\begin{array}{c}\xi _1=\mathrm{exp}\left(N_{\mathrm{m},0}\sigma _\mathrm{m}(p_1+11.778Kp_2)\right)\hfill \\ p_1=c^1\mathrm{sec}z(1\mathrm{exp}(cu\mathrm{cos}z)+\frac{16p_3\mathrm{tan}^2z}{9\pi 2cR_T}\hfill \\ p_2=a^1\mathrm{sec}z(1\mathrm{exp}(au\mathrm{cos}z)+\frac{16p_4\mathrm{tan}^2z}{9\pi 2aR_T}\hfill \\ p_3=(c^2u^2\mathrm{cos}^2z+2cu\mathrm{cos}z+2)\mathrm{exp}(cu\mathrm{cos}z)2\hfill \\ p_4=(a^2u^2\mathrm{cos}^2z+2au\mathrm{cos}z+2)\mathrm{exp}(au\mathrm{cos}z)2\hfill \end{array}$$
(32)
$$\begin{array}{c}\xi _2=\mathrm{exp}\left(N_{\mathrm{m},0}\sigma _\mathrm{m}(f_1+11.778Kf_2)\right)\hfill \\ f_1=c^1\mathrm{sec}\psi (1\mathrm{exp}(cs\mathrm{cos}\psi )+\frac{16f_3\mathrm{tan}^2\psi }{9\pi 2cR_T}\hfill \\ f_2=a^1\mathrm{sec}\psi (1\mathrm{exp}(as\mathrm{cos}\psi )+\frac{16f_4\mathrm{tan}^2\psi }{9\pi 2aR_T}\hfill \\ f_3=(c^2s^2\mathrm{cos}^2\psi +2cs\mathrm{cos}\psi +2)\mathrm{exp}(cs\mathrm{cos}\psi )2\hfill \\ f_4=(a^2s^2\mathrm{cos}^2\psi +2as\mathrm{cos}\psi +2)\mathrm{exp}(as\mathrm{cos}\psi )2\hfill \end{array}$$
(33)
For $`\psi =90`$° the reader is referred to the cited paper.

## Appendix B Used DMSP orbits

Orbits of the DMSP satellite F12 used in the composite radiance calibrated image.

| 199603191827 | 199702061858 | 199701091752 |
| --- | --- | --- |
| 199603192009 | 199702062040 | 199701091934 |
| 199603192151 | 199702062222 | 199701092116 |
| 199603211803 | 199702081833 | 199701111727 |
| 199603211945 | 199702082015 | 199701111909 |
| 199603212127 | 199702082157 | 199701112051 |
| 199603231738 | 199702101809 | 199701131703 |
| 199603231920 | 199702101951 | 199701131845 |
| 199603232102 | 199702102133 | 199701132027 |
| 199701061828 | 199603171709 | 199701132209 |
| 199701062010 | 199603171851 | 199702031752 |
| 199701062152 | 199603172215 | 199702031934 |
| 199701081804 | 199603181839 | 199702032116 |
| 199701081946 | 199603182021 | 199702032258 |
| 199701082128 | 199603182203 | 199702051728 |
| 199701101740 | 199603201815 | 199702051910 |
| 199701101922 | 199603201957 | 199702052052 |
| 199701102104 | 199603202139 | 199702071703 |
| 199701121715 | 199603221750 | 199702071845 |
| 199701121857 | 199603221932 | 199702072027 |
| 199701122039 | 199603222114 | 199702072209 |
| 199701122221 | 199701051841 | 199702091821 |
| 199702041740 | 199701052023 | 199702092003 |
| 199702041922 | 199701052205 | 199702092145 |
| 199702042104 | 199701071816 | 199701071958 |
| 199702061716 | 199701072140 | 199701091752 |


format: text

text:
# BARYON NUMBER ASYMMETRY INDUCED BY COHERENT MOTIONS OF A COSMOLOGICAL AXION-LIKE PSEUDOSCALAR

## 1 Introduction.

The origin of the baryon asymmetry of the universe remains one of the most fundamental open questions in high energy physics and cosmology. In 1967 Sakharov noticed  that three conditions are essential for the creation of a net baryon number in a previously symmetric universe: 1) baryon number non-conservation; 2) C and CP violation; 3) out of equilibrium dynamics. Since then many different hypothetical scenarios for baryogenesis have been proposed. A dramatic conclusion emerged from the studies of these scenarios: new physics beyond the Standard Model (SM) is required to explain baryogenesis .

It has been recently realized  that topologically non-trivial configurations of hypercharge gauge fields can be relevant players in the electroweak (EW) scenario for baryogenesis. Hypercharge fields couple anomalously to fermionic number densities in the symmetric phase of the EW plasma, while their surviving long-range projections onto usual electromagnetic fields in the broken phase of the plasma do not. As a consequence, the hypercharge Chern-Simons (CS) number stored in the symmetric phase just before the transition can be converted into a fermionic asymmetry along the direction $`BL=0`$ when the EW symmetry is spontaneously broken.

A hypothetical axion-like pseudoscalar field coupled to hypercharge topological number density can amplify hyperelectric and hypermagnetic fields in the unbroken phase of the EW plasma, while coherently rolling or oscillating around the minimum of its potential. The coherent motion provides the three Sakharov’s conditions and is capable of generating a net CS number that can survive until the phase transition and then be converted into baryonic asymmetry , as first noted in . The mechanism could explain the origin of the baryon number of the universe, if the EW phase transition is strong enough such that the generated asymmetry is not erased by B-violating processes in thermal equilibrium in the broken phase of the plasma.

Pseudoscalar fields with axion-like coupling to CS number densities appear in several possible extensions of the SM. They were originally proposed as an elegant solution to the strong CP-problem. In models with an extended higgs-sector the physical pseudoscalar can couple to hypercharge topological density through quantum effects. In supergravity or superstring models axions that couple to extended gauge groups are common.

Experimental signatures of this hypothetical particle could appear in present and/or future colliders. Since the hypercharge photon is a linear combination of the ordinary photon and the $`Z`$, the particle can be produced in association with a photon or a $`Z`$, and detected through its decay into a pair of neutral gauge bosons, if these signatures are not overshadowed by their SM backgrounds .

## 2 The model

We will assume that the universe is homogeneous and isotropic, and can be described by a conformally flat metric. In addition to the SM fields we consider a time-dependent pseudoscalar field $`\varphi (\eta )`$ with coupling $`\frac{\lambda }{4}\varphi Y\stackrel{~}{Y}`$ to the $`U(1)_Y`$ hypercharge field strength and a potential $`V_0^4V(\varphi /f)`$ generated by processes at energies higher than the EW scale. We will also assume that the universe is radiation dominated at some early time before the scalar dynamics becomes relevant. The coupling constant $`\lambda =1/M`$ has units of mass-1. For a typical axion coupled to a non-abelian gauge topological density the potential is generated by non-perturbative effects at the confinement scale $`V_0`$, and $`fM`$. In general, this is not always the case, and we will allow $`\mathrm{\Lambda }f/M>1`$, but keeping the pseudoscalar mass $`mV_0^2/f`$ much smaller than the scale $`M`$.

Maxwell’s equations describing hyper EM fields in the symmetric phase of the highly conducting EW plasma, coupled to the heavy pseudoscalar are the following,

$`(i)`$ $`\stackrel{}{E}=0(ii)\stackrel{}{B}=0`$ (1)
$`(iii)`$ $`\stackrel{}{J}=\sigma \stackrel{}{E}(iv){\displaystyle \frac{\stackrel{}{B}}{\eta }}=\times \stackrel{}{E}`$ (2)
$`(v)`$ $`{\displaystyle \frac{\stackrel{}{E}}{\eta }}=\times \stackrel{}{B}\lambda {\displaystyle \frac{d\varphi }{d\eta }}\stackrel{}{B}\stackrel{}{J}.`$ (3)
Equations $`(i)`$ and $`(iii)`$ are valid for wavelengths larger than the typical collision length $`r_{Debye}T^1`$ in the hot plasma, so that individual charges are screened by collective effects. The description of short wavelength modes should take into account charge separation in the medium.

We have rescaled the electric and magnetic fields $`\stackrel{}{E}=a^2(\eta )\stackrel{}{}`$, $`\stackrel{}{B}=a^2(\eta )\stackrel{}{}`$ and the physical conductivity $`\sigma =a(\eta )\sigma _c`$, where $`a(\eta )`$ is the scale factor of the universe and $`\eta `$ is conformal time. In the EW plasma $`\sigma _c10T`$. The fields $`\stackrel{}{}`$, $`\stackrel{}{}`$ are the flat space EM fields, and we have assumed for simplicity vanishing bulk velocity $`\stackrel{}{v}`$ of the plasma and zero chemical potentials for all species.

The equation for the pseudoscalar $`\varphi `$ is the following

$$\frac{d^2\varphi }{d\eta ^2}+2aH\frac{d\varphi }{d\eta }+a^2m^2\varphi =\lambda a^2\stackrel{}{E}\stackrel{}{B},$$
(4)
where $`H=\frac{1}{a^2}\frac{da}{d\eta }`$ is the Hubble parameter. We will neglect the backreaction of the electromagnetic fields on the scalar field since is irrelevant for most of the physics we would like to explore . We therefore solve eq.(4) with vanishing r.h.s., and substitute the resulting $`\varphi (\eta )`$ into eq.(3).

## 3 Amplification of Primordial Hypermagnetic Fields

We will describe solutions to eq.(3) of the form $`\stackrel{}{E}(\stackrel{}{x},\eta )=d^3\stackrel{}{k}e^{i\stackrel{}{k}\stackrel{}{x}}\stackrel{}{e}_\stackrel{}{k}ϵ_\stackrel{}{k}(\eta )`$, $`\stackrel{}{B}(\stackrel{}{x},\eta )=d^3\stackrel{}{k}e^{i\stackrel{}{k}\stackrel{}{x}}\stackrel{}{b}_\stackrel{}{k}\beta _\stackrel{}{k}(\eta )`$, for which the electric and magnetic modes are parallel to each other. We find that the Fourier modes $`\stackrel{}{e}_\stackrel{}{k}`$ and $`\stackrel{}{b}_\stackrel{}{k}`$ are related, $`\stackrel{}{b}_\stackrel{}{k}{}_{}{}^{\pm }=b_k^\pm (\widehat{e}_1\pm i\widehat{e}_2),k\stackrel{}{e}_\stackrel{}{k}{}_{}{}^{\pm }ϵ_{\stackrel{}{k}}^{\pm }(\eta )=\pm \stackrel{}{b}_\stackrel{}{k}{}_{}{}^{\pm }\frac{\beta _\stackrel{}{k}}{\eta }_{}^{\pm },`$ where $`\widehat{e}_1`$, $`\widehat{e}_2`$ are unit vectors in the plane perpendicular to $`\stackrel{}{k}`$ such that $`(\widehat{e}_1,\widehat{e}_2,\widehat{k})`$ is a right-handed system. The function $`\beta _\stackrel{}{k}(\eta )`$ obeys the following equation,

$$\frac{^2\beta _k}{\eta ^2}^\pm +\sigma \frac{\beta _k}{\eta }^\pm +\left(k^2\pm \lambda \frac{d\varphi }{d\eta }k\right)\beta _k^\pm (\eta )=0.$$
(5)
Some qualitative behaviour of the solutions of eq.(5), can be inferred from the simple case of a constant $`\frac{d\varphi }{d\eta }`$. Then the solutions are simply linear superpositions of two exponentials. Only if $`\lambda \left|\frac{d\varphi }{d\eta }\right|>k`$, one of the two modes ($`\pm `$) is exponentially growing. Otherwise both of them are either oscillating or damped, as in ordinary magnetohydrodynamics. To obtain significant amplification, coherent scalar field velocities $`\frac{d\varphi }{d\eta }`$ over a duration are necessary, larger velocities leading to larger amplification. The amplified mode is determined by the sign of $`\frac{d\varphi }{d\eta }`$. Modes with wavenumber $`k_{max}=\frac{1}{2}\lambda \left|\frac{d\varphi }{d\eta }\right|`$ get maximally amplified.

For values of the scalar mass in the TeV range and temperature above $`100`$ GeV, the cosmic friction term $`2aH\frac{d\varphi }{d\eta }`$ in eq.(5) is negligible compared to the mass term. The scalar field oscillates and its velocity changes sign periodically over a time scale much shorter than the Hubble time at the epoch, so that both modes can be amplified. Each mode is amplified during one part of the cycle and damped during the other part of the cycle. Net amplification results when amplification overcomes damping. It occurs for a limited range of Fourier modes, peaked around $`k/m\frac{1}{2}\lambda f`$. The modes of the EM fields are oscillating with (sometimes complicated) periodic time dependence and an exponentially growing amplitude. Total amplification is exponential in the number of cycles. For the range of parameters in which fields are amplified, the amount of amplification per cycle for each of the two modes is very well approximated by the same constant $`\mathrm{\Gamma }(k/m,\lambda f,\sigma )`$. In Fig. 1 we show: a) an example of the time dependence for a specific mode and a selected set of parameters, and b) amplification rates as a function of the wave number $`k`$.

Another interesting approximate solution can be obtained for $`mT`$. In this extreme limit, it would take a time interval of the order of the characteristic cosmic expansion time for the scalar velocity to change its value significatively. We say that the scalar rolls. Eq.(5) can be approximated by a first order equation $`\sigma \frac{\beta }{\eta }^\pm +\left(k^2\pm \lambda \frac{d\varphi }{d\eta }k\right)\beta ^\pm (\eta )=0`$, which can be solved exactly. The amplified mode is determined by the sign of $`\mathrm{\Delta }\varphi (\eta )=\varphi (\eta )\varphi (0)`$. Looking at $`\eta \eta _{EW}`$ the amplification is maximal for $`k_{max,EW}10^7T_{EW}`$.

A detailed discussion of the end of oscillations or rolling is beyond the scope of this paper. However, we do know that once the oscillations or rolling stop, the fields are no longer amplified and obey a diffusion equation. Modes with wave number below the diffusion value $`k<k_\sigma T10^8T`$, where $`\frac{k_\sigma ^2}{\sigma }\frac{1}{\eta _{EW}}=1`$, remain almost constant until the EW transition, their amplitude goes down as $`T^2`$, and energy density as $`T^4`$, maintaining a constant ratio with the environment radiation. Modes with $`k/T>k_\sigma `$ decay quickly, washing out the results of amplification. We have seen that for oscillating fields the range of momenta $`k^1>r_{Debye}`$ that get amplified is not too different than $`T`$, therefore scalar field oscillations have to occur just before, or during the EW transition. In that case, the amplified fields do not have enough time to be damped by diffusion. If the field is rolling, momenta $`kT`$ can be amplified, and therefore the rolling can end sometime before the EW transition.

## 4 Implications for Electroweak Baryogenesis

To obtain the average magnetic energy density in amplified fields we have to average over the initial conditions, which may result, for example, from initial thermal or quantum fluctuations. We assume translation and rotation invariance, $`B^\pm (\stackrel{}{k},\eta )B^\pm (\stackrel{}{\mathrm{}},\eta )^{}_{|\eta =0}=\delta ^3(\stackrel{}{k}\stackrel{}{\mathrm{}})\delta _+f^\pm (k).`$ Using the definition for magnetic power spectrum $`P_B^\pm (k)=4\pi k^3f^\pm (k),`$ the average energy density in amplified fields $`\rho _B=\frac{1}{8\pi }<B^2>`$, can be expressed as $`\rho _B=\frac{1}{8\pi }d\mathrm{ln}k`$ $`\left\{P_B^{}(k)\left|𝒜^{}(k,\eta )\right|^2+P_B^+(k)\left|𝒜^+(k,\eta )\right|^2\right\}`$, where $`𝒜^\pm (k,\eta )`$ are the amplification factors for each one of the two magnetic modes.

Since in our case $`E_\stackrel{}{k}^\pm (\eta )=\pm \frac{1}{k}_\eta B_\stackrel{}{k}^\pm (\eta )`$, and therefore $`E^\pm (\stackrel{}{k},\eta )B^\pm (\stackrel{}{\mathrm{}},\eta )^{}=\pm \frac{1}{2}\frac{1}{k}_\eta B^\pm (\stackrel{}{k},\eta )B^\pm (\stackrel{}{\mathrm{}},\eta )^{},`$ the CS number density stored in the amplified magnetic modes $`\mathrm{\Delta }n_{CS}=\frac{g^2}{4\pi ^2}\underset{0}{\overset{\eta }{}}𝑑\stackrel{~}{\eta }EB`$, where $`g^{}`$ is the hypercharge gauge coupling, can be related to $`\rho _B`$

$$n_{CS}=\frac{g^2}{\pi }d\mathrm{ln}k\frac{1}{k}\rho _B(k)\frac{\gamma _{AS}^B+\gamma _{AS}^\varphi }{1+\gamma _{AS}^B\gamma _{AS}^\varphi },$$
(6)
introducing the asymmetry parameters, $`\gamma _{AS}^B=\frac{P_B^{}(k)P_B^+(k)}{P_B^{}(k)+P_B^+(k)},`$ and $`\gamma _{AS}^\varphi =\frac{\left|𝒜^{}(k,\eta )\right|^2\left|𝒜^+(k,\eta )\right|^2}{\left|𝒜^{}(k,\eta )\right|^2+\left|𝒜^+(k,\eta )\right|^2}.`$ Further, we can relate the fractional energy density in coherent magnetic field configurations $`\mathrm{\Omega }_B(k)=\rho _B(k)/\rho _c`$ to the CS fractional number density $`n_{CS}/s`$, assuming the universe is radiation dominated and the amplification factors $`𝒜_{}^{\pm }{}_{k}{}^{}`$ are, as we have seen, sharp functions of $`k`$, peaked at $`k_{max}`$222As we have pointed out our description is valid for short wavenumber modes $`k<k_Dr_D^1`$. $`k_{max}`$ is understood to be the maximally amplified mode in the range $`0<k<k_D`$.,

$$\frac{n_{CS}}{s}0.01\frac{T}{k_{max}}\mathrm{\Omega }_B(k_{max})\frac{\gamma _{AS}^B+\gamma _{AS}^\varphi }{1+\gamma _{AS}^B\gamma _{AS}^\varphi }(k_{max}).$$
(7)
This Chern-Simons number will be released in the form of fermions which will not be erased if the EW transition is strongly first order, and will generate a baryon asymmetry , $`\frac{n_B}{s}=\frac{3}{2}\frac{n_{CS}}{s}`$. An equal lepton number would also be generated by the same mechanism so that $`BL`$ is conserved. Note that the fact that baryon number asymmetry is generated at $`k_{max}0`$ does not mean that baryon density is actually inhomogeneous on this short length scale $`L_{max}1/k_{max}`$. Comoving neutron diffusion distance at the beginning of nucleosynthesis is much longer than $`L_{max}`$, so that by that time inhomogeneities would have been erased by free streaming . If $`T/k_{max}`$ is not too different than unity, as we have seen for the case of oscillating field, and $`\gamma _{AS}^B`$ and $`\gamma _{AS}^\varphi `$ are small, it is possible to obtain $`\frac{n_B}{s}10^{10}`$ and have strong magnetic fields $`\mathrm{\Omega }_B1`$ present during the EW transition. If $`T/k_{max}`$ is large and $`\gamma _{AS}^\varphi `$ is order unity as we have seen in the rolling case, it is not possible to have strong magnetic fields without producing too many baryons.

## 5 Experimental Signatures of HyperAxions in Colliders

A pseudoscalar field that couples to hypercharge topological number density as we have described, could be produced and detected in colliders . Since the hypercharge photon is a linear combination of the ordinary photon and the $`Z`$, the hypothetical particle could be produced, in hadronic or leptonic colliders, in association with photons or $`Z`$’s through the channels $`f\overline{f}Z^{},\gamma ^{}ZX`$ and $`f\overline{f}Z^{},\gamma ^{}\gamma X`$, (here $`Z^{}`$ and $`\gamma ^{}`$ denote a virtual $`Z`$ and a virtual photon, and $`f`$ is a charged fermion). The produced particle would decay into two neutral gauge bosons, $`\gamma \gamma `$, $`\gamma Z`$ or $`ZZ`$. The experimental signature of these processes would , then, be a triplet of neutral gauge bosons produced in well defined ratios. The almost isotropic angular distribution of momenta of the outgoing bosons could help in separating it from the QED background, that is strongly peaked in the forward/backward directions.

In the simple case of a singlet elementary pseudoscalar whose only coupling to SM fields is to hypercharge topological density, the extra field can be described by two parameters with units of mass: the mass, $`m`$, and the inverse coupling, $`M=1/\lambda `$. The plots in Fig. 2 show our estimation of the regions of the $`(m,M)`$ parameter space where the particle signature could be separated from the background, and so the particle can be detected in present or future colliders. Below the curves, a sufficient number of events to allow detection are expected. Since the particle has not been detected yet, the first plot, that corresponds to LEPII and RunI of the Tevatron, can be used to rule out an interesting region of parameters. In future colliders detection capabilities will be increased significatively in the range of parameters relevant for baryogenesis, as shown in Fig 2.

## References