Léo Bosse
email: leo.bosse@univ-grenoble-alpes.fr

Jean Lilensten

Nicolas Gillet

Colette Brogniez

Olivier Pujol

Sylvain Rochat

Alain Delboulbé

Stéphane Curaba

Magnar G Johnsen

Night light polarization: modeling and observations of light pollution in the presence of aerosols and background skylight or airglow

' First polarized radiative transfer code for night light and atmospheric studies including Lorenz-Mie and Rayleigh polarized scattering.

' Application to aerosols, light pollution and thermospheric nightglow studies.

Introduction

The main auroral emissions are due to the atomic oxygen and molecular ion nitrogen N 2 . The former produces the red (630 nm) and green (557.7 nm) light at the altitudes of about 220 and 110 km respectively. The latter emits in a large band amongst which the most prominent emissions are the blue (427.8 nm) and the purple (391.4 nm) radiation, around an altitude of 85 km. They mainly originate in the collisions between precipitated electrons and the ambient atmosphere (see for example [START_REF] Banks | Aeronomy[END_REF].

These emissions are pronounced at high latitudes, in the auroral ovals. At other latitudes, emissions of the upper atmosphere are called nightglow [START_REF] Leinert | The 1997 reference of diffuse night sky brightness[END_REF]. The literature is abundant and here we only refer to some recent works. They exist at all latitudes [START_REF] Parihar | An investigation of the ionospheric F region near the EIA crest in India using OI 777.4 and 630.0 nm nightglow observations[END_REF]. Their origin is due to chemical reactions [START_REF] Plane | On the sodium D line emission in the terrestrial nightglow[END_REF] and collisions, either between neutral molecules or atoms, or between ambient ionospheric electrons and gases [START_REF] Tashchilin | Modeling nightglow in atomic oxygen red and green lines under moderate disturbed geomagnetic conditions at midlatitudes[END_REF]. The drivers may be gravity waves [START_REF] Vargas | Traveling Ionosphere Disturbance Signatures on Ground-Based Observations of the O( 1 D) Nightglow Inferred From 1-D Modeling[END_REF] or ionospheric currents [START_REF] Dymond | A Comparison of Electron Densities Derived by Tomographic Inversion of the 135.6-nm Ionospheric Nightglow Emission to Incoherent Scatter Radar Measurements[END_REF].

For the last decade, a serie of experiments have shown that the red emission, when measured from the ground, is polarized [START_REF] Lilensten | The thermospheric auroral red line Angle of Linear Polarization: Auroral angle of linear polarization[END_REF], and references herein). The direction of polarization for this red line was shown theoretically to be parallel to the magnetic field [START_REF] Bommier | The theoretical impact polarisation of the O I 6300 Å red line of Earth Aurorae[END_REF]. Correlations between variations in the magnetic field and in the angle of linear polarization (AoLP) had indeed been observed experimentally, but not systematically, with a first prototype of auroral polarimeter [START_REF] Lilensten | Polarization in aurorae: A new dimension for space environments studies[END_REF].

Recently, a new nightglow polarimeter has been developed in order to observe low radiant flux emissions (i.e. not limited to bright aurorae). It also allows to target simultaneously several atmospheric emissions (for a full description, see [START_REF] Appendix A Bosse | On the nightglow polarisation for space weather exploration[END_REF]. Three findings by Bosse et al. (2020) lead us to reconsider our first understanding of the upper atmosphere polarization:

' All of the four auroral emission lines appear polarized when measured from the ground, ' Although in several circumstances this polarization appears to be linked to the local magnetic activity and to the state of the ionosphere, it is far from being systematically aligned on the apparent direction of the magnetic field, as foreseen theoretically, ' Light pollution from nearby cities significantly impacts, via scattering, the polarization measurements.

These series of observations questioned the geophysical origin of the polarization: how much is it affected by light pollution scattering in the lower atmosphere? Are auroral lights polarized at the emission or during their propagation toward to the instrument? To answer these questions, it is necessary to develop a polarized radiative transfer code able to account for sources potentially spread all over both the sky and the ground.

In this paper, we describe the first code developed to this purpose, called 'POMEROL' (standing for 'POlarisation par Mie Et Rayleigh des Objets Lumineux').. It takes into account the light emitted in the upper atmosphere, as well as all possible pollution sources on the ground. We present some examples of applications, with comparisons to geophysical data obtained in the French Alps. In this paper the code is restricted to a single scattering approximation (SSA). The goal here is to prove the feasibility of night-time polarisation measurements. Thus, the model is currently kept as simple as possible. We show that even in this relatively simple configuration, it is possible to reproduce quite well the measurements and to deduce physical parameters behind the polarization. However, we also discuss the limitations of this first approach, and what improvements could be brought by considering a polarized multiple scattering radiative scheme. In this first paper about night-time polarization, we think it is useful to know what the single scattering approximation can give before diving into the complexity of multiple scattering. We will see that the SSA can already catch the essence of the physics at hand. Our study paves the way to future investigations of the night-light polarization, and calls for evolution of the model, including in particular multiple scattering. Such a code is not designed only to better understand the upper atmosphere. It could also help to characterize aerosols in a passive experimental way in the absence of the moon or the sun. It may also lead to more accurate monitoring of light pollution, which represents a growing concern over the last decades. The area covered by light from human origin is spreading, along with its impacts on life [START_REF] Grubisic | Insect declines and agroecosystems: Does light pollution matter[END_REF], energy consumption [START_REF] Kyba | Artificially lit surface of Earth at night increasing in radiance and extent[END_REF] or astronomy. Among the literature displaying a varied list of negative effects, we can cite the insect population decline and the "Ecological Armageddon" [START_REF] Grubisic | Insect declines and agroecosystems: Does light pollution matter[END_REF] or health issues (Garcia-Saenz et al., 2018; Zielinska-Dabkowska, 2018). Most of the studies concentrates on the exposure rate and magnitude. However, we lack studies on light pollution polarization [START_REF] Horvath | Polarized light pollution: A new kind of ecological photopollution[END_REF][START_REF] Kyba | Lunar skylight polarization signal polluted by urban lighting[END_REF]). Yet it plays an important role for vast groups of insects that use it to navigate. Below, we briefly describe the instrument in Section 2. We then present the principles of the radiative transfer model in Section 3.2 (the full details are given in Appendix A). Next we provide in Section 4 a series of synthetic experiments in order to assess the influence of the input parameters: effect of a localized source on the ground, of atmospheric properties (ozone, aerosols) and of multiple scattering. In Section 5, we compare the measurements from an experimental campaign at mid-latitude to the model outputs. Finally, we discuss our findings in Section 6.

Description of the polarimeter

Here is a brief overview of the experimental set-up. The polarimeter used in this study has been fully described in Bosse et al. (2020). We therefore only recall here its basics.

The incoming light along the line of sight is filtered through a narrow optical filter (of 2 nm width for the red line and 10 nm width for all other lines). Behind it a polarizing lens rotates at 2 Hz. The light passing the lens then hits a photomultiplier and is converted into an electrical current with a 1 kHz sampling rate. Data are smoothed over a given time window during which the polarization is assumed constant (10 seconds for all data presented here). A lock-in analysis is performed in real time. This powerful method allows a fast and accurate computation of the polarization. However, when the degree of linear polarization (DoLP) becomes too small, typically below 0.5%, the AoLP can hardly be computed and becomes very noisy. Both the DoLP and AoLP have been calibrated, but not the radiant flux.

We note F 0 the incident radiant flux received at a given wavelength. The DoLP ranges between 0 and 1 (or, in the figures, 0 to 100%), and we define the AoLP with respect to the vertical (0 ˝being upward, ˘90 ˝horizontal as it is π-periodic). ϕ t is the angle of the polarizing filter with the vertical at time t. We suppose F 0 , DoLP and AoLP to not change during one rotation of the polarizing filter. From basic optics, the radiant flux passing through the polarizing filter can be decomposed in two parts: a polarized one that varies as cos 2 pϕ t ´AoLPq, and an unpolarized part, assumed to be constant over one rotation. For an incident radiant flux F 0 , the polarized and unpolarized fluxes are

" F pola 0 " DoLP ˆF0 F unpola 0 " p1 ´DoLPq ˆF0 . (1)
Following Malus Law, after the polarizing filter they become

