F André
email: frederic.andre@insa-lyon.fr

Radiative transfer in the O 2 A-band -a fast and accurate forward model based on the -distribution approach

Keywords: remote sensing, O 2 A-band, scattering atmosphere, -distribution, k-distribution, LBL

An efficient method based on accurate transmissivity calculations to account for gaseous absorption in Monte Carlo / quasi-Monte Carlo radiative transfer codes is presented. The modeling approach is based on an improved nonuniform transmission formulation, called the -distribution method. The technique is founded on two components: 1/ an original uniform method to tabulate band averaged transmissivities of gaseous paths and 2/ its extension to non-uniform paths, based on the Godson-Weinreb-Neuendorffer's (GWN) approximation that consists of the definition of effective scaling factors to relate gaseous spectra in distinct states. As the GWN method is known to provide results dependent on the ordering of the gas layers, two path reordering strategies are introduced and compared. One of them is founded on results from statistical theory and requires the introduction of a new coefficient (Kendall's Ke). This coefficient together with its role on the path reordering strategy is introduced, detailed and analyzed. The two schemes are then assessed against Line-By-Line calculations in line-of-sight geometries. Two configurations representative of radiative transfer in the O 2 A-band for scattering atmospheres are then studied: the first one only considers molecular (Rayleigh) scattering and the second one involves various cloud configurations (single or bi-layer clouds located at

Introduction

Satellite and ground based imaging devices are widely used in remote sensing applications for the inference of atmospheric parameters from radiance (reflected and/or emitted) measurements. The reliability of the inferred atmospheric profiles (temperature, species concentrations, etc) strongly depends on the accuracy of the forward radiation model used in combination with the inversion method.

The most accurate technique to account for gaseous absorption in the simulation of radiative transfer processes in the atmosphere is the Line-By-Line (LBL) method [START_REF] Clough | Line-by-line calculation of atmospheric fluxes and cooling rates: Application to water vapor[END_REF][START_REF] Clough | Atmospheric radiative transfer modeling: a summary of the AER codes[END_REF]. But this technique, even if some optimized versions such as 4A/OP (Automatized Atmospheric Absorption Atlas) exist [START_REF] Scott | A fast line-by-line method for atmospheric absorption computations: The automatized atmospheric absorption atlas[END_REF], is widely recognized as time-consuming which prevents its use in many practical problems, especially for operational calculus. Efficient though accurate parameterizations for the treatment of gaseous absorption remain of primary importance in atmosphere sciences, both in an operational context at a time where the amount of data produced by remote sensing devices increases exponentially, or for other applications, such as weather forecasting or radiative budget studies.

The high computational cost of LBL calculations is caused by the double integration scheme required to evaluate band averaged transmissivities of nonuniform paths: 1/ a path integration is first needed to account for local contributions of the gaseous species encountered along the radiation path. Path integrals need to be evaluated for each wavenumber ν inside the spectral band of interest and, 2/ a spectral integration needs to be performed in order to average the contributions of all spectral values inside the considered band. This double integration scheme, and more critically the first one as it needs to be performed for each wavenumber, is the main source of computational cost in LBL calculations. It should be noticed that inversion of measured spectral data requires iterative schemes, associated with the minimization process of a cost function. This makes LBL calculation cost even more critical in this particular application context, since all path integrals need to be evaluated many times during the exploration of the parameter search space before to reach the optimal solution. Presently, k-distribution methods [START_REF] Lacis | A description of the correlated k distribution method for modeling nongray gaseous absorption, thermal emission, and multiple scattering in vertically inhomogeneous atmospheres[END_REF] are among the most widely used approximate techniques to treat gaseous absorption. These approaches allow reducing significantly the amount of path integrals required to estimate band averaged radiances from thousands in the LBL approach to a few dozens in k-distributions. k-distribution methods require assumptions to treat path non-uniformities. These assumptions are their main source of errors in radiative transfer applications, with orders of magnitude of errors usually reported on the literature of a few percents in non-uniform cases [START_REF] Drummond | A novel k-distribution parameter development system and its application to MAS/SUCCESS channels[END_REF].

In order to circumvent this problem of multiple integration, encountered in LBL but also to a less extent in k-distribution models, the -distribution method was proposed recently [START_REF] André | The -distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media[END_REF]. The approach, which was initially developed and applied in high temperature combustion applications, is the main focus of the present paper. This method is highly computationally efficient because it neither requires spectral or spatial integration schemes to estimate band averaged transmission functions of non-uniform gaseous paths. This results in a very fast direct model that can nevertheless achieve LBL accuracy both in uniform and non-uniform situations.

Up to now, most of the application cases reported in the literature for this method were dedicated to high temperature applications, including standard line-of-sight benchmarks and turbulence-radiation interactions calculations in [START_REF] André | The -distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media[END_REF], and to radiative heat transfer simulations in benchmark [START_REF] André | Accuracy of engineering methods for radiative transfer in CO2-H2O mixtures at high temperature[END_REF] and actual tridimensional industrial (glass furnace geometries) configurations [START_REF] Galtier | Assessment of engineering gas radiation methods in an industrial glass furnace configuration[END_REF]. The aim of the present work is to evaluate this technique for atmospheric calculations by simulating reflectance at the top of the atmosphere, for different incident and satellite viewing angles, and for typical atmospheres composed of gaseous molecules and / or clouds.

The paper is structured as follows. The second section is dedicated to a description of the treatment of uniform gaseous paths by the -distribution approach. The third section focuses on the extension of the uniform approach to non-uniform paths. This adaptation of the method is founded on an application of the Godson-Weinreb-Neuendorffer's (GWN) method which consists in the definition of effective scaling factors. This method is recasted within the frame of Archimedean copula's theory to analyse some of its propertie.s The link between these mathematical properties and their physical interpretation is provided. The use of recent results from copula theory to improve the GWN method is the main originality of the present work. Section 4 is dedicated to applications of the method in clear sky and scattering atmosphere configurations. Absorption by a single molecular specie (O 2 ) is considered all along the paper.

The -distribution method in uniform gaseous layers

In this section, we introduce the main equations required to apply thedistribution method for the calculation of transmissivities of uniform atmospheric paths averaged over the spectral response function of an optical filter φ (ν). Only the main results are given. Interested readers will find more details and explanations about the method in Refs. [START_REF] André | The -distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media[END_REF][START_REF] André | An analysis of the symmetry issue in the -distribution method of gas radiation in non-uniform gaseous media[END_REF].

We start with the gaseous transmittance of a single gas (as considered all along the present paper) in a uniform layer defined as:

τ ∆ν φ (L) = 1 Φ ∆ν φ (ν) exp (-κ ν L) dν (1)
where τ ∆ν φ (L) is the band averaged gaseous transmittance weighted by the filter reponse function φ (ν), L (in cm) the total length of the gas path inside the uniform layer, κ ν (in cm -1 ) the spectral absorption coefficient of the gas inside the uniform layer at wavenumber ν (in cm -1 ) and Φ the integral of the filter response function over the band ∆ν, viz. Φ = ∆ν φ (ν) dν.

The -distribution approximation of the transmittance given by Eq. ( 1) can be written as:

τ ∆ν φ (L) = Gr τ ∆ν 0,φ (L) (2)
where τ ∆ν 0,φ (L) is called the "germ" model and Gr is a mapping function that associates values of the germ model for a given L with the true LBL transmittance at the same length. The germ model needs to be chosen in such a way that analytical mathematical expressions exist for both the transmissivity τ ∆ν 0,φ (L) and its inverse Λ ∆ν 0,φ (X), X ∈ [0, 1] defined as Λ ∆ν 0,φ τ ∆ν 0,φ (L) = L. A simple choice for the germ is thus the Statistical Narrow Band model for Lorentz lines with the Malkmus' distribution of linestrengths (written from now on SNB-LM), viz.:

τ ∆ν 0,φ (L) = exp - β π 1 + 2πk P L β -1 (3)
where k P is the mean absorption coefficient of the gas weighted by the filter response function:

k P = 1 Φ ∆ν φ (ν) κ ν dν (4)
Parameter β, that represents an overlapping parameter in the SNB-LM [START_REF] Young | Band model theory of radiation transport[END_REF],

only plays here the role of a coefficient that adjusts the curve-of-growth of the germ model at the optically thick limit. This parameter β is defined as:

β = 1 k P /k R -1 (5)
where we have introduced a band averaged Rosseland mean absorption coefficient k R as:

1 k R = 1 Φ ∆ν φ (ν) κ ν dν (6)
The inverse Λ ∆ν 0,φ (X) of the transmission function ( 3) is given analytically as:

Λ ∆ν 0,φ (X) = β 2πk P 1 - π β ln X 2 -1 (7)
From these definitions, one can immediatly notice that the mapping function,

Gr, is simply given as:

Gr(X) = 1 Φ ∆ν φ (ν) exp -κ ν Λ ∆ν 0,φ (X) dν (8)
The role of function Gr is to associate values of the germ model at a given value of the gas path length L, which provides some rough estimate of the true LBL transmissivity, to its actual LBL value. Gr maps the unit interval [0, 1] into itself and can be readily tabulated at any desired accuracy by: 1/ discretizing the interval [0, 1], 2/ estimating the function at the corresponding X ∈ [0, 1] by a direct application of Eq. ( 8) to the LBL dataset. Then, as soon as the One interesting feature of this formulation in terms of mapping function Gr is that as Gr is strictly increasing, its inverse Gr -1 exists. This function can be obtained easilly by reverting the look-up table of the Gr function. This allows obtaining the inverse of function τ ∆ν φ (L), written and defined through the relationship • τ ∆ν φ (L) = L as:

look-up table {X i , Gr(X i )}, i = 1, ..,
(X) = Λ ∆ν 0,φ Gr -1 (X) (9)
The -distribution method is probably the first approach to allow simultane-ously estimating the direct model (band averaged transmissivity) and its inverse.

This inverse can be used to construct a large range of possible models for the extension of the -distribution technique from uniform to non-uniform situations.

Indeed, in the case of n distinct thermophysical states represented by a sequence κ 1 ν , .., κ n ν of spectral absorption coefficients and for a set L 1 , .., L n of pathlengths, the band averaged transmissivity of the non-uniform path τ ∆ν 1..n (L 1 , .., L n ) (from now on, we will simplify the general notation, which considers a filter response function, for legibility but application cases described later in this paper will account for its effect) defined as:

τ ∆ν 1..n (L 1 , .., L n ) = 1 ∆ν ∆ν exp(-κ 1 ν L 1 + .. + κ n ν L n ) dν (10)
can be written under the following mathematical form:

τ ∆ν 1..n (L 1 , .., L n ) = C 1..n τ ∆ν 1 (L 1 ), .., τ ∆ν n (L n ) (11)
where C 1..n is a n-dimensional copula given in terms of functions i , i = 1, .., n as:

C 1..n (X 1 , .., X n ) = 1 ∆ν ∆ν exp -κ 1 ν 1 (X 1 ) + .. + κ n ν n (X n ) dν (12)
This formulation was successfully applied in high temperature configurations in Ref. [START_REF] André | A polynomial chaos approach to narrow band modeling of radiative heat transfer in non-uniform gaseous media[END_REF] for dimension up to n = 4 based on multivariate polynomial representations of function C. For higher dimensions, polynomial expressions can quickly require a large amount of coefficients and become computationally inefficient both in terms of calculation cost (which increases with the number and orders of the polynomials) and memory space (to store the polynomial coefficients): the use of neural networks is then recommended. Such developments are scheduled as future work but are not considered further in the following.

Another possible use of function is its application to the GWN (Godson-Weinreb-Neuendorffer, see [START_REF] Young | Band model theory of radiation transport[END_REF]) method which is detailed in the next section.

Adaptation of -distributions to non-uniform gaseous paths through the GWN method

The GWN method is among the oldest non-uniform techniques. It was first proposed by Godson in 1953 [START_REF] Godson | The evaluation of infra-red radiative fluxes due to atmospheric water vapour[END_REF], two years before the widely used Curtis-

Godson approximation [START_REF] Young | Nonisothermal band model theory[END_REF][START_REF] Young | Band model theory of radiation transport[END_REF]. The same method was "rediscovered" independently in 1973 by Weinreb and Neuendorffer [START_REF] Weinreb | Method to apply homogeneous-path transmittance models to inhomogeneous atmospheres[END_REF] and, in 1981, by Gordley and

Russel [START_REF] Gordley | Rapid inversion of limb radiance data using an emissivity growth approximation[END_REF] who called it the Emissivity Growth Approximation (EGA). More recently, Modest used the same technique within the frame of the k-distribution method and called it the scaled-k approximation [START_REF] Modest | Narrow-band and full-spectrum k-distributions for radiative heat transfer-correlated-k vs. scaling approximation[END_REF]. Recent applications of the GWN method (called in this reference the EGA technique) in atmospheric configurations can be found in [START_REF] Marshall | BANDPAK: Algorithms for modeling broadband transmission and radiance[END_REF].

The GWN method can be rather naturally introduced as a generalization of scaled models. This is the aim of the next section. Connections with the theory of Archimedean copulas are also emphasized.

Scaled models and Archimedean copulas

We consider two spectra κ 1 ν and κ 2 ν associated with two gas layers in distinct thermophysical states. Layers are represented by exponents 1 and 2 respectively.

The length of the gas path in state 1, at temperature T 1 , is L 1 and that in state 2, at temperature T 2 , is L 2 . The transmissivity of the non-uniform path L 1 + L 2 is defined as:

τ ∆ν 12 (L 1 , L 2 ) = 1 ∆ν ∆ν exp(-κ 1 ν L 1 -κ 2 ν L 2 ) dν (13)
As the transmissivities of the uniform layers τ ∆ν 1 (L 1 ) and τ ∆ν 2 (L 2 ) are strictly decreasing with respect to L 1 and L 2 , their inverses 1 and 2 defined as solu- 13) can be rewritten in terms of these inverses to provide:

tions of i • τ ∆ν i (L i ) = L i , i = 1, 2 , exist. Equation Eq. (
τ ∆ν 12 (L 1 , L 2 ) = C 12 τ ∆ν 1 (L 1 ), τ ∆ν 2 (L 2 ) ( 14
)
where function C 12 is defined for X, Y ∈ [0, 1] as:

C 12 (X, Y ) = 1 ∆ν ∆ν exp -κ 1 ν 1 (X) -κ 2 ν 2 (Y ) dν (15)
It can be readily checked from its definition that function C 12 has the following properties:

C 12 (X = 1, Y ) = Y, C 12 (X, Y = 1) = X (16)
C 12 (X = 0, Y ) = C 12 (X, Y = 0) = 0 (17) ∂ 2 C 12 (X, Y ) ∂X∂Y ≥ 0 (18)
Accordingly, function C 12 is mathematically called a copula [START_REF] Nelsen | An introduction to copulas[END_REF]. It is important to notice here that the introduction of the concept of copula is only useful to clarify the type of function that relates the transmissivity of uniform paths to those of non-uniform layers. A detailed knowledge of copula's theory is not required to follow the next developments. Some useful definitions and theorems related to copula theory are provided in Reference [START_REF] Nelsen | An introduction to copulas[END_REF].

The set of equations Eqs. ( 16) to ( 18) can be interpreted physically. Indeed, a value of X or Y equal to 1 indicates that the length of one of the two layers is zero ( i (1) = 0). In this case, the property given by Eq. ( 16) simply means that the application of function C 12 to the couple of transmissivities provides the transmissivity of the layer with non-nul length. In a similar way, if one of the values of X or Y is equal to zero, this means that one of the gas paths has an infinite length ( i (0) = +∞) in which case the transmissivity of the non-uniform paths needs to be zero: this is actually what is meant by Eq. ( 17). In order to interpret physically the inequality [START_REF] Marshall | BANDPAK: Algorithms for modeling broadband transmission and radiance[END_REF], we need to introduce the next exchange ϕ δL1↔δL2 between two small path elements δL 1 and δL 2 defined as:

ϕ δL1↔δL2 = ∂ 2 τ ∆ν 12 (L 1 , L 2 ) ∂L 1 ∂L 2 δL 1 δL 2 [I b,ν (T 2 ) -I b,ν (T 1 )] (19)
where:

∂ 2 τ ∆ν 12 (L 1 , L 2 ) ∂L 1 ∂L 2 = 1 ∆ν ∆ν κ 1 ν κ 2 ν exp -κ 1 ν L 1 -κ 2 ν L 2 dν ≥ 0 ( 20
)
The sign of ϕ δL1↔δL2 thus only depends on the difference

[I b,ν (T 2 ) -I b,ν (T 1 )].
Consequently, it is positive if T 2 ≥ T 1 indicating that the transfer of radiative energy arises in the direction of decreasing temperatures. Eq. ( 20) can be rewritten in terms of function C 12 :

∂ 2 τ ∆ν 12 (L 1 , L 2 ) ∂L 1 ∂L 2 = ∂τ ∆ν 1 (L 1 ) ∂L 1 ∂τ ∆ν 2 (L 2 ) ∂L 2 ∂ 2 C 12 X = τ ∆ν 1 (L 1 ), Y = τ ∆ν 2 (L 2 ) ∂X∂Y (21)
which is positive if and only if the joint derivative of C 12 with respect to X and Y is strictly positive i.e. if Eq. ( 18) is true. Consequently, Eq. ( 18) can be interpreted as a condition that ensures the sign of net exchanges between gaseous cells to be properly evaluated when formulated in terms of function

C 12 .
It should be noticed at this level that if the "copula" property of function C 12 can be easilly justified physically in the bi-variate case, the more general situation of n layers cannot be, unfortunately, interpreted in the same simple way.

If gas spectra are scaled, then the ratio u =

κ 2 ν κ 1 ν
does not depend on wavenumber, ν. Then, the previous equation simplifies into:

τ ∆ν 12 (L 1 , L 2 ) = 1 ∆ν ∆ν exp -κ 1 ν (L 1 + uL 2 ) dν = τ ∆ν 1 (L 1 + uL 2 ) ( 22
)
and the transmissivity of the non-uniform path can thus be treated exactly in the same way as the transmissivity of the gas in state 1 with a total length L 1 + uL 2 (a symetrical relationship can also be written in terms of state 2). In this situation, it thus suffices to know how to calculate average gas radiative transmissivities of uniform layers to treat any non-uniform situation.

Equation ( 22) can be rewritten in terms of the inverse 1 of the transmission function of the gas in state 1 defined by the relationship 1 • τ ∆ν 1 (L 1 ) = L 1 providing:

τ ∆ν 12 (L 1 , L 2 ) = C 11 τ ∆ν 1 (L 1 ), τ ∆ν 2 (L 2 ) (23)
If gas spectra are scaled, the copula C 12 = C 11 thus only depends on a single thermophysical state (state 1 in the considered case) and the copula is called

Archimedean. Archimedean copulas are among the most widely used copula models because of their relatively simple extension to high dimensions.

In practice there exist two main physical reasons that explain the departure of gas spectra from scaling:

gas spectra are made of many thin spectral lines whose profiles depend non-linearly on the local gas properties (temperature, pressure, species concentrations).

the content of gas spectra (amount of single lines that participate to the value of the absorption coefficient at some given wavenumber ν) varies significantly with temperature. This is because linestrengths (surface of a single line calculated over its spectral profile) vary exponentially with respect to the temperature of the gas through a Boltzmann law.

Accordingly, the definition of a constant scaling coefficient is not physically realistic in the case of real gaseous media. Nevertheless, the use of Equation Eq.

(23) remains possible, as a mean to relate the transmissivity of the non-uniform path to those of the uniform sub-layers. This is the principle of the GWN method. Indeed, this technique mostly consists of the definition of effective scaling coefficients that depend on the length of the gas paths through the [START_REF] André | The -distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media[END_REF][START_REF] André | Effective scaling factors in non-uniform gas radiation modeling[END_REF]. A detailed physical analysis of the concept of effective scalling factor is given in Ref. [START_REF] André | Effective scaling factors in non-uniform gas radiation modeling[END_REF]. With this notation, it is in fact easy to check that:

relationship u(L 2 ) = 1•τ ∆ν 2 (L2) L2
τ ∆ν 12 (L 1 , L 2 ) = C 11 τ ∆ν 1 (L 1 ), τ ∆ν 2 (L 2 ) = τ ∆ν 1 [L 1 + u(L 2 ) L 2 ] ( 24
)
which is similar to Eq. [START_REF] Górecki | An approach to structure determination and estimation of hierarchical Archimedean copulas and its application to bayesian classification[END_REF]. But, instead of a constant scaling coefficient, an effective scaling factor that can be defined as solution of

τ ∆ν 2 (L 2 ) = τ ∆ν 1 [u(L 2 ) L 2 ]
appears.

In general, the two copulas C 11 and C 22 are different and their application to the same couple of transmissivities τ ∆ν 1 (L 1 ) and τ ∆ν 2 (L 2 ) has no reason to provide the same result. Criteria thus needs to be defined in order to choose the "best direction" viz. the index i = 1 or 2 that provides the most physically realistic approximation of C 12 by the Archimedean copula C ii . The choice of such a criterion is the objective of the next section.

Physical requirements for a reordering of the path

Let us consider an assembly of spectral lines inside a band ∆ν whose centers:

1/ are statistically independent from each other and 2/ do not depend on the thermophysical state of the gas (line shift is thus not considered). For a nonuniform path made of n uniform layers of lengths L i , i = 1, .., n over each of which absorption coefficients constructed with the help of the preceeding set of lines can be found, the derivative of the transmissivity of the non-uniform path with respect to any L i , i = 1, .., n follows the following inequality:

- ∂τ ∆ν 1..n (L 1 , .., L n ) ∂L i ≤ k i P τ ∆ν 1..n (L 1 , .., L n ) ( 25
)
This inequality on the spatial derivative of the mean equivalent width, given as -ln τ ∆ν (L) , is rather common in the litterature of band models based on transmissivities [START_REF] Young | Band model theory of radiation transport[END_REF]. However, the treatment proposed here introduces two subtleties. The first one is that instead of treating the problem in approximate form, viz. using an estimate of this mean equivalent width as done is most of the litterature on transmissivity based models, the problem is studied here in its general form. This requirement is related to the -distribution formulation.

The second one is based on the observation that as soon as an idea of path reordering is considered, the transmissivity of any non-uniform path involving the same set of absorption coefficients κ i , i = 1, .., n needs to be represented by the same ordering of the indices i = 1, .., n. In some configurations, the emission point will be associated with index 1, in others with index 2, etc. This means that in reordered models, it is not sufficient to consider inequalities of the type of Eq. ( 25) at the last point encountered along the radiation path, but all possible configurations must be considered at once. This is why we have formulated the problem in terms of the system of inequalities Eq. ( 25) and not as:

- ∂τ ∆ν 1..j (L 1 , .., L j ) ∂L j ≤ k j P τ ∆ν 1..j (L 1 , .., L j ), j = 1, .., n (26)
which is the more usual formulation obtained by considering the propagation of radiation emitted at the j-th location along the path up to the layer with index 1. The main difference between the two systems of inequalities is that in the set Eq. ( 26) only transmissivities up to an index j are considered (all indices higher than j + 1 are skipped) whereas in the set Eq. ( 25), all indices up to the highest possible one, n, are treated simultaneously. Notice that if Eq. ( 25) is correct then Eq. ( 26) is valid too: this is due to the fact that Eq. ( 26) can be written as a particular case of Eq. ( 25) for which L j+1 = 0, .., L n = 0. However, the reciprocal may be wrong in the general frame.

In the GWN method, the non-uniform transmissivity τ ∆ν 1..n (L 1 , .., L n ) is modeled as the transmissivity of the gas in state 1 with an equivalent length L 1..n , viz. τ ∆ν 1..n (L 1 , .., L n ) = τ ∆ν 1 (L 1..n ) where the length L 1..n is defined recursively as:

L nn = L n (27)
L i..n = L i + i • τ ∆ν i+1 (L i+1..n ), i = 1, .., n -1 (28)
The previous equation can be also written in terms of effective scaling factors, as defined in the previous section, as:

L i..n = L i + u(L i+1..n )L i+1..n , i = 1, .., n -1 ( 29
)
This process is equivalent to the recursive definition of function C 1..n as:

C nn (X n , 1) = X n (30)
C i..n (X i , .., X n ) = C ii [X i , C i+1,..,n (X i+1 , .., X n )] , i = 1, .., n -1 (31)
For our model to be physically realistic, we may impose the transmissivity

τ ∆ν 1 (L 1.
.n ) to follow the inequalities observed for the true non-uniform path transmissivities, i.e, Eqs. 25. In the following, only the cases of two and three cells are treated. The more general case follows the same steps and its treatment is not informative, but is more tedious.

Our objective is thus to find in which condition on the choice of indices 1, 2 the following inequalities are simultaneously verified:

- ∂τ ∆ν 1 (L 12 ) ∂L 1 ≤ k 1 P τ ∆ν 1 (L 12 ) (32)
- ∂τ ∆ν 1 (L 12 ) ∂L 2 ≤ k 2 P τ ∆ν 1 (L 12 ) (33)
The first inequality is always verified, as a direct application of Chebyshev inequality for integrals and using the fact that κ ν and exp (-κ ν L 1..n ) have opposite monotonicity with respect to wavenumber ν, i.e.:

-

∂τ ∆ν 1 (L 12 ) ∂L 1 = 1 ∆ν ∆ν κ 1 ν exp -κ 1 ν L 1..n dν ≤ 1 ∆ν ∆ν κ 1 ν dν k 1 P × 1 ∆ν ∆ν exp -κ 1 ν L 1..n dν τ ∆ν 1 (L1..n) (34)
The second inequality mostly follows the same steps. Indeed:

∂τ ∆ν 1 (L 12 ) ∂L 2 = ∂τ ∆ν 1 (L 12 ) ∂L 12 ∂L 12 ∂L 2 = ∂τ ∆ν 1 (L 12 ) ∂L 12 ∂ 1 • τ ∆ν 2 (L 2 ) ∂L 2 (35)
Application of Chebyshev inequality to the first derivatives at the RHS then provides:

- ∂τ ∆ν 1 (L 12 ) ∂L 2 ≤ k 1 P τ ∆ν 1 (L 12 ) ∂ 1 • τ ∆ν 2 (L 2 ) ∂L 2 (36)
It follows that a sufficient condition for Eqs. ( 25) and [START_REF] Desmons | Improved information about the vertical location and extent of monolayer clouds from POLDER3 measurements in the oxygen A-band[END_REF] to coincide is:

∂ 1 • τ ∆ν 2 (L 2 ) ∂L 2 ≤ k 2 P k 1 P = ∂ 1 • τ ∆ν 2 (L 2 = 0) ∂L 2 (37)
The second equality arises directly from the definition of function

1 •τ ∆ν 2 (L 2 ).
Consequently, it is sufficient that 1 • τ ∆ν 2 (L 2 ) is concave (its second derivative with respect to L 2 is in this case negative and thus

∂ 1•τ ∆ν 2 (L2) ∂L2
decreases and is lower than its value in L 2 = 0) to ensure that the system of inequalities Eqs.

(32,33) is actually verified. The case of three layers (and more generally of n uniform layers, by induction) follows the same steps for i = 1, 2. Only the case of the third layer differs and is detailed below.

The derivative of the transmissivity of the non-uniform path with respect to the length of the third path is given as:

∂τ ∆ν 1 (L 123 ) ∂L 3 = ∂τ ∆ν 1 (L 123 ) ∂L 123 ∂L 123 ∂L 3 = ∂τ ∆ν 1 (L 123 ) ∂L 123 ∂ 1 • τ ∆ν 2 (L 23 ) ∂L 3 (38)
which is the same as Eq. ( 35). However, following the definition of L 23 as

L 23 = L 2 + 2 • τ ∆ν 3 (L 3
) one has (by application of the chain rule):

∂ 1 • τ ∆ν 2 (L 23 ) ∂L 3 = ∂ 1 • τ ∆ν 2 (L 23 ) ∂L 23 ∂L 23 ∂L 3 = ∂ 1 • τ ∆ν 2 (L 23 ) ∂L 23 ∂ 2 • τ ∆ν 3 (L 3 ) ∂L 3 (39)
Once reported inside Eq. ( 36), one then obtains:

- ∂τ ∆ν 1 (L 123 ) ∂L 3 ≤ k 1 P τ ∆ν 1 (L 123 ) ∂ 1 • τ ∆ν 2 (L 23 ) ∂L 23 ∂ 2 • τ ∆ν 3 (L 3 ) ∂L 3 (40)
A sufficient condition for this inequality to follow Eq. ( 25) with i = 3, and assuming that levels 1 and 2 are organized in such a way that inequality Eq. ( 37) is actually verified, in which case

∂ 1•τ ∆ν 2 (L23) ∂L23 ≤ k 2 P k 1 P
, is:

∂ 2 • τ ∆ν 3 (L 3 ) ∂L 3 ≤ k 3 P k 2 P = ∂ 2 • τ ∆ν 3 (L 3 = 0) ∂L 3 (41)
This follows directly from the fact that if the layers are organized this way, we have:

- ∂τ ∆ν 1 (L 123 ) ∂L 3 ≤ k 1 P τ ∆ν 1 (L 123 ) k 2 P k 1 P k 3 P k 2 P = k 3 P τ ∆ν 1 (L 123 ) (42)
The problem of path reordering can thus be reformulated in the following terms:

find a permutation σ : {1, .., n} → {1, .., n} such that:

∂ 2 σ(i) • τ ∆ν σ(i+1) (L) ∂L 2 ≤ 0 (43)
This problem can benefit from results from copulas theory, as explained in the next section.

Useful mathematical results

As emphasized previously, the GWN method can be naturaly formulated in terms of Archimedean copulas. It is out of the scope of the present paper to provide a full description of this mathematical framework, but only to use some elements of this theory, useful for the purpose of obtaining a reorganization of the gaseous layers in such a way that the constraints provided in the previous section are actually verified or can be at least considered as such. More details on copulas can be found for instance in Ref. [START_REF] Nelsen | An introduction to copulas[END_REF][START_REF] Joe | Multivariate models and multivariate dependence concepts[END_REF].

Following Archimedean copula's theory, the structure determination of hierarchical (or nested) Archimedean copulas (HAC) can be made by analysis of Kendall's Ke correlation coefficient [START_REF] Górecki | An approach to structure determination and estimation of hierarchical Archimedean copulas and its application to bayesian classification[END_REF][START_REF] Górecki | Kendall's tau and agglomerative clustering for structure determination of hierarchical Archimedean copulas[END_REF]. In the case of a bivariate Archimedean copula associated with a twice differentiable generator φ with φ > 0 for all t ∈ [0, +∞), this coefficient is defined as [START_REF] Górecki | An approach to structure determination and estimation of hierarchical Archimedean copulas and its application to bayesian classification[END_REF]:

Ke(φ) = 1 -4 +∞ 0 t [φ (t)] 2 dt (44)
which provides, all calculations done with φ = τ ∆ν , i.e., t ← L and φ ← ∂τ ∆ν (L) ∂L

:

Ke(τ ∆ν ) = 1 ∆ν 2 ∆ν ∆ν κ ν -κ ν κ ν + κ ν 2 dν dν (45)
From this formula, it is obvious that two scaled spectra share the same value of Kendall's coefficient.

Kendall's coefficient is a so-called mesure of concordance [START_REF] Nelsen | An introduction to copulas[END_REF][START_REF] Joe | Multivariate models and multivariate dependence concepts[END_REF]. It can be shown (see Appendix A for details) that given two transmissivity curves τ ∆ν 1 and τ ∆ν 2 :

Ke(τ ∆ν 1 ) ≥ Ke(τ ∆ν 2 ) ⇒ one can find L ≥ 0 such that

∂ 2 2 • τ ∆ν 1 (L) ∂L 2 > 0 (46)
The analysis of Kendall's coefficient thus allows determining the "wrong" direction and thus its opposite can be used as the "good" one. However, it should be noticed that the ordering of the layers in terms of the Kendall's coefficients does not ensure that Eq. ( 43) is actually verified because the two sides of Eq. ( 46) are not rigorously equivalent. But from the analysis of Kendall's coefficient, we can however find one direction for which it is sure that our constraint is not verified and thus choose the opposite direction to construct our iterative scheme.

The physical meaning of Kendall's coefficient can be explored further by considering two spectra, κ 1 ν and κ 2 ν made of the same spectral lines with distinct profiles and linestrengths and such that, for instance, Ke(τ ∆ν 2 ) > Ke(τ ∆ν 1 ). We can construct a third absorption spectrum by κ 3 ν =

k 1 P k 2 P κ 2 ν .
As κ 2 ν and κ 3 ν are scaled by definition, they share the same Kendall's coefficient, viz., Ke(τ ∆ν 2 ) = Ke(τ ∆ν 3 ). Combining this equality with Ke(τ ∆ν 2 ) > Ke(τ ∆ν 1 ), we obtain that Ke(τ ∆ν 3 ) > Ke(τ ∆ν 1 ). The derivative with respect to the length L of the band averaged transmissivities for the spectra κ 1 ν and κ 3 ν are the same at the optically thin limit (they are equal to the Planck mean absorption coefficients which are obviously the same for these spectra) but, following the previous inequality, are associated with distinct values of Kendall's coefficients. This means that values of Kendall's coefficients are more affected by the wings of the spectral lines than by their centers. This idea can be explored further in the case of the SNB-LM model, for which the transmissivity of a gas path of total length L is given as:

τ SNB-LM (L) = exp - β π 1 + 2πk P L β -1 ( 47
)
where k P is the mean absorption coefficient (Planck mean) over the band and β characterizes the overlapping between spectral lines. The corresponding Kendall's coefficient Ke is:

Ke(τ SNB-LM ) = 1 2 + 2β 2 π 2 exp 2β π Ei 2β π - β π ( 48
)
where Ei represents the Exponential integral:

Ei(p) = +∞ 1 exp(-px) x dx (49)
As noticed in Eq. ( 48), in the case of the SNB-LM, the Kendall's coefficient only depends on the overlapping parameter β. This result is in accordance with our previous analysis on the physical interpretation of this coefficient as a characteristic of line wings. Function Ke(τ SNB ) is depicted as a function of the overlapping parameter β in Figure 2. Its maximum value is 1 2 when β → 0 and admits the limit 0 when the gas spectrum is gray, i.e., κ ν = constant. Ke(τ SNB ) decreases with respect to the parameter β. For two thermophysical states for which the gas spectra follow rigorously the assumptions of the SNB-LM, it was shown in Ref. [START_REF] André | An analysis of the symmetry issue in the -distribution method of gas radiation in non-uniform gaseous media[END_REF] that:

1 • τ ∆ν 2 (L 2 ) = k 2 P β 2 k 1 P β 1 L 2 + k 2 P k 1 P 1 - β 2 β 1 β 2 πk 2 P   1 + 2πk 2 P L 2 β 2 -1   (50) = k 2 P β 2 k 1 P β 1 L 2 + +∞ 0 f (s) [1 -exp(-sL 2 )] ds
where:

f (s) = 1 √ 2π β 2 πk 2 P 1 2 k 2 P k 1 P 1 - β 2 β 1 s -3 2 exp - β 2 s 2πk 2 P ( 51
)
From Eqs. ( 48) and ( 50), one can observe that a decrease in the Ke coefficient is associated with an increase of the overlapping parameter β. In this case of increasing β parameter (from cell 2 to cell 1, viz.

β 1 ≥ β 2 ), Eq. (50) shows that function 1 • τ ∆ν 2 (L 2
) is then concave (its second derivative with respect to L 2 is in this case negative). A calculation made in the direction of decreasing Ke provides in this situation a so-called Levy-subordinated Archimedean copula [START_REF] Hering | Constructing hierarchical Archimedean copulas with Lévy subordinators[END_REF][START_REF] André | An analysis of the symmetry issue in the -distribution method of gas radiation in non-uniform gaseous media[END_REF] (this name is related to the fact that 1 • τ ∆ν 2 (L 2 ) is then a Bernstein function, i.e., its derivative is a Laplace transform, as shown in Eqs. (50,51), more usually called in the context of statistical studies, a Laplace exponent of a Levy subordinator [START_REF] Schilling | Bernstein functions: theory and applications[END_REF]).

Application and results

In this section, the accuracy and computational cost of the -distribution method based on two reordering schemes is assessed against LBL calculations.

The first reordering scheme, called β reordering, uses the set of indices of the various atmospheric layers obtained by a reordering of the layers with respect to decreasing values of the parameter β as defined by Eq. ( 5). It thus corresponds to the reordering scheme introduced in Ref. [START_REF] André | An analysis of the symmetry issue in the -distribution method of gas radiation in non-uniform gaseous media[END_REF] and mostly consists of neglecting the effect of the mapping functions Gr since this path ordering strategy is the same as that of the germ models. The second reordering scheme, called Ke reordering, uses a reordering of the gas path with respect to decreasing values of the Kendall coefficient defined by Eq. ( 44). This second reordering scheme thus accounts for possible effects of the Gr functions. It can be noticed here that both reordering schemes suggest a treatment of the gaseous layers from the top of the atmosphere down to the ground. The only difference is that in the case of a reordering with respect to the Ke coefficient, the 11th layer starting from the ground needs to be placed after all other layers. Physical reasons that could explain this result have not been found yet. However, following the analysis provided in Ref. [START_REF] Gordley | Rapid inversion of limb radiance data using an emissivity growth approximation[END_REF] this result suggests that the spectral content of full paths for the TOA down to the ground is closer (more linearly scaled) with the spectrum at the 11th layer in the atmosphere than that close to the ground surface. Results of a Correlated k-distribution (Ck) [START_REF] Lacis | A description of the correlated k distribution method for modeling nongray gaseous absorption, thermal emission, and multiple scattering in vertically inhomogeneous atmospheres[END_REF] model with various number of gray gases based on Gauss-Legendre quadratures at orders from 16 up to 256 are also provided for completeness. For the -distribution method, high resolution mapping functions made of 50,000 values were used. This large number of values has no influence on the calculation cost of the technique but requires more memory space for the application of the method. Using high resolution mapping functions allows focusing on the effect of the reordering scheme only, which is the main aim of the present work.

Two types of conditions are studied. The first one considers clear sky nonscattering situations, for which transmissivities of non-uniform paths are compared (this is the situation for which a perfect miror is put at some given altitude). The second type of comparison involves scattering events (molecular and/or by clouds). This more realistic situation was treated by incorporating the -distribution model inside the atmospheric code 3DMCPOL [START_REF] Cornet | Three-dimensional polarized monte carlo atmospheric radiative transfer model (3DMCPOL): 3D effects of polarized visible reflectances of cirrus cloud[END_REF] based on a Monte Carlo RTE solver, developed at the Laboratory of Atmospheric Optics in Lille, France. All calculations were made under the standard Mid-Latitude Summer (MLS) atmospheric profile [START_REF] Mcclatchey | Optical properties of the atmosphere[END_REF]. The transmission in the A-band of oxygen and the filter response function used in this study are shown in Figure 3 in the case of an Air Mass Factor of 1 (one single Nadir path in the atmosphere).

This spectral response is characteristic of instruments that use the A-band of oxygen, with measurements both in the O 2 absorption and outside the absorption band or with a weak O 2 absorption (such as MERIS [START_REF] Rast | The ESA medium resolution imaging spectrometer MERIS: A review of the instrument and its mission[END_REF], POLDER [START_REF] Buriez | Cloud detection and derivation of cloud properties from POLDER[END_REF] or the future instrument 3MI [START_REF] Manolis | The MetOp second generation 3MI mission[END_REF]). The O 2 A-band is widely used in remote sensing (see [START_REF] Barton | Remote measurement of surface pressure using absorption in the oxygen A-band[END_REF][START_REF] Dubuisson | Surface pressure estimates from satellite data in the oxygen A-band: Applications to the MOS sensor over land[END_REF][START_REF] Tran | Line mixing and collision-induced absorption by oxygen in the A band: Laboratory mesurements, model, and tools for atmospheric spectra computations[END_REF] and references therein) for the determination of the atmospheric pressure near the ground or above clouds, as well as for the determination of the altitude at which scattering aerosols can be found. The LBL code used for the high resolution calculations is described in [START_REF] Dubuisson | High spectral resolution solar radiative transfer in absorbing and scattering media: application to the satellite simulation[END_REF]. Results, assuming a surface albedo of 0.2 fixed in all simulations involving scattering events, are provided in terms of the Air Mass Factor (AMF) defined as:

AMF = 1 cos(θ 0 ) + 1 cos(θ v ) (52)
where θ 0 and θ v are the solar and satellite viewing angles respectively. Relationships between values of AMF used in the various figures and the corresponding couples of angles θ 0 and θ v are given in Table 1. The largest relative errors for the -distribution method are below 0.8 % in all the cases considered. The mean value of relative errors evaluated over the full range of pathlengths considered in the case of Figure 6 (240 values) is 10.2 10 -4 for Ck vs LBL whereas it is 7.54 10 -4 /beta (7.4766 10 -4 / Kendall) fordistribution vs LBL. This means that at a global scale the -distribution method is more accurate than the correlated k-distribution. However, the -distribution method provides in some cases (for some values of pathlengths) results less accurate than correlated k-distributions.The effect of β vs Ke reordering scheme is noticeable in the lower layers of the atmosphere, where the orders of the paths are different. Nevertheless, the Ke reordering scheme is found to provide slightly more accurate estimates of the transmissivities of the non-uniform paths than the β reodering scheme in all the cases treated. The two methods share in this case the same computational cost as only the orders of treatment of the set of lenghts L 1 , ...L n differs between the two calculations.

AMF

Scattering conditions -case of cloudy atmospheres

All calculations in the following scattering situations were made with the code 3DMCPOL. This Monte Carlo radiative transfer simulator was compared recently to several other atmospheric codes in Refs. ( [START_REF] Emde | IPRT polarized radiative transfer model intercomparison project -phase A[END_REF], [START_REF] Emde | IPRT polarized radiative transfer model intercomparison project -three-dimensional test case (phase B)[END_REF]).

These cases treat radiative transfer in the atmosphere considering scattering events. They were simulated by separating the treatment of scattering from molecular absorption. The method used in combination with the -distribution approach is the following: 1/ a photon is launched from space and propagates in the atmosphere until a first scattering event occurs, 2/ one part of the radiative energy associated with the photon is sent to the imaging device, including the process is iterated (back to step 2) until the remaining energy of the photon reaches some cut-off or leaves the atmosphere. This method is applied to a significant amount of (10 6 ) photon paths in order to provide a proper statistical estimate of the simulated radiance.

Pure Molecular Scattering

The first scattering configuration only involves molecular scattering (Rayleigh).

Results are depicted in Figures 7 (calculation cost) and 8 (accuracy). In these cases, both -distribution methods (β and Ke reordering schemes) provide the highest improvements in terms of computational cost (higher than 1,000) compared to LBL calculations. On the same cases, the CPU time ratio for the Ck model with 16 coefficients remains lower than 800. Furthermore, thedistribution models based on the Ke reordering scheme appears to be the most accurate, confirming the results provided in the previous section (cases without scattering).

Scattering by clouds -mono-and multi-layers

O 2 A-Band information is widely used to retrieve cloud [START_REF] Desmons | Improved information about the vertical location and extent of monolayer clouds from POLDER3 measurements in the oxygen A-band[END_REF][START_REF] Merlin | Cloud information content analysis of multi-angular measurements in the oxygen A-band: application to 3MI and MSPI[END_REF] and aerosol [START_REF] Zeng | Constraining aerosol vertical profile in the boundary layer using hyperspectral measurements of oxygen absorption[END_REF] geometrical properties. The second series of tests was consequently made in cloudy atmosphere composed of uniform layers with cloud droplet size distribution of effective radius of 10µm. As cloud are optically dense media, scattering processes are numerous and the photons path lengths can be high in the cloud layers. In order to evaluate this effects on the -distribution method at different In addition, the accuracy is higher for -distribution and in particular for the Kendall reordering in the low cloud case (Figure 9). For the high cloud case, the accuracy of the -distribution method is in most cases, higher than the Finally, -distribution was applied to a multi-layer cloud made of the two mono-layer situations studied previously. In this case, the reflections between the two clouds increase the lengths of the path of the photons in the most absorbing parts of the atmosphere. This configuration is thus a complete test for the -distribution method. Results, depicted in Figure 13, show again a better accuracy than the k-distribution for the airmass factors of 2 and 4 but, driven certainly by the high cloud results, the accuracy is lower for the airmass factor of 16. However, as for other cases, the computation cost is reduced by 30% for the airmass factor 2 and about 60% for the airmass factors 4 and 16, as shown in figure 14. Archimedean copulas. This allows using directly elements from this theory to propose a reorganization of the paths in such a way that some physical constraints, discussed in the paper, are actually verified.

Once reordered, it is shown that the method is capable of reproducing mean radiance calculated LBL with an accuracy higher than 1% at a tiny fraction of the computational cost of a high resolution model (up to 1,000), both in clear sky and cloudy simulations. In addition, it is shown that the present - distribution method can compete and in some cases outperform the more usual Ck approach both in terms of accuracy and computational cost.

The main defect of the present technique is its formulation in transmission form which restricts its application to methods based on the integral formulation of the RTE such as ray-tracing, but do not allow its use together with a discrete ordinate [START_REF] Stamnes | Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media[END_REF], Matrix operator or adding-doubling solver of the RTE, for instance. However, due to its high efficiency, the technique can be a good candidate for applications that can accomodate with such a formulation, including Monte-Carlo (as used here) or quasi-Monte Carlo methods. For these kinds of RTE solvers, the -distribution approach can provide accurate results at small computational cost and thus find applications in near real time sensing problems.

It should be noticed that the present paper only focuses on the application of the -distribution approach in the visible range. Evaluations of the technique in the thermal infrared region has always been conducted in high temperature applications and found to provide more accurate results than the Ck method for the treatment of gas mixtures: applications of the method for thermal infrared radiative transfer in the atmosphere are scheduled as a continuation of the present work. The work presented here is for moderate resolution imagers.

A preliminary study shows that the developed method works for high resolution spectrometers, such as IASI-NG, too: high resolution calculations will be explored further and these preliminary results will be confirmed in future works.

Appendix A. Theoretical justification of relationship (46)

The aim of this appendix is to provide a justification of the relationship (46). It is based on the following definitions and theorems taken from Nelsen [START_REF] Nelsen | An introduction to copulas[END_REF]. The notation < c represents the concordance ordering, i.e., C 1 < c C 2 if C 1 (x, y) ≤ C 2 (x, y) for any couple x, y inside the unit square. The concept of 1. µ is defined for every pair of random variables Kendall's coefficients as defined by Eq.( 44) are measures of concordance. A proof can be found for instance in Nelsen [START_REF] Nelsen | An introduction to copulas[END_REF]. The following theorem can be also found on the same reference.

Theorem 1. Let C 1 and C 2 be Archimedean copulas generated respectively by 1 (inverse of τ ∆ν 1 ) and 2 (inverse of τ ∆ν 2 ) in Ω. If 1 • τ ∆ν 2 is concave then

C 1 < c C 2 .
Accordingly, the combination of theorem 1 with property 5. from the definition of measures of concordance shows that if This result is typical of one parameter copula families as generated by SNB-LB models which yield copulas that depend on a single parameter (the β coefficient in the present case). It is generalized here but the principle is founded on the same kind of mathematical analysis.

Acronyms

  at an altitude of 120 km viz. a total pressure of 2.10 -5 hPa and a temperature of 380 K) for high values of the transmissivities (X close to 1) are far from the identity line. This is because at this optically thin limit, large errors can be observed while trying to estimate transmissivities in the Doppler-dominant regime by a model for Lorentz lines.

