C Eyraud

H Saleh

J.-M Geffrin

Influence of the description of the scattering matrix on the permittivity reconstruction with a quantitative imaging procedure: polarization effects

This paper focuses on the role of polarization -and more specifically, the effect of its selection -in 3D quantitative imaging obtained from scattered field measurements. Although polarization is now commonly used in linear imaging procedures (when unknowns are linked by a linear relationship to the measured signal), the influence of polarization choice is generally ignored in non-linear imaging problems. In this paper, we propose a formulation to obtain the 3D permittivity map, by a non-linear imaging procedure, from the scattering matrix. This allows, from the same data set, to select the desired polarization case as input data for the imaging algorithm. We present a study of the influence of the input data polarization choice on the reconstructed permittivity map. This work shows that a suitable basis choice for the description of the scattering matrix and an appropriate selection of the element of this scattering matrix can greatly improve imaging results.

I. INTRODUCTION

Electromagnetic wave probing is a powerful tool to reach the physical features of a scene in a non destructive way and has a wide area of potential applications in several domains. Indeed, after interaction with the scene, the radiated field contains the characteristics of the targets (position, shape, electromagnetic properties). Considering real measurements, the field can not obviously be measured on an entire surface enclosing the targets, only some directions of the transmitter and of the receiver are possible and the data are thus spatially truncated in such problems. It is crucial to exploit any opportunity to have information on the scene, and an appropriate selection of polarization is one of them. The polarimetric information of the electromagnetic wave is of great interest for imaging purposes and has been used for several years in radar [1] for radar remote sensing applications [START_REF] Boener | Lectures notes POLSARPRO[END_REF], [START_REF] Tyo | [END_REF], and in light scattering [START_REF] Mishchenko | Scattering, Absorption, and Emission of Light by Small Particles[END_REF], [START_REF] Garcia | [END_REF] such as, for example, for optical biomedical applications [6]. In these last cases, the recorded signal often allows a direct access to an image of the scene. To extract as much information as possible, the recorded signal can be acquired with different polarizations for both the transmitter and the receiver. Indeed, the polarimetric contrast is an interesting way to detect targets [7], [8]. Different methods have thus been developed to optimize the polarimetric contrast between two regions of the observed scene or between the target and the embedding medium [9]. When the recorded signal does not directly provide an image of the scene and when an imaging procedure is required, polarization must be integrated into the imaging process. In the case where unknowns are linearly related, under certain assumptions, to the measured signal, polarization has been used for several years to improve images. This is the case, for example, with synthetic aperture radar (SAR) [START_REF] Zyl | Synthetic Aperture Radar Polarilmetry[END_REF], [START_REF] Ma | [END_REF] scattering tomography [12] or Optical Coherence Tomography [13]. These procedures provide "qualitative" maps (except in cases where the targets are weak scatterers), allowing the targets to be detected and their positions and shapes to be obtained.

When one want, in addition, a permittivity map of the imaged area, "quantitative" imaging procedures must be used. In such procedures, a significant difference from the processes described above is that the unknowns (3D permittivity map) are related to the measured data by a non-linear relationship. In these procedures, the influence of polarization is very rarely taken into account. When different polarizations are taken into account, they are generally considered to contain the same amount of information and the processing of each polarization state is performed in the same way [14], [15], [16], [17]. In this contribution, we are interested in quantitative imaging and we propose a formulation to perform this quantitative imaging from the scattering matrix (which can be written in different basis) in order to take polarization into account easily. For 3D objects that have a high degree of freedom, it is all the more necessary to make the most of all the available information and polarization is important information. We propose here a study of the influence of the polarization selection on the reconstructed permittivity map. To focus only on the influence of polarization, no a priori information (neither on the object sought nor on the noise disturbing the measurements) is incorporated into the inversion procedure and no regulation term is added to the cost function. Images are obtained from both calculations and measurements. Linear and circular polarizations are studied. For the different cases, we are interested in the consequence of choosing a given scattering element of the scattering matrix over another and in the choice of the basis for writing the scattering matrix. We show that appropriate choices can significantly improve the reconstructed map. Part II is devoted to the scattering matrix construction in the various basis and to the presentation of the different polarization cases. The imaging algorithm is explained in part III. The details of the experimental setup and the measurement of the fields are explained in part IV. The imaging results are presented in part V and a discussion on these results follows in part VI. Some concluding remarks are finally given in part VII.

II. SCATTERING MATRIX

In far field, for a non-absorbing medium, the complex amplitude scattering matrix [S B ] describes the relationship between the incident field vector [E i B ] and the scattered field vector

[E s B ] [4]. [E s B ] = [S B ] [E i B ] (1)
This complex amplitude scattering matrix depends on the polarization type (linear, circular, elliptical) and on the choice of the basis B for this polarization description. Since far-field conditions are assumed, the field vector has no radial component and is then totally determined by its components along the two basis vectors. In a multistatic configuration, this scattering matrix [S B ] is defined as follows for each direction of the incident wave vector k i and each direction of the scattered wave vector k s :

E s p E s q = S pp S qp S pq S qq E i p E i q (2)
In this paper, we focused on two linear polarization representations and one circular polarization (Figure 2).

A. Linear polarization

When a linear polarization is selected, two main conventions are commonly used to express the scattering matrix : spherical convention -which is widespread in the microwave community [18] -and Bohren Huffman convention -commonly used in the light scattering community - [START_REF] Mishchenko | Scattering, Absorption, and Emission of Light by Small Particles[END_REF].

Spherical convention

For a given transmitter/receiver pair, the scattering matrix can be constructed with a spherical convention (Sp), using the vectors (e i θ , e i φ ) as basis for the incident wave and the vectors (e s θ , e s φ ) as basis for the scattered wave (Figure 1). Therefore, the element S p,q of the scattering matrix (with p, q ∈ {φ, θ}) corresponds to a scattered field measured along the vector e s q for an incident field E i polarized along the vector e i p (Figure 2 (c) Circular.

Bohren and Huffman convention

In the Bohren and Huffman (BH) definition, the scattering elements are linked to the scattering plane, i.e, the plane containing both the incident wave vector k i and the scattered wave vector k s [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF] (see Figure 2(b)). In the far-field region, the scattered field components (E s and E s ⊥ ) are then linked to the incident ones (E i and E i ⊥ ) by the scattering matrix. E i (resp. E i ⊥ ) corresponds to the parallel (resp. perpendicular) incident field component to the scattering plane and E s (resp. E s ⊥ ) corresponds to the parallel (resp. perpendicular) scattered field component to the scattering plane. The element S p,q of the scattering matrix (with p, q ∈ {⊥, }) corresponds to a scattered field measured along the e s q vector for an incident field E i polarized along the e i p vector (Figure 2(b)). The scattered field E s ⊥ along e s ⊥ (resp. the scattered field E s along e s ) is also called in polarization S (resp. in polarization P ). The scattering matrix described in Bohren Huffman convention can be obtained by a simple basis change from the scattering matrix described in spherical convention and vice versa.

B. Circular polarization

The left-circular polarized incident (resp. scattered) wave E i l (resp. E s l ) and right-circular polarized wave E i r (resp. E s r ) can be obtained from a linear combination of the linear polarized incident wave in the Bohren Huffman convention. In the FIG. 3: Definition of the left and right circular polarization with the IEEE convention [START_REF] Toh | [END_REF].

far field region and with the convention of the Figure 3 (outgoing wavevector for emission and incoming wavevector for reception) [START_REF] Toh | [END_REF], [21], the scattering elements S p,q of the scattering matrix (with p, q ∈ {l, r}) can thus be expressed from the linear polarized elements in the Bohren Huffman conven-tion:

S rr = 1 2 [(S -S ⊥⊥ ) -j (S ⊥ + S ⊥ )] (3)
S rl = 1 2 [(S + S ⊥⊥ ) -j (S ⊥ -S ⊥ )] S lr = 1 2 [(S + S ⊥⊥ ) + j (S ⊥ -S ⊥ )] S ll = 1 2 [(S -S ⊥⊥ ) + j (S ⊥ + S ⊥ )] III. IMAGING PRINCIPLE
The complex 3D permittivity map of the investigation zone is extracted from the fields measured in the receiving area by solving a non-linear inverse problem. In this work we started from an algorithm already implemented for the different cases of polarizations expressed in a spherical linear basis [22], [23]. For the current work, our algorithm has been adapted for the case where the scattering matrix is expressed in other basis. To do this, the cost functional, which estimates the difference between the measured scattering matrix element S m B,p,q and the simulated one S c B,p,q for a given permittivity map ε B,p,q (n) r , both expressed in the basis B, are written as follows

F (n) B,p,q = 1 2Ns Nr Ns s=1 Nr r=1 |S c B,p,q (ε B,p,q (n) r ;rs,rr)-S m B ,p,q (rs,rr)| 2 (4)
where N s (resp. N r ) is the number of source (resp. receiver) positions and ε

B,p,q (n) r
is the (n th ) 3D reconstructed permittivity map with ε B,p,q r = ε B,p,q r + jε B,p,q r . During the iterative process, this cost functional is minimized by an iterative procedure by modifying the permittivity map. Indeed, the permittivity of the voxels is updated through a gradient calculation:

ε B,p,q (n+1) r (r) = ε B,p,q (n) r (r) + α (n) B,p,q d (n) B,p,q (5)
with α (n) B,p,q the optimal step, is determined by a classical line search to minimize the cost functional at each iteration. d

(n) B,p,q
is the descent direction derived from the gradients of the cost functional ∇ εr F B,p,q (n) which are calculated with the adjoint problem [24]. At each step of the inversion procedure, the scattered field is calculated for the permittivity map ε

B,p,q (n) r .
The calculated scattered field is computed for the entire scattering matrix for the linear polarization in the spherical basis with a homemade forward problem based on a volume integral formulation and resolved with a biconjugated gradient stabilized FFT method [25] based on a 1D FFT to improve the calculation speed and to reduce the memory requirement [26] (see [27] for more details on this algorithm). The construction of the scattering matrix in the other basis is then deduced from the elements of the scattering matrix in the spherical basis. Initially, the permittivity of the voxels into the investigation domain is chosen to be very close to one (because of the constraint of positivity, a very low constant value is added to this value in order to avoid the case where the gradients disappear). The final permittivity map is obtained when the cost functional reaches a minimum. Very weak a priori information is introduced on the permittivity map during the inversion procedure, i.e., we use only positivity constraints to ensure physically realistic values (ε B,p,q r ≥ 1 and ε B,p,q r ≥ 0) and, in order to concentrate only on the polarization effects, no regularization has been made.

IV. EXPERIMENTS

The measurements, used for the 3D reconstructions, were made in the anechoic chamber of the CCRM (Centre Commun de Ressources en Microonde) (Figure 4 (left)). A spherical exprimental setup is included in this anechoic chamber, a sketch of its positioning devices is presented in Figure 4 (right)). Thanks to the different possible azimuthal and elevation displacements of the antennas in the setup, it is possible to measure the scattered field strength over a large part of a sphere surrounding the target. The target under test is placed at the center of a sphere of around 1.7 m diameter on a vertical expanded polystyrene mast considered transparent to microwaves. Two horn antennas, which transmit alternately, are placed on a circular vertical arch and each can have a total angular range of 180 • . The two transmitting antennas are in e φ and e θ polarizations respectively. The polarization of each antenna with respect to the spherical setup is selected by the antenna rotation as they are linearly polarized. A receiving horn antenna can also move circularly around the target in the azimuthal plane to capture the scattered signal, with a total angular range of 260 • . The target can make a full rotation (360 • ) around the vertical axis. More details on this experimental setup can be found in [28], [29]. In this study, in order to acquire the measured scattered fields corresponding to the four cases of co-polarization and cross polarization, two measurements are made. A first measurement with the receiving antenna polarization fixed at e φ allows to achieve the φφ and θφ measurement cases, and a second measurement with the receiving antenna polarization fixed at e θ allows to achieve the φθ and θθ measurement cases. With the experimental setup, the signal is measured in both magnitude and phase which allows to obtain the complex value of the electric field. The scattered electric field E s p,q cannot be obtained in a unique step, but it is rather obtained from the complex subtraction of two different fields: the measured electric field with the target under test is at its place in the chamber (total field) and measured electric field with no target (incident field). After obtaining the scattered field, a drift-correction procedure is applied [30] to remove the drift errors that could happen due to the time delay between the total and incident field. At the end, the drift corrected scattered field is calibrated, using a measurement of a reference target, to refer the incident field to a magnitude equal to one and a null phase at the target's center.