" F pola t " F pola 0 cos 2 pϕ t ´AoLPq F unpola t " F unpola 0 {2 . (2)
The 1{2 factor on the unpolarized radiant flux comes from the averaging of Malus law on all

AoLP. The radiant flux measured at time t can therefore be written as:

F t " F pola t `Funpola t " F 0 ˆDoLP ˆcos 2 pϕ t ´AoLPq `1 ´DoLP 2 ˙. (3)
Over one rotation of period T r , this allows computing the Stokes parameters in spherical coordinates as:

$ ' ' ' ' ' ' & ' ' ' ' ' ' % I " 1 T r ż T r 0 F t dt Q " 1 T r ż T r 0 F t cos 2ϕ t dt U " ´1 T r ż T r 0 F t sin 2ϕ t dt . ( 4
)
Note that we do not consider circularly polarized light, such that the last Stokes parameter V " 0. Injecting (3) in ( 4), one deduces the polarization parameters:

$ ' ' ' & ' ' ' % F 0 " 2I DoLP " 2 I a Q 2 `U2 AoLP " 1 2 arctan ˆU Q ˙(5)
The data may be smoothed over time by averaging I, Q and U over the desired number of rotations and then calculating the corresponding polarization values.

The radiative transfer model

In order to interpret the data from the instrument described in Section 2, we need to solve the polarized radiative transfer equations. This is the goal of the POMEROL code. We describe its inputs in Section 3.1 while the polarized radiative transfer equations are summarized in Section 3.2 and fully described in Appendix A. For the purpose of the present study, our analysis is based on single scattering. The limits of this approximation are discussed later on in Section 4.4. We aim at modeling the polarization observations under different configurations. We thus consider potentially several sources of light: direct light (from the nightglow, or the star-light), single scattered light (e.g. from cities, or auroral lights at high latitude).

Inputs of the model

Instrument related entries

The experimental characteristics are the first inputs, in particular the surface Σ " 20 cm 2 of the detector and its half aperture angle ε " 1 ˝, as well as its geographical position (latitude, longitude and altitude) and its pointing direction. This latter is defined by the elevation e (angle between the horizontal and the line of sight) and the azimuth a, reckoned with respect to the North, positive Eastward, i.e. clockwise rotation. One can specify discrete azimuths and elevations, or span over an almucantar (i.e. a full rotations in azimuth at a constant elevation) automatically.

Cities and pollution map

The model also takes as input a light pollution map. These are ground images of Earth at night produced by the NOAA Earth Observations Group, using satellite data from the Visible Infrared Imaging Radiometer Suite Day/Night Band. These images sum the emissions for wavelengths from 500 to 900 nm [START_REF] Mills | VIIRS day/night band (DNB) stray light characterization and correction[END_REF][START_REF] Appendix | Illuminating the Capabilities of the Suomi National Polar-Orbiting Partnership (NPP) Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band[END_REF]. The maps are processed by their authors to remove ephemeral light and are averaged over one year. They are provided by their authors in two modes, with the minimum emission set to zero or not. At high latitudes, auroral light are not always removed with other ephemeral lights, thus emission maps are overestimated. To better suppress the auroral contribution to the ground emission map, we use the second mode. At mid-latitude, the difference is minimal and we use the same mode for consistency. We use the most recent release from 2016. The ground emission are in units of nW/m 2 /sr.

The map is centered on the instrument in polar coordinates (see Figure 8). The maximum range from the instrument (typically up to 200 km) and the number of bins (the best resolution being that of the map, namely 46 m) are adjustable parameters. To reduce artifacts due to digitization, the size of the bins increases with the square root of the distance to the instrument.

We can also consider instead synthetic emission maps. For example, these may consist of a point source of given radiance at a given distance, azimuth and elevation with respect to the instrument, or a uniform emission map of given radiance.

Natural background

In the following, the natural background light designates any source of light from the sky that we approximate as constant and isotropic. It includes two main contributions: the nightglow [START_REF] Leinert | The 1997 reference of diffuse night sky brightness[END_REF]) and the integrated starlight (e.g. [START_REF] Staude | Scattering in the earth's atmosphere: Calculations for milky way and zodiacal light as extended sources[END_REF]. The nightglow is specific as it is well defined in wavelengths (see Table 2). Its emissions can change with time and it can be highly structured. All-sky camera images, which are unfortunately not available in the present study, could be used to model the nightglow more precisely. The integrated starlight covers a wide spectrum over all observed wavelengths and is unpolarized. It varies with the time and place and could be recovered with astronomical tools in up-coming studies. For now on, since the absolute values of these contributions cannot be estimated, and also because the absolute radiant flux measured by the instrument is not calibrated, we cannot disentangle the several background contributions from the confrontation of our measurements with the model outputs. Thus, the background contribution may be considered to account for the model uncertainties. These are different for each wavelengths, time and location of observation.

Atmospheric properties

To compute the Rayleigh scattering in the atmosphere, we need different atmospheric parameters. We use the 2001 MIPAS Model Atmospheres [START_REF] Remedios | Mipas reference atmospheres and comparisons to v4. 61/v4. 62 mipas level 2 geophysical data sets[END_REF] up to the lower thermosphere, at 120 km of altitude. It provides the temperature T pzq and pressure Ppzq vertical profiles (z is the altitude) as well as an ozone vertical profile. For the present purpose, we use their standard night-time mid-latitude profiles, displayed in Figure 1.

To account for the influence of the aerosols, we consider a Lorenz-Mie scattering model [START_REF] Lorenz | Lysbevaegelsen i og uden for en af plane lysbolger belyst kugle". det kongelige danske videnskabernes selskabs skrifter[END_REF][START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF][START_REF] Born | Principles of optics : electromagnetic theory of propagation, interference and diffraction of light[END_REF][START_REF] Van De Hulst | Light scattering by small particles[END_REF]. This implies a wide range of input parameters, such as the complex refractive index, the aerosol sizes and their vertical profiles. Three aerosol models are considered in this paper, named 1-low, 2-high and 3-mid, whose parameters are listed in Table 1. We use complex optical indices from [START_REF] Dubovik | Accuracy assessments of aerosol optical properties retrieved from Aerosol Robotic Network (AERONET) Sun and sky radiance measurements[END_REF] assumed to be the same for all wavelengths considered in this paper. The size distribution nplnprqq of the aerosol is supposed to be log-normal,

nplnprqq " dN d lnprq " N ? 2π lnpσq exp ˆ´ln 2 pr{Rq 2 ln 2 pσq ˙, (6)
with r the aerosol radius (in µm), N the total number of aerosols, R the mode radius (the radius where the distribution is maximal) and lnpσq controls the dispersion of the aerosol sizes around R. The vertical number density distribution of aerosol npzq is given by

npzq " n 0 " exp ˆ´z H ˙`ˆn B n 0 ˙ , (7)
with n 0 the number density at the surface (in cm ´3), n B the minimum number density when parameters are taken from Jaenicke (1993). In the context of the present paper, we consider these standard aerosol profiles to be sufficient in order to illustrate our purpose. However, in future studies they could be recovered from independent measurements (e.g. LIDAR) or using an inversion scheme in order to best fit polarization data.

z

Topographic map

Our model also incorporates the topography around the instrument. It is used for the computation of occultation: when a light beam hits the ground between the emission and the scattering point, its contribution to the light received by the instrument is set to zero. The topographic altitude data are taken from ALOS GLobal Digital Surface Model AW3D30 DSM of the Japan Aerospace Exploration Agency [START_REF] Tadono | Generation of the 30 m-mesh global digital surface model by alos prism[END_REF]. It has a resolution of 30 m, which can be downgraded in order to reduce the computation time. In a first step, we describe the properties of the light (at a given wavelength) reaching our modeled instrument located at point A (hereafter called virtual instrument), coming from a single emission point source E, and scattered at a single point R along the line of sight (see Figure 2).

In this scenario, we shortly discuss here the computation of the radiant flux F A measured by the instrument, as well as the DoLP and AoLP of the light scattered towards it.

Firstly, the radiant flux F A measured by the virtual instrument depends on many physical quantities described in detail in Appendix

DoLPpθq is always taken as positive. The sign of DoLP aer gives information on the direction of polarization of the light. This is not useful here because the direction of polarization is given by the AoLP in the following work. To keep consistency with the real instrument, we set the zero of the AoLP along the vertical direction, and thus AoLP " ˘90 ˝in the horizontal direction, when the virtual instrument points horizontally. The AoLP is positive in the trigonometric direction (anti-clockwise) when looking towards the source. The AoLP is π periodic, such that we define it in the interval r´90 ˝, 90 ˝s.

Integration over all sources and along the entire line of sight

Once the light polarization characteristics for a single path are computed (the radiant flux F using equation A17, the DoLP and the AoLP using equation 9), the model repeats this computation for every point R along the line of sight and every point source E given in input.

We then integrate along the line of sight and over all point sources by summing all these contributions. However, this integration is not straightforward, since the DoLP and AoLP are not linear quantities. For example, the resulting DoLP from the mix of two light rays with different DoLPs is not their average if the AoLPs are different. In order to maintain the same formalism between the virtual and the real instruments, we use the Stokes parameters in spherical coordinates [START_REF] Mcmaster | Polarization and the Stokes Parameters[END_REF][START_REF] Van De Hulst | Light scattering by small particles[END_REF], which are equivalent to equation (4):

I " F 2 , Q " F DoLP 4 cosp2 AoLPq , (10)
U " F DoLP 4 sinp2 AoLPq .

When integrating over the line of sight (l.o.s.) and over all the sources given in input, we can add up the Stokes parameter and get back to the resulting DoLP and AoLP via the relations given in equation ( 5). In the end, the final polarisation parameters are:

I total " ÿ sources ÿ l.o.s. I , Q total " ÿ sources ÿ l.o.s. Q , (11)
U total " ÿ sources ÿ l.o.s.
U .

Model validation and influence of the parameters

In this section we study:

i) The influence of a localized source on the ground (such as a cities), at several distances from the instrument, and in the absence of any source in the sky ;

ii) The influence of the atmospheric parameters (absorption by O 3 , aerosols) for a localized source on the ground;

iii) The effect of a uniform source in the sky, polarized or not, in the absence of any source on the ground.

iv) The domain of validity of the single scattering approximation (in comparison with multiple scattering) by comparing POMEROL outputs with modeling and observations provided by [START_REF] Pust | Comparison of skylight polarization measurements and MODTRAN-P calculations[END_REF].

These reduced configurations have been chosen so as to illustrate the most important factors that impact ground-based polarization measurements. In the following, we present results by means of clockwise almucantars of elevation 45 ˝, the starting direction being the North. The virtual instrument parameters (Σ and ε) correspond to that of the real polarimeter (see Section 3.1.1). We show here results at λ " 557.7 nm (auroral green line, Table 2). In the following series of tests, we use a polarized Lorenz-Mie scattering model for scattering by aerosols (see Section 3.1.4). The comparison between the aerosols models and real data is shown later on in Section 5.1.2. The atmospheric profile used for the tests is the MIPAS standard night-time mid-latitude profile (see Section 3.1).

4.1 Influence of a localized source on the ground in absence of any source in the sky.

We consider an isotropic point source on the ground of radiance 100 nW/m 2 /sr, at different distances d AE away South of the virtual instrument. The purpose of this arbitrary setup is only to showcase the influence of a point source on the ground at different distances from the instrument. Thus, the radiance is arbitrary and the distances are chosen to be representative of a real environment (between 1 and 100 km). There is no mountain obstruction. The ground surface takes the Earth curvature into account. Only Rayleigh scattering and ozone absorption are taken into account, no aerosols are present in the atmosphere. The ozone is taken into account through its number density vertical profile (Figure 1) and its absorption cross sections as a function of the wavelength (Figure A1). For a source located 100 km away from the observation point, the ozone decreases the measured radiant flux by less than 3%.

There is no effect on the AoLP since this angle depends only on the scattering plane (defined by the emission and scattering direction). The effect on the DoLP is smaller than 0.4%, and the closer the source, the lower the effects. More complex sources than a point may slightly increase these values. From now on, we take into account the impact of O 3 .

In Figure 3, we show the polarization results for d AE " 1, 5, 10, 50 and 100 km. In all such figures, the upper panel shows the radiant flux measured by the virtual instrument, representing the energy per unit time on the collector, the middle panel shows the DoLP and the lower panel the AoLP.

The further the source, the lower the measured radiant flux. However, the effect is not merely a decrease as 1{ ? d AE , because the source illuminates points at all altitudes along the line of sight that are integrated on the virtual instrument. Moreover, the light parameters change during its crossing into the atmosphere from the source to the line of sight, and along the line of sight. As seen in Figure 1, the concentration in ozone peaks at about 20 km. As an example, the decreasing factor between a source lying at 1 and 10 km is about 25 (and not 10 2 if only the effects of the distance were considered). It becomes approximately 5 000 (and not 10 4 ) with a source at 100 km compared to 1 km. Because of the exponential decrease of the atmospheric pressure with the altitude, this effect is less important when the source moves away from the polarimeter. However, in all cases, the radiant flux is maximum in the Southern direction, i.e. toward the source.

Geometrically, for a single point source and a given line of sight, the AoLP of every scattering point is the same. In this particular case, we do not have to use the I, Q and U notation (equations 11). The total DoLP integrated along the line of sight is an average of the DoLPpθ i q of each scattering points i weighted by the scattered radiant flux F i as DoLP "

ř i F i DoLP i pθ i q ř i F i , (12)
where θ i is the scattering angle at point i, F i follows equations (A17) and DoLP i equation (9).

In the limit of a source infinitely far away from the instrument (on a flat surface), all paths from the source to the instrument are parallel. Thus, the scattering angle is the same everywhere on a given line of sight. The total DoLP measured by the virtual instrument should then follow equation (8) (dashed line in Figure 3, middle) since there is no aerosols in this case. When looking eastwards and westwards, the value of the DoLP reaches 100% as expected from the theory. When the source gets closer to the virtual instrument, each point along the line of sight has a different scattering angle, which smooths the variations along the almucantar. On a flat surface, the AoLP is the same whatever the distance from the source (Figure 3, bottom).

Influence of the aerosols for a localized source on the ground

The aerosols are taken into account through their number density vertical profiles and their cross sections (Table 1). The aerosol cross section depends on the wavelength as computed using Lorenz-Mie theory (Lorenz, 1890; [START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF][START_REF] Born | Principles of optics : electromagnetic theory of propagation, interference and diffraction of light[END_REF]. For the sake of simplicity, we consider the aerosol refractive index to be the same for all wavelengths. A variation of the refractive index with wavelength could be taken into account, but is out of the scope of the present study.

In order to illustrate the effect of the aerosols, we consider three different profiles listed in Table 1: South of the instrument. The continuous line is the case without aerosols (similar to Figure 3). The dashed lines correspond to three different aerosol models listed in Table 1: 1-low (black), 2-high (blue) and 3-mid (red).

For an increasing aerosol contribution, the DoLP decreases in the direction of the source.

With the 1-low model, it decreases by about 50%, with the two maxima along the almucantar still present. The AoLP does not change since the aerosol radius is small compared to the wavelength, and the polarization direction is the same as for Rayleigh scattering (AoLP ray ).

However, when the aerosol size is large compared to the wavelength (2-high and 3-mid), the DoLP behavior changes drastically, with a maximum at about 25% when pointing away from the source. The AoLP is shifted by 90 ˝when pointing away from the source due to larger aerosol size, while it is the same as AoLP ray when pointing towards the source. The impact of aerosols on the virtual instrument observations is complex. They may for instance either increase or decrease the intensity depending on their size or on the scattering angle. They do therefore play a major role. Spanning all possible aerosol models is not a crucial point for this study, and we limit ourselves to the three models presented here.

Influence of the skylight with different polarization parameters

We document here the effect of a simple skylight (nightglow and/or integrated star light).

It is modeled as an infinitely thin uniform sky emission of arbitrary intensity at 110 km of altitude. The model computes the sum of the non scattered light and the single Rayleigh scattered light measured by the instrument. In order to prepare for further auroral studies (out of the scope of this article), we allow this emission to be polarized [START_REF] Bommier | The theoretical impact polarisation of the O I 6300 Å red line of Earth Aurorae[END_REF]. In cases where the nightglow is polarized, the DoLP depends on the angle α between the line of sight and the polarization direction, as DoLP " p 0 `p1 sin α. For the EW and NS cases, p 0 « 0.85p source and p 1 « 0.1p source . When the nightglow is polarized along the magnetic field lines, the behaviour is much different as the lines are close to vertical, with p 0 « 0.6p source and p 1 « 0.3p source . For a uniform sky emission and in the absence of other sources, Figure 5 shows that the polarization from Rayleigh scattering in the atmosphere is relatively small compared to the polarization of the non scattered light. Thus, in this simple configuration, any polarized source in the sky with a DoLP higher than the instrument noise level should be detectable.

The AoLP (lower panel) for the unpolarized nightglow emission (continuous black line)

is undefined (it is set to zero by default, with no physical meaning). The EW and NS cases show an AoLP regularly rotating over all 360 ˝with the same behavior, but shifted in azimuth by 90 ˝. As expected from the DoLP definition, the maximum DoLP corresponds to a 90 AoLP. When the nightglow polarization direction is aligned along the magnetic field, the magnitude of the AoLP variations along an almucantar is significantly weaker, within ˘40 ˝, which led Bosse et al. (2020) to reject this single source as the origin of the measured polarization of auroral lights.

On the impact of multiple scattering

The model described here does not consider the effects of multiple scattering. We recall that the main goal of this approach is to test the feasibility of night-time polarisation measurements and modelisation, and not to develop a new state of the art radiative transfer code.

However, these effects can not be neglected, in particular at the shortest wavelengths, and we document here their impact. The literature on the atmospheric polarization is already large, both in term of observation and modeling. The reader may refer to the seminal work of Hovenier (1971) and [START_REF] Hansen | The doubling method applied to multiple scattering of polarized light[END_REF] on the modelisation of terrestrial and planetary are not erased by MS [START_REF] Hansen | Multiple scattering of polarized light in planetary atmospheres part ii. sunlight reflected by terrestrial water clouds[END_REF]. Generally, in comparison with single scattering, MS tends to enhance the radiant flux and reduce the DoLP, with relative effects that depend on the optical depth, the solar elevation and the wavelength [START_REF] Hovenier | Multiple Scattering of Polarized Light in Planetary Atmospheres[END_REF]). An order of magnitude of the MS effect is proposed for instance by [START_REF] Staude | Scattering in the earth's atmosphere: Calculations for milky way and zodiacal light as extended sources[END_REF] who mentioned an increase in the radiant flux ranging from 10 to 45%, and a drop in the DoLP from 5 to 20%, depending on the optical depth. It is acknowledged that simply using a scaling factor to accurately match DoLP and radiant flux observations with single scattering is over-simplistic [START_REF] Hansen | The doubling method applied to multiple scattering of polarized light[END_REF]). In atmospheric optics, the SSA is more of an issue towards short wavelengths (e.g. following a forest fire with high aerosol content and finally shortly after the fire, a day with moderate aerosol content. For each day, they measured the DoLP of the sky with an all-sky camera at different wavelengths, and compared the maximum DoLP found in the sky with their model prediction, as a function of solar elevation. They found that single scattering overestimates the DoLP, while it is under-predicted by unpolarized MS. However, in low aerosol conditions, the SSA is sufficient to reproduce their data for long wavelength (ą 630 nm).

We report their measurements in Figure 6 at three representative wavelengths (630, 530 and 450 nm) and for the two extreme days with low and high aerosol content. We see a decrease of the maximum DoLP with increasing solar elevation. This effect is more important for high aerosols conditions and stronger for large elevations and short wavelengths. To reproduce the above configuration with POMEROL, we set a point source at a given elevation representing the Sun, and compute the DoLP for each observation direction in the sky. We then take the maximum DoLP, which is always found in a direction around 90 ˝of the Sun for low aerosol contents. No ground reflections are taken into account. We use the 1-low aerosol profile for the day before the fire, and the 2-high aerosol profile for the time just after the fire. Note that the vertical concentration profiles have been tweaked so that the optical depth matches that of the aerosol profiles considered by PS11.

In the first case with low aerosol content, POMEROL overestimates the maximum DoLP for short wavelengths: at 450 nm, the modeled DoLP is " 20% too high (in relative values).

At 630 nm however, single scattering predictions from POMEROL are much closer to the observations, although our model does not reproduce the " 10% relative increase seen in the maximum DoLP for low elevation angles. This result is overall consistent with the behavior of MODTRAN-P by PS11 in similar conditions: single scattering model with a low aerosol content (see Figure 7). The several low aerosol models considered by PS11 give a hint of the sensitivity to this parameter, with a spread of the predicted maximum DoLP that ranges from 5% to 8% from 450 to 630 nm (i.e. less than " 10% in relative values).