Figure 1 :

 1 Figure 1: Anatomy of mapping functions Gr -parameter X represents the value of the germ function, that contains the information on the gas path length.

Figure 2 :

 2 Figure 2: Kendall coefficient as a function of overlapping parameter β

Figure 3 :

 3 Figure 3: Filter response function and gas path transmissivity accounting for the filter response for AMF=1

Figures 4 to 6 .

 6 Figures 4 to 6. In these figures, the accuracy of both -distribution methods is found higher than the Ck model with 256 gray gases. It was observed, but the results are out of the scope of the present work, that this Ck method provides the same results as a true correlated LBL model obtained by associating at the spectral scale values of absorption coefficient inside different layers through the same relationship, based on equality of the corresponding cumulative distribution functions, as the gray gases Ck model. Results of this technique can thus be considered as the limit in terms of accuracy of a model based on the socalled correlation assumption for this specific problem (O 2 A-band in the MLS configuration).

Figure 4 :

 4 Figure 4: Non-uniform transmission curves for AMF=2

Figure 5 :

 5 Figure 5: Non-uniform transmission curves for AMF=4

Figure 6 :

 6 Figure 6: Non-uniform transmission curves for AMF=16

Figure 7 :

 7 Figure 7: Gain in terms of CPU time -RTE solver: 3DMCPOL Monte Carlo code -Pure molecular Scattering