FIG. 4: Experimental setup V. IMAGING RESULTS

A. Targets

We consider two targets, i.e., the first is a canonical object with a relatively low permittivity value. It is a small cube with permittivity equal to ε r = 1.45. Its side is equal to 20 mm (1.2 λ at 18 GHz). The results for the first target may, in some cases, be interpreted with approximations on the induced current (see subsection VI B). The second object has a more complicated shape, i.e., two dimensions are very small compared to the wavelength (< λ 3 at the work frequency). Unlike the previous object, with this second object, the multiple scattering effect plays an important role, that cannot be neglected in the interpretation. This second target is a circular helix, made by the Centre de Transfert de Technologies du Mans through additive manufacturing (Figure 5). It can be included in a sphere of diameter 70 mm (4.2 λ at 18 GHz). Its length is 45 mm and its diameter 40 mm. The wire diameter is 5 mm (< λ 3 at 18 GHz). The aspect ratio of the equivalent cylinder, defined as the ratio between radius and length, is equal to 0.135. This target was built in acrylate and its permittivity was measured with a method based on the bistatic far-field scattering patterns of a sphere made with the same material [31] and was found equal to ε r = (2.97 ± 0.02) + j (0.085 ± 0.001) at 18 GHz [START_REF] Saleh | Analogie microonde appliquée à l'étude de la diffraction par des arbres, par des particules atmosphérqiues et des micros-organismes[END_REF].

FIG. 5: Circular helix B. Reconstructed 3D permittivity maps

This paragraph is devoted to the presentation of the 3D reconstructed permittivity maps. All the inversions are performed at 18 GHz, with 99 source positions (φ s ∈ [30 • : 10

• : 130 • ], θ s ∈ [110 • : 40 • : 430 • ]) and 53 receivers (φ r = 90 • , θ r ∈ [50 • + θ s : 5 • : 310 • + θ s ]
). These angles refer to the spherical coordinates of our experimental setup (Figure 4). The investigation zone is a cubic zone with a side of 30 mm (6 λ 3 at 18 GHz) for the cube and a parallelepipedic area (49 × 95 × 82) mm 3 (83 λ 3 at 18 GHz) for the circular helix.

To compare the various reconstructions in a quantitative way, some criteria are introduced to compare the N voxels of the reconstructed map with the ones of the theoretical map. It can be noticed that misalignment errors in the measurements can introduce a bias in the value of these criteria (but misalignment errors are low, as shown in [18]). Five criteria are considered, plus one C B that mixes the different effects. The first one is the classical quadratic norm C L2 B . As we shall see, this criterion is often "good" even "very good" and is too weakly sensitive to differences between maps. Two other criteria are also used, C L2 B,with and C L2 B,empty , that are also quadratic norms, but calculated by considering only the voxels with material and only the empty voxels (taking as reference the theoretical map). The fourth criterion is a crosscorrelation between the reconstructed and theoretical ε r -1 maps, obtained for the correct geometrical offset d * between the two maps. C c B is the correlation between the two maps.

C cc εr-1,B (d * ) = N (ε true r (r + d * ) -1)(ε reco r (d * ) -1) N (ε true r (r + 0) -1)(ε true r (0) -1) (6)
C c B = N (ε reco r,B (r) -ε reco r,B )(ε true r (r) -ε true r ) N (ε reco r,B (r) -ε reco r,B ) 2 N (ε true r (r) -ε true r ) 2 (7)
C B = (1 -C L2 B ) × C cc εr-1,B × C c B (8)
It can be noticed that for these criteria (1

-C L2 B , 1 -C L2 B,with , 1 -C L2 B,empty , C cc εr-1,B , C c B , C B
) the best value is 1.

Cube reconstructions

For the cube, we perform the inversion from calculated scattered fields. These scattered fields were calculated with the same home-made forward software used for the inversion procedure, based on a volume integral formulation [27]. To avoid the "inverse crime", the discretization of the target zone is not similar in this forward problem and in the inversion procedure. As it is explained in the section VI, we present in this section only the co-polarization and the well-adapted polarization elements. In figures 6, 7 and 8, the permittivity maps in the three planes x = 0 mm, y = 0 mm and z = 0 mm are plotted in the various polarization cases. The criteria values are summarized in table I.

Circular helix reconstructions

For this target, the inversion was done both from measured scattered fields and from calculated scattered fields. The criteria values, summarized in the tables II and III, were calculated for the reconstructions from the two cases. We present only the reconstructed maps obtained from the measurements (see Figure 9 for the 3D views and Figures 10,[START_REF] Ma | [END_REF]

Comparison of the different polarization cases

As it can be seen in Tables II andIII, the maps reconstructed from measurements for the circular helix are close to those obtained from calculations in the different polarization cases. However, the maps reconstructed from experimental fields have higher discrepancy (in particular the θθ case). This is logical, especially since no information on the experimental noise is taken into account in the inversion procedure, i.e., no regularization is carried out in this inversion concerning this noise [22]. In tables I, II and III, we can see that the quadratic standard norm calculated on all voxels (1 -C L2 B ) is very good (val- ues close to one) in all the polarization cases for the two objects. It is important to note that this criterion is the most common. Nevertheless, it gives only a global information and is not very discriminating. Indeed, it mainly provides information on empty cells, especially when they are in majority. Then, to compare the results in details, other criteria have been used. One solution is to separate the cells with material (according to the theoretical map) from the empty cells and to calculate the quadratic norm on these two sets (1 -C L2 B,empty , 1 -C L2 B,with ). Considering that, 1 -C L2 B,empty , empty voxels are well reconstructed for the two objects in all polarization cases except in the ⊥⊥ case. Indeed, in this polarization case, there are artefacts on the reconstructed maps, which results in FIG. 9: 3D view of the reconstructed maps in the different polarization cases from measured scattered fields for the circular helix (the threshold is chosen at 1.1). The real object is superimposed in cyan.

a deterioration of this criterion. Considering C L2 B,with , the voxels that should contain material are better reconstructed for the φφ, θθ for the cube case. For the circular helix, the four cases φφ, θθ , lr and rl are almost equivalent with respect to this criterion and are better than the two polarizations in the BH basis (⊥⊥ and ). The cross correlation C cc εr-1,B drops (i) if artifacts exist on the reconstruction (voxels that should be empty but are not on the reconstructed map) or (ii) if a part of the target is not reconstructed or (iii) if the object's permittivity value is not well recovered. C cc εr-1,B is very low for the cases ⊥⊥ and for the two targets. For these two polarization cases, this is because many parts of the object are not reconstructed (see Figures 9101112) and because, specifically for the ⊥⊥ case, there are many artifacts (see Figures 678). For the circular helix, in φφ and θθ cases, this criterion has an average value. For φφ, there are artifacts on the reconstructed map and for θθ, the value of the reconstructed permittivity is too low. The two best cases are lr and rl, where the reconstituted value of the permittivity is well recovered and there are few artifacts. For the cube case, the best values of this criterion are for the φφ and θθ cases, the ones where the entire shape of the object is better reconstructed. Finally, the C c correlation, that takes also into account the empty voxels, goes in the same direction, i.e., the best reconstruction cases are the ones in circular polarization for the circular helix and φφ and θθ for the cube. The global criterion C B , which is a simple multiplication of the three cri- teria, confirms that.

VI. ANALYSIS OF THE POLARIZATION EFFECTS

This section is devoted to the analysis of differences in reconstructions from the various elements of the scattering matrix.

A. Singular Value Decomposition

To analyze the information content in each scattering matrix element [S B, p q ] of polarization (p q) in the basis B, a Singular Value Decomposition of this multistatic element has been made.

[S B, p q ] = [U B, p q ] [Λ B, p q ] [V B, p q ] (9)
with [Λ B, p q ], the matrix that contains the singular values, [U B, p q ] and [V B, p q ], the two matrix that contain the associated bases singular vectors for (p, q) ∈ [φ, θ, ⊥, , l, r]. All these quantities expressed in the basis B for a scattered electric field measured along e s q and for an incident electric field polarized along e i p .

Figures 13 and14 show, for the two objects, the sum of the spectrum in each polarization case considering the spherical basis, the BH basis and the circular basis. The cases that contain most information are the co-polarization elements (S φφ , S θθ , S ⊥⊥ , S ) and the well-adapted circular polarization elements (S lr , S rl ). To take a closer look at the spectra, the items with most information are detailed in Figures 15 and16.

The size of the signal space for all polarization cases appears to be the same. This is also in agreement with the Efficient Bandwidth approach [START_REF] Bucci | [END_REF]. Indeed, with such an approach, for a target of this size, illuminated at this frequency, the width of this spectrum should be, given our limited circular arc, equal to 16 for the cube case and 27 for the circular helix case.

As might be expected, the co-polarized elements (S φφ , S θθ , S θθ , S ⊥⊥ , S ) and the well-adapted circular polarized elements (S rl , S lr ) contain more information than the crosspolarized elements (S φθ , S θφ , S ⊥ , S ⊥ ) or than the notadapted circular polarized elements (S rr , S ll ). Differences between the co-polarized/well-adapted circular polarized elements (S φφ , S θθ , S θθ , S ⊥⊥ , S , S rl , S lr ) are small and appear to be relative to the target (below, an analysis for the two targets under studied). This multistatic matrix analysis does not allow to discriminate between these elements. In particular, it does not allow to choose one of the elements rather than another. For the cube case, the sum of the singular values of the elements S lr , S rl , S φφ and S ⊥⊥ is a bit higher than the elements S θθ and S ) (Figure 13). These elements correspond also to the ones having the higher values in the signal space (Figure 15). Within these elements, the differences are too small to FIG. 12: 2D cross section in the (y0z) plane at x = 0 mm of the reconstructed maps in the principal polarization cases from measured scattered fields for the circular helix.

B, p q for the circular helix: full spectra (left) and spectra corresponding to the signal space (right).

values of the S φφ elements are higher than the those of the other elements in the signal space (only the first three singular values are almost identical for all the polarization cases) (Figure 16). Then, the spectra of the elements S ⊥⊥ , S lr , S rl are almost similar and are higher than the spectra of the elements S and S θθ .

B. Polarization vector analysis