In the case with high aerosol content the observed maximum DoLP decreases from around 45% at low solar elevation down to " 30% at 630 nm and it decreases frrom 35% to 15% at 450 nm for solar elevation above 20 ˝. This trend is overall recovered with POMEROL that shows a net decrease in the DoLP at solar elevation above 20 ˝in all wavelength (see Figure 6). At 450 nm, the maximum DoLP is well reproduced at high and low solar elevation, but is overestimated in between. At 630 nm and 530 nm, the DoLP is overestimated for all angles, with larger differences around 20 ˝solar elevation. In comparison, the single-scattering model predictions by PS11 overestimate the maximum DoLP by a factor as large as " 2 at 450 nm (see Figure 7). Here again the difference with our model might come from distinct aerosol profiles used in both studies.

As shown with the above example, POMEROL tends to overestimate the DoLP for a day-like environment, i.e. a single point source at infinity. This effect is sensitive to the aerosol profile, the wavelength and whether or not MS is considered in PS11. The benefits that we could gain by implementing MS is expected to be larger at short wavelengths. Meanwhile, we notice that the above test situation differs from our experimental environment in several ways. During the day, the Sun is approximated as a single point source because it dominates other light sources. By night on the contrary, light pollution from urban lightning and spread emissions from the integrated star light and airglow are mixed all together. In this case, the measured polarization depends heavily on the relative contributions of the different sources rather than on the absolute contribution of each one. We thus expect MS impact on the DoLP and AoLP to decrease for night-time scenarios.