Figures 10 and 12

 12 Figures 10 and 12 show that, as previously, the computation cost is reduced significantly with the -distribution method in comparison with the LBL and Ck methods. The computation time is decreased by a factor comprised between 500 and 1000 compared to the reference LBL calculation. It corresponds to a gain of about 14 to 84 % compared to k-distribution computation with 16 coefficients.

Figure 8 :

 8 Figure 8: Relative Error / LBL -Pure molecular Scattering

Figure 9 :

 9 Figure 9: Relative Error / LBL -cloud between 1 and 2 km with COT = 10

Figure 10 :

 10 Figure 10: Gain in terms of CPU time -RTE solver: 3DMCPOL Monte Carlo code -cloud between 1 and 2 km with COT = 10

Figure 11 :

 11 Figure 11: Relative Error / LBL -cloud between 10 and 11 km with COT = 2

Figure 12 : 1 .

 121 Figure 12: Gain in terms of CPU time -RTE solver: 3DMCPOL Monte Carlo code -cloud between 10 and 11 km with COT = 2

Figure 13 :

 13 Figure 13: Relative Error / LBL -double layer of clouds between 1 and 2 km with COT = 10 and between 10 and 11 km with COT = 2

Figure 14 :

 14 Figure 14: Gain in terms of CPU time -RTE solver: 3DMCPOL Monte Carlo code -double layer of clouds between 1 and 2 km with COT = 10 and between 10 and 11 km with COT = 2