In this section, we analyze, in the case of linear polarization, the filter effect created by the scalar product between the transmitter and the receiver polarization directions. It should be noticed that in the case of the circular base, this effect does not exist. As it is explained in the following, we restrict our analysis to the cases where, multiple scattering effects are weak as a first approximation.

For a source position at r s and a receiver position at r r , the scattered field can be written as

E s (r s , r r ) = Ω G(r r , r ) χ(r ) E(r )dr (10)
with G the dyadic Green function, χ(r ) = k 2 (r ) -k 2 o , the contrast term (k o being the wavenumber in the free space) and E, the electric field vector. The scattered field measured by the receiver, polarized along e s B,q , is:

E s (r s , r r ).e s B,q = ( Ω G(r r , r ) χ(r ) E(r )dr ).e s B,q (11)
B Sp Sp Sp Sp BH BH BH BH p, q φφ φθ θφ θθ ⊥⊥ ⊥ ⊥ QB 0.84 0.30 0 0.55 1 0 0 0.46 The induced current vector in the target zone Ω, for a permittivity map ε r (r), obtained for an incident field vector E i B,p = E i e i B,p is:

∀r ∈ Ω, J B,p (r) = χ(r) E B,p (r) (12)
This analysis is restricted to scattering with small multiple scattering effects and we assume that the induced current J B,p = χ E B,p is collinear to the incident field

E s (r s , r r ).e s B,q ≈ [ Ω G(r r , r ) J B,p (r )dr e i B,p ].e s B,q (13)
In the far field conditions, i.e. when each receiver position r r in the receiver zone Γ, satisfies the following conditions:

∀r r ∈ Γ, ∀r ∈ Ω, k o r r 1, r r r et k o r 2 r r 1 (14)
G can be approximated by [START_REF] Mishchenko | Scattering, Absorption, and Emission of Light by Small Particles[END_REF]:

G(r r , r ) ≈ e jkorr 4πr r e -jkoe s B,q .r [e θ ⊗ e θ + e φ ⊗ e φ ] (15)
The the scattered field measured by the receiver can thus be written as follows:

E s (r s , r r ).e s B,q ≈ ( Ω J B,p (r )e -jkoe s B,q .r dr ) (16) g(r r ) [e θ ⊗ e θ + e φ ⊗ e φ ]e i B,p .e s B,q with g(r r ) = e jko rr 4πrr .

Under the previous assumptions, we were interested in the part of the field which is not filtered by the scalar product between the polarization of the transmitting antenna and the polarization of the receiving antenna. Note that this quantity does not depend on the object, it only depends on the geometry of the configuration, i.e., it changes depending on the position of the source and receiver. This provides information on the main directions examined in the target. The quantity Q B is constructed as follows:

Q B = 1 N s N r Ns Nr q B (φ s , θ s , θ s -θ r ) (17)
with q B (φ s , θ s , θ s -θ r ) = |[(e θ ⊗ e θ + e φ ⊗ e φ )e i B,p ].e s B,q | and its are provided in Table IV.

According to Table IV, the ⊥⊥ seems to be the best case: Q B is equal to one for this polarization case and none of the measurement points is filtered in this polarization. Then, the φφ case is an interesting case because Q B is high and all

: Q ≈ 0.5. With
case, only measurement points in the forward scattering zone are contributing. For the θθ, a few more measurement points contribute because there is no filtering effect along the φ s angle.

The analysis made in this section can be used only if the induced current vector J B, p can be considered as collinear with the incident field vector E i B, p . In the following, we use this analysis to interpret the imaging results in the cube case. First of all, to determine if it is in agreement with the assumption of this analysis, we calculated the quantity A B, p for the two linear basis (Table V

): Spherical basis, A B,p = A Sp,p = 1 N 1 N s s | J B,p .e i B,p | ||J B,p || (18)
BH basis,

A B,p = A BH,p = 1 N 1 N s N r s r | J B,p .e i B,p | ||J B,p || . ( 19
)
If this value is close to 1, the image obtained from an element with such an illumination can be interpreted with the following analysis. As it can be seen in Table V, for the cube case, the results from the elements S φφ , S φθ , S θφ and S θθ can be analyzed with the reasoning proposed above, but it is not possible to use it to interpret the results from the elements of the scattering matrix in the BH basis. With this analysis, the best case should be the S φφ element and then, the S θθ element (see table IV). It's consistent with the imaging results obtained from these elements (see Figures 6,7, 8 and table I).

C. Discussion

If one considers the information contained in the multistatic element of the scattering matrix (see subsection VI A), the copolarized/well-adapted elements seem roughly equivalent. As we have seen in section V, the imaging results from the different elements are not equivalent at all. As a first intuition one might think that the BH basis should be the best approach and in particular that the S ⊥⊥ should give the best results as the component of the scattered field is measured colinearly to the incident field vector (see subsection VI B). On the contrary, all the reconstructions from this element failed. This surprising result may be understandable. First, as it can be seen in subsection VI B, even for a small low-contrasted object (the cube), the induced current vectors created in the target are not collinear to the incident field. This can be understood keeping in mind that in the BH basis, the basis vectors are defined in relation to the (k i , k s ) plane. This implies that with this basis, the polarization selection for source and receiver also plays the role of a filter and its influence can be strong. Thus, for a given illumination, the incident field vector direction is different for each receiver. It implies that for each receiver the induced currents in the target are different. This means that the first step in the inverse algorithm, which is equivalent to a backpropagation process, does not provide a first good estimate of the induced current map. This result is potentially a bad starting point for the algorithm and implies that the algorithm can more easily converge towards a local minimum that is not the global minimum. On the contrary, the co-polarized elements in the spherical basis and the well-adapted elements in the circular basis allow to provide a good estimate of the induced current map at the first step as, for a given illumination, the induced current map is the same for all the receivers. Looking to the differences between the reconstructions with the co-polarized element in linear polarization in the spherical basis (S φφ , S θθ ) and the well-adapted elements in the circular polarization (S lr , S rl ), in the case of a relatively simple target, like the cube, all these elements allow to reconstruct the target with a correct value of the permittivity (Figures 678). It can be noticed that the reconstruction with the S θθ is anyways worst along the z-axis. This makes sense as, in this polarization case, the induced current vectors are nearly collinear to the incident field polarization and the choice of the receiver polarization (along e θ ) plays the role of a filter (see section VIB) and also because the polarization vectors of the illuminated field and of the scattered field have no component along the z-axis. Now considering the circular helix, with the SVD analysis of the multistatic element (see section VIA), the elements S φφ , S lr and S rl should be better than the S θθ , that is indeed the case (see . Based on this analysis, the element S φφ contains a bit more information than the elements S lr and S rl , nevertheless the reconstruction with the elements S lr and S rl are a bit better. This can be explain by the geometry of this object, particularly adapted to a circular polarization.

VII. CONCLUSION

In this study, we have adapted our quantitative imaging procedure from scattering matrix elements in various polarization cases. To focus only on the influence of the polarization, we chose to introduce no a priori information into the inversion procedure (on the searched object or on the noise disturbing the measurements) and no regularization term was added in the cost functional. The scattering matrix was measured for a multistatic configuration in the spherical basis and then, obtained in the other basis by a simple basis change. In this paper, we considered the scattering matrix in different kinds of polarization (linear or circular) and in the linear case, with the two usual basis choices. In each case, the various elements of the scattering matrix were taken as input of our imaging procedure. Even if the different polarization cases are generally considered as equivalent inputs in such non-linear inverse procedures, we have shown that the polarization has a great impact on the reconstructed permittivity maps. The first result is not surprising, i.e., for a given choice of polarization (linear or circular), in a given basis, the co-polarization elements or the well-adapted polarization elements are the best ones. This is in agreement with the analysis of the information content in the multistatic elements of the scattering matrix with a Singular Value Decomposition. Moreover, we have seen that the scattering matrix elements having a similar information content (considering the SVD) can lead to very different reconstructed permittivity maps and thus, a SVD is not sufficient to understand which elements are the optimal ones. We have also shown that a good choice of the basis to describe the polarization is important. For the linear polarization case, results obtained with scattering matrix expressed in the spherical basis are not equivalent to those obtained with Bohren and Huffman basis. This can be understood considering the induced current created inside the target. It can be also noticed that, when one is interested in relatively simple targets, like the cube, even if the imaging results are not similar, almost all the scattering elements in co-polarization (for the linear polarization) or well-adapted polarization (for the circular polarization) allow to reconstruct this object. Looking to more complex objects, as the circular helix which in addition of its shape presents a wire diameter smaller than a third of the wavelength, that is no longer true. For these object, the choice of the polarization (linear or circular), of the basis and of the element of the scattering matrix have an important consequence on the imaging results. In the future, possible combinations of the most relevant scattering matrix elements can help to improve the imaging results according to the studied problem.

FIG. 1 :

 1 FIG. 1: Definition of the spherical convention: coordinate (a) and the polarization basis (b).

  FIG. 2: Definition of the three basis: (a) Linear -Spherical convention, (b) Linear -Bohren and Huffman convention and (c) Circular.

FIG. 6 :

 6 FIG.6: 2D cross section in the (y0z) plane at x = 0 mm of the reconstructed maps in the principal polarization cases from calculated scattered fields for the cube case. The edges of the target are superimposed in cyan.

  and 12 for the permittivity maps).

FIG. 7 :

 7 FIG. 7: 2D cross section in the (x0z) plane at y = 0 mm of the reconstructed maps in the principal polarization cases from calculated scattered fields for the cube case. The edges of the target are superimposed in cyan.

FIG. 8 :

 8 FIG. 8: 2D cross section in the (x0y) plane at z = 0 mm of the reconstructed maps in the principal polarization cases from calculated scattered fields for the cube case. The edges of the target are superimposed in cyan.

FIG. 10 :

 10 FIG. 10: 2D cross section in the (x0y) plane at z = 0 mm of the reconstructed maps in the principal polarization cases from measured scattered fields for the circular helix.

FIG. 11 :

 11 FIG. 11: 2D cross section in the (x0z) plane at y = 0 mm of the reconstructed maps in the principal polarization cases from measured scattered fields for the circular helix.

FIG. 13 :FIG. 14 :FIG. 16 :

 131416 FIG. 13: Sum of all the singular values for the Cube (left: in spherical basis, middle: in BH basis and right: in circular basis) ( Nr n=1 λ (n) B p q ).

TABLE I :

 I Comparison of the different reconstruction cases from the calculated scattered field. The target is the cube at the frequency of 18 GHz. It can be noticed that these criteria are better when there are 1.

TABLE II :

 II Comparison of the different reconstruction cases from calculated scattered field. The target is the circular helix at the frequency of 18 GHz.

TABLE III :

 III Comparison of the different reconstruction cases from the measured scattered field. The target is the circular helix at the frequency of 18 GHz.

TABLE IV :

 IV Values of the quantity Q B for different polarization cases.

TABLE V :

 V Values of the quantity A B, p for the different illuminations of the two linear basis in the cube case. measurements contribute. The part of the measurement that has the lowest weight corresponds to the smallest values of the angle φ s . θθ and correspond to low values of Q B

VIII. ACKNOWLEDGMENTS

The authors acknowledge the opportunity provided by the Centre Commun de Ressources en Microonde to use its fully equipped anechoic chamber. The authors want also to thank the Centre de Transfert de Technologie du Mans for the realization of the circular helix.


Paola Goatin

Elena Rossi