MS can not be neglected in most cases, particularly at short wavelengths («ă 600 nm).

The conditions of our experiment corresponds to a low aerosol content. In this case, only a moderate sensitivity to the source angle has been reported, and we can expect MS to decrease the DoLP by at most about 20% at 450 nm, and less for longer wavelengths. MS is thus expected to result in a more isotropic contribution than single scattering, with a lower DoLP.

For the sake of simplicity, we then incorporate its effect together with that of the uniform, isotropic and unpolarized emission in the sky (both tend to reduce the DoLP and increase the radiant flux). One consequence of this coarse approximation is that we cannot isolate unambiguously the signature of the airglow and/or starlights. Still, this parsimonious approach appears to be enough to account for a significant part of our observations for wavelengths down to the blue line at 427.8 nm, as we shall see later.

Comparison with geophysical measurements at mid-latitude

In order to validate our model, we performed a series of observations in the French Alps.

In the following, we focus on the night from 19 to 20 January, 2021 when the moon was down (below ´10 ˝elevation during the whole observation). The latitude is 45.2123 ˝and longitude 5.9369 ˝. The altitude is 770 m. The nearest city is Grenoble 15 km away with downtown at an azimuth of 260 ˝. However, the valleys around this bright city produce also some light pollution. Figure 8 shows the geographical configuration with lines of constant elevation at 500 m, 1000 m and 2000 m.