Table 1 :

 1 Relationship between viewing angles, in degrees, and Air Mass Factor (AMF)viewing azimuth angle is 180 degrees from the sun azimuth in the specular direction4.1. Clear sky non-scattering conditionsThese cases are mostly used to evaluate the accuracy of the various approximate methods in situations of pure gaseous absorption. Results are depicted in

		θ 0	θ v
	2	0	0
	4	60	60
	16	84.4 80

  LM Statistical Narrow Band model (Lorentz lines, Makmus' distribution of linestrengths)

	4A/OP Automatized Atmospheric Absorption Atlas
	RTE Radiative Transfer Equation
	Ck Correlated k-distribution method
	GWN Godson-Weinreb-Neuendorffer's method
	SNB-

Acknowledgments

This work has been supported by the Programme National de Télédétection Spatiale (PNTS, http://www.insu.cnrs.fr/pnts grant PNTS-2017-04).


Xingran Xu

Huawen Xu

S Mandal

R Banerjee

Sanjib Ghosh

T C H Liew

Interaction induced bi-skin effect in an exciton-polariton system

The non-Hermitian skin effect can be realized through asymmetric hopping between forward and backward directions, where all the modes of the system are localized at one edge of a finite 1D lattice. However, achieving such an asymmetric hopping in optical systems is far from trivial. Here we show theoretically that in a finite chain of 1D exciton-polariton micropillars with symmetric hopping, the inherent non-linearity of the system can exhibit a bi-skin effect, where the modes of the system are localized at the two edges of the system. To show the topological origin of such modes, we calculate the winding number.

Introduction.-Exciton-polaritons are bosonic quasiparticles arising from the strong coupling between excitons and photons in microcavities [START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Deng | Exciton-polariton bose-einstein condensation[END_REF][START_REF] Keeling | Excitonpolariton condensation[END_REF][START_REF] Byrnes | Excitonpolariton condensates[END_REF]. The combination of photons and excitons allows polaritons to possess the characteristics of both its root parts, for example, low effective mass from the photonic part and strong nonlinearity from the excitonic part. Moreover, due to the photons escaping from the microcavity, exciton-polaritons are noticeably an open quantum (non-Hermitian) system, which requires additional pumping to maintain a steady state in the system [START_REF] Savvidis | Off-branch polaritons and multiple scattering in semiconductor microcavities[END_REF]. This character makes the exciton-polariton system an ideal platform to study non-Hermitian transitions with gain and loss [START_REF] Nalitov | Optically trapped polariton condensates as semiclassical time crystals[END_REF][START_REF] Xu | Nonequilibrium landau-zener tunneling in excitonpolariton condensates[END_REF]. Although various research has been done in the Hermitian regime to realize polariton topological phases with the interplay between Zeeman shift resulting from the application of magnetic field and the transverse electric-transverse magnetic (TE-TM) splitting of the photonic modes [START_REF] Bardyn | Topological polaritons and excitons in garden-variety systems[END_REF][START_REF] Nalitov | Polariton Z topological insulator[END_REF][START_REF] Li | Lieb polariton topological insulators[END_REF][START_REF] Gulevich | Kagome lattice from an exciton-polariton perspective[END_REF][START_REF] Klembt | Exciton-polariton topological insulator[END_REF], accounting for nonlinearity [START_REF] Bardyn | Chiral bogoliubov excitations in nonlinear bosonic systems[END_REF][START_REF] Sigurdsson | Spontaneous and superfluid chiral edge states in exciton-polariton condensates[END_REF], and using the polarization splitting of elliptical micropillars [START_REF] Banerjee | Realization of hofstadter's butterfly and a one-way edge mode in a polaritonic system[END_REF], several works have been reported in recent years to explore the non-Hermitian physics in the exciton-polariton system. Exceptional points (EPs) in the exciton-polariton system have been realized in Refs. [START_REF] Gao | Observation of nonhermitian degeneracies in a chaotic exciton-polariton billiard[END_REF][START_REF] Gao | Chiral modes at exceptional points in excitonpolariton quantum fluids[END_REF][START_REF] Gao | Continuous transition between weak and ultrastrong coupling through exceptional points in carbon nanotube microcavity exciton-polaritons[END_REF]. Recently, the measurement of non-Hermitian topological invariants [START_REF] Su | Direct measurement of a non-hermitian topological invariant in a hybrid light-matter system[END_REF], topological end mode lasing [START_REF] Comaron | Non-hermitian topological endmode lasing in polariton systems[END_REF], and non-reciprocal polaritons [START_REF] Mandal | Nonreciprocal transport of exciton polaritons in a non-hermitian chain[END_REF] have been reported.

Non-Hermitian physics has attracted a tremendous attention in recent years due to the discovery of the non-Hermitian skin effect, where all the modes are localized to one end of a lattice . This is different to the case of Hermitian systems, where modes within a topological bandgap are localized at the edges of a finite sample due to bulk-boundary correspondence. Due to non-Hermiticity, the Bloch theory does not even hold approximately for finite sized non-Hermitian systems and one must turn to the generalized Brillouin zone (GBZ) theory based on a complex momentum to explain the non-Hermitian skin effect [START_REF] Okuma | Topological origin of non-hermitian skin effects[END_REF][START_REF] Kunst | Biorthogonal bulk-boundary correspondence in non-hermitian systems[END_REF][START_REF] Kawabata | Symmetry and topology in nonhermitian physics[END_REF][START_REF] Gong | Topological phases of non-hermitian systems[END_REF][START_REF] Yao | Edge states and topological invariants of non-hermitian systems[END_REF][START_REF] Lieu | Topological phases in the non-hermitian suschrieffer-heeger model[END_REF]. To realize the skin effect, the simplest model is the so-called Hatano-Nelson(HN) model [START_REF] Hatano | Localization transitions in non-hermitian quantum mechanics[END_REF][START_REF] Hatano | Vortex pinning and non-hermitian quantum mechanics[END_REF] without disorder where the hopping in a lattice is different in different directions.

Along with the same idea of the HN model, the Su-Schrieffer-Heeger(SSH) model with asymmetric hopping has been proposed or realized in many systems like waveguides [START_REF] Zeuner | Observation of a topological transition in the bulk of a non-hermitian system[END_REF][START_REF] Yokomizo | Non-bloch band theory and bulk-edge correspondence in nonhermitian systems[END_REF], photonic lattices [START_REF] Xiao | Non-hermitian bulk-boundary correspondence in quantum dynamics[END_REF][START_REF] El-Ganainy | The dawn of non-hermitian optics[END_REF], circuits [START_REF] Helbig | Generalized bulk-boundary correspondence in non-hermitian topolectrical circuits[END_REF] etc. The non-Hermitian topological bulk-boundary correspondence can be obtained by using the GBZ theory and by solving the boundary equations with complex momentum [35-37, 49, 57-60]. Due to the chiral symmetry, the system will have an energy pair of (E,-E) and they will collapse at zero energy modes. The winding number can be well defined by the calculation with the periodic boundary energies with the GBZ. Meanwhile, more and more work is going to explain the topological origin of the skin modes and the GBZ is not necessary for line gap non-Hermitian Hamiltonians where the basis can be enlarged to construct a Hermitian Hamiltonian [START_REF] Hua | Anatomy of skin modes and topology in non-hermitian systems[END_REF][START_REF] Okuma | Topological origin of non-hermitian skin effects[END_REF][START_REF] Zhang | Correspondence between winding numbers and skin modes in nonhermitian systems[END_REF][START_REF] Gong | Topological phases of non-hermitian systems[END_REF][START_REF] Liu | A new way to construct topological invariants of non-hermitian systems with the non-hermitian skin effect[END_REF][START_REF] Jong Yeon Lee | Topological correspondence between hermitian and non-hermitian systems: Anomalous dynamics[END_REF]. Although the winding in the whole Bloch region is still zero for zero energy pairs, the system can still have a topological non-trivial phase if energy pairs can be made to collapse at non-zero (complex) energy.

In this letter, we use the particle-hole symmetry of the fluctuation theory of exciton-polariton systems to realize "bi-skin" modes where all the wavefunctions are localized at both sides of a one-dimensional lattice. Being compatible with recent experiments, we consider a non-equilibrium polariton condensate corresponding to a plane-wave with non-zero momentum. By the Bogoliubov-de-Gennes(BdG) transformation [START_REF] Xu | Spinor polariton condensates under nonresonant pumping: Steady states and elementary excitations[END_REF][START_REF] Banerjee | Coupling between exciton-polariton corner modes through edge states[END_REF], the fluctuations of the condensate will have two modes with particle-hole symmetry. By considering a suitable BdG lattice by tight-binding theory, the bi-skin effect will appear with symmetric hopping between sites. To classify the symmetry and the topologies of the system, we calculate the winding number of the system without assuming the GBZ theory. Because of the particle-hole symme- try [START_REF] Okuma | Topological phase transition driven by infinitesimal instability: Majorana fermions in non-hermitian spintronics[END_REF], the energy spectrum will not collapse at zero modes but some purely imaginary energies. The different topological phases can be well defined by the winding of the EPs and the localization of the wavefunction. Finally, we analyze the stability of the system and calculate the energy dispersion, which can be observed via photoluminescence [START_REF] Utsunomiya | Observation of bogoliubov excitations in exciton-polariton condensates[END_REF][START_REF] Pieczarka | Ghost branch photoluminescence from a polariton fluid under nonresonant excitation[END_REF] or four-wave mixing spectroscopy [START_REF] Kohnle | From single particle to superfluid excitations in a dissipative polariton gas[END_REF].

Model.-We consider a simple chain with symmetric hopping J between nearest neighbours, formed by exciton-polariton micropillars, as shown in Fig. 1(a). Such chains have been realized experimentally [START_REF] Winkler | Collective state transitions of exciton-polaritons loaded into a periodic potential[END_REF]. Typically the behavior of a polariton condensate is studied, which itself has no topological features in such a system. The system can be described by the driven-dissipative nonlinear Gross-Pitaevskii (GP) equation

i ∂ψ n ∂t = (∆ -iγ C ) ψ n + J(ψ n-1 + ψ n+1 ) + g |ψ n | 2 ψ n + F n , (1)
Here, ∆ = E 0 -ω p is the onsite detuning between the bare polariton mode having energy E 0 and a coherent drive F n with frequency ω p . g represents the nonlinear interaction strength and γ C is the effective dissipation. In order to find the optimum parameter range we move to the dimensionless units by making the following substitutions:

t → t /γ C , J → J/γ C , ∆ → ∆/γ C , ψ n → ψ n γ C /g, and F n → F n g/γ 3
C . One of the main ingredients of our scheme is a plane-wave like stationary state ψ 0,n = Ae ikpn (as shown in Fig. 2(a) ) realized by spatially modulating a coherent field [START_REF] Ohadi | Synchronization crossover of polariton condensates in weakly disordered lattices[END_REF][START_REF] Alyatkin | Optical control of couplings in polariton condensate lattices[END_REF][START_REF] Pickup | Synthetic band-structure engineering in polariton crystals with non-hermitian topological phases[END_REF][START_REF] Carusotto | Probing microcavity polariton superfluidity through resonant rayleigh scattering[END_REF] and the appropriate choice of F n is obtained by setting the time derivative in Eq. (1) to zero:

F n = -∆ + i -2J cos(k p ) -|A| 2 Ae ikpn , ( 2
)
where A is the square root of the density of the condensate and k p is the pumping momentum and n = 1 or N . Because of the boundary, we should set the first and the last site pumping by replacing 2J cos(k p ) in Eq. ( 2) with Je ikp and Je -ikp , respectively. To show the bi-skin effect we consider the fluctuations, by introducing the fluctuation amplitudes u n and v n , and writing ψ n = ψ (0,n) + u n e iEt/ + v * n e -iE * t/ . The BdG modes as shown in Fig. 1(b) are derived by substituting ψ n into the mean-field GP equation, which yields the time-independent equations:

Eu n = (∆ -i) u n + J (u n+1 + u n-1 ) + ψ 2 (0,n) v n , (3) Ev n = -(i + ∆ ) v n -J (v n+1 + v n-1 ) -ψ * 2 (0,n) u n .( 4
)
Here, ∆ = ∆ + 2 ψ (0,n) 2 is an effective detuning. The BdG transformation ensures that u n and v n modes have the particle-hole symmetry, so the hopping energies for u n and v n are ±J and the inter hopping energies in one unit cell are ±ψ ( * )2 (0,n) . We need to mention that the hopping strength between different pairs of neighboring sites have the same magnitude but differ in their phases because of the symmetry. The Bogoliubov matrix H n with N sites obtained by the Eqs. ( 3)-( 4) can be diagonalized by

V -1 Ĥn V = diag E N ×N , -E * N ×N
with V the matrix formed by eigenvectors, meanwhile,

S x V -1 Ĥn V S x = diag E N ×N , -E * N ×N
where S x V can also be the eigenvectors of the system [START_REF] Xu | Spinor polariton condensates under nonresonant pumping: Steady states and elementary excitations[END_REF][START_REF] Banerjee | Coupling between exciton-polariton corner modes through edge states[END_REF], here, S x,y,z = ÎN,N ⊗ σ x,y,z and σ x,y,z are the Pauli matrices. On the other hand, the system has pseudo-Hermiticity [START_REF] Kawabata | Symmetry and topology in nonhermitian physics[END_REF][START_REF] Yokomizo | Non-bloch band theory in bosonic bogoliubov-de gennes systems[END_REF] defined as S z Ĥn S z = Ĥn † . The particle-hole symmetry and pseudo-Hermiticity guarantee that the system can have a 1D topological transition [START_REF] Okuma | Topological origin of non-hermitian skin effects[END_REF][START_REF] Yokomizo | Non-bloch band theory in bosonic bogoliubov-de gennes systems[END_REF] and the energy spectrum in real space is highly symmetrical.

The bi-skin modes.-The stationary states of the bulk polaritons can have non-zero momentum induced by the specific F n pumping laser and the nonlinearity makes the fluctuations modes u n and v n have a momentum shift. The momentum shift will just influence the phases of the energy spectrum under the periodic boundary condition (PBC), however, the energy spectrum under the open boundary condition (OBC) will completely change for the fluctuations at the edge can not go further and then become localized.

The model is equivalent to say that the hopping energies between the u n and v n sites are different, but we stress that this is an effect induced by nonlinear interactions, where the fluctuations are affected by the non-zero momentum of the considered polariton condensate (stationary state). The underlying physical system is still none other than a regular polariton lattice with symmetric Hermitian hopping between lattice sites. As is shown in Figs. 2(b) and (c), all eigenstates for u n and v n sites are localized at both sides. The energy spectrum of u n sites are (E, E * ) -i pairs, while the energy spectrum of v n sites are (-E * , -E)-i pairs. All eigenvectors of E -i and -E * -i will be localized at one side and wavefunctions of the E * -i and -E -i will be localized at the opposite side. Remarkably, the localization here does not mean the polariton's density localization but the fluctuations on the top of the polaritons. Stability.-Our model is based on the fluctuations of exciton-polaritons with a non-zero momentum stationary state. Because of the symmetry of the system, the imaginary parts of the eigenenergies can be beyond zero. However, with sufficient dissipation, there will be a total energy shift iγ C that keeps all imaginary energy components below zero. 3)-( 4). Other parameters are the same as in Fig. 2.

ized at the edges of the system, while we find that those with least negative imaginary part are localized on the left side.

Because the imaginary parts of the energy spectrum are larger than the real parts, a fundamental problem is whether the spectrum can be observed experimentally.

To prove this, we numerically calculate the BdG lattice's energy dispersion in Fig. 3(c). We take u n and v n as perturbations with small random number compared with ψ (0,n) and solve the BdG equations. The dashed fitting curves have the same energy dispersion with diagonalization results but doubly shrinks the period for two inner sites here. The imaginary parts of the energy dispersion are different for the linewidth of the energy dispersion are different as a function of momentum. Remarkably, we solve the Eqs. ( 3)-( 4) without coherent pumping here, which reveals our model is universal and can be extended to other systems.

Non-Hermitian topological invariants-The effective Hamiltonian can be written into the momentum space with use of the Fourier transformation

u n → 1 √ N k âk e ikn and v n → 1 √ N
k bk e ikn with âk and bk the annihilation operators for the momentum k, giving

H k = ∆ + 2J cos (k) A 2 -A 2 -∆ -2J cos (k -2k p ) -i Î, ( 5
)
where Î is the 2 × 2 identity matrix and the periodic bounary condition is considered and the basis vector is (â k , bk-2kp ). The imaginary energy shift i can be ignored for the theoretical study, however, it can determine the stability of two localization and is necessary to be observed by the energy dispersion. The non-Hermitian topology of H k with respect to the reference point E R is equivalent to the Hermitian topology of the following doubled Hamiltonian [START_REF] Hua | Hybrid higher-order skin-topological modes in nonreciprocal systems[END_REF][START_REF] Gong | Topological phases of non-hermitian systems[END_REF][START_REF] Okugawa | Topological exceptional surfaces in non-hermitian systems with paritytime and parity-particle-hole symmetries[END_REF] with the basis vector ṽ = (â k , bk , âk+2kp , bk-2kp ) can be obtained in the supplementary materials where the open boundary condition is considered.

Hk = [∆ + 2J cos(k)] σz igA 2 σy igA 2 σy ∆ σz + h(k, k p ) , (6)
where h(k, k p ) = 2J cos(k) cos(2k p )σ z -sin(k) sin(2k p ) Î and Î is the 2 × 2 identity matrix. The absolute value of the fluctuations (on top of the condensate) hopping from u n to v n will move to the -k p direction (â k b † k-2kp term ), meanwhile the polariton fluctuations hopping from v n to u n will move to the +k p direction ( bk â † k+2kp term). Therefore, the absolute squre of the eigenstates of the BdG lattice are localized at both sides of the lattice as shown in Figs. 2(b)-(c).

The construction of the Hermitian Hamiltonian in the enlarged Hilbert space not only allows interpretation with the Bloch theory [START_REF] Zhang | Correspondence between winding numbers and skin modes in nonhermitian systems[END_REF][START_REF] Hua | Hybrid higher-order skin-topological modes in nonreciprocal systems[END_REF][START_REF] Gong | Topological phases of non-hermitian systems[END_REF][START_REF] Okugawa | Topological exceptional surfaces in non-hermitian systems with paritytime and parity-particle-hole symmetries[END_REF], but also highlights the differences of the four basis operators under the open boundary condition. The boundary will let the momentum be cut off at the different values and cause the fluctuations to move to the edge and then be localized. For example, the range of k is (-π, π), however, the range of k ± 2k p is (-π ± 2k p , π ± 2k p ). Some related work has been done in the SSH model with asymmetric hopping and interpreted with the GBZ theory [START_REF] Liu | Helical damping and dynamical critical skin effect in open quantum systems[END_REF][START_REF] Yokomizo | Non-bloch band theory in bosonic bogoliubov-de gennes systems[END_REF][START_REF] Yang | Non-perturbative breakdown of bloch's theorem and hermitian skin effects[END_REF].

There are two EPs in the energy spectrum, so we need to calculate each winding separately. By unitary transition, the enlarged Hamiltonian can be written into a two blocks off-diagonal matrix and for each one we can define the winding around the reference energy E R [START_REF] Zhang | Correspondence between winding numbers and skin modes in nonhermitian systems[END_REF]

W E = 1 2πi π -π d log det [H(k) -E R I 2×2 ] dk dk. (7)
The total winding number is zero for the vanishing of E and E * because the winding direction is completely opposite. However, we can get the integer winding number ±1 or 0 in different topological regions. In our model, the reference energy is purely imaginary energies for the gapless E and -E * . The energies as a function of Hopping energies in the real space with OBC are shown in Figs. [START_REF] Byrnes | Excitonpolariton condensates[END_REF] (a) and (c). The critical points are labeled for the callapse of the imaginary energies. The different topological phases can be distinguished by the eigenenergy spectrum of Eq.( 5) and how many OBC energies in real space are encircled by the energy band in momentum space with PBC (see supplementary material):

• Nontrivial phase 0 < |J| < (A 2 +∆ ) sec(kp) 2
: All eigenstates are localized at both sides and all OBC energies are circled by the PBC energies. or no eigenstates are localized and no EPs are encircled by the periodic boundary spectrum.

• Intermediate phase (A 2 +∆ ) sec(kp) 2 |J| < (A 2 -∆ ) sec(kp
Here, the critical points are calculated by the phases of the eigenvalues of Eq. ( 5) and this result is corresponding with the winding number obtained by Eq. ( 7) (see the supplementary materials). If the PBC energy dispersion has different linewidth (imaginary parts of the energies) as a function of momentum, the skin modes will appear.

Although the skin modes can appear in the above three phases, the topological invariants are only ±1 for two EPs in the nontrivial phase and the intermediate phase.

The energies profile under the PBC can not encircle all the energies under the OBC in the intermediate phase.

Therefore, all eigenstates are localized only in non-trivial phase and the quantities of the localized eigenstates are corresponding with the quantities of OBC energies encircled by the PBC energies.

In the considered model all the calculations are done using the dimensionless variables. To get an idea about the physical parameters, we set γ C = 0.1 meV, which corresponds to a polariton lifetime around 3 ps. The coupling between the neighbouring micropillars becomes J=40µeV, which can be tuned by adjusting the overlap between the micropillars [START_REF] Michaelis De Vasconcellos | Spatial, spectral, and polarization properties of coupled micropillar cavities[END_REF]. The non-linearity induced blueshift becomes 0.1 meV which is a routine observation in GaAs based samples [START_REF] Estrecho | Direct measurement of polariton-polariton interaction strength in the thomas-fermi regime of exciton-polariton condensation[END_REF]. Given the room temperature polariton micropillar chain has been demonstrated experimentally [START_REF] Su | Observation of exciton polariton condensation in a perovskite lattice at room temperature[END_REF][START_REF] Dusel | Room temperature organic exciton-polariton condensate in a lattice[END_REF], our scheme is also compatible with perovskite and organic samples. To make it more experimental friendly, the fluctuations can be made to have different polarization from that of the steady state by proper choice of the resonant driving and excitonphoton detuning [START_REF] Sigurdsson | Spontaneous and superfluid chiral edge states in exciton-polariton condensates[END_REF][START_REF] Mandal | Antichiral edge states in an exciton polariton strip[END_REF].

Conclusion-We consider the behaviour of fluctuations on a polariton mean-field stationary state in a onedimensional lattice of coupled micropillars. By taking the stationary state as a plane wave, which can be resonantly injected with a suitably patterned optical field, we find that interactions in the system allow the presence of a bi-skin effect where fluctuations are localized at both edges of the lattice. Even though the underlying hopping of the system is Hermitian, the interactions allow an effective non-reciprocal coupling between particlehole fluctuations. This results in a non-trivial topology, confirmed by calculations of the winding number.

STATIONARY AND EXCITED STATES OF THE EXCITON POLARITONS UNDER RESONANT PUMP

In this section, we consider the behaviour of the Gross-Pitaevskii equation directly in the presence of fluctuations, which are modelled numerically. This will allow us to study the stability of the polariton mean-field solutions and access the fluctuation directly. We begin with the Gross-Pitaevskii equation for a driven-dissipative one dimensional micropillar chain, identical to Eq. ( 1) of the main text but accounting for an additional noise term [START_REF] Winkler | Collective state transitions of exciton-polaritons loaded into a periodic potential[END_REF]

i ∂ψ n ∂t = (∆ -iγ C ) ψ n + <m> J m,n ψ m + g |ψ n | 2 ψ n + F n + i dψ st dt . ( 8
)
Here, ψ st = 2DdW accounts for the fluctuations induced by white noise, where dW is a Gaussian random variable and D is the strength of the noise. The sum over < m > represents summation over nearest neighbour sites. Other parameters are mentioned in the main text. As is shown in Fig. 6, we plot the time evolution of the wavefunctions of a finite-sized lattice with drive field in Eq. ( 2) and the boundary is considered in the main text. If there is no dissipation in the system, all fluctuations with different energies will be excited (see Fig. 6(e)). However, if the dissipation gets larger, we can observe a time crystal like density distribution, where the intensity has a self-induced periodic oscillation and the parts of the energies get excited [START_REF] Nalitov | Optically trapped polariton condensates as semiclassical time crystals[END_REF] (see Figs. 6(b),(c),(f),(g)). If the dissipation is larger enough that can shift all imaginary parts of the energies below zero, the system can have a stable plane-wave solution and only the zero energies modes.

As is shown in Fig. 7(a)-(c), the energy dispersion of the free particle, in the absence of interactions, is symmetric at k = 0, and the coherent pumping does not make the momentum shift but excites a strong signal of k p and -k p (because of the scattering from the end of the chain). However, if we take the nonlinearity and the dissipation into account, the original energy dispersion is not stable as illustrated in Fig. 7(d). If the fluctuations are added in Eq. ( 8), we can see some obscure energy profile in Fig. 7(e). In the third column of Fig. 7, we plot ψ n (k, ω) -ψ (0,n) (k, ω) , which shows the energy dispersion of the fluctuations. The energy spectrum of the fluctuations under the periodic potential is still symmetric at k = 0 and the linewidth is almost the same. The nonlinearity and the dissipation make energy dispersion of the fluctuations have a non-zero momentum ground state.

THE TOPOLOGICAL INVARIANTS OF THE SYSTEM

In this section, we will calculate the eigenvalues and the winding number in detail. The OBC energies have two exceptional points (EPs) as is shown in Fig. 8(b). We need to calculate each winding of them separately. The Hamiltonian in the real space is k bk e ikn with âk and bk the annihilation operators for the momentum k, we can obtain the Hamiltonian with

H n = (∆ -iγ C ) n u * n u n -(∆ + iγ C ) n v * n v n +J n u * n+1 u n -v * n+1 v n + h.c. +gA 2 n e 2ikp u * n v n -e -2ikp v * n u n , (9)
H = ∆ k â † k âk -b † k bk + gA 2 k â † k bk-2kp -b † k âk+2kp + 2J k â † k âk -b † k bk cos k -iγ C . (10)
When the system is large enough and the periodic boundary condition is taken, we can get the 2 × 2 Hamiltonian with the vector (â k , bk-2kp ) as illustrated in the main text. To prove this, we can build a unit cell that contains l = 4π/k p inner sites

H(k) = H diag A 2 , • • • , A 2 e 2ikp(l-1) l×l diag -A 2 , • • • , -A 2 e -2ikp(l-1) l×l -H * , (11)
with H = (∆ -iγ C ) Îl×l + diag (J, -1) l×l + diag (J, 1) l×l . The matrix diag (J, ±1) are the upper and lower shift diagonal matrix with the matrix element J and the Bloch theory gives H 1,l = Je ik and H l,1 = Je -ik . The energy profile is corresponding with the 2 × 2 Hamiltonian in the main text. However, if we consider the finite size and assume the momentum k changes continuously, the effective Hamiltonian in the momentum space with vectors v = (â k , bk , âk+2kp , bk-2kp ) ignoring the energy shift γ C can be written as

H (k, k p ) =     ∆ + 2J cos(k) 0 0 gA 2 0 -∆ -2J cos(k) -gA 2 0 0 gA 2 ∆ + 2J cos(k + 2k p ) 0 -gA 2 0 0 -∆ -2J cos(k -2k p )     . (12)
The Hamiltonian in the momentum space has the symmetry of H(-k, -k p ) = H(k, k p ), S z H(k, k p )S z = H(k, k p ) † , and S x H(k, k p )S x = H(k, k p ) with S x,y,z = I 2×2 σ x,y,z with σ i is the Pauli matrix. The profile of the energy dispersion is highly symmetric and having two EPs if the Hamiltonian is non-Hermitian with A is a purely real number. By diagonalization, the four eigenenergies can be obtained

E 1,2 (k) = ± f (k) -2J sin(k p ) sin(k -k p ), (13)
with

f (k) = -g 2 A 4 + ∆ 2 + J 2 cos(2k p ) cos(2k -2k p ) + J 2 + 2∆ J cos(k) + 2J cos(k -2k p )(∆ + J cos(k)), (14)
and E 3,4 (k)=E 1,2 (k + 2k p ). Meanwhile, the related wavefunctions are

|ψ 1(2) =      - ∆ -2J cos kp cos(k-kp)± √ f (k) √ 2A 2 0 0 1/ √ 2      , |ψ 3(4) =      0 - ∆ -2J cos kp cos(k+kp)± √ f (k+2kp) √ 2A 2 1/ √ 2 0      , (15)
where, the function f (k) determine all the topological properties of the system. The real energy will collapse at k = ±π/2 for E 1,2 (k) = E 3,4 (k), but the imaginary energy will have a jump across the EPs. To define each winding of the EPs, Eq.( 12) can be written into two blocks off-diagonal form

H(k) = U H(k)U -1 =     0 0 -∆ -2J cos(k) -gA 2 0 0 gA 2 ∆ + 2J cos(k + 2k p ) -∆ -2J cos(k -2k p ) -gA 2 0 0 gA 2 ∆ + 2J cos(k) 0 0     , (16)
with the unitary transition U =     0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0     , and the winding number around two EPs can be calculated by

W 1 = - dk 2πi ∂ k log Det -∆ -2J cos(k) -gA 2 gA 2 ∆ + 2J cos(k + 2k p ) -E R Î , (17)
W 2 = - dk 2πi ∂ k log Det -∆ -2J cos(k -2k p ) -gA 2 gA 2 ∆ + 2J cos(k) -E * R Î , (18)
where E R and E * R are two reference energies with open boundary condition(OBC) in real space. If we take the loop from the whole Bloch zone, the total winding number is zero, but for each winding, we have W 1,2 = ±1. E and -E * will give a positive winding number while the E * and -E give the negative one. As is marked in Figs.

CRITICAL POINTS OF THE TOPOLOGICAL TRANSITION

In this section, we will calculate the critical points of the topological phase transition. As is defined by Eqs. ( 17 The HN model [START_REF] Hatano | Localization transitions in non-hermitian quantum mechanics[END_REF][START_REF] Hatano | Vortex pinning and non-hermitian quantum mechanics[END_REF] can be easily theoretically realized by H = J 1 e ikx + J 2 e -ikx with J 1 and J 2 are asymmetric hoppings between two nearest sites and the energy spectrum can be obtained by: H = (J 1 + J 2 ) cos (kx) -i(J 2 -J 1 ) sin (kx) .

(

) 19
The skin modes will vanish for J 2 = J 1 and appear for J 2 = J 1 . Actually, if the system has the same energy dispersion like Eq. ( 19) where different directions of the momenta have gain and loss at the same time, the skin modes will appear even with the same hopping by inducing other terms into the Hamiltonian.

To realize the skin effect, the complex dispersion will be like

E k = C 1 sin(k + C 2 ), (20)
with C i must be complex numbers. The complex number of C i ensures the condensates have gain in one direction and loss in another direction forming skin modes. The derivative of Eq. ( 14) is

f (k) = -4J cos(k p ) sin(k -k p ) [∆ + 2J cos(k p ) cos(k -k p )] , (21)
with the zero points k = k p for the maximum and k = arcos(-∆ 2J cos kp ) + k p for the minimum. If Eq. ( 14) is a purely imaginary number in the whole Bloch zone, all eigenstates are localized. The critical points for the topological transition are J 1 = gA 2 + ∆ sec(k p )/2 and J 2 = (gA 2 -∆ ) sec(k p )/2 which is corresponding with the numerical results J 1 ≈ 0.46 in Fig. 8 (a) and (c) and J 2 ≈ 0.69 in Figs. 8(a) and (e).

The phase can de defined in two ways, one is the appearance of the skin modes and the other one is the winding of the EPs energies. In the non-trivial phase, 0 < |J| < (gA 2 +∆ ) sec(kp) 2 all eigenstates are localized and the PBC energies can encircle all the OBC energies including EPs as is shown in Fig. 8 (e) andf). At the critical points J = ±J 1 or J = ±J 2 , the PBC energy will cross at E R = 0, and the modes -iγ C arise as extended eigenstates that can't be encircled by the OBC energies.

Figure 1 .

 1 Figure 1. Scheme of polariton lattice (a) and the related Bougoliubov lattice (b) with un and vn representing the different Bogoliubov fluctuation amplitudes on each site.

Figure 2 .

 2 Figure 2. The time evolution of the absolute square of the wavefunction |ψn| 2 (a) described by Eq. (1) and the spatial profiles of Bogoliubov modes un and vn are plotted in (b) and (c), which are obtained by Eqs.(3)-(4). The parameters are used ∆ /γC =-0.2, J=0.4, A=1, and kp=π/6.

Figure 3 .

 3 Figure 3. Profiles of the real (a) and imaginary (b) eigenstates as a function of un sites and eigenenergies. (c) Numerical results of the energy spectrums solved by Eqs. (3)-(4) with continuous pumping and the dashed lines are the fitting curves of the energy dispersion. Parameters are the same as in Fig. 2.

  Figs 3(a) and (b) show the distribution of intensity of the real and imaginary parts of the fluctuation eigenstates in real space. All fluctuation eigenstates are local-

Figure 4 .

 4 Figure 4. The real (a) and the imaginary (b) energy spectrum obtained by Eqs.(3)-(4). Other parameters are the same as in Fig.2.

) 2 :

 2 Parts of eigenstates are localized and parts of OBC energies including EPs are encircled by the PBC energies. • Trivial phase J = 0 or J (A 2 -∆ ) sec(kp) 2 : Parts

Figure 5 .

 5 Figure 5. The real space energy profile with OBC ('circles') and the momentum space energy profile with PBC ('dots'). J=0.4 and other parameters are the same as in Fig. 2.

Figure 6 .

 6 Figure 6. The time evolution of the absolute square of the wavefunctions as a function of sites and time ψ (n,t) 2 (the first row) or energies log 10 ψ (n,ω) 2 /max( ψ (n,ω) 2 ) (the second row). The parameters are used ∆=-0.22meV, J=40µeV, kp=π/6, gA 2 =0.1meV, and γC =0.00meV, 0.05meV, 0.08meV, 0.10meV for different columns.

Figure 7 .

 7 Figure 7. Energy dispersions for the mean-field (the first column), the mean-field with fluctuations (the second column), and only the fluctuations (the third column) obtained from Eq.(8) under the coherent pumping. Parameters are used: ∆=-0.22meV, J=40µeV, kp=π/6, and without the nonlinearity the dissipation in the first row and gA 2 =0.1meV, γC =0.1meV in the second row.

  where u n and v n are different Bogoliubov modes. If we just consider the nearest neighbouring hopping and use the Fourier transformation u n → 1 √ N k âk e ikn and v n → 1 √ N

Figure 8 .

 8 Figure 8. (a) Phase transition diagram along with the change of the hopping energy J and (b)-(f) the energy profiles obtained by OBC 'dots' and PBC 'circles' . Parameters are used: ∆ =-20µeV, γC =0.1meV, gA 2 =0.1meV, kp=π/6, and J=40µeV, 46µeV, 60µeV, 69µeV, 80µeV for (b)-(f).

  8(b)-(f) with red circles, all EPs are occurring at the pure imaginary energies. The eigenenergies that go around the left EPs can give the winding 1 and -1 for the right EPs in Figs. 8(b)-(d). The total winding is always 0 for the vanishing of E and E * .

  )-[START_REF] Gao | Continuous transition between weak and ultrastrong coupling through exceptional points in carbon nanotube microcavity exciton-polaritons[END_REF], two EPs have opposite windings. Based on the winding number and the localization of the wavefunction, we can define three phases: topological, intermediate, and trivial phase.

Figure 9 .

 9 Figure 9. The wavefuntions of un (the first row) and vn (the second row) sites obtained by OBC. The parameters are used ∆ =-20µeV, gA 2 =0.1meV, γC =0.1meV, kp=π/6, and J=60µeV (the first column), 80µeV(the second column), 0.2meV(the third column).

2 < 2 ,

 22 (b). In the intermediate phase, (gA 2 +∆ ) sec(kp) |J| < (gA 2 -∆ ) sec(kp) 2 parts of eigenstates are localized (see Figs.9 (a) and (d)) and the PBC energies can encircle parts of the OBC energies concluding the EPs (see Fig. 8 (d)). In the trivial phase, J = 0 or J (gA 2 -∆ ) sec(kp) parts of (see Figs. 9 (b) and (e)) or no eigenstates (see Figs. 9 (c) and (f)) are localized and the PBC energies can encircle parts of the OBC energies except the EPs (see Figs. 8

Acknowledgements-This work was supported by the Singaporean Ministry of Education, via the Tier 2 Academic Research Fund project MOE2018-T2-2-068.

𝜓 (",$%&) * ) 𝜓 (",$%&) ) -𝜓 (",$) * ' 𝜓 (",$) ' -𝜓 (",$()