Comparative study of macroscopic traffic flow models at road junctions

Keywords: 2010 Mathematics Subject Classification: 35L65, 90B20, 82B21 macroscopic traffic flow models on networks, multilane models, multi-population models, finite volume schemes

We qualitatively compare the solutions of a multilane model with those produced by the classical Lighthill-Whitham-Richards equation with suitable coupling conditions at simple road junctions. The numerical simulations are based on the Godunov and upwind schemes. Several tests illustrate the models' behaviour in different realistic situations.

Introduction

Starting form the well-known Lighthill-Whitham-Richards (LWR) model [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF], a variety of macroscopic traffic flow models based on hyperbolic balance laws have been proposed to capture traffic behaviour in different situations. In this paper, we focus on the description of traffic dynamics on road networks. More precisely, we investigate the role of coupling conditions at road junctions. To this aim, we compare the behaviour of the multilane junction model introduced in [START_REF] Goatin | A multilane macroscopic traffic flow model for simple networks[END_REF], and the classical LWR model at junctions [START_REF] Garavello | Traffic flow on networks[END_REF]Chapter 5].

The multilane junction model [START_REF] Goatin | A multilane macroscopic traffic flow model for simple networks[END_REF] allows to handle in detail various realistic cases of road junctions, with the major exception of the diverging ones, which requires some additional information on drivers' routing preferences upstream. To avoid cumbersome notation, we recap here the main features of the multilane model in the case of a junction consisting of 2 incoming roads and 1 outgoing road, the number of lanes in each road to be specified later. Let M = M 1 + M 2 be the number of lanes on the half line x < 0, M j being the number of lanes on the j-th incoming road, j = 1, 2. Similarly, let M r be the number of lanes on x > 0, corresponding to the outgoing road. Each lane has maximal possible density normalized to 1, but each road may have different maximal speed. The detailed model is described in Section 3 below. Clearly, the case of a 1-to-1 junction is included, but not the 1-to-2 junction. This specific situation will be addressed separately in Section 5. 1 Inria Sophia Antipolis -Méditerranée, Université Côte d'Azur, Inria, CNRS, LJAD, 2004 route des Lucioles -BP 93, 06902 Sophia Antipolis Cedex, France. E-mail: paola.goatin@inria.fr 2 Università degli Studi di Milano-Bicocca, Dipartimento di Matematica e Applicazioni, via R. Cozzi 55, 20126 Milano, Italy. E-mail: elena.rossi@unimib.it The author is a member of INdAM-GNAMPA (Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni).

In order to draw a comparison with the LWR model on networks, we sum the vehicle densities on the various lanes, and we compare the total density profile with the solution given by the LWR model with corresponding maximal density. Notice that, in this way, we are led to consider the LWR model at a 2-to-1 junction with maximal density M j on the j-th incoming road and M r on the outgoing road. This choice is driven by the fundamental principle of conservation of the number of vehicles across the junction. The same approach is used on diverging junctions.

The numerical tests described in Sections 4 and 5 point out similarities and differences between the multilane model and its LWR counterpart. Interestingly, in the case of a 1-to-1 junction, the multilane description captures a different behaviour than those commonly described by the classical LWR approach. Merging junctions display pretty the same dynamics with both approaches, while diverging junctions show a more complex behaviour, acting as FIFO (first-in-first-out) when one outgoing road is fully congested, and as non-FIFO otherwise. Further assessment on the validity of the models necessitate the comparison with suitable real data.

LWR on n-to-m junctions

We recall the LWR model on n to m junctions, that is on junctions with n incoming roads and m outgoing roads. More details can be found in [START_REF] Garavello | Models for vehicular traffic on networks[END_REF][START_REF] Garavello | Traffic flow on networks[END_REF].

A road network is modelled by a finite collection of edges and vertices, representing respectively unidirectional roads and junctions. On each edge I, and thus on each road, traffic flow is described through the Lighthill-Whitham-Richards (LWR) model: the density of vehicles ρ(t, x) satisfies

     ∂ t ρ + ∂ x f I (ρ) = 0, ρ(0, x) = ρ o (x), (2.1)
where f I is assumed to be strictly concave and of class C 2 , such that f (0) = f (R I ) = 0, R I being the maximal possible density on road I. In particular, we may choose

f I (u) = V I u 1 - u R I , (2.2)
where V I is the maximal speed on road I. Note that, with this particular choice of flux, f I attains its unique maximum at ϑ I = R I /2. At network junctions, the conservation of the number of vehicles has to be ensured and, moreover, the usual choice is to maximise the flux through the junction. We introduce therefore the demand and supply functions (see for instance [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF]) and express the flux at the junction through an ad hoc minimisation between the two:

D I (u) = f I (u) if u ≤ ϑ I , f I (ϑ I ) if u ≥ ϑ I ,
and

S I (u) = f I (ϑ I ) if u ≤ ϑ I , f I (u) if u ≥ ϑ I . (2.
3)

The demand and supply functions introduced above are displayed in Figure 1.

Notice that, as in [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF], the flux function f I , as well as the corresponding demand and supply functions, may depend also on time, if traffic regulations are applied.

u D I (u) R I 2 R I f I (u) u S I (u) R I 2 R I f I (u)
Figure 1: The demand (left) and supply (right) function for f I as in (2.2). In both pictures, the dashed line represents f I (u).

At each junction, the model includes also of the preferences of the drivers, which are prescribed and known a priori. These preferences describe the distribution of traffic from incoming to outgoing roads and are expressed as elements of a matrix A, called traffic distribution matrix : for a vertex J with n incoming edges I 1 , . . . , I n and m outgoing edges I n+1 , . . . , I n+m , the matrix A is given as follows A = (a j,i ) ∈ R m×n with 0 ≤ a j,i ≤ 1 for every i ∈ {1, . . . , n} and j ∈ {n + 1, . . . , n + m} and n+m j=n+1 a j,i = 1 for every i ∈ {1, . . . , n} .

In other words, each element a j,i of the matrix A gives the percentage of traffic from road I i going into road I j at the junction.

In view of the examples considered in this paper, we focus on the cases of one-to-one junctions (modeling speed limit or lane number changes), two-to-one (merging) junctions and one-to-two (diverging) junctions, see Figure 2. In the case of merging junctions, where the One-to-one Two-to-one One-to-two number of incoming roads is greater than the number of outgoing roads, it is necessary to introduce a priority parameter for each incoming road to single out a unique solution. In particular, for two-to-one junctions, the two incoming roads are characterised by the priority parameters P and (1 -P ), with P ∈ ]0, 1[.

1-to-1 junction

Consider a one-to-one junction and call I 1 , I 2 respectively the incoming and outgoing road. Denote by ρ i the density of vehicles in road I i , i = 1, 2. Observe that the LWR model (2.1) on this particular junction is equivalent to a Cauchy problem for a conservation law with discontinuous flux (see, for instance, [START_REF] Holden | Front tracking for hyperbolic conservation laws[END_REF]Chapter 8] for a brief introduction to the topic):

     ∂ t ρ + ∂ x f (x, ρ) = 0, ρ(0, x) = ρ o (x), (2.4) with v(x, u) = H(x) v r (u) + (1 -H(x)) v (u), (2.5
)

f (u) = u v (u) f r (u) = u v r (u), (2.6
)

f (x, u) = u v(x, u) = H(x) f r (u) + (1 -H(x)) f (u), (2.7)
H being the Heaviside function. The functions v and v r are strictly decreasing, positive and such that v d (R d ) = 0 for d = , r. One admissible choice is

v d (u) = V d 1 - u R d , (2.8)
for a suitable positive constant V d , leading to the flux function (2.2).

Concerning the initial datum ρ o , it satisfies the following constraints:

ρ o (x) ∈ [0, R ] for x ∈ ] -∞, 0[, ρ o (x) ∈ [0, R r ] for x ∈ ]0, +∞[.
Hyperbolic conservation laws with discontinuous flux such as (2.4) arise in the modelling of several phenomena, such as, for instance, two phase flow in heterogeneous media and traffic flow with rough road conditions. Among the rich literature on the subject, we refer the interested reader to [START_REF] Adimurthi | Godunov-type methods for conservation laws with a flux function discontinuous in space[END_REF][START_REF] Adimurthi | Optimal entropy solutions for conservation laws with discontinuous flux-functions[END_REF][START_REF] Diehl | Scalar conservation laws with discontinuous flux function. I. The viscous profile condition[END_REF][START_REF] Karlsen | L 1 stability for entropy solutions of nonlinear degenerate parabolic convection-diffusion equations with discontinuous coefficients[END_REF][START_REF] Towers | Convergence of a difference scheme for conservation laws with a discontinuous flux[END_REF] and the references therein. The main issue when dealing with conservation laws with discontinuous flux function is the lack of a unique solution, since the classical theory by Kružkov does not apply. Clearly, Kružkov entropy conditions are valid away from the point(s) of discontinuity of the flux, but they are not enough to provide the uniqueness of solution. With the aim of proving the well-posedness of the problem, various notions of solutions have been introduced in the literature, based of different admissibility conditions. We recall for example the minimal variation criterion introduced by Gimse and Risebro [START_REF] Gimse | Solution of the Cauchy problem for a conservation law with a discontinuous flux function[END_REF] and the Γ-condition described by Diehl [START_REF] Diehl | On scalar conservation laws with point source and discontinuous flux function[END_REF]. We refer to [START_REF] Adimurthi | Optimal entropy solutions for conservation laws with discontinuous flux-functions[END_REF] for a more complete overview.

In the present study, we will focus on the solutions given by the supply-demand flux maximizing criterion classically used in traffic flow modeling (satisfying the entropy criterion of [START_REF] Adimurthi | Conservation law with discontinuous flux[END_REF]) and by the vanishing viscosity limit [START_REF] Diehl | Scalar conservation laws with discontinuous flux function. I. The viscous profile condition[END_REF].

2-to-1 junction

In the case of merging junctions, a priority parameter has to be assigned to each incoming road, see [START_REF] Garavello | Traffic flow on networks[END_REF]Section 5.2.2]. Therefore, in the specific situation of a two-to-one junction, the two incoming roads have priority P and 1 -P respectively, with P ∈ ]0, 1[. The idea is that, in general, not all vehicles can pass at the junction: denoting by γ the amount of vehicles that can enter the outgoing road in the supply-limited case, P γ vehicles comes from the first incoming road, and (1 -P ) γ cars from the second incoming road.

Clearly, since there is only one outgoing road, the traffic distribution matrix A reduces to the two-dimensional vector [START_REF] Adimurthi | Conservation law with discontinuous flux[END_REF][START_REF] Adimurthi | Conservation law with discontinuous flux[END_REF].

The initial datum satisfies the following constraints: [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF], denoting by ρ i and γ i the density and the flux on the incoming roads (i = 1, 2) and by ρ r , γr the density and the flux on the outgoing one, for P ∈ ]0, 1[ the fluxes at the junctions are defined as

ρ o, 1 (x) ∈ [0, R 1 ] for x ∈ ] -∞, 0[, ρ o,r (x) ∈ [0, R r ] for x ∈ ]0, +∞[. ρ o, 2 (x) ∈ [0, R 2 ] for x ∈ ] -∞, 0[, Following
γ 1 (ρ 1 , ρ 2 , ρ r ) := min D 1 (ρ 1 ), max P S r (ρ r ), S r (ρ r ) -D 2 (ρ 2 ) , γ 2 (ρ 1 , ρ 2 , ρ r ) := min D 2 (ρ 2 ), max (1 -P ) S r (ρ r ), S r (ρ r ) -D 1 (ρ 1 ) , (2.9) γr (ρ 1 , ρ 2 , ρ r ) := γ 1 (ρ 1 , ρ 2 , ρ r ) + γ 2 (ρ 1 , ρ 2 , ρ r ).

1-to-2 junction

The one-to-two junction is a particular diverging junction, see [10, Section 5.2.1] for more details. (Since there is only one incoming road, there is no need to prescribe a right of way.) In this case, the traffic distribution matrix A is given by:

A =    α 1 -α    , with α ∈ ]0, 1[.
In [START_REF] Garavello | Traffic flow on networks[END_REF] a FIFO (first-in-first-out) rule is applied at the junction. In terms of the demand and supply functions (2.3), the FIFO rule amounts to the following: denoting by ρ the density on the incoming road, by ρ r 1 and ρ r 2 the densities on the outgoing roads, by γ , γ r 1 and γ r 2 the fluxes on the incoming and outgoing roads respectively, we have

γ (ρ , ρ r 1 , ρ r 2 ) := min D (ρ ), S r 1 (ρ r 1 ) α , S r 2 (ρ r 2 ) 1 -α , γ r 1 (ρ , ρ r 1 , ρ r 2 ) := α γ (ρ , ρ r 1 , ρ r 2 ), (2.10) γ r 2 (ρ , ρ r 1 , ρ r 2 ) := (1 -α) γ (ρ , ρ r 1 , ρ r 2 ).
As a consequence, if one of the outgoing roads is fully congested (i.e. min{S r 1 (ρ r 1 ), S r 2 (ρ r 2 )} = 0), the vehicles are stuck at the junction (γ = γ r 1 = γ r 2 = 0), even if some of them could proceed in the other outgoing road. In order to overcome this problem, non-FIFO models have been developed, see for instance [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF]. Non-FIFO model allows some flow through the junction even if one of the outgoing road is fully congested. Using the same notation as before, the non-FIFO rule reads as follows, in terms of the demand and supply functions (2.3):

γ r 1 (ρ , ρ r 1 , ρ r 2 ) := min α D (ρ ), S r 1 (ρ r 1 ) , γ r 2 (ρ , ρ r 1 , ρ r 2 ) := min (1 -α) D (ρ ), S r 2 (ρ r 2 ) , (2.11)
γ (ρ , ρ r 1 , ρ r 2 ) := γ r 1 (ρ , ρ r 1 , ρ r 2 ) + γ r 2 (ρ , ρ r 1 , ρ r 2 ).
In the case of a diverging junction, the initial datum satisfies the following constraints:

ρ o,r 1 (x) ∈ [0, R r 1 ] for x ∈ ]0, +∞[, ρ o, (x) ∈ [0, R ] for x ∈ ] -∞, 0[, ρ o,r 2 (x) ∈ [0, R r 2 ] for x ∈ ]0, +∞[. 3 A multilane model for n-to-1 junctions (n = 1, 2)
We recall the main features of the multilane junction model introduced in [START_REF] Goatin | A multilane macroscopic traffic flow model for simple networks[END_REF], based on the multilane traffic flow model proposed in [START_REF] Holden | Models for dense multilane vehicular traffic[END_REF]. The model provides a description of traffic on road networks with several lanes, allowing for lane changes and overtaking, as well as change in the speed laws and in the number of lanes along the road. In particular, M ⊂ N + represents the set of indexes of the active lanes on ] -∞, 0[, with cardinality M = |M | ≥ 1; on the other hand, M r ⊂ N + is the set of indexes of the active lanes on ]0, +∞[, with cardinality M r = |M r | ≥ 1. We let M ≥ max {M , M r }, its choice depending on the specific situation considered, and we assume, for technical needs, that there are M -M and M -M r additional phantom lanes on ] -∞, 0[ and ]0, +∞[ respectively. We add moreover a condition (see (3.8) below), which prevents vehicles from passing from the active to the fictive lanes. The model reads then as follows: for x ∈ R and t > 0, the vehicle density on lane j, ρ j = ρ j (t, x), solves the Cauchy problem

     ∂ t ρ j + ∂ x f j (x, ρ j ) = G j-1 (x, ρ j-1 , ρ j ) -G j (x, ρ j , ρ j+1 ), j = . . . , M, ρ j (0, x) = ρ o,j (x), j = . . . , M, (3.1)
with, for j = 1, . . . , M ,

v j (x, u) = H(x) v r,j (u) + (1 -H(x)) v ,j (u), (3.2)
f ,j (u) = u v ,j (u), f r,j (u) = u v r,j (u), (3.3)
f j (x, u) = u v j (x, u) = H(x) f r,j (u) + (1 -H(x)) f ,j (u), (3.4)
H being the Heaviside function. The velocities v d,j , for d = , r and j = 1, . . . , M , are strictly decreasing positive functions such that v d,j (1) = 0. We assume that each map f d,j (u) = u v d,j (u) admits a unique global maximum attained at u = ϑ j d . We set ρ o,j : R → [0, 1] for j = 1, . . . , M and

ρ o,j (x) = 0 for x ∈ ] -∞, 0[ and j / ∈ M , (3.5)
ρ o,j (x) = 1 for x ∈ ]0, +∞[ and j / ∈ M r . (3.6)
The source terms, which account for the flow rates across lanes, are defined as in [START_REF] Holden | Models for dense multilane vehicular traffic[END_REF]:

G d,j (ρ j , ρ j+1 ) = K v d,j+1 (ρ j+1 ) -v d,j (ρ j ) + ρ j -v d,j+1 (ρ j+1 ) -v d,j (ρ j ) -ρ j+1 = K v d,j+1 (ρ j+1 ) -v d,j (ρ j ) •      ρ j if v d,j+1 (ρ j+1 ) ≥ v d,j (ρ j ), ρ j+1 if v d,j+1 (ρ j+1 ) < v d,j (ρ j ), (3.7)
for d = , r and j = 1, . . . , M -1, where (a) + = max {a, 0} and a -= -min{a, 0} and K is a dimensional constant (1/m), which can be assumed equal to 1 after rescaling space and time.

To account for separate lanes, such as different roads or fictive lanes, we set

G d,j d (u, w) = 0 for some j d ∈ {1, . . . , M -1} , d = , r. (3.8)
The functions appearing in the source term of (3.1) are then defined as follows

G j (x, u, w) = H(x) G r,j (u, w) + (1 -H(x)) G ,j (u, w) for j = 1, . . . , M -1, (3.9)
G 0 (x, u, w) = G M (x, u, w) = 0. (3.10)
For the sake of brevity, we introduce the notation ρ = (ρ 1 , . . . , ρ M ), so that the initial data associated to problem (3.1)

-(3.5)-(3.6) read ρ(0, x) = ρ o (x).
For simplicity, and with slight abuse of notation, we consider ρ = ρ(t, x) for t > 0, x ∈ R. However, as shown in [START_REF] Goatin | A multilane macroscopic traffic flow model for simple networks[END_REF], by (3.5), (3.6) and (3.8), there holds ρ j (t, x) = 0 for all t > 0,

x ∈ ] -∞, 0[ and j / ∈ M , respectively ρ j (t, x) = 1 for all t > 0, x ∈ ]0, +∞[ and j / ∈ M r .

In order to compare the multilane model to the LWR model at junctions, we make the following choices. In all the numerical experiments, for the multilane model we choose

v I,j (u) = V I (1 -u) for j = 1, . . . , M,
with V I > 0. Therefore, the maximal speed is the same for all lanes before, respectively after, x = 0. This corresponds to the following velocity for the LWR model:

v I (u) = V I 1 - u M I ,
where M I is the cardinality of the active lanes. Through our numerical integrations, we show that the outcome may be different when not considering the number of lanes involved.

Comparison study

We compare now numerically the solution given by the multilane n-to-1 model with the corresponding LWR solutions. We first describe the numerical schemes used to compute approximate solutions of the considered models. Then we report on some numerical tests, which illustrates the behaviour of solutions in different cases.

Numerical schemes

We provide in this section the details on the numerical schemes used throughout the paper.

We introduce a uniform mesh in space, of width ∆x, and a time step ∆t, subject to a suitable CFL condition. For k ∈ Z set

x k = k + 1 2 ∆x, x k-1/2 = k ∆x,
where x k is the centre of the cell and x k±1/2 its interfaces. Notice that nonnegative integers denote the cells on the positive part of the x-axis. Set λ = ∆t/∆x. For the multilane model, we approximate the initial data as follows: for j = 1, . . . , M ρ 0 j,k = 1 ∆x

x k+1/2 x k-1/2 ρ o,j (x) dx , (4.1)
recalling that (3.5) and (3.6) hold. The approximated initial data for the LWR model is given consequently as the sum of the initial data on various lanes, depending on the configuration under consideration:

• 1-to-1 junction (Section 4.2): the initial data for the LWR model is merely given by the sum, thus

ρ 0 k = M j=1 ρ 0 j,k for k ≤ -1; ρ 0 k = Mr j=1 ρ 0 j,k
for k ≥ 0.

• 2-to-1 junction (Section 4.3): for the merging junction, there are two incoming and one outgoing roads:

for k ≤ -1 : ρ 0 i ,k = M i j=1 ρ 0 j,k for i = 1, 2;
for k ≥ 0 :

ρ 0 r,k = Mr j=1 ρ 0 j,k .
• 1-to-2 junction (Section 5): for the diverging junction, there are one incoming and two outgoing roads:

for k ≤ -1 : ρ 0 ,k = M j=1 ρ 0 j,k ; for k ≥ 0 : ρ 0 r i ,k = Mr i j=1 ρ 0 j,k for i = 1, 2.
The solution to the multilane model

ρ n+1/2 j,k = ρ n j,k -λ F j (x k+1/2 , ρ n j,k , ρ n j,k+1 ) -F j (x k-1/2 , ρ n j,k-1 , ρ n j,k ) , ρ n+1 j,k = ρ n+1/2 j,k + ∆t G j-1 (x k , ρ n+1/2 j-1,k , ρ n+1/2 j,k ) -∆t G j (x k , ρ n+1/2 j,k , ρ n+1/2 j+1,k ),
where

F j (x, u, w) = min D j (x, u), S j (x, w) if x = 0, min D ,j (u), S r,j (w) if x = 0, (4.2)
with The LWR model is numerically integrated through a Godunov type scheme, which provides the solution to the Cauchy problem (2.1) that maximises the flux through the junction:

D j (x, u) = H(x)D j,r (u) + (1 -H(x))D j, (u), S j (x, u) = H(x)S j,r (u) + (1 -H(x))S j, (u)
• 1-to-1 junction (Section 4.2):

ρ n+1 k = ρ n k -λ F (x k+1/2 , ρ n k , ρ n k+1 ) -F (x k-1/2 , ρ n k-1 , ρ n k ) , (4.3)
with F defined as in (4.2), omitting index j.

In Section 4.2, we apply also a different strategy for the numerical integration of the LWR model: we make use of an upwind scheme, which provides the solution to the Cauchy problem (2.1) coming from the vanishing viscosity approach. The numerical scheme reads as in (4.3), where the numerical flux is now chosen as F (x, u, w) = u v(x, w).

• 2-to-1 junction (Section 4.3): we follow [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF] so that, exploiting the notation introduced in (2.9), setting γn

i = γi (ρ n 1 ,-1 , ρ n 2 ,-1 , ρ n r,0 ) for i = 1 , 2 , r, we have for k < -1, = 1 , 2 : ρ n+1 ,k = ρ n ,k -λ F (ρ n ,k , ρ n ,k+1 ) -F (ρ n ,k-1 , ρ n ,k ) , for k > 1 : ρ n+1 r,k = ρ n r,k -λ F r (ρ n r,k , ρ n r,k+1 ) -F r (ρ n r,k-1 , ρ n r,k ) , for k = -1, = 1 , 2 : ρ n+1 ,-1 = ρ n ,-1 -λ γn -F (ρ n ,-2 , ρ n ,-1 ) ,
for k = 0 : ρ n+1 r,0 = ρ n r,0 -λ F r (ρ n r,0 , ρ n r,1 ) -γn r .

• 1-to-2 junction (Section 5): as already recalled in Section 2.3, we use two different schemes, satisfying different rules at the junction.

FIFO rule: we follow [START_REF] Garavello | Traffic flow on networks[END_REF], and exploit the notation already introduced in (2.10), with γ n i = γ i (ρ n ,-1 , ρ n r 1 ,0 , ρ n r 2 ,0 ) for i = , r 1 , r 2 . Non-FIFO rule: we follow [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF] and exploit the notation already introduced in (2.11), with

γ n i = γ i (ρ n ,-1 , ρ n r 1 ,0 , ρ n r 2 ,0 ) for i = , r 1 , r 2 .
For both rules, the scheme amounts to the following:

for k < -1 : ρ n+1 ,k = ρ n ,k -λ F (ρ n ,k , ρ n ,k+1 ) -F (ρ n ,k-1 , ρ n ,k ) , for k > 1, r = r 1 , r 2 : ρ n+1 r,k = ρ n r,k -λ F r (ρ n r,k , ρ n r,k+1 ) -F r (ρ n r,k-1 , ρ n r,k ) , for k = -1 : ρ n+1 ,-1 = ρ n ,-1 -λ γ n -F (ρ n ,-2 , ρ n ,-1 ) , for k = 0, r = r 1 , r 2 : ρ n+1 r,0 = ρ n r,0 -λ F r (ρ n r,0 , ρ n r,1 ) -γ n r ,
with the above choices of γ n i , i = , r 1 , r 2 .

1-to-1 junction: from 2 to 3 lanes

We consider the case of a junction with one incoming road with 2 lanes and one outgoing road with 3 lanes. Case 1. We choose V = 1.5, V r = 1. The initial datum for the multilane model is

ρ o,1 (x) = 0.6, ρ o,2 (x) = 0.4, ρ o,3 (x) = 0.5 * χ [0,+∞[ (x), (4.4)
while the initial datum for the LWR model is given by the sum of the above functions, thus

ρ o (x) = 1 * χ ]-∞,0[ (x) + 1.5 * χ [0,+∞[ (x),
corresponding to the critical densities of the flux functions f and f r respectively, see Figure 3, i.e. the points where the flux functions attain their unique global maximum.

u f (u) 1 2 3 f (u) f r (u) ǔ f (ǔ) = f r (ǔ) f (M /2) = f r (M r /2) U f (U )
Supply-demand Vanishing viscosity multilane Figure 4, left, displays the solutions of the considered models at time t = 1: in particular, it allows to compare the multilane model, namely the sum of the densities on the various lanes, to the LWR model. For the latter, there are many admissible solutions: indeed, depending on the physics of the problem, the solution may be different, see [START_REF] Adimurthi | Optimal entropy solutions for conservation laws with discontinuous flux-functions[END_REF] and references therein for a thorough discussion on the choice of the "right" solution and on the physical models behind. In particular, the solutions obtained through the Godunov type scheme and through the upwind scheme are different. The former, obtained via the Godunov type scheme and thus through an ad hoc minimisation between the demand and supply function, is the solution that maximises the flux at x = 0, i.e. at the point of discontinuity of the flux: the solution connects f (M l /2) and f r (M r /2), see also Figure 3. On the other hand, the intermediate value attained around x = 0 by the solution coming from the upwind scheme corresponds to the point ǔ in Figure 3: this point is such that f (ǔ) = f r (ǔ). When higher order effects are taken into account for choosing the "right" solution, the choice of the vanishing viscosity approximation leads to this particular solution. Concerning the solution of the multilane model, the low value of the sum of the densities right downstream x = 0 is due to the fact that, before x = 0, the third lane does not exist. Thus, even though vehicles pass from the second to the third lane, the solution in the third lane is given by a shock with positive speed. On the other hand, the queue right upstream x = 0 is caused by the congestion of lanes 1 and 2, since after x = 0 the maximal speed decreases. Moreover, notice that the solution of the multilane model maximises the flux on each lane at the junction, see also Figure 5, right.

We make a second comparison between the LWR and the multilane model, namely we take the latter in the form of the average of the densities on the various lanes. In this way, the maximal density for the LWR model equals 1 everywhere. The initial datum is therefore given by the average of (4.4):

ρ o (x) = 0.5, which corresponds to the critical density of the flux function on both incoming and outgoing lanes. Figure 4, right, displays the solutions at time t = 1. We notice how the solutions resemble each other on the incoming lane: the queue upstream x = 0 is well captured by all models. On the other hand, also in this case the behaviour of the multilane model on the outgoing lane strongly differs from the LWR model, for which in this case the supply-demand and the vanishing viscosity solutions coincide.

Recall that, by [START_REF] Goatin | A multilane macroscopic traffic flow model for simple networks[END_REF]Lemma 2.3], the multilane model (3.1)-(3.5)-(3.8) preserves the total number of vehicles over time. This property is clearly valid when considering the multilane model in the form of the sum of densities on the various lanes. However, this second type of comparison considering the average of the densities on the various lanes implies that the total number of vehicles is not conserved anymore through the junction: indeed, the total mass is divided by the number of lanes, which is equal to 2 upstream x = 0 and 3 downstream. Case 2. We choose now V = 1 and V r = 1.5 and the same initial data of Case 1. The flux functions for the LWR model are displayed in Figure 6. In this case, the solution obtained with the Godunov type scheme is identical to that obtained through the upwind method, see Figure 7. Notice that the stationary shock at x = 0 is slightly smoothed out when using the upwind scheme, due to numerical viscosity. The small slope right before x = 0 in the sum of the solutions to the multilane model is caused by the higher flux of vehicles in lane 2 through x = 0, since right after x = 0 the speed is higher and vehicles can pass also in lane 3.

u f (u) 1 f (u) f r (u) 0.5 u f (u) 1 f (u) f r (u)
u f (u) 1 2 3 f (u) f r (u) ũ f (M /2) = f r (ũ)

2-to-1 junction: from 2+2 to 2 lanes

We consider the case of a junction consisting of two incoming roads with 2 lanes each and one outgoing road with 2 lanes, that is problem (3.1)-(3.6)-(3.8), with M = {1, 2, 3, 4}, M r = {2, 3} and initial data

ρ o,1 (x) = 0.6 χ ]-∞,0] (x) + 1 χ ]0,+∞[ (x), ρ o,2 (x) = 0.4, ρ o,3 (x) = 0.6χ ]-∞,0] (x) + 0.8 χ ]0,+∞[ (x), ρ o,4 (x) = 0.4 χ ]-∞,0] (x) + 1 χ ]0,+∞[ (x), (4.5)
with the additional assumption that there is no flow of vehicles between the second and the third lane on ] -∞, 0[, i.e. S ,2 (u, w) = 0 (we also impose S r,1 (u, w) = S r,3 (u, w) = 0). The situation under consideration looks as follows: Concerning the LWR model, we denote by 1 and 2 the two incoming roads, and by r the outgoing one. Then, the initial data for problem (2.1) corresponding to (4.5) read as follows:

ρ o, 1 (x) = 1, ρ o, 2 (x) = 1, ρ o,r (x) = 1.2, (4.6)
since road 1 is given by lanes 1 and 2 for x ∈ ] -∞, 0[, road 2 is given by lanes 3 and 4 for

x ∈ ] -∞, 0[ and road r is given by lanes 2 and 3 for x ∈ ]0, +∞[. Moreover, we prescribe the priority P = 1 -P = 1/2 on each incoming road. Observe that the maximal density on each road is equal to 2. Figure 8 displays the solution on each road at time t = 1: in each picture we can see both the LWR model and the multilane model, in the form of the sum of the densities on the involved lanes. Notice the asymmetric behaviour of the multilane model: this is due to the asymmetry of the initial data (4.5) with respect to the lanes that are left after the junction. Indeed, the initial datum in lane 2 has a lower value than in lane 1, while the opposite is true for lane 3 with respect to lane 4. In the former case, this implies that more vehicles have to pass form lane 1 to lane 2 in order to go through the junction, causing the formation of a longer queue, which become visible when summing the two densities, see Figure 8. Figure 9 provides a detail of the multilane model at the same time t = 1, allowing a better insight on the queue formation.

A multilane model for a diverge

In this section, we aim to compare the LWR model on a one-to-two (diverging) junction with a multilane model different from the one considered in Section 3. Indeed, we need to account for the drivers' routing preferences upstream, a feature that is not included in the multilane model (3.1)-(3.7).

The specific multilane junction looks as follows: there are one incoming road with two lanes and two outgoing roads, with one lane each. lane 2

x = 0 lane 1

The core idea is to consider a multi-population model [START_REF] Cristiani | A destination-preserving model for simulating Wardrop equilibria in traffic flow on networks[END_REF][START_REF] Garavello | Source-destination flow on a road network[END_REF], where each population is identified by its desired final destination. In the particular case under consideration, we have two populations, which coexist upstream x = 0. We assume that, in each lane of the incoming road, a fraction α of the drivers belong to the first population, i.e. wants to go in lane 1 downstream x = 0, while the rest is directed towards lane 2. Overtaking is possible only upstream x = 0, and, differently from the multilane model considered in Section 3, only vehicles in lane 1 (respectively lane 2) targeting lane 2 (respectively lane 1) change lane. The rate of vehicles changing lane depends on the total density in the target lane: the more the target lane is crowded, the less vehicles are allowed to move into that lane. Moreover, the desire to go in the target lane increases as drivers approach the junction (x = 0). Since we focus on the dynamics at the junction, we do not consider the classical overtaking dynamics, which can occur before x = 0. This could be added in the source term with minor changes in the overall approach. Our multilane multi-population model is the following: given the initial data ρ o,1 and ρ o,2 on R, the initial data on x < 0 read

ρ 1 o,1 (x) = α ρ o,1 (x) χ ]-∞,0] (x), ρ 2 o,1 (x) = (1 -α) ρ o,1 (x) χ ]-∞,0] (x), ρ 1 o,2 (x) = α ρ o,2 (x) χ ]-∞,0] (x), ρ 2 o,2 (x) = (1 -α) ρ o,2 (x) χ ]-∞,0] (x), (5.1)
where ρ i o,j , for i, j ∈ {1, 2}, denotes the amount of drivers in lane j targeting lane i. Setting ρ 1 = ρ 1 1 + ρ 2 1 and ρ 2 = ρ 1 2 + ρ 2 2 , the densities ρ i j , i, j = 1, 2, satisfy the following system:

                           ∂ t ρ 1 1 + ∂ x ρ 1 1 v 1 (ρ 1 ) = K ρ 1 2 max {1 -ρ 1 , 0} |x| χ ]-∞,0[ (x)
∂ t ρ 2 1 + ∂ x ρ 2 1 v 1 (ρ 1 ) = -K ρ 2 1 max {1 -ρ 2 , 0} |x| χ ]-∞,0[ (x)
∂ t ρ 1 2 + ∂ x ρ 1 2 v 2 (ρ 2 ) = -K ρ 1 2 max {1 -ρ 1 , 0} |x| χ ]-∞,0[ (x)
∂ t ρ 2 2 + ∂ x ρ 2 2 v 2 (ρ 2 ) = K ρ 2 1 max {1 -ρ 2 , 0} |x| χ ]-∞,0[ (x), (5.2)
where K is a dimensional constant (m/s), which can be assumed equal to 1 after rescaling space and time. In the source terms, the factor |x| -1 represents the increase of the urgency to change lane as the junction approaches.

We observe that a close model which describes traffic on a multilane highway under the hypotheses that traffic is neither perfectly FIFO nor perfectly non-FIFO has been introduced in [START_REF] Schnetzler | A multilane junction model[END_REF]. Here, we compare the multilane multi-population model (5.2) to the LWR model for a diverging junction, both in the case of a FIFO rule and of a non-FIFO rule at the junction, see Section 2.3. The incoming roads are paired together, and denoted by the subscript , while so that

ρ 1,n+1 1,k = ρ 1,n+1/2 1,k + ∆t G(x k , ρ n+1/2 1,k , ρ 1,n+1/2 2,k ), ρ 2,n+1 1,k = ρ 2,n+1/2 1,k -∆t G(x k , ρ n+1/2 2,k , ρ 2,n+1/2 1,k ), ρ 1,n+1 2,k = ρ 1,n+1/2 2,k -∆t G(x k , ρ n+1/2 1,k , ρ 1,n+1/2 2,k ), ρ 2,n+1 2,k = ρ 2,n+1/2 2,k + ∆t G(x k , ρ n+1/2 2,k , ρ 2,n+1/2 1,k
).

Example 1: fully congested outgoing road

The first example we take into account considers a fully congested outgoing road. We choose α = 0.4, the maximal speed on the incoming road is V = 1.5 and on both outgoing roads it is V r = 2. The initial data for the multilane multi-population model (5. Figure 10 displays the solution on each road at time t = 0.5: in each picture we see the LWR model with both FIFO and non-FIFO rule and the multilane multi-population model, in the form of the sum of the densities for x < 0. Lane 2 is initially fully congested: no vehicles can enter, and the solution on that lane is clearly the same with all the three models. In this situation, the solution of the multilane multi-population model looks very close to that of the LWR model with FIFO rule. Indeed, non-FIFO rule allows vehicles targeting lane 1 to overcome the junction: they are not queuing waiting for their turn to pass the junction.

Example 2

In this second example, none of the outgoing lane is fully congested. We keep the same parameters as before, thus α = 0.4, V = 1.5, V r = 2. The initial data for the multilane multi-population model (5. Figure 11 displays the solution on each road at time t = 0.5: in each picture we see the LWR model with both FIFO and non-FIFO rule and the multilane multi-population model, in the form of the sum of the densities for x < 0. Even though lane 2 is not fully congested downstream x = 0, all the three models produce the same solution there. However, on the outgoing road corresponding to lane 1 the solution to the multilane multi-population model (5.2) is very close to the solution to the LWR model (2.1) with non-FIFO rule. In the incoming roads, the similarity is reduced, but nevertheless the solution to the LWR model with FIFO rule looks much different than the others, since it reaches a greater value upstream x = 0.

Figure 2 :

 2 Figure 2: The types of junction considered in this work for the LWR model.

  (3.1)-(3.5)-(3.6)-(3.8) is obtained through a Godunov type scheme, with fractional step to take into account the source terms, see [13, Algorithm 2.1]:

  , D j,d and S j,d , d = , r, defined as in (2.3).

  For the multilane model, this corresponds to problem (3.1)-(3.5)-(3.8) with M = {1, 2}, M r = {1, 2, 3} and S ,2 (u, w) = 0. For the LWR model, this is a 1-to-1 junction with maximal density on ] -∞, 0[ equal to M = 2, while on ]0, +∞[ it is M r = 3.

Figure 3 :

 3 Figure 3: Flux functions f and f r related to the LWR model, for M = 2, M r = 3, V = 1.5 and V r = 1. The point ǔ is such that f (ǔ) = f r (ǔ). The value U corresponds to the left trace at x = 0 of the sum of the solutions of the multilane model on the two incoming lanes, and f (U ) is the corresponding value of the flux.

Figure 4 :

 4 Figure 4: The dashed blue line corresponds to the multilane model (3.1)-(3.5)-(3.8): on the left, it is the sum of its solutions, on the right it is the average. The dash-dotted orange line corresponds to the solution to the LWR model (2.1) obtained via the Godunov type scheme; the dotted green line is the solution to the LWR model (2.1) obtained through the upwind scheme. Here: V = 1.5, V r = 1, initial datum (4.4).

Figure 5 :

 5 Figure 5: Left: Flux functions f and f r related to the LWR model when comparing it to the multilane model in the form of the average of the densities on the various lanes: V = 1.5 and V r = 1, in both cases the maximal density is 1. The orange line represents the solution to the Riemann problem with initial datum ρ o (x) = 0.5. Right: flux functions related to the multilane model for the incoming (f ) and the outgoing (f r ) lanes. The magenta line represents the solution on lane 1, with initial datum ρ o,1 = 0.6; the dotted blue line corresponds to the solution on lane 2, with initial datum ρ o,2 = 0.4.

Figure 6 :

 6 Figure 6: Flux functions f and f r related to the LWR model, for M = 2, M r = 3, V = 1 and V r = 1.5. The point ũ is such that f r (ũ) = f (1).

Figure 7 :

 7 Figure 7: The dashed blue line is the sum of the solutions to the multilane model (3.1)-(3.5)-(3.8); the dash-dotted orange line corresponds to the solution to the LWR model (2.1) obtained via the Godunov type scheme; the dotted green line is the solution to the LWR model (2.1) obtained through the upwind scheme. Here: V = 1, V r = 1.5, initial datum (4.4).

Figure 8 :

 8 Figure 8: In each picture, the dashed blue line is the sum of the solutions to the multilane model (3.1)-(3.5)-(3.8): from left to right, lanes 1 and 2; lanes 3 and 4; lanes 2 and 3. The dash-dotted orange line corresponds to the solution to the LWR model (2.1), obtained through a Godunov type scheme, with priorities 1/2, 1/2 : from left to right, incoming roads a, b, outgoing road c. The initial data are given in (4.5) and (4.6) respectively. In each lane we set V = 1.5.

Figure 9 :

 9 Figure 9: Solution to the multilane model (3.1)-(3.5)-(3.8) at time t = 1, with initial data (4.5) and V = 1.5 on each lane.

  2) areρ o,1 (x) = 0.6 χ ]-∞,0] (x) + 0.4 χ ]0,+∞[ (x), ρ o,2 (x) = 0.7 χ ]-∞,0] (x) + 1 χ ]0,+∞[ (x),(5.4)and through (5.3) we recover the initial data for the corresponding LWR model (2.1).

Figure 10 :

 10 Figure 10: In each picture, the dashed blue line corresponds to the solutions to the multilane multi-population model (5.2), and in particular to their sum on x < 0. The dotted orange line corresponds to the solution of the LWR model (2.1) with non-FIFO rule, while the dashdotted green line is the solution to the LWR model (2.1) with FIFO rule. The initial data is (5.4), while V = 1.5, V r = 2 and α = 0.4.

  from(5.4) only in the density on the second road downstream x = 0. Through (5.3) we recover the initial data for the LWR model (2.1).

Figure 11 :

 11 Figure 11: In each picture, the dashed blue line corresponds to the solutions to the multilane multi-population model (5.2), and in particular to their sum on x < 0. The dotted orange line corresponds to the solution of the LWR model (2.1) with non-FIFO rule, while the dashdotted green line is the solution to the LWR model (2.1) with FIFO rule. The initial data is (5.5), while V = 1.5, V r = 2 and α = 0.4.

the two distinct outgoing roads are denoted by subscripts r 1 and r 2 respectively. The initial data for problem (2.1) is related to the initial data ρ o,1 , ρ o,2 for the multilane multi-population model as follows:

.

Numerical scheme

We detail here the numerical scheme exploited for the integration of the multilane multipopulation model (5.2). The scheme is inspired by that presented in [START_REF] Festa | Modeling the impact of on-line navigation devices in traffic flows[END_REF][START_REF] Samaranayake | Discretetime system optimal dynamic traffic assignment (SO-DTA) with partial control for physical queuing networks[END_REF], and the source terms are treated through fractional step.

The space and time mesh are defined as in Section 4.1. The initial data on each lane, ρ o,1 and ρ o,2 , are approximated by ρ 0 1,k and ρ 0 2,k , for k ∈ Z, as in (4.1), then distinguished on the negative part of the x-axis according to their target lane: for k ≤ -1

The solution to the multilane multi-population model (5.2) is obtained through a Godunov type scheme, with fractional step to account for the source terms. In particular, set ρ n j,k = ρ 1,n j,k + ρ 2,n j,k and introduce the quantity γ n j = min Dj (ρ j,n j,-1 , ρ i,n j,-1 ), S j (ρ n j,0 ) , for j = 1, 2, i = 3 -j, with S j as in (2.3) and Dj defined as let ϑ j (w) the point of maximum of the function u → uv j (u + w), then

Then, for i, j = 1, 2, we set


Isaure De Kernier

Stéphanie Bressieux

Nelly Rongeat

Anaïs Ali-Cherif

Sophie Morales

Serge Monneret

Pierre Blandin

Statistical study of blood cells population by very wide-field bimodal phase/ fluorescence imaging

Keywords: optics, statistics, fluorescence, phase, in-line holography, microscopy

Very wide-field of view imaging can provide statistical data on large cell populations in a single acquisition. It can be used as an alternative to standard fluorescence flow cytometry methods for applications in hematology. We will describe a multimodal imaging system combining brightfield, phase and fluorescence contrasts. Its greater simplicity and lower cost compared to flow cytometry make it suitable for Point-Of-Care applications. The system's resolution was characterized on calibrated beads and targets. We illustrate the potential of the single-shot imaging approach in hematology by showing classification of white blood cells morphologies.

INTRODUCTION

The study of cellular populations requires acquiring information on a significant number of cells to ensure statistically reliable results. In imaging flow cytometry, the statistics comes from the flow of objects in front of a high magnification imaging system. Such a solution typically ensures a throughput of 500 -1000 cells per second [START_REF] Zuba-Surma | The ImageStream System: a key step to a new era in imaging[END_REF] . This is much lower than the statistics achieved using flow cytometry, that tops 5 000 cells per seconds [START_REF] Adan | Flow cytometry: basic principles and applications[END_REF] . This affects the period of time over which the acquisition is done. Recent research work have examined the method and applications of very-wide field imaging. If very large objects such as embryos [START_REF] Mcconnell | A novel optical microscope for imaging large embryos and tissue volumes with sub-cellular resolution throughout[END_REF] or large numbers of small objects [START_REF] Hatch | 1-Million droplet array with wide-field fluorescence imaging for digital PCR[END_REF] such as blood cells need to be studied, then acquiring a large Field Of View (FOV) in a single-shot minimizes the total number of images required. The figure below illustrates that the FOV of a standard microscope with a 10x objective is 100 times smaller than that of a 1x objective. In this work we aim at developing a high throughput single-shot imaging system. Such a system enables statistical studies of cell populations without requiring any flow of cells as in imaging flow cytometry or scanning as in laser-scanning cytometry [START_REF] Darzynkiewicz | Laser-Scanning Cytometry: A New Instrumentation with Many Applications[END_REF] . We demonstrate that we are able to acquire transmission, phase and fluorescence images of a 30 mm 2 FOV. Furthermore, we illustrate the potential of our system for statistical analysis of blood cell populations. The morphologies of labeled nuclei of purified White Blood Cells (WBC) sub-types are considered in fluorescence. This proof of concept experiment suggests that statistical analysis of blood cell populations is possible.

DESIGN OF AN IN-LINE IMAGING SYSTEM

Combining contrasts in an imaging systems requires registration of the images from each modality. We chose to design a system with a single imaging pathway in order to obtain hardware registrated images, such that no additional image postprocessing is necessary to merge information. The common imaging path set up with a single objective (referred to as Macro Lens in the figure) and a CMOS sensor is illustrated in Figure 2. Advantageously this solution is cost-effective. The aim of the system is to perform single-shot statistical imaging of micrometric objects. This implies that both FOV and resolution should be optimized. To obtain a large FOV, we chose a SRL camera lens from macrophotography series of Canon Inc. (Canon EF 100 mm f/2.8 Macro USM, Canon Inc., USA). The numerical aperture of this lens is measured to be 0.113 ± 0.003. The PSF diameter is theoretically determined to be 2.81 µm at 520 nm. Hence, according to the Shannon-Nyquist criteria, the optimal pixel pitch is 1.4 µm. However, as the development of sensors has been driven by low-cost high resolution cameras for phones, small pixel pitch down to 1.1 µm are available but at the detriment of the chip size. For our system, we make the compromise of choosing a monochromatic CMOS 11 Mpx 6.4 by 4.6 mm 2 sensor with a 1.67 µm pixel pitch (UI-1492LE-M, IDS GmbH, Germany) resulting in a 29.4 mm 2 FOV. This solution is found to be a good resolution-FOV compromise. However, a larger chip with smaller pixels would enhance the performances of the set-up.

The chosen macro lens is designed to achieve a 1:1 magnification ratio. In this configuration, the object plane is 145 mm away from the front lens of the objective. The sensor is mounted on a manual micrometric z-axis translation stage (SM1Z, Thorlabs Inc., USA) with a 25 mm range to enable precise positioning of the image plane.

Brightfield images can be acquired by illuminating the object in an in-line configuration, as illustrated in Figure 2. The illumination module comprises a 520 nm low-cost LED (XLamp® XM-L™ Color LED, Cree Inc., USA) coupled into an optical fiber (FG200UEA, Thorlabs Inc., USA) and shaped into a convergent beam by an additional lens (LAT075, Thorlabs Inc., USA).

Several methods enable to obtain a phase contrast from out-of-focus brightfield images. With the specifications mentioned, our system's brightfield illumination is coherent enough to enable in-line holography [START_REF] Singh | Label-free, high-throughput holographic screening and enumeration of tumor cells in blood[END_REF] acquisitions. State of the art algorithms [START_REF] Hervé | Method for observing a sample, by calculation of a complex image[END_REF] are used to reconstruct an amplitude and a phase map from the acquired hologram.

Both brightfield and phase imaging allows to obtain morphological information on the studied object. However, this is not sufficient to discriminate sub-types of cells whose morphology is similar. For this reason, it is interesting to introduce a fluorescence contrast that is more specific. In this work, we chose to work with the Thiazole Orange fluorescence dye [START_REF] Lee | Thiazole orange: A new dye for reticulocyte analysis[END_REF] . The excitation and emission specifications for this dye are given in Figure 3. A laser (LBX-505, Oxxius, France) was combined to a bandpass excitation filter (FF01-504/12, Semrock, USA) and used as the excitation source. The excitation beam was tilted with an angle ~45° from the optical axis to avoid direct excitation to be incident on the sensor, as illustrated in Figure 2. This reduced the constraints on the emission filter. The emission filter (FF01-536/40-25) was added in the imaging path between the objective and the sensor to filter out the excitation wavelength. Because the brightfield illumination wavelength was chosen higher than the filter's cut-off wavelength, brightfield and holographic imaging could be performed without removing the fluorescence emission filter.

CHARACTERIZATION OF THE SYSTEM

The system combines a very wide field of view of 29.4 mm 2 with a micrometric resolution, as illustrated in Figure 4. To characterize the resolution, we imaged an amplitude 1951 USAF resolution test chart. Here resolution is defined as the width of the last group of bars whose contrast is over 10%. Under blue LED illumination, we measured a resolution of 2.76 ± 0.15 µm in the center of the FOV. Furthermore, the fluorescence resolution was experimentally determined. A solution of 1 µm GFP-fluorescent beads was evaporated on a standard microscope glass slide. In the acquired fluorescence image, individual beads resulted in Point Spread Functions (PSF). Resolution was defined as the FWHM of the experimental discrete data points of a line profile of a PSF. An average resolution of 3.72 ± 0.25 µm was obtained in the center of the FOV, as illustrated in Figure 5. These experimental evaluations of the multi-modal resolution indicate that resolved imaging of small biological objects should be possible.

STATISTICAL STUDY OF WHITE BLOOD CELL MORPHOLOGIES

The classification of the three sub-populations of White Blood Cells (WBC) is an important basic information, usually provided in routine blood tests. WBC can be separated into three distinct populations: monocytes, lymphocytes and granulocytes. The proportion of these different types of blood cells is used for medical diagnosis, such as Dengue fever [START_REF] Forcucci | All-plastic, miniature, digital fluorescence microscope for three part white blood cell differential measurements at the point of care[END_REF][START_REF] Buttarello | Automated Blood Cell Counts: State of the Art[END_REF] .

To demonstrate the imaging capacity of our system, we acquired bi-modal images of blood smears of WBC purified populations. We show that the fluorescence image can be used to discriminate the sub-types of WBC.

The study was achieved on three whole blood samples. For each 2mL whole blood sample, a specific sub-population of WBC was purified. Purification of lymphocyte was done using the RosetteSep Human total Lymphocyte enrichment kit (15223, StemCell). The purified solution was diluted in 500 µL of PBS + 2% FBS (07905, StemCell, USA). Purification of granulocytes was done using the EasySep Direct Human Pan-Granulocytes Isolation kit (19659, StemCell, USA). Finally, purification of monocytes was done in two steps. First the platelets were disposed of by three consecutive centrifugations for 15 min at 200g. After each centrifugation the supernatant was eliminated and 2 mL PBS + 2% FBS and 1 mM EDTA (ED2SS, Sigma, USA) were added. Then, the monocytes were purified using the RosetteSep Human Monocyte Enrichment Cocktail kit (15028, StemCell, USA). The purified fraction was diluted in 500 µL of PBS + 2% FBS and 1 mM EDTA. The ABX Pentra DX 120 (Horiba Medical, Japan) hematology analyzer was used to verify the purification process. We established that the sub-population were purified with a rate above 98% for granulocytes, 95% for lymphocytes and 74% for monocytes.

Then, the purified WBC samples were specifically labeled with Thiazole Orange dye (390063, Sigma-Aldrich Inc., USA). This nucleic acids marker is known to have a high quantum yield [START_REF] Lee | Thiazole orange: A new dye for reticulocyte analysis[END_REF] . A marking solution was obtained by dissolving at 1 mg/mL Thiazole Orange in methanol, and by then diluting it to 18 µg/mL in PBS (D8537, Sigma-Aldrich Inc., USA). The samples were incubated for 30 minutes with 1 mL of the marking solution.

Finally, blood smears were made on strandard microscope slides using the automated blood smearing device Hemaprep (Horiba Medical, Japan).

For a healthy person, WBC account for roughly 0.1% of the blood cells. The wide field of view of the systems allowed to detect a large number of labelled cells. In this experiment, we acquired sequentially a brightfield and a fluorescence image. As pointed out previously, the bimodal images required no registration thanks to the common imaging path. On a single image, we were able to detect 1 000 -3 000 WBC, as illustrated in Figure 6. To demonstrate the performance of our system, we compared the fluorescence images from our set up to those acquired using a standard microscope (BX51, Olympus, USA) with a 40x objective. In particular, the morphology of the nuclei of WBC were examined. Figure 7 shows similar morphologies of WBC nuclei imaged using the two systems. Although the images were acquired on the same slide and on the same day, it is not certain that the exact same cells were imaged. The registration of the images from the reference microscope and our system is difficult given the size of the FOVs.

Granulocytes are the most common WBC type (40 to 75%), they are roughly 12 µm in diameter and their nucleus has a characteristic multi-lobe or irregular-shape. As illustrated on a bi-lobes nucleus in Figure 7, our system's resolution was sufficient to distinguish the morphological features of granulocytes. Lymphocytes account for 20 to 40 % of WBC, are regular and round-shaped cells of 7 µm diameter with a nucleus occupying almost the entire cell. Monocytes represent 2 to 10% of WBC, are large 15-20 µm cells with a granular horse shoe shaped nucleus. As illustrated in Figure 7, lymphocytes and monocytes could be differentiated. The size difference of the nuclei of two subtypes were discernable, and the fluorescence signal of monocytes were more diffuse. These images indicate that our system has a sufficient resolution to point out the key features of each WBC sub-type. Thus, these results suggest that it might be possible to perform statistical determination of the ratios of WBC populations within a whole blood sample.

CONCLUSION

We described an imaging system combining in-line holography, brightfield and fluorescence over a 30 mm 2 FOV, achieving a 3-4 µm resolution. This system takes advantage of a SLR camera lens with a front lens wider than those usually used in microscopy. Such a system can target many biological applications. To illustrate its potential, we chose to show its performances in imaging WBC sub-types. Results suggest that the system does combine a resolution sufficient to distinguish the morphological features of the nuclei, while the FOV is sufficiently large to be able to infer statistically relevant information. The next step of this work would be to image a blood smear of a whole i.e. not purified blood sample and show that it is possible to classify the detected WBC into distinct populations.

Figure 1 -

 1 Figure 1 -Theoretical comparison of the FOV of objectives depending on their magnification.

Figure 2 -

 2 Figure 2 -Wide FOV multi-modal imaging system

Figure 3 -

 3 Figure 3 -Excitation and emission spectra of Thiazole Orange dye. The pass-band filters used for excitation and emission are also shown.

Figure 4 -

 4 Figure 4 -USAF resolution target illustrating the wide FOV and the resolution.

Figure 5 -

 5 Figure 5 -PSF. Image of a fluorescent GFP -1 µm bead. Line profile pointing out the FWHM.

Figure 6 -

 6 Figure 6 -Superposed transmission (in grey scale) and fluorescent (artificial green color) images of labelled WBC.

Figure 7

 7 Figure 7 Acquisitions of the fluorescence of granulocytes, lymphocytes and monocytes (respectively) with a standard epifluorescence microscope with a 40x objective and with our system. The field of views are cropped. Scale bar is 10 µm.