The light pollution was clearly visible with the naked eyes. At the time of the observations, the snow was covering the ground above about 500 m elevation. The altitudes below h « 1500 m are mostly covered with forests in the mountains which lowers the albedo significantly, even though some snow remains on the trees. The air temperature was below 0 ˝C all over the surrounding area, and reached about ´5˝C at the point of observation. The modeling below takes the relief into account (see Section 3.1). We could observe simultaneously five wavelengths, summarized in Table 2: three of them concern emissions bands emitted in the upper atmosphere, either in the thermosphere (green, due to the O 1S excitation state) or in the thermosphere (blue and purple, due to the 1 st N 2 negative band). The "orange" band around 620 nm is not present in the thermosphere, but is emitted in the mesosphere by OH(5-0) and OH(9-3), though with a relatively lower radiant flux [START_REF] Broadfoot | The airglow spectrum, 3100-10,000 A[END_REF] We show below the results of an almucantar at a constant elevation e " 45 ˝from 1:18 UT to 3:35 UT taken with the real instruments and compared to the virtual one (i.e. the modeling). We rotate clockwise from North to the East and back to North with 10 ˝increments in azimuth. At each step in azimuth, we record the radiant flux, DoLP and AoLP during about 1'30". We recall that our instrument provides the relative (uncalibrated) measurements of the radiant flux, and the calibrated DoLP and AoLP. The orientation angles e and a are set using stars as references, together with a compass (to help on the increments in azimuth) and an inclinometer (for the elevation). We consider the accuracy to be of the order of ˘2˝. We focus first on a single wavelength (the green line, see Section 5.1) for which we detail our analysis of the observations. We next discuss the full picture in Section 5.2.

Model predictions versus observations in the green line

Figure 9 shows an almucantar in the green line. The radiant flux maximizes above downtown Grenoble (a " 260 ˝, this is clearer in Figure 10). In the theoretical case of a point source, the DoLP maximizes at ˘90 ˝of the direction of this point (see Section 4 and Figure 3). Here, the DoLP shows a single maximum near 240 ˝in azimuth, illustrating the influence of an extended pollution source. The value of its measured maximum is around 12%. The AoLP rotates regularly with a behaviour very similar to that of a point source on the ground westward with respect to the instrument.

In the absence of aerosols and background skylight

We first run our POMEROL model with the light pollution as single input (i.e. with neither aerosol nor natural background skylight). We use a ground emission map composed of 1000 pixels mapping an area of 100 km radius around the instrument. The elevation map 1): without aerosols nor background lights (black ‹); without aerosols and with background lights adjusted to fit the radiant flux variations (red `); without background lights and with aerosols model 1-low (blue ˆ), 2-high (magenta ) and 3-mid (purple ). Abscissa are given as the pointing direction under the bottom panel (azimuth where 0 is North, 90 is East) and the corresponding time in UTC format above the top panel.

covering the same area is used to model the mountain obstructions, with a 30 m resolution.

Note that the emission maps are not intercalibrated with our instrument, so that the units for the radiant fluxes are arbitrary. The ozone is taken into account, although its effect remains small. The results are shown in Figure 9 (black stars). Both the modeled and measured radiant fluxes peak around the same azimuth (260 ˝). However, the radiant flux variations along one almucantar are significantly larger for the model: the maximum radiant flux is 110% higher than the minimum for the model, with respect to only 43% for the measurement.

The modeled DoLP reaches 47% at its maximum around 200 ˝azimuth, about 4 times higher than the highest measured DoLP, and is offset by 40 ˝in azimuth. A second local maximum at 300 ˝azimuth is also present in the model, and absent from the measurements.

Two maxima in the DoLP are indeed expected for a point source (see Figure 3 This first simple approach does not take into account one or more additional sources.

We will now show the effects of introducing either a natural background light or the aerosols.

Impacts of aerosols and background skylight

The fact that the model overestimates the radiant flux variations and the DoLP in the case above is an indication that non-polarized (or little-polarized) additional sources must be taken into account. Such sources exist naturally: they are the nightglow and the integrated star light. This trend can also be due to the lack of multiple scattering in the model (see Section 4.4). The exact contributions of these effects is delicate to evaluate at this point, so we fit them as a whole called a background skylight in such a way that the relative variations of the modeled radiant flux match the measured ones. This method constitutes therefore an indirect estimate of this background skylight (see red `in Figure 9). Here, it is considered constant, isotropic and unpolarized (these approximations could be reconsidered in the future by estimating each contribution from independent measurements or more accurate models, see Section 3.1.3). Averaged over an almucantar, this results in approximately doubling the total radiant flux received by the instrument. This is not visible in Figure 9 because of the Arbitrary Unit radiant flux scale. As expected, the DoLP is decreased in any direction by about 50%, peaking at 25% around 220 ˝in azimuth. The two maxima in the modeled DoLP over one almucantar are still present, although the secondary maximum appears attenuated. The AoLP is unchanged in comparison with the previous case with no background sky, which is expected since the natural background is not polarized.

We now turn to a case with no background skylight, but with the aerosol polarized models included (see Section 3.1.4). The amplitude of the modeled radiant flux variation along the almucantar become larger, whatever the aerosol model, as seen with the synthetic case in Figure 4. The higher the aerosol contribution, the higher the radiant flux variations.

The modeled DoLP decreases when the aerosol contribution increases, so that it is possible to find a model (2-high: magenta ) that matches the measured DoLP. However, this makes the modeled AoLP depart from the measurements, so that it is not possible to find a set of aerosol parameters that allows fitting the measured radiant flux, AoLP and DoLP at once.

Comparison of the observations and the modeling in five wavelengths.

From the results presented above, we deduce that even though both additional sources (background and aerosols) significantly affect the polarization parameters, none of them reproduce the measurements alone. We consider below the model predictions for a combination of those two contributions. We first have to choose an aerosol profile. We use the 1-low aerosol model, as it fits best the measurements when combining with background skylight. It could already be intuited from Figure 9 as the 1-low profile affects less the AoLP than the other ones.

This choice is also coherent as we observe in the Alps at an altitude where the air is very clean.

We show hereafter model outputs with this aerosol profile together with a natural background skylight. At all wavelengths, the modeled radiant flux is the highest in the direction of the city of Grenoble, and the model predictions for the AoLP convincingly fit the observed values (see

Figures 10 to 14).
The magnitude of the signal from the ground is likely very small in the purple (391.4 nm) and turquoise (413 nm), because in these wavelengths the radiant flux of city lights is in average very weak (see for instance the spectral analysis of Fig. 6 in Bosse et al. (2020) and references herein). The DoLP observed for almucantars in these two lines are therefore associated with large error bars. This in turn also affects the uncertainty level for the AoLP.

The changes of their DoLPs with the azimuth appear rather flat, but are also less tightly constrained. Combined with the stronger limitation of the single scattering approach towards short wavelengths, it explains that the model fails at predicting the low values of the DoLP observed in the direction of Grenoble at these two wavelengths. Furthermore, ground emissions maps are less sensitive in these wavelengths [START_REF] Appendix | Illuminating the Capabilities of the Suomi National Polar-Orbiting Partnership (NPP) Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band[END_REF], so that the map used in the model might also differ from the real emissions.

Examination of Figure 8 shows that, without any contribution from the background or aerosols, significantly larger values of the DoLP are modeled in the green line (557.7 nm). The blue (427 nm) and the orange lines (620 nm) behave similarly. For these three wavelengths, the measured radiant fluxes are fitted when adding the 1-low aerosol model and respectively 114% (Figure 10), 213% (Figure 11) and 43% (Figure 12) of natural background radiant flux (integrated star light and/or nightglow, plus the effect of MS). The percentages are calculated with respect to the maximum radiant flux measured by the instrument around azimuth a " 260 ˝. Interestingly, for these three colours we now reproduce also the minima and maxima of the observed DoLP along one almucantar, without adding any further complexity to the model.

If slight discrepancies remain for the shape of DoLPpaq between the model and measurements for the blue line, this characteristic is very convincingly reproduced by the model in the green and orange lines.

It is striking that we manage to replicate these observations with our model in three different wavelengths simultaneously, as they are affected differently by light pollution, the natural background from the sky, not speaking of MS. Indeed, there exist some orange emissions in the star light and in the mesospheric nightglow [START_REF] Bellisario | O2 and OH night Airglow emissions derived from GOMOS-EnvisatInstrument[END_REF](Bellisario et al., , 2020)). The green and the blue lines exist both in the stellar light and the natural nightglow originating in thermospheric / ionospheric emissions. From Broadfoot & Kendall (1968), the orange is the tiniest of the three. The blue in the urban light spectrum is much dimmer than the green or the orange. These considerations are compatible with our results.

Finally, the misfit in the DoLP observed for shorter wavelengths (turquoise and purple) might require adjusted aerosol profiles (see e.g. [START_REF] Bergstrom | Estimates of the spectral aerosol single scattering albedo and aerosol radiative effects during safari 2000[END_REF], Figure 3). For example a different size distribution which would influence differently the shorter and the longer wavelengths. This would demand a thorough study of aerosols contribution, which is out of the scope of this paper. It is also likely related to the stronger impact of MS towards short wavelengths (see Section 4.4).

Discussion

We have presented in this study a polarized radiative transfer model in the single scattering approximation, which can account for a spread light pollution from the ground (including the surrounding topography), natural skylights (from stars or the nightglow) and the effect of aerosols. To our knowledge, it is the first time that such a radiative transfer code, that combines Rayleigh and Lorenz-Mie scattering, is applied to answer the question of the nightlight and its polarization. We have confronted this model to measurements performed at mid-latitudes in the French Alps in five different wavelengths, either within or outside the lines of emissions for the natural nightglow. We obtain a convincing comparison between the model outputs and the observations of the relative radiant flux, degree and angle of the linear polarization in the three wavelengths with the highest signal-to-noise ratio (green, orange and blue lines).

The results presented in this paper show the feasibility of nocturnal polarization measurements and modelisation. Our model, although currently limited to a single scattering scheme, is doomed to evolve, and constitutes a proof of concept that paves the way to further investigations of the night-sky light polarization.

Our results show that several contributions must be considered in order to explain observations: extended sources on the ground (light pollution from the nearby cities) and sources in the sky (the background) via two possible mechanisms (starlight and nightglow), plus a correction for multiple scattering. Furthermore, the atmospheric model must include Rayleigh scattering and Lorenz-Mie scattering by aerosols. We discuss below several potential applications of this model.

where wide areas are covered with snow, the reflection of auroras on the ground might be of importance. These effects will be accounted for in an upcoming study.

In the above interpretation of mid-latitude observations, considering an unpolarized background source and the 1-low aerosol model seem compulsory to explain most of our measurements. It is difficult at this stage to reach a definitive conclusion on its origin, because we lack an absolute measure of the radiant flux. However, in the absence of any other source, we have shown that our code can also handle a polarized background source, and that its signature (if any) should be detectable despite the Rayleigh and Lorenz-Mie scattering. We We have shown that the relative contribution of background lights to scattering strongly varies with the wavelength: from `43% for the longest one (orange), to `114% in the green and `219% in the shortest one (blue). Measuring how much this means in absolute units would require the knowledge of the city light spectrum (that is not provided with the input maps) as well as the instrument transfer function (such a calibration would require a dedicated study, as it depends on the photo-sensors and on the opacity of the filters used in the instrument). This higher background contribution at shorter wavelengths can be a sign that the uniform, isotropic and unpolarized background partly corrects for the effect of multiple scattering. Indeed, its contribution is expected to be higher at short wavelengths. This result is also compatible with the expected relative contributions from the city lamps, the nightglow and the starlight. On the one hand, the city contribution is likely higher in the green and orange than in the blue (i.e. the same sky emission for all three lines should lead to a relatively higher contribution for the blue). On the other hand, if we were to detect a contribution from the nightglow, we would expect relatively higher radiant fluxes in the blue and green than in the orange (which is not emitted in the ionosphere), for a constant contribution of the cities at the observation point. Thus our results are fully compatible with the nightglow emissions,

showing that the method may be used for such studies. Moreover, adjusting the modeling to the observations requires a contribution from the integrated starlight (at least that for the orange). The combined contributions of the nightglow and the starlights can therefore be detected with this method even in the presence of light pollution, although we cannot distinguish yet between the two sources.

Another major result is the influence of the aerosols. The good fit between the model and the measurements at several wavelengths comes with a determination of the size distribution, profile distribution and refractive index of the aerosols. Any attempt to change drastically these parameters degrades the fit. The precise determination of the aerosol parameters is out of scope of this preliminary study. It would require for example using independent measurements (such as LIDAR) or dedicated inversion schemes in order to optimize the input model parameters.

Nevertheless, we showed that the measurement of the light polarization constitutes an original way for aerosol studies at night.

Of course our model suffers from limitations, the first of which being the omission of multiple scattering. However, we have shown that this can be partly accounted for via the unpolarized isotropic background for the blue and longer wavelengths, without ruling out our model. Then, pollution maps are not provided for all wavelengths (but these might be calibrated with in situ measurements). Also, although we found that an isotropic sky was enough to fit our observations in the green and orange lines, a more complex model of starlight and nightglow that takes into account variations in time and space could further improve the model.

-25-manuscript submitted to Journal of Geophysical Research Nevertheless, our results suggest that our single-scattering polarized radiative transfer model together with our polarimeter may be used during night-time to detect light pollution even in remote areas, and to diagnose the aerosol content with for instance potential applications concerning the quality of air. The experimental technique is passive, contrary to lidars, low power consumption (less than 5 W), fully transportable in a suitcase, and can use any kind of visible source (artificial or natural such as the moon, the skylight, the airglow). These polarimeters may be deployed over a large area, thus allowing in the future a global coverage.

Coefficient 0.2 ă λ pµmq ă 0.5 0.5 ă λ pµmq [START_REF] Bucholtz | Rayleigh-scattering calculations for the terrestrial atmosphere[END_REF] centered around R (see Fig. 2), is a truncated cone defined by ε the half aperture angle of the virtual instrument (similar to that of the real polarimeter), d AR and l the half-height of the cone.

A
It is given by

Γ " π 3 tan 2 pεq l `3d 2 AR ´l2 {4 ˘. (A8)
Finally, the radiant flux in nW measured by the virtual instrument in A is:

F A " I R Ω A expp´τ AR q (A9) " I R Σ d 2 AR expp´τ AR q (A10)
Where Ω A is the solid angle (in sr) of the detector as seen from the scattering volume in R, Σ

is the surface area of the detector (in m 2 ), d AR is the distance from the scattering point to the detector (in km) and τ AR is the effective optical depth along the path from R to A. Developing I R into the initial parameters gives us:

F A " I R Σ d 2

AR

expp´τ AR q (A11)

" The total radiant flux measured by our virtual instrument is then the sum of both contributions as:

F A " L E 4π A E d 2 ER Σ d 2

AR

Γ expp´τ ER ´τAR q rβΦ ray pθq `ωC ext Φ aer pθqs (A17)

Figure 1 :

 1 Figure 1: Mid-Latitude night time atmospheric profiles from the 2001 MIPAS model (see text for details). From left to right: temperature [K], pressure [Pa], Ozone number density [cm ´3].

Figure 2 :

 2 Figure 2: Geometry of the problem for a light source in E of area A E , a scattering volume (Γ) centered at R and the instrument in A of area Σ and half aperture angle . The instrument pointing direction is defined by its azimuth a and elevation e. θ is the scattering angle.

manuscript submitted to Journal of Geophysical Research 3 . 2 One-dimensional radiative transfer model 3 . 2 . 1

 32321 Computations for a single point source over a single scattering volume

  figuration, the AoLP is AoLP aer determined by D aer , possibly 90 ˝off with respect to AoLP ray . In all cases, the DoLP for the sum of the Rayleigh and aerosol scattered radiant fluxes is DoLPpθq " |F ray A DoLP ray pθq `Faer A DoLP aer pθq| F A . (9)

  1-low, 2-high, 3-mid. The impact of the aerosol profiles is drastic, but not straightforward to interpret. It is illustrated in Figure4, for a point source located 5 km South from the observation point. With the 1-low model (lowest aerosol contribution), the radiant flux is increased by 600% (with respect to the case with no aerosols) when pointing above the source, and by 40% in the opposite direction. Using the 2-high profile (highest aerosols contribution), the radiant flux is increased by 150% when pointing above the source, and decreased by 40% in the opposite direction. With the 2-mid model, the radiant flux above the source increases by 120%, against 10% in the opposite direction. The ratio of the maximum to minimum radiant flux along the almucantar is in all cases amplified by the presence of aerosols.

Figure 3 :

 3 Figure 3: Almucantar for e " 45 ˝, for a point source on ground, on a flat Earth surface, South of the instrument at varying distances d AE (with no aerosols). Upper panel: measured radiant flux [nW] for a source of radiance 100 nW{m 2 {sr. Middle panel : DoLP [%]. Lower panel: AoLP [ ˝]. From blue to red: d AE " 1, 5, 10, 50, 100 km. In the upper panel, the highest value (in blue) corresponds to the closest source position and the lowest value (in red) to the most distant point source. Similarly, in the middle panel, the most distant point source has two very marked maxima (in red) while the closest is the flattest (in blue). The dashed line in the middle panel corresponds to a theoretical case for a source at infinite distance.

Figure 4 :

 4 Figure 4: Almucantar for e " 45 ˝, for different aerosol profiles. A point source is located 5 km South of the instrument. The continuous line is the case without aerosols (similar to Figure3). The dashed lines correspond to three different aerosol models listed in Table1: 1-low (black), 2-high (blue) and 3-mid (red).

Figure 5

 5 Figure5shows the model outputs over an almucantar (with no aerosols), for a nightglow either unpolarized, or polarized with a DoLP " p source " 1% (arbitrary choice) along several directions: East-West (EW), North-South (NS) or along the magnetic field (B). For this last case, the instrument is positioned at latitude 45.2123 ˝N and longitude 5.9369 ˝E and we use the CHAOS-6 internal magnetic field model[START_REF] Finlay | Recent geomagnetic secular variation from Swarm and ground observatories as estimated in the CHAOS-6 geomagnetic field model[END_REF]) evaluated in 2019. The atmospheric profile described in Section 3.1 is used. As the background skylight is isotropic, the measured radiant flux does not change over an almucantar.

Figure 5 :

 5 Figure 5: Polarization parameters for an almucantar at e " 45 ˝in the presence of an isotropic skylight of radiance taken at 100 nW/m 2 /sr at 557.7 nm. No aerosols are present in the model. The direct radiant flux is taken into account, with different polarization parameters. No polarization of the skylights (black line), with a 1% polarization along the East-West (red) and North-South (green) directions, and along the magnetic field direction (blue).

Figure 6 :

 6 Figure 6: Maximum DoLP found over the sky as a function of solar elevation angle at three wavelengths: 630 nm (red), 530 nm (green) and 450 nm (blue). Observations from PS11 are shown in plain lines (resp. dashed lines) for a low (resp. high) aerosols content. Predictions from POMEROL are shown in dotted lines with circled (resp. crosses) for low (resp. high) aerosols content.

Figure 7 :

 7 Figure 7: Comparison of POMEROL and MODTRAN-P maximum DoLP with the PS11 observations at maximum solar elevation, at three wavelengths: 630 nm (red), 530 nm (green) and 450 nm (blue). The group of points with DoLP ă 40% correspond to a high aerosol content, while symbols in the upper right quadran correspond to a low aerosol content. Crosses represent the result of POMEROL. Single (resp. multiple) scattering predictions obtained with MODTRAN-P by PS11 are indicated with squares (resp. circles). The spread within the several low aerosol contents investigated by the authors is indicated by errorbars.

  ; Bellisario et al., 2014, 2020). Finally, the turquoise line can be produced by Herzberg and Chamberlain O 2 lines and/or ground light pollution. The former indeed induce spectral lines of weak intensity (Broadfoot & Kendall, 1968; Leinert et al., 1998), while the latter, which we account for, might be related to the presence of Hg in city lights. The two contribution are likely weak, and the main Hg line falls outside the filter bandwidth.

Figure 8 :

 8 Figure 8: Input emission map and contour of the elevation map for altitudes 500 m (white), 1000 m (yellow) and 2000 m (orange). The instrument is at the center of the map (red cross). Grenoble is the bright emission west to the instrument. The map covers 100 km around the instrument.

Figure 9 :

 9 Figure 9: Comparison between model outputs and observational data. Measured data are shown in green, the width of the lines indicating the errorbars (for their computation, we refer to Bosse et al., 2020). Errors on the azimuth due to the pointing direction are of the order of a few degrees (not shown). The other symbols indicate the model output for a ground light pollution map (see text for details). Several aerosol profiles are considered. (see Table1): without aerosols nor background lights (black ‹); without aerosols and with background lights adjusted to fit the radiant flux variations (red `); without background lights and with aerosols model 1-low (blue ˆ), 2-high (magenta ) and 3-mid (purple ). Abscissa are given as the pointing direction under the bottom panel (azimuth where 0 is North, 90 is East) and the corresponding time in UTC format above the top panel.

  ) and are due to the scattering in the first few kilometers of atmosphere. Here it is smoothed out because of the multiple light sources and the occultation from the mountains. This feature is not present in the data and may hint that the model overestimates the scattering below 10 km. The modeled AoLP is very similar to the measurement, showing the same rotation pattern at almost the same angles, except around 90 ˝in azimuth (i.e. pointing away from Grenoble). This difference can be explained by the poor quality of the data in this direction where the radiant flux is minimal and the DoLP close to zero, inducing large uncertainties on the AoLP measurements (see Bosse et al., 2020).

Figure 10 :

 10 Figure 10: Almucantar in the green line. Measured data are shown in green, the width of the lines indicating the errorbars. Errors on the azimuth due to the pointing direction are of the order of a few degrees (not shown). Model predictions (black stars) are obtained with the 1-low aerosol profile plus an unpolarized, isotropic background skylight adjusted so as to fit the radiant flux (see text).

Figure 11 :

 11 Figure 11: Same as Figure 10 for the blue line.

Figure 12 :

 12 Figure 12: Same as Figure 10 for the orange line.

Figure 13 :

 13 Figure 13: Same as Figure 10 for the purple line.

Figure 14 :

 14 Figure 14: Same as Figure 10 for the turquoise line.

  have shown that a horizontally polarized uniform sky source produces AoLP changes along an almucantar that resemble to the ones produced by scattering within the atmosphere. It might thus be insightful to consider the possibility of the polarization by horizontal ionospheric electrical currents. These can be partially reconstructed from ground-based magnetometers (e.g.Amm & Viljanen, 1999;[START_REF] Pulkkinen | Ionospheric equivalent current distributions determined with the method of spherical elementary current systems[END_REF]. This possibility is motivated by the observation of coherent fluctuations in time series for both the polarization parameters of auroral lights, and the electron density content in the ionosphere (Bosse et al., 2020). This may have further consequences for the study of the equatorial airglow, below the electrojet.

Table 1 :

 1 tends to infinity (in cm ´3) and H the scale height (in m). The values chosen for the above Parameters used to define the aerosol model (see text for details).

		Complex refractive index Size distribution	Vertical number density profile
	Aerosol
	profile	Real	Imaginary	R (µm) lnpσq H (m) n 0 (cm ´3) n B (cm ´3)
	name
	1-low	1.45	0.0035	0.15	0.29	440	4000	10
	2-high 3-mid	1.61	0.03	0.557	0.266	500	1000 500	1

  However, the configuration of our experiment significantly differs with most of the already documented work. PS11 for instance evaluate the polarization of the solar light scattered in the atmosphere. They do not study the AoLP, and only provide measurements of the maximum DoLP, i.e. in a direction at approximately 90 ˝from the Sun direction. Many of such studies also consider the rising Sun only. In our case, we consider a mixture of sources, either extended or localised, during the night-time, and observed over a multiplicity of angles.

	Then, because our study constitutes a first approach to a fully new principle of night-
	time atmospheric observations, we wish to keep our set-up as simple as possible, although
	realistic, and we do not account for MS in this prototype version of the POMEROL model. Its
	effect is out of the scope of the present work, and should constitute a dedicated future study.
	As we will see later, already a large part of our night-time measurements can be interpreted in
	the single scattering approximation, which means that the SSA an already catch a great part of
	the phenomenon we are dealing with. This is particularly true not only for the DoLP profiles
	along one almucantar, but also for the AoLP data, a quantity that is almost never considered
	in the literature. Meanwhile, we provide below an estimate of the uncertainty level associated
	with this limitation of our model, by presenting comparisons of POMEROL with previously
	published observations and MS model predictions from PS11.
	Several studies addressed the comparison between single and multiple scattering of po-
	larized light in the atmosphere. One can cite Evans & Stephens (1991), but also a series of
	balloon experiments (Herman et al., 1986; Santer et al., 1988), which interpretation was first
	based on a single scattering approach, and later on refined with MS. Ougolnikov & Maslov
	(2002), followed by Ugolnikov et al. (2004) proposed a thorough study on the impact of MS as
	a function of the wavelength in the visible range. They conclude to a radiant flux contribution
	from single scattering ranging from about 55% in the blue to over 80% in the red (slightly

[START_REF] Hovenier | Multiple Scattering of Polarized Light in Planetary Atmospheres[END_REF]

, while it can give decent predictions for red and longer wavelengths. There are several approaches to polarized MS, such as adding-doubling methods with several levels of refinement, Monte-Carlo methods, etc. (de Haan et al., 1987;

[START_REF] King | Comparative accuracy of selected multiple scattering approximations[END_REF][START_REF] Ramella-Roman | Three monte carlo programs of polarized light transport into scattering media: part i[END_REF]

. For a review of multiple scattering of waves from the theoretical perspective, see for example

[START_REF] Lax | Multiple Scattering of Waves[END_REF]

. less than theoretical estimates). A comprehensive study over a wide range of wavelengths and zenith angles can also be found in Pust & Shaw (2011) (mentioned as PS11 hereafter), in various aerosols conditions. These authors studied the output of the MODTRAN-P model (e.g.

[START_REF] Berk | Modtran® 6: A major upgrade of the modtran® radiative transfer code[END_REF]

, accounting or not for the effect of MS. Their model was tested against measurements corresponding to various aerosol contents: a clear day with low aerosol content, another day

Table 2 :

 2 List of emission lines observed during the January, 2021 campaign (see text). All emissions are present in the light pollution spectrum.

	Name	Wavelength (nm)	Instrumental width (nm)	Atmospheric source	Layer
	Orange	620	2	OH	High Mesosphere (ă 80 km)
	Green	557.7	10	O	Thermosphere (110 km)
	Blue	427.8	10	N 2	Ionosphere (90 km)
	Turquoise	413	10	O 2	Mesosphere (ă 80 km)
	Purple	391.4	10	N 2	Ionosphere (85 km)

Table A1 :

 A1 Parameters used in equation (A20).

	ray	7.68246 ˆ10 ´4	10.21675 ˆ10	´4
	A op	6.49997 ˆ10 ´3	8.64145 ˆ10 ´3
	B	3.55212	3.99668
	C	1.35579	0.00110298
	D	0.11563	0.0271393

  At this point, we replace the scattering coefficient by the corresponding expression, whether we consider the Rayleigh scattering on air molecules, or the Mie scattering on aerosols. This gives: AR expp´τ ER ´τAR qΓ ωC ext Φ aer pθq (A16)

					E E σ	Φpθq 4π	Σ AR d 2	expp´τ AR q	(A12)
			" L E		A E d 2 ER	expp´τ ER ´τAR qσ	Φpθq 4π	Σ AR d 2	(A13)
			"	L E 4π	A E d 2 ER	Σ AR d 2	expp´τ ER ´τAR qσΦpθq	(A14)
		F ray A "	L E 4π	A E d 2 ER	Σ AR d 2	expp´τ ER ´τAR qΓβΦ ray pθq	(A15)
	and:	F aer A "	L E 4π	A E d 2 ER	Σ d 2

First, it will be used to interpret properties of the polarized lights measured at high latitudes under auroral conditions (Bosse et al., 2020), to validate or not a polarization when these are emitted in the ionosphere (Lilensten et al., 2016). The present study shows that the use of wavelengths outside the ionospheric emissions lines is crucial to determine a correct profile of aerosols, as this latter is of prime importance to understand the measured radiant flux and DoLP. Furthermore, in auroral conditions, ionospheric lights will strongly dominate over the integrated starlight. As a consequence, the complex picture of these emissions (2D maps) will be needed as an entry of the model. We can also imagine that in these countries -24-manuscript submitted to Journal of Geophysical Research Acknowledgments This work was partly funded by the French polar institute IPEV under program POLARLIS 3, number 1026. It was also funded by the Prematuration CNRS program and by the Maturation program CM180023, project PTCU Number 180018M from the SATT company Linksium.

The polarimeter is patented by the CNRS and valorised by Linksium (international patent number 1873378). We thank the Programme National Soleil-Terre for its financial support.

NG was partially supported by the French Centre National d'Etudes Spatiales (CNES) for the study of Earth's core dynamics in the context of the Swarm mission of ESA.

We thank Pierre Simoneau (ONERA) for usefull discussions on the mesopheric nightglow.

The source code will be freely available online as opensource as soon as it is documented and portable. Ground emission maps are available here: https://eogdata.mines.edu/ products/vnl/ Atmospheric profiles are available here: http://eodg.atm.ox.ac.uk/RFM/atm/ Elevation data are available here: https://www.eorc.jaxa.jp/ALOS/en/aw3d30/ index.htm

Appendix A Computation of the radiant flux measured by our virtual instrument

In this appendix, we compute the radiant flux measured by a virtual instrument in A when the light is emitted by a single point source E, and scattered at a single point R along the line of sight (see Figure 2). We describe each intermediate step in order for the reader to understand and reproduce our results.

We use here the units used by POMEROL, with the scale factor implemented in the code. The scale factor has no effect on the dimensional analysis, but are kept here for more transparency on the code.

First, the source is fully described by its radiance L E in nW/m 2 /sr as given by the input map (see section 3.1), and its surface area A E in km 2 . It is considered small enough to be approximated as a point source of radiant intensity I E " L E A E in nW/sr. Its emission is considered isotropic. From this, the irradiance reaching the scattering volume in R is (in nW/m 2 ):

where Ω E " A E {d 2 ER is the solid angle of the emission surface as seen from the scattering point at a distance d ER (in km), and τ ER is the effective optical depth of the atmosphere between the emission in E and the scattering in R (see equations A21 to A24).

The radiant intensity I R (in nW/sr) scattered by the volume Γ of air or aerosol around R is:

where σ (in km 2 ) is the total scattering coefficient, θ is the scattering angle and Φpθq is the scattering phase function of the molecules or particles (see equation A5 and following text).

By definition the phase function Φ is dimensionless and normalized across a sphere (e.g. 1. For the air molecules, the total scattering coefficient is

and for the aerosols, it is:

with β the Rayleigh scattering volume coefficient in km ´1 as described in [START_REF] Bucholtz | Rayleigh-scattering calculations for the terrestrial atmosphere[END_REF] (see equations A18 to A20), ω the single scattering albedo of the aerosol, C ext the extinction coefficient of the aerosol in km ´1 and Γ the scattering volume in km 3 . Γ, the scattering volume We describe below in more detail the computations of all the physical quantities used to compute F A . Classically, the Rayleigh scattering cross section per molecule is [START_REF] Rayleigh | On the scattering of light by small particles[END_REF][START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF][START_REF] Bucholtz | Rayleigh-scattering calculations for the terrestrial atmosphere[END_REF]:

with λ the wavelength, N s the molecular number density for standard air, m s the refractive index of standard air at λ and ρ n the depolarization factor. In our model, this is accounted for via the Rayleigh volume scattering coefficient βpλ, zq (in km ´1). Following [START_REF] Bucholtz | Rayleigh-scattering calculations for the terrestrial atmosphere[END_REF], it is approximated as βpλ, zq " β 0 pλq Ppzq P 0

with Ppzq and T pzq the atmospheric pressure and temperature profiles. P 0 " 101 325 Pa and T 0 " 288.15 K are the pressure and temperature at sea level, and

with λ the wavelength (in µm). This approximation takes into account the depolarization (or King) factor. Values for coefficients A ray , B, C and D are given in Table A1, and βpλ, zq |z"0 is shown in Figure A1.

The effective optical depth τ of the atmosphere between E and R (and similarly between R and A), is the sum of three contributions from Rayleigh scattering, ozone absorption and aerosol extinction:

τ " τ ray `τO 3 `τaer . (A21)

The optical depth from Rayleigh scattering between E and R is

where τ 0 pλq follows equation (A20) with A ray replaced by A op (value given in Table A1). The optical depth for ozone absorption is

where N O 3 pzq is the ozone number density at altitude z (provided with the atmospheric profile, see Figure 1) and σ O 3 pλq (shown in Figure A1) is the ozone absorption cross section at wavelength λ and a fixed temperature of 273 K (Burrows et al., 1999). This is an approximation which could be refined in future versions with a dependence on temperature. The optical depth for aerosols is τ aer pz, λq "

ż R E C aer ext pz 1 , λqdz 1 . (A24)