Rémi Colom
email: remi.colom@fresnel.fr

Ross Mcphedran

Brian Stout

Nicolas Bonod
email: nicolas.bonod@fresnel.fr

Modal Analysis of Anapoles, Internal Fields and Fano Resonances in Dielectric Particles

All-dielectric nanostructures have aroused strong interest because of their potential to trap light at subwavelength scales and to yield strong internal electric or magnetic field enhancements. Optimizing the internal fields appears as a crucial challenge for enhancing light matter interactions in all-dielectric nanostructures. Mie resonators host radiationless states called anapoles associated to a pronounced minimum of light scattering. However, the question is to know whether these radiationless states maximize the internal field intensities. Here we use a modal expansion of the internal and external fields to demonstrate that anapoles in dielectric Mie resonators result from a Fano resonance produced by the interference between two mode contributions with low and high quality factors plus an additional non-resonant term. A modal expansion of the internal field enhancement averaged over the whole scatterer volume shows that the maximum of the internal field enhancement does not occur at the anapole frequency but at the real part of the eigen-frequency associated with the high quality factor. This analysis is carried out for both electric and magnetic modes of the dielectric scatterer and evinces that a larger internal field enhancement is found for the first magnetic dipole resonance even though this resonance does not feature an anapole behavior.

Introduction

The field of resonant nanophotonics aims at enhancing light matter interactions at nanometer scales due to the resonant coupling of light with subwavelength sized particles [START_REF] Koenderink | Nanophotonics: Shrinking light-based technology[END_REF][START_REF] Kuznetsov | Optically resonant dielectric nanostructures[END_REF][START_REF] Krasnok | Spectroscopy and biosensing with optically resonant dielectric nanostructures[END_REF]. The excitation of plasmonic and Mie resonances in such particles permits the enhancement of light scattering as well as the near field in the direct vicinity of the particle [START_REF] Devilez | Mimicking localized surface plasmons with dielectric particles[END_REF]. In the case of dielectric particles, the internal field can be strongly enhanced without loss, which opens unique opportunities in non-linear optics to generate second or third harmonics with semi-conductor particles [START_REF] Shcherbakov | Enhanced third-harmonic generation in silicon nanoparticles driven by magnetic response[END_REF][START_REF] Tribelsky | Giant in-particle field concentration and fano resonances at light scattering by high-refractive-index particles[END_REF][START_REF] Gili | Monolithic algaas second-harmonic nanoantennas[END_REF][START_REF] Kapitanova | Giant field enhancement in high-index dielectric subwavelength particles[END_REF][START_REF] Colom | Enhanced four-wave mixing in doubly resonant si nanoresonators[END_REF]. It was shown in particular that the magnetic modes yield stronger confinement of the electric field intensity inside the resonant particles than the electric modes [START_REF] Shcherbakov | Enhanced third-harmonic generation in silicon nanoparticles driven by magnetic response[END_REF][START_REF] Shcherbakov | Ultrafast all-optical switching with magnetic resonances in nonlinear dielectric nanostructures[END_REF]. Recently, a very interesting phenomenon was also identified in the field of all-dielectric nanophotonics with anapoles [START_REF] Miroshnichenko | Nonradiating anapole modes in dielectric nanoparticles[END_REF][START_REF] Basharin | Dielectric metamaterials with toroidal dipolar response[END_REF][START_REF] Luk'yanchuk | Hybrid anapole modes of high-index dielectric nanoparticles[END_REF][START_REF] Luk'yanchuk | Suppression of scattering for small dielectric particles: anapole mode and invisibility[END_REF][START_REF] Svyakhovskiy | Transient effects at resonant light scattering by particles: Anapole as a storage for" frozen light[END_REF][START_REF] Gongora | Fundamental and high-order anapoles in all-dielectric metamaterials via fano-feshbach modes competition[END_REF][START_REF] Baryshnikova | Optical anapoles: Concepts and applications[END_REF]. Anapoles are classically defined in physics as non-radiating sources of current [START_REF] Zel'dovich | Electromagnetic interaction with parity violation[END_REF][START_REF] Afanasiev | Some remarkable charge-current configurations[END_REF][START_REF] Papasimakis | Electromagnetic toroidal excitations in matter and free space[END_REF]. They have aroused a keen interest on account of their potential ability to enhance the light absorption in the context of dielectric resonators with promising applications in non-linear optics, light absorption and enhanced Raman scattering [START_REF] Grinblat | Enhanced third harmonic generation in single germanium nanodisks excited at the anapole mode[END_REF][START_REF] Grinblat | Efficient third harmonic generation and nonlinear subwavelength imaging at a higher-order anapole mode in a single germanium nanodisk[END_REF][START_REF] Wei | Excitation of the radiationless anapole mode[END_REF][START_REF] Wang | Engineering non-radiative anapole modes for broadband absorption enhancement of light[END_REF][START_REF] Timofeeva | Anapoles in free-standing iii-v nanodisks enhancing second-harmonic generation[END_REF][START_REF] Rocco | Tuning the second-harmonic generation in algaas nanodimers via non-radiative state optimization[END_REF][START_REF] Yang | Anapole-assisted strong field enhancement in individual all-dielectric nanostructures[END_REF][START_REF] Baranov | Anapole-enhanced intrinsic raman scattering from silicon nanodisks[END_REF][START_REF] Xu | Boosting third-harmonic generation by a mirror-enhanced anapole resonator[END_REF].

In many studies, the choice has been made to forsake the multipolar theory in spherical coordinates, a theory that is widely used for explaining Mie resonances in high refractive index materials, in favor of the multipolar theory in cartesian coordinates including the toroidal multipoles introduced by Zel'dovitch in 1958 [START_REF] Zel'dovich | Electromagnetic interaction with parity violation[END_REF][START_REF] Afanasiev | Some remarkable charge-current configurations[END_REF][START_REF] Papasimakis | Electromagnetic toroidal excitations in matter and free space[END_REF]. In fact, the toroidal-type currents are the higher-order terms appearing in the power series expansion of the spherical electric multipoles with respect to the size of the current source. [START_REF] Fernandez-Corbaton | On the dynamic toroidal multipoles from localized electric current distributions[END_REF][START_REF] Alaee | Exact multipolar decompositions with applications in nanophotonics[END_REF]. In this framework, anapoles occur when the electromagnetic field radiated by these toroidal multipoles cancels the field radiated by the electric cartesian multipoles. This phenomenon can also be understood by using a Quasi-Normal Mode (QNM) expansion. A QNM analysis associated with a multipolar decomposition pointed out the destructive interference between the first two modes of the electric dipolar coefficient that drives to a dip in the scattering spectrum [START_REF] Powell | Interference between the modes of an all-dielectric meta-atom[END_REF]. The role of the first two electric modes of the electric dipolar coefficient in the excitation of anapoles in silicon scatterers was confirmed in the framework of a multipolar analysis [START_REF] Stout | Egocentric physics: just about mie[END_REF]. Lastly, a toy model was developed on the analogy between optical resonators and the dynamic of a forced pendulum to study the transient internal field yielded by anapoles [START_REF] Svyakhovskiy | Transient effects at resonant light scattering by particles: Anapole as a storage for" frozen light[END_REF].

Here we use a modal expansion of the internal and scattered fields to demonstrate that anapoles in spherical Mie resonators arise from the fact that the scattered field features Fano-type resonances while the internal field features Lorentzian responses only. We perform the modal analysis for both electric and magnetic cases and explain why anapoles are not observed for the first resonance of the magnetic dipole. It is also pointed out that the internal electric field enhancement is larger for the first magnetic dipole resonance than for its electric dipole counterpart. The fact that the strongest field enhancement is observed in absence of anapole condition is properly predicted by our modal analysis. Results are obtained by carrying out a modal analysis of the T-matrix scattering operator in the multipolar Mie basis [START_REF] Stout | Egocentric physics: just about mie[END_REF][START_REF] Grigoriev | Optimization of resonant effects in nanostructures via weierstrass factorization[END_REF][START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF]. We believe that determining the anapole conditions in high refractive index scatterers in a widely used multipolar theory facilitates the understanding of the origin of this phenomenon.

Pole expansion of the T-matrix coefficients

T-matrix coefficients

Here, we are studying the anapole condition for a spherical scatterer located in vacuum with a radius R and a dielectric perimittivity ε s . Since only non-magnetic materials are considered, the refractive index possesses the following definition n s = √ ε s . When studying light-scattering in 3D, in particular for spherically-symmetric scatterers, it is convenient to make use of the set of Vector Partial Waves (VPWs) as a basis for the fields of the problem [START_REF] Hulst | Light scattering by small particles, Structure of matter series[END_REF][START_REF] Bohren | Absorption and scattering of light by small particles[END_REF][START_REF] Chew | Waves and fields in inhomogeneous media[END_REF]. In this section, we shall consequently start by introducing the definitions used throughout this paper for the Vector Partial Waves. The outgoing VPWs are defined as follows:

M (+) n,m (kr) = h (+) n (kr)X n,m (θ, φ), N (+) n,m (kr) = 1 kr n(n + 1)h (+) n (kr)Y n,m (θ, φ) + ξ (+) n (kr)Z n,m (θ, φ) , (1)
where X n,m (θ, φ), Y n,m (θ, φ) and Z n,m (θ, φ) are the Vector Spherical Harmonics [START_REF] Stout | Multipole methods for nanoantennas design: applications to yagi-uda configurations[END_REF], k = 2π λ , λ being the wavelength in the embedding medium and ξ (+) n (kr) is the derivative with respect to kr. n and m are the total and projected angular momentum of the VPW, n is classically called the multipolar order. h (+) n (x) and ξ (+) n (x) are respectively the spherical Hankel and the associated Ricatti-Hankel functions of the first kind. The regular VPWs are defined as follows:

M (1) n,m (kr) = j n (kr)X n,m (θ, φ), N (1) n,m (kr

) = 1 kr n(n + 1) j n (kr)Y n,m (θ, φ) + ψ n (kr)Z n,m (θ, φ) , (2)
where j n (x) and ψ n (x) are respectively the spherical Bessel and associated Ricatti-Bessel functions. In the T-matrix framework, the total field outside the scatterer is split into the excitation and scattered field contributions. The expansions of the excitation and scattered electric fields read as follows:

E exc (kr) = E 0 ∞ n=1 n m=-n e (h)
n,m M (1) n,m (kr) + e (e) n,m N (1) n,m (kr)

E scat (kr) = E 0 ∞ n=1 n m=-n f (h) n,m M (+) n,m (kr) + f (e) n,m N (+) n,m (kr) (3)
The internal field takes the following expression:

E int (k s r) = E 0 ∞ n=1 n m=-n u (h) n,m M (1) n,m (k s r) + u (e) n,m N (1) n,m (k s r) (4)
where k s = n s k. The T-matrix coefficients are defined as follows

T (i) n (z) = f (i) n, m (z) e (i) n, m
where i = e, h and z = k R. For spherical scatterers, the T-matrix coefficients can be computed using the continuity conditions leading to the following expression

T (i) n (z) = - N (i) T, n (z) D (i) n (z)
in terms of the following numerator and denominator functions [START_REF] Stout | Egocentric physics: just about mie[END_REF]:

N (e) T,n (z) = z ε s j n (n s z)ψ n (z) -ψ n (n s z) j n (z) in s , D (e) n (z) = z ε s ξ (+) n (z) j n (n s z) -ψ n (n s z)h (+) n (z) in s , N (h) T,n (z) = -iz ψ n (z) j n (n s z) -ψ n (n s z) j n (z) , D (h) n (z) = -iz ξ (+) n (z) j n (n s z) -ψ n (n s z)h (+) n (z) , (5)
It is also necessary to introduce the Ω-matrix coefficients linking the internal field to the incident field,

Ω (i) n (z) = u (i) n, m (z) e (i) n, m
. These coefficients can also be simply expressed by use of the denominator functions

Ω (i) n (z) = 1 D (i) n (z)
.

Modal expansion of the scattered efficiency and internal field enhancements

QNM analysis is common in wave physics [START_REF] Kokkotas | Quasi-normal modes of stars and black holes[END_REF][START_REF] Chang | Optical processes in microcavities[END_REF] and has attracted a strong interest in photonics when studying resonant light interactions with nanocavities [START_REF] Grigoriev | Optimization of resonant effects in nanostructures via weierstrass factorization[END_REF][START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF][START_REF] Sauvan | Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators[END_REF][START_REF] Vial | Quasimodal expansion of electromagnetic fields in open two-dimensional structures[END_REF][START_REF] Kristensen | Normalization of quasinormal modes in leaky optical cavities and plasmonic resonators[END_REF][START_REF] Zambrana-Puyalto | Purcell factor of spherical mie resonators[END_REF][START_REF] Muljarov | Exact mode volume and purcell factor of open optical systems[END_REF][START_REF] Lassalle | Interplay between spontaneous decay rates and lamb shifts in open photonic systems[END_REF][START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF][START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF]. QNMs are the eigen-vector solutions of the Maxwell equations without sources and are associated with the complex eigenvalue frequencies p (i) n,α that also correspond to the poles of the scattering operators introduced in the previous section. For each multipolar order n, an infinite number of poles α exists. Poles associated with positive α have a positive real part while poles associated with negative α have a negative real part and are associated via the relation, p n,-α = -p * n,α . In this study, we determine the eigen-frequencies p (i) n,α by calculating the poles of the scattering matrix. We use root-finding methods to locate the position of these poles starting from the asymptotic formulas of these poles provided in ref. [START_REF] Stout | Egocentric physics: just about mie[END_REF]. We plot the real and imaginary parts of the first eigen-frequencies α = [-5; 5] of the electric dipole (Fig. 1a) and magnetic dipole (Fig. 1b) in the case of a spherical particle made of dielectric permittivity ε s = 16. The fundamental mode is α = 0 for the magnetic case and α = ±1 for the electric case. It is characterized by a much larger leakage, i.e. a larger imaginary part in absolute value since poles in passive media have negative imaginary parts, than higher order modes.

Following the derivation of the scattering operator by Grigoriev et al. [START_REF] Grigoriev | Optimization of resonant effects in nanostructures via weierstrass factorization[END_REF], the modal expansions of the T (i) n and Ω (i) n coefficients were recently derived [START_REF] Stout | Egocentric physics: just about mie[END_REF][START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF]:

T (i) n (z) = S (i) nr,n e -2iz -1 2 + e -2iz 2 +M α=-M r (i) n,α z -p (i) n,α
,

Ω (i) n (z) = 1 2 +M α=-M r (i) Ξ,n,α z -p (i) n,α , (6)
where r (i) n,α is the residue of T (i) n at the poles p (i) n,α , r (i) Ξ,n,α is the residue of Ω (i) n (z) and

S (i) nr,n = 1 + +M α=-M r (i) n, α p (i) n, α .
It can be shown that the convergence of this truncated series can be improved for a given truncation order M by approximating the higher order tem α > M [START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF]. At the first order, the term that is removed when truncating the sum can be approximated by [START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF]. The contribution from higher order poles |α| > M is then included in the term S (i) nr,n in order to improve the accuracy of the truncated expansion. The notation

-∞ α=M+1 r (i) n, α p (i) n, α --∞ α=-(M+1) r (i) n, α p (i) n, α . This term is actually equal to -∞ α=M+1 r (e) n, α p (e) n, α + r (e) n,-α p (e) n,-α = -(-1) n+1 -1 + M α=-M r (e) n, α p (e) n, α and
-∞ α=M+1 r (h) n, α p (h) n, α + r (h) n,-α p (h) n,-α = -((-1) n -1) + M α=-M r (h) n, α p (h) n, α [
A (i) n (z) = S (i) nr, n e -2i z -1 2
will be used in what follows.

Pole expansions in Eqs. ( 6) thus show the existence of a non-resonant term for the T-matrix coefficient only. The existence of a non-resonant contribution is also predicted by the Temporal Coupled Mode Theory (TCMT) [START_REF] Fan | Temporal coupled-mode theory for the fano resonance in optical resonators[END_REF][START_REF] Ruan | Temporal coupled-mode theory for fano resonance in light scattering by a single obstacle[END_REF][START_REF] Ruan | Temporal coupled-mode theory for light scattering by an arbitrarily shaped object supporting a single resonance[END_REF]. However, TCMT is a phenomenological model that requires fitting parameters while Eqs. ( 6) derived here are free from fitting parameters and only rely on energy conservation principle [START_REF] Grigoriev | Optimization of resonant effects in nanostructures via weierstrass factorization[END_REF][START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF]. This pole expansion permits us to derive the modal expansion of the multipolar scattering coefficients and therefore also the modal expansion of the scattering efficiency Q (i) scat,n and internal field intensities averaged over the volume of the particle |E (i) int,n | 2 . The analytical expression of the residues for spherical scatterers can be obtained from the analytical expressions of the Ω and T matrices [START_REF] Stout | Egocentric physics: just about mie[END_REF][START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF]:

r (e,h) Ξ,n,α = 1 d dz D (e,h) n (z)| z=p (e , h) n, α (7)
that yields the following expression for r (e,h) Ξ,n,α :

r (e) Ξ,n,α = 1 ε s -1 × 2in s ξ (+) n p (e) n,α ψ n n s p (e) n,α + n(n + 1)h (+) n p (e) n,α j n n s p (e) n,α r (h) Ξ,n,α = 2i (ε s -1) p (h) α 2 j n n s p (h) α h (+) n p (h) α
where the following identities of the Ricatti-Hankel and Ricatti-Bessel functions have been used:

ξ (+) n (z) = n(n+1)-z 2 z h (+) n (z) and ψ (+) n (z) = n(n+1)-z 2 z j n (z). r (e,h)
n,α can also be analytically determined for spherically-symmetric scatterers from: r (e,h) n,α = -

e 2ip (e ,h) n, α N (e,h) n p (e,h) n,α d dz D (e,h) n (z)| z=p (e ,h) n, α (8)
The set of r (e,h) n,α consequently takes the following expression:

r (e) n,α = j n n s p (e) n,α h (+) n p (e) n,α 1 ε s -1 × 2iε s e 2ip (e) n, α ξ (+) n p (e) n,α ψ n n s p (e) n,α + n(n + 1)h (+) n p (e) n,α j n n s p (e) n,α r (h) n,α = j n n s p (h) n,α h (+) n p (h) n,α 2ie 2ip (h) n, α (ε s -1) p (h) n,α 2 j n n s p (h) n,α h (+) n p (h) n,α
The scattering efficiency Q scat of a scatterer can be calculated with respect to the electric and magnetic scattering Mie coefficients a n and b n that are simply the opposite of the T-matrix coefficients a n = -T (e) n and b n = -T (h) n :

Q (e) scat,n = 2(2n + 1) z 2 |a n | 2 = 2(2n + 1) z 2 |T (e) n | 2 Q (h) scat,n = 2(2n + 1) z 2 |b n | 2 = 2(2n + 1) z 2 |T (h) n | 2 (9)
A similar approach can be followed with the internal Mie coefficients d n and c n that can be derived from the Ω matrix elements:

d n = Ω (e) n , c n = Ω (h) n .
The average of the electric field intensity over the volume inside the spherical particle |E int | 2 can be calculated with respect to the multipolar components of the electric field E (i) int,n [START_REF] Stout | Egocentric physics: just about mie[END_REF]:

|E int | 2 = n |E (e) int,n | 2 + |E (h) int,n | 2 (10)
where :

|E (e) i,n | 2 |E 0 | 2 = 3(2n + 1) 2z 2 I n (ε s , z)|Ω (e) n | 2 , |E (h) i,n | 2 |E 0 | 2 = 3(2n + 1) 2z 2 J n (ε s , z)|Ω (h) n | 2 , (11)
with:

I n (ε s , z) = z n 2 s × ψ n (n s z) ψ n (n s z) + j n (n s z) + (n s z) 2 -n(n + 1) j 2 n (n s z) J n (ε s , z) = z n 2 s × ψ n (n s z) ψ n (n s z) -j n (n s z) + (n s z) 2 -n(n + 1) j 2 n (n s z).
The modal expansion of the internal field enhancement in the whole scatterer allows to determine the mode of each multipolar order that yields the strongest field enhancement inside the particle.

Combined with the modal expressions of the scattered efficiency (Eq. 9), we can compare the spectra of the internal and external fields and analyze the origin of the scattering anomalies in terms of QNMs.

Results and discussions

Let us first compare the spectra of the averaged field enhancement with the scattering efficiencies in the case of a spherical particle of dielectric permittivity ε = 16, value closed to the dielectric permittivity of silicon in visible and near infrared spectrum (see Fig. 2). For that purpose, we plot in Fig. 2(a-b) the scattering efficiency associated with the electric Q (e) scat,1 and magnetic

Q (h) scat,1
dipolar contributions, and of the averaged field intensity inside the particle |E (i) int,1 | 2 , with respect to the size parameter z = k R when the particle is illuminated by a plane wave.

Q (i)
scat,1 and |E (i) int,1 | 2 both exhibit a resonant behavior with the presence of several peaks when increasing the parameter size z. These plots are obtained using conventional multipolar formalism and will be compared later in Figs. 3,4 with results calculated with our modal expansion. As can be seen in Fig. 2(a-b), the peaks associated with

Q (i)
scat,1 and |E (i) int,1 | 2 , for both the electric i = e (Fig. 2a) and magnetic i = h (Fig. 2b) cases occur at nearby frequencies. A careful analysis reveals that the peak frequencies are slightly different for the electric case, while they are almost equal for the magnetic case.

A more fundamental difference can be noticed when observing the shape of the resonances: while |E (e) int,1 | 2 and |E (h) int,1 | 2 feature pure Lorentzian shapes, the scattering efficiencies

Q (i)
scat,1 In particular, for the electric case, the maximum of Q (e) scat,1 is closely followed by a minimum of the scattering efficiency. For the magnetic case, the first peak at z ≈ 0.75 features an almost perfect Lorentzian shape while the second peak at z ≈ 1.5 is associated with a strong minimum.

(i = e,
As the different multipolar contributions are well separated, these minima of Q (e) scat,1 and

Q (h)
scat,1 correspond to low total scattered fields as can be seen in Fig. 2(c). Importantly, position of the anapole condition, corresponding to a minimum of Q (e) scat,1 and Q (h) scat,1 and indicated by vertical dashed lines in Fig. 2(c), can be well predicted when considering only dipolar contributions. The minimum is more pronounced for the electric dipole coefficient than for the magnetic dipole due to the presence of non-negligible contributions from higher order modes. The higher order modes are at the origin of the narrow peaks observed in Fig. 2c between z = [1; 1.5]. It is also observed that at these minimum frequencies, the average field intensity remains enhanced by the resonances which means that, at certain frequencies, the scattering efficiency can be strongly minimized while the internal field is enhanced by a resonant behavior of the internal scattered field coefficients: this condition is commonly referred as an anapole condition [START_REF] Miroshnichenko | Nonradiating anapole modes in dielectric nanoparticles[END_REF][START_REF] Grinblat | Enhanced third harmonic generation in single germanium nanodisks excited at the anapole mode[END_REF][START_REF] Wei | Excitation of the radiationless anapole mode[END_REF][START_REF] Wang | Engineering non-radiative anapole modes for broadband absorption enhancement of light[END_REF][START_REF] Timofeeva | Anapoles in free-standing iii-v nanodisks enhancing second-harmonic generation[END_REF].

The different terms of the modal expansion of T (e) 1 (z) (top line) and T (h) 1 (z) (bottom line) (Eq. 6) are plotted in Fig. 3. Two types of terms can be distinguished: there is (i) a non-resonant term associated with

A (i) n (z) = S (i) nr, n e -2i z -1 2
that corresponds to a slowly varying background (see Fig. 3a,d) [START_REF] Grigoriev | Optimization of resonant effects in nanostructures via weierstrass factorization[END_REF][START_REF] Grigoriev | Singular analysis of fano resonances in plasmonic nanostructures[END_REF] and (ii) a resonant term e -2i z (see Fig. 3b,e). The plot of the non-resonant term in Fig. 3(a,d) highlights the π phase shift between the electric and magnetic cases: while this non-resonant term is maximum at z = 0 and null at z = 1.5 for the magnetic case, it is null at z = 0 and maximum at z = 1.5 for the electric case. Therefore the non-resonant contribution plays a larger role for the electric low order resonances than for the magnetic ones.

Let us now investigate the role of the resonant terms by plotting in Fig. 3(b,e) the contributions of the first three modes separately: α = [1; 3] for the electric case and α = [0; 2] for the magnetic case. In all cases, it can be seen that they correspond to pure Lorentzian shapes. For both electric and magnetic cases, the first mode (α = 1 for electric, α = 0 for the magnetic) has a very poor quality factor which means that it overlaps with the higher order modes. An important question arises now: can a minimum of the far-field spectrum arise when taking into account the sum of the mode contributions only, i.e. without taking into account the A (i) 1 (z) term? To answer this question, we plot in grey line the resonant term 4b,e (see Eq.6) while considering a large number of poles (e.g. M = 100), but neglecting the non-resonant term

| e -2i z 2 +M α=-M r (i) 1, α z-p (i) 1, α | 2 in Fig.
A (i)
1 (z). The sum of the different pole contributions does not lead to Lorentzian responses but to asymmetric Fano-like shapes featuring strong minima. This result is a direct consequence of the fact that the different modes interfere. In fact, the QNM are not orthogonal in the sense of the usual conjugated scalar product. This means that energy can be exchanged between different modes leading to a phenomenon of interference. That is the reason why interference between modes may appear in the scattering efficiency spectrum. However, let us point out that an unconjugated scalar product can also be defined for which QNM become orthogonal even though they can exchange energy [START_REF] Muljarov | Exact mode volume and purcell factor of open optical systems[END_REF][START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF][START_REF] Stout | Eigenstate normalization for open and dispersive systems[END_REF]. In Fig. 3b, we see that, in the electric case, the first

minimum of | e -2i z 2 +M α=-M r (e) 1, α z-p (e) 1, α
| 2 occurs when the modes α = 1 and α = 2 that have the same real frequency are in phase opposition. The difference of phase between the modes involved in the destructive interference are plotted in Fig. 4e-f, i.e. φ 1 -φ 2 in the electric case and φ 0 -φ 2 in the magnetic case where

φ α = Arg e -2i z 2 r (i) 1, α z-p (i) 1, α
. We can see in Fig. 4e-f that the difference of phase turns out to be equal to π near the destructive interference between the modes. For these two modes to destructively interfere, not only should they be in phase opposition but they also need to have the same amplitude. The amplitude is quantified by the expression |

r (i) α, n z-p (i) α, n | 2 , which
shows that a resonant contribution is proportional to the residue r (i) α,n and requires poles p (i) α,n with small imaginary parts to maximize 1 z-p | (residues are plotted with respect to α in Fig. 1). This condition is fulfilled near p (e) 1,1 ≈ p (e) 1,2 . In the electric case, the sum of the resonant terms | r (h) 1,α /(zp 1,α ) (h) | displayed in Fig. 3(e) shows that a minimum is associated with both α = 1 and α = 2 modes since they both interfere destructively with the low quality α = 0 mode. Minima always occur at the vicinity of the real part of one pole,

p (i) 1,α , since the function e -2i z 2 r (i) 1, α z-p (i) 1, α
experiences an abrupt change of phase.

The resonant contribution from this mode can thus become out of phase with the contribution of the low quality mode. A similar explanation holds for the appearance of minimum scattering in the magnetic case: it results from the destructive interference between the wide resonance associated with pole α = 0 with a narrow resonance associated with poles α = 1 and α = 2. However, in this case, the first minimum close to the pole α = 1 disappears when the non-resonant contribution is taken into account. Moreover, concerning the minimum associated with α = 2, the shift due to the non-resonant contribution is minor as A (h) 1 (z) is almost null in the vicinity of the peak due to the pole α = 2 (see Fig. 3d).

At this stage, it has been established that minima in the scattering spectra result from the interference between a low and a high quality factor mode. We now study the influence of the non-resonant contribution on the spectral location of such minima. For that purpose, the spectrum

φ α = Arg e -2i z 2 r (i) 1, α z-p (i) 1, α .
of the resonant contribution only is compared with exact calculations and modal calculations that include both resonant and non-resonant terms in Fig. 3(c,f). A spectral shift between these spectra with or without non-resonant contribution can be observed, especially for the electric case. The spectral shift is less pronounced for the magnetic case (Fig. 3(f)) as the magnetic non-resonant term A (h) n (z) is almost null when z ≈ 1.6, i.e. close to the frequency where the modes destructively interfere (Fig. 3(c)). In all cases, it is necessary to add up the non-resonant contribution A (i) n (z) in order to accurately retrieve the position of the peaks and dips of

|T (i) 1 | 2 .

Let us now explain why the first resonance of |T (h)

1 | at z ≈ 0.75 is not associated with an anapole. Indeed, a strong minimum is observed near p (h) 1,2 when the resonant contributions only are summed. This minimum completely disappears when the non-resonant contribution is included as shown in Fig. 3(f). This shows that the non-resonant contribution destructively interferes with the contribution from α = 0. For the first peak of |T (h) 1 | 2 , the contribution from the mode α = 1 is predominant and a peak with almost a Lorentzian shape is observed in Fig. 3(f). In Fig. 3(c,f), it is shown that a near perfect agreement with exact calculations (Mie theory) is found while taking into account both resonant (M = 100) and non-resonant contributions. Minima of the scattering spectrum therefore result not only from (i) the interference between the mode contributions (the eigenfields associated with complex frequencies are not orthogonal) but also from (ii) the presence of a non-resonant term.

The internal field can also be studied in the framework of the modal expansion in Eqs. 10-11. The accuracy and the convergence of this pole expansion of the internal field are illustrated in Fig. 4 where the exact value of |E (i) int,1 | 2 is compared to the pole expansions in Eqs. 10-11 when considering one or two poles, i.e. M = 5 (top line: electric case (i = e); bottom line: magnetic case (i = h)). Each of the Lorentzian-type resonances appearing in the spectrum of the internal field is associated with one complex resonance frequency. Remarkably, it can be observed that the modal expansion fits perfectly the rigorous calculations with a very limited number of poles. While a large number of poles was required to match the rigorous calculations of the scattered field, the convergence is much more rapidly achieved in the case of the internal field.

In Fig. 4(a,c), we focus on the first resonance of the electric dipolar scattering coefficient that exhibits an anapole behavior. The explanation of the dip in the spectrum of the scattering coefficient was earlier explained at length. Just as in the electric case, the resonance of the internal field is mostly due to the resonant contribution associated with the mode α = 2. The mode α = 1 only contributes weakly to the internal field enhancement; its contribution is too weak to be visible in Fig. 4 whereas it played a key role in the appearance of a dip in the scattering coefficient spectrum. This is the main difference between external and internal fields: the fundamental modes, α = 0 in the magnetic case and α = 1 in the electric case, contribute only for the external field, and not for the internal field. This is due to the fact that the residues of the fundamental modes in the internal fields are much smaller than those of higher order modes. The first poles of the dipolar electric Mie coefficients and the modulus of the associated residues r (i) Ξ,1,α and r (i) 1,α for an ε = 16 particle in air are plotted in Fig 1.

As the quality factor associated with the fundamental mode is also much weaker than that of higher order modes (see Fig. 1), we get the relation:

r (e) Ξ, 1, α=1 z-p (e) 1, α=1 r (e) Ξ,1, α>1 z-p (e) 1, α>1 and r (h) Ξ,1, α=0 z-p (h) 1, α=0 r (h) Ξ,1, α>0 z-p (h) 1, α>0
.

Physically, this means that in the external field, the fundamental mode with low quality factor can therefore interfere with the higher order modes of higher quality factor leading to Fano resonances, while in the internal field, this fundamental mode does not interfere with higher order modes and the averaged internal field only exhibits Lorentzian-like resonances driven by higher order modes that are spectrally well separated. Fig. 4(a,c) also show that the position of the maximum internal field enhancement would correspond to the position of the dip of the scattering coefficients without the additional blue-shift due to the non-resonant term. However, this minimum shifts when the non-resonant contribution is added so that the minimum of the scattered field does not correspond in general to the maximum of the internal field.

In Fig. 4(b,d), we focus on the second resonance of T (h) 1 . We do not concentrate on the first resonance since it does not exhibit an anapole behavior as explained earlier. The mode α = 0 that interfered with the mode α = 2 to yield the dip in the spectrum of T (h) 1 (z) contributes only weakly to the internal field. The resonance of the averaged internal field is then almost entirely due to the mode α = 2 and thus possesses a Lorentzian shape.

The distributions of the electric field intensity in the dielectric Mie scatterer are plotted in Fig. 5 at frequencies that maximize the average internal field intensity (z = 1.065) and minimize the scattering efficiency (z = 1.160). The field distributions are significantly different and we remark (i) the anapole-like field distribution at the minimum of the scattering efficiency and (ii) the stronger field enhancement at a frequency different from the anapole frequency. These results highlight the fact that anapoles do not necessarily lead to an optimal internal field enhancement in average inside the scatterer. Also, when comparing the averaged field intensity enhancement between the electric dipole, electric anapole and magnetic dipole, it is clear that the maximum of field enhancement is achieved for the first magnetic dipolar mode (see also Fig. 2).

After having thoroughly analyzed the anapole conditions for the important case of ε = 16, one now aims at generalizing this analysis for arbitrary dielectric permittivities. As anapole conditions can be predicted by calculating the poles of the resonant scatterer, the trajectories of the first 5 poles of the electric and magnetic dipolar Mie coefficients p (e) 1,α , p (h) 1,α , α = [1; 5] are plotted in the complex frequency planes when increasing the dielectric permittivity of the resonant scatterer. For both magnetic and electric cases, it turns out that the real and imaginary parts of the eigen-frequencies p 1,α of the first poles (α = 0 for the magnetic case and α = 1 for the electric case) are almost insensitive to the dielectric permittivity and that the imaginary part is much higher in magnitude than all the other poles. p (h) 1,0 features the peculiarity of having a real part equal to 0. On the other hand, all the other poles depend strongly on the dielectric permittivity of the scatterer. Their real parts decrease, meaning that resonances occur for smaller size parameters z, and their imaginary parts decrease in magnitude and tend toward zero. The high quality factor of the higher order modes will prohibit any significant interference effect between them but the low quality factor of the fundamental mode leads to interference effects with higher order modes. The anapole location on the spectrum is therefore driven by the high quality factor resonance, i.e. the scattering spectrum features sharp minima close to the real parts of the high quality modes. For the electric case, the blue-shift of the anapole due to the non-resonant term decreases with increasing ε since the non-resonant term A (e) 1 (z) → 0 when z → 0. The contrary happens for the magnetic mode since the electric and magnetic non-resonant terms are in phase opposition as illustrated in Fig. 3(a,d). Regarding the internal field, it should be first noticed that the internal field enhancement depends on the quality factor that describes the ability of the cavity to trap light. Therefore, the internal field associated with the fundamental modes characterized by large values of (p (i) n,α ) do not lead to significant internal field enhancements. As all the higher order modes (α > 0 for the magnetic case and α > 1 for the electric case) do not overlap due to their narrowness, and due to the absence of a non-resonant term for the internal coefficients, the averaged field inside the scatterer features almost pure Lorentzian shapes.

Conclusion

To conclude, anapoles in Mie resonators result from the existence of electric and magnetic fundamental modes characterized with a weak quality factor that interfere with higher order modes associated with much higher quality factors. The interference is also affected by a non-resonant term that can shift the spectral position of the scattering minimum. The destructive interference leads to a Fano resonance associated with a minimum of the scattering efficiency. This modal analysis also clarifies the link between Fano resonances, anapole conditions and internal field enhancements. We used a modal expansion of the light intensity enhancement over the whole volume of the scatterer to identify the conditions that yield the maximum of the internal field enhancement in spherical Mie scatterers. The maximum of the field enhancement does not occur at the anapole condition but at the real part of the eigen-frequencies of high quality factor modes. The field inside the scatterer does not feature any minimum since it is only weakly affected by the fundamental modes with low quality factor and the modal expansion of the internal field coefficients does not possess a non-resonant term. Consequently, the internal field spectrum shows peaks associated with eigenmodes (except the fundamental low Q mode) and features almost pure Lorentzian shapes. This also means that the internal field associated with an anapole signature in the far field does not yield a particular field confinement inside the scatterer since it results simply from the excitation of a given mode in the dielectric particle. This modal analysis is carried out in the framework of a multipolar expansion of the field and could be further extended to non-spherical objects using appropriate multipolar decompositions.

Fig. 1 .

 1 Fig. 1. a) Position of the first poles of the dipolar electric and magnetic scattering coefficients in the complex frequency plane. b) Modulus of the associated residues of T (h) n :|r (e) 1,α | and |r (h) 1,α |. c) Modulus of the associated residues of Ω (h) n : |r (e) Ξ,1,α | and |r (h) Ξ,1,α |.

Fig. 2 .

 2 Fig. 2. (a,b) Multipolar partial scattering efficiencies versus internal field enhancement. Left scale, black full line : dipolar scattering efficiencies calculated with Eqs. 9, right scale, blue dashed line: averaged internal field intensity |E int,1 | 2 calculated with Eqs. 10 with respect to the parameter size z. Vertical dashed lines indicate the minimum of the scattering efficiency. (a) Electric dipole, (b) Magnetic dipole. (c) Scattering spectrum. Blue line: total scattering (with N=10); grey dashed line: electric dipole (N=1,e); green dashed line: magnetic dipole (N=1,h).

  h) feature resonances with asymmetric line shapes referred as Fano resonances [6, 53-59].

Fig. 3 . 2 r| 2 2 Mwith

 3222 Fig. 3. Analysis of the different terms in the modal expansion of the T (i) 1 (z)-matrix element.

.

  Here the pole α = 1 is characterized by a much larger imaginary part than the pole α = 2 but this weak quality factor is compensated by a larger residue r (e) α=1,1 leading to | r

Fig. 4 .

 4 Fig. 4. Spectra of the average internal field compared with the spectra of the scattering efficiencies for the magnetic and electric cases. (a) and (b) are taken from Fig. 3(b,c,e,f). Grey line: resonant contribution only in Eq. 6 with M = 100, red dotted line: resonant + non resonant contributions with M = 100. (c) and (d): Square of the norm of the internal electric field averaged over the volume of the particle < |E (i) int,n | 2 >, i = e for (c) and i = h for (d). Orange and green lines+shadow: modal expansion calculated with Eqs.10-11 when taking into account a single pole: α = 2 (orange), α = 3 (green) or when summing 2M poles (blue dashed line with M = 5). Full red line: full calculations. The vertical dashed lines indicate the spectral position of the anapole condition. (e) and (f): Difference of phase between the modes φ 1 -φ 2 in the electric case and φ 0 -φ 2 in the magnetic case with

Fig. 5 .

 5 Fig. 5. Internal field enhancement log|E (e) int,1 | 2 in the plane containing the incident electric field plotted inside the spherical scatterer at the electric dipole resonance, z = (p (i) 1,α ) = 1.065 (left), and at the anapole condition z = 1.160 (right). The field enhancement is stronger at the electric dipole resonance (left).

Fig. 6 .

 6 Fig. 6. Trajectories of the 5 first dipolar eigen-frequencies in the complex frequency plane when increasing the dielectric permittivity ε of the scatterer. Real and imaginary parts of the eigen-frequencies p (i) 1,α of the electric (i = e, a and b) and magnetic (i = h, c and d) dipolar modes n = 1 with respect to the dielectric permittivity of the scatterer embedded in air. The 5 first dipolar modes are plotted: α = [1; 5] for the electric case (a and b); α = [0; 4] for the magnetic case (c and d). The imaginary part of the first mode (α = 0 for i = (h) and α = 1 for i = (e)) is plotted in inset.

Acknowledgments

Research conducted within the context of the International Associated Laboratory ALPhFA: Associated Laboratory for Photonics between France and Australia. This work has been carried out thanks to the support of the A*MIDEX project (no. ANR-11-IDEX-0001-02) funded by the Investissements d'Avenir French Government program, managed by the French National Research Agency (ANR). We acknowledge support by the Australian Research Council.


G Galati

N Petrochilos
email: nicolas.petrochilos@gmail.com

E G Piracci
email: piracci@ing.uniroma2.it

Degarbling Mode S replies received in single channel stations with a digital incremental improvement

Multilateration (MLAT) and Automatic

Dependent Surveillance -Broadcast (ADS-B) systems are used in air traffic control to detect, locate and identify cooperating aircraft using signals emitted by airborne transponders and received by dedicated ground stations. In areas with a high traffic density, these stations may receive simultaneously several superimposed signals. Present operational systems use only one receiving channel connected to a nondirectional antenna. When the received replies are superimposed, i.e. "garbled", their detection and/or decoding are severely affected in nowadays equipment. The aim of this paper is to transform the single channel problem into a multiple channels problem in order to solve it using specific knowledge about the Mode S signals. In fact, the multiple channels problem is a typical signals separation problem applied to Mode S mixture for which several algorithms already exist. Our algorithm, named PASA, is based on the existing Projection Algorithm [13]

and can be easily implemented on existing receiving stations. The effectiveness of our method is demonstrated using real data collected from our experimental receiver.

Introduction 1.Operational frame

In the air traffic control (ATC) the modern cooperating surveillance systems are the well-known Automatic Dependent Surveillance -Broadcast systems (ADS-B) and multilateration (MLAT and wide area MLAT, WAM) systems. ADS-B, MLAT and WAM use the secondary surveillance radar (SSR) Mode A/C and Mode S signals [START_REF] Stevens | Secondary Surveillance Radar[END_REF]- [START_REF]DO-260A MOPS for 1090 MHz ES ADS-B and TIS-B[END_REF], for the detection, the tracking and the identification of aircraft transmitting messages at 1090 MHz with their standard transponder as well as of ground vehicles equipped with a 1090 MHz "non-transponder device".

The SSR operates by exchanging interrogations-replies between a ground station and the targets: the ground station emits, at 1030 MHz, interrogations addressed to the targets, which receive the interrogation, and then transmit to the ground station (at 1090 MHz) a reply containing the requested data (identity or flight level). The Mode S transponder transmits either short (64 μs) or long (120 μs) replies at the standard, nominal carrier frequency of 1090 MHz [START_REF] Stevens | Secondary Surveillance Radar[END_REF]- [START_REF]Annex 10 to the convention on international civil aviation, Aeronautical Telecommunications, Vol. IV Surveillance Radar and Collision Avoidance systems[END_REF], encoding the information by pulse position modulation (PPM). The transponders also transmit with a periodical timing, non-elicited signals called "squitters" in which several information about the aircraft are contained (i.e. identification code, 3D position, altitude, speed vector etc.). The "non-transponder devices" installed on a ground vehicle transmits squitters signals using the standard 1090 MHz Mode S signals, containing the identification code and the ground position of the vehicle. ADS-B, MLAT and WAM systems use the SSR and the Mode S 1090 MHz signals for surveillance purposes. In MLAT and WAM, the localization is independent of the information content of the received signals, and is achieved via triangulation by measuring the signal time difference of arrival (TDOA) at a number of ground receivers at different fixed locations. The accuracy of the target localization depends on the number of receiver stations and their layout on the operational area. The MLAT/WAM ground segment is composed by a network of distributed stations (at least four measurements are needed for a 3D localization) receiving the signals and estimating their time of arrival (by timestamp), and a central processor station collecting and processing the measurements to localize the emitters. The ADS-B uses the spontaneous Mode S messages (squitters) discarding all the other 1090 MHz traffic: by decoding the messages, the ADS-B receiver (part of the ADS-B station) extracts identification, position and velocity of the emitter. Finally, the ground ADS-B station sends the surveillance information to the users.

State of the art and research frame

Due to the fact that ADS-B, MLAT and WAM receiver stations use a single channel receiver with a non-directional antenna, the probability of receiving Mode S signals overlapped in time is not negligible and grows up with the traffic density increasing. When replies or squitters from different sources overlap, likely the received messages are corrupted and cannot be decoded, and the emitters (either airborne or vehicular) cannot be located and identified. As remarked in [START_REF]1090 MHz Capacity studyfinal report[END_REF] the interfering signals rate for a ground receiver placed in the Core Europe (Brussels) is likely from 25000 msg s -1 up to 95000 msg s -1 depending on different operating conditions (time of the year, number of operational interrogators etc.). With these values of the FRUIT (False Replies Unsynchronized In Time) rate, the detection probability, depending on the target range, can be less than 70% for target ranges greater than 20 NM [START_REF]1090 MHz Capacity studyfinal report[END_REF].

In order to solve the problem of overlapping, source separation using array processing has been deeply investigated. There is a large literature also on the specific applications to SSR signals [START_REF] Chaumette | An ICA-based technique for radiating sources estimation: application to airport surveillance[END_REF]- [START_REF] Zhou | Improved blind separation algorithm for overlapping secondary surveillance radar replies[END_REF]. In [START_REF] Chaumette | An ICA-based technique for radiating sources estimation: application to airport surveillance[END_REF] and [START_REF] Tol | SSR reply separation using array signal processing methods[END_REF] an adaptation of existing array processing algorithms was proposed for SSR signals. In [START_REF] Petrochilos | Algorithms for separation of SSR replies[END_REF]- [START_REF] Petrochilos | Algebraic algorithms to separate overlapping secondary surveillance radar replies[END_REF] the authors presented several novel algebraic methods based on the sources property. In [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF]- [START_REF] Petrochilos | Array processing of SSR signals in the multilateration context, a decade survey[END_REF] different methods exploiting the time-domain sparsity property of the sources are proposed (paragraph 2.2 contains a detailed description). Finally, in [START_REF] Zhou | Improved blind separation algorithm for overlapping secondary surveillance radar replies[END_REF] an improvement of previous algebraic methods is proposed.

In the last decade different single channel Blind source Separation (BSS) methods were proposed but none of them address specifically the SSR signals. The single-channel ICA proposed in [START_REF] Davies | Source separation using single channel ICA[END_REF] needs independent sources with disjoint spectral support. Some methods are useful for sound and speech signals [START_REF] Jang | Single-Channel Signal separation using time-domain basis functions[END_REF]- [START_REF] Hopgood | Single Channel Nonstationary Stochastic Signal Separation Using Linear Time-Varying Filters[END_REF], other are suitable in the radar and communications context [START_REF] Warner | Single channel blind signal separation of filtered MPSK signals[END_REF]- [START_REF] Rui | Blind separability for single-channel mixtures in satellite channels[END_REF] and in the bio-medical field [START_REF] Mijović | Source separation from single-channel recordings by combining empirical-mode decomposition and independent component analysis[END_REF]- [START_REF] Gao | Blind separation of fetal ECG from single mixture using SVD and ICA[END_REF]. A number of papers in [START_REF]Source Separation and Applications[END_REF] provide a wide overview of many aspects and recent advances on source separation and applications.

In this frame, we present here a single channel separation method for SSR signals suited to any single receiving station with one non-directional antenna. The algorithm is based on a data adaptation needed to use the Projection Algorithm, PA [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF], and leads to a simpler and more effective method, the Projection Algorithm -Single Antenna, PASA. As explained before, the aim is to mitigate the garbling problem and to improve the channel capacity even in high traffic density We first describe the single antenna data model. Next, we present its transformation to obtain an equivalent array antenna data model. The link between this model and the signals separation method is described, and after a short reminder of the PA [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF] algorithm, the PASA is presented.

The Data Model

We denote scalars by italic lowercase or uppercase letters (e.g. a, A), vectors by lower case boldface letters (e.g. x, x[n]), vector elements with lower case letter with in-bracket position reference (i.e.

x[n]), and matrices by uppercase boldface letters (e.g. X).

There are two formats for the 1090 MHz signals: the conventional mode (including mode A/C and military signals) and the Mode S. This work is focused on Mode S signals (reply, squitter) and on their mixtures. A squitter is a randomly self-triggered reply that follow Mode S format and contains 56 (short) or 112 (long) binary symbols bn such as a reply. The symbol period is 1 μs, and each symbol is made up by two 0.5 μs chips. Figure 2 shows a schematic view of the time intervals for the preamble and the data block, and also shows the concept of superimposed signals. with a total length of 128 (short reply) or 240 (long reply) elements (64 or 120 µs).

The Mode S reply/squitter has the form:

𝐛(𝑡) = ∑ 𝐛[𝑛] 𝒑(𝑡 -𝑛𝑇) 127/239 𝑛=0 ( 1 )
where p(t) is a rectangular pulse of width T= 0.5 μs. This signal is up-converted to the nominal transmission frequency fo= 1090 MHz, with a ± 1 MHz tolerance as set by ICAO standards [START_REF]Annex 10 to the convention on international civil aviation, Aeronautical Telecommunications, Vol. IV Surveillance Radar and Collision Avoidance systems[END_REF]. Due to this tolerance, after down-conversion to base band, a residual frequency fr remains, adding a progressive phase rotation to the transmitted symbols. The received base band signal is:

𝐬 [𝑘] = 𝐛 [𝑘] exp(𝑗 2 𝜋 𝑘 𝑓 𝑟 𝑇 𝑠 ) = 𝐛[𝑘] 𝜑 𝑘 , 𝑘 = 1, 2, … , 𝑁 ( 2 )
where 𝜑 = 𝑒𝑥𝑝(𝑗2𝜋𝑓 𝑟 𝑇 𝑠 ) is the elementary phase shift due to the residual carrier frequency, TS the sampling period, and j is the imaginary unit.

Considering d signals emitted by independent sources and impinging on a single antenna, the general data model for the N samples of the received signal is:

𝐱 = 𝐦 𝐒 + 𝐧 ( 3
)
where 𝐱 is the [1×N] data vector, 𝐒 is the [d×N] sources matrix:

𝐒 = [ s 1 [1] s 1 [2] … s 1 [N] ⋮ ⋮ ⋮ s 𝑑 [1] s 𝑑 [2] … s 𝑑 [N] ] ( 4 )
𝐦 is the mixing vector, (in this case, the sources mixing is a weighted sum) i.e. a [1×d] vector containing the complex gain of the antenna in the directions of the sources; finally, 𝐧 is the [1×N] additive noise vector. In the present discussion we assume the two sources case, i.e. d=2, and N to be large enough so that the vector x contains both signals, with the vectors s1[n] and s2[n] zeropadded at their trailing and leading edges, respectively. The aim of the remaining of this paper is to describe and evaluate the method to extract the sources signals s1[n] and s2[n] from the observed vector x.

The sources (i.e. the signals) share the frequency band, are not synchronised, and have the same encoding; therefore, most of existing single-channel algorithms in the literature are not suited to the Mode S replies separation problem in the MLAT-ADS-B context. In this frame, if the leading has a free interference preamble it can be detected and decoded, but the trailing reply does suffer much more, and cannot be detected (its preamble is garbled with the data block of the leading signal), hence the need for a different solution.

Transforming the single channel problem into a multiple channel problem

Let us consider the signal emitted by one source as a stream, i.e. a vector s[n]. Such a single source is described by the vector s which differs from the received vector x by the (complex) antenna gain and by the additive (receiver) noise. In this paragraph the noiseless case is described, therefore, being immaterial the constant antenna gain, the vector x, s (resp. the matrices X, S) are equivalent, and may be exchanged each other (when the noise is neglected) in the following formulas (see also where l is the rounded part of N/m, and contains in the first column the first m samples of the vector s, in the second column the second m samples of s and so forth. Let s be the source stream:

𝐬 = [ s[1] s[2] ⋮ ⋮ s[𝑁]]
( 5.a ) and let us construct the derived data matrix SA:

𝐒 𝐀 = [ s[1] ⋮ s[𝑚 + 1] ⋮ ⋮ s[𝑚] ⋮ s[2𝑚] ⋯ s[𝑙 • 𝑚 -𝑚 + 1] ⋮ ⋱ ⋯ ⋮ s[𝑙 • 𝑚] ] ( 5.b )
Note that we lose Nl•m samples in the transformation, which is smaller than m by the definition of l. From now on, we call m the reshaping factor. With the correct reshaping factor, the matrix S has a low-rank factorization; let us confirm that fact empirically, and then see why. Figure 3 shows a case of a received signal in field experiments sampled at 50 Megasamples per seconds that has been transformed using a reshaping factor m from 1 to 125. For each value of m, the matrix SA was constructed, and its singular values were estimated. Then the eigenvalues of SA were normalized by division with the value of the greatest one. Figure 3 shows these eigenvector grey-coded by amplitude; the bright white corresponds to 1, the maximum, whereas black represents values below or equal to 0.01 (-40 dB). The first and the last column are the same, up to noise, and the second column is nearly orthogonal to the first and to the fourth, and finally the third column contains only noise. Even if the reply and the reshaping are not synchronised together, their time lengths are commensurate, therefore the cut induced by the reshaping always happen at the same position on the chip. Considering figure 4 and

denoting 𝐩 u = [ 0 0 0.25 1 1 ]
, the vector from slot 1, and

𝐩 d = [ 1 1 0.75 0 0 ]
, the vector from slot 2, all subsequent pulses can be described as a combination of both vectors, neglecting the noise. In this case the noiseless matrix SA is written as:

𝐒 𝐀 = [ 𝐩 𝐮 𝐩 𝐝 ]. [ 1 0 0 1 0 1 0 0 ] = 𝐏 𝐀 𝐒 𝐀 ̃ (6)
Although pu and pd are not orthogonal, the matrix PA is full rank, so is 𝐒 𝐀 ̃. This means that, adopting m1= ½ fs Tb, the rank of the noiseless matrix is two, and for the noisy version, two large Tb Tb singular values (above noise threshold) are expected. Note that the received reply does not have to be synchronised, as every pulse will be sliced at the same position. There are two non-trivial vectors: first half, pu, and second half, pd. Similarly, when considering values of reshaping factor integer multiples of m1, for each additional length m1 there are two more options: first half up or second half up; thus for mk=k m1, the rank is 2•k.

Revisiting the case on Figure 4 with a non-zero residual frequency, we note that pu and pd are changed to reflect the phase change:

𝐩 u = [ 0 0 0.25𝜑 2 𝜑 3 𝜑 4 ]
and

𝐩 d = [ 1 𝜑 0.75𝜑 2 0 0 ]
. So equation ( 6) is modified to:

𝐒 𝐀 = [ 𝐩 𝐮 𝐩 𝐝 ]. [ 1 0 0 𝜑 15 0 𝜑 5 0 0 ] = 𝐏 𝐀 𝐒 𝐀 ̃ (7)
Note that the px's have a residual frequency shift included in their model, therefore two different replies with different residual frequency will have a set of px's that are independent although not orthogonal.

Link with Blind Source Separation problems

Now that we have a transformation from the initial single channel problem into a low rank multiple channels problem, let see how to transform equation 3. In a sum form, we get:

𝐱 = ∑ m[i] • 𝐬 𝑖 𝑑 𝑖=1 + 𝐧 (8)
by linearity of the reshaping transformation, we have:

𝐗 = ∑ m[i] • 𝐒 𝑖 𝑑 𝑖=1 + 𝐍 (9)
where the obtained X is constructed in the same fashion as SA in equation (5.b):

𝐗 = [ x[1] x[𝑚 + 1] ⋮ x[𝑚] ⋮ x[2𝑚] ⋯ x[𝑙 • 𝑚 -𝑚 + 1] ⋯ ⋮ x[𝑙 • 𝑚] ] (10)
from equation ( 9) and ( 7), X is then equal to:

𝐗 = [m[1] • 𝐏 𝟏 ⋯ m[d] • 𝐏 𝐝 ] • [ 𝐒 𝟏 ⋮ 𝐒 𝟐 ̃] + 𝐍 (11)
and in reduced form:

𝐗 = 𝐌 𝐒 ̃+ 𝐍 (12)
M is a [mk × 2kd] matrix, and 𝐒 ̃ is [2kd × l] matrix, so by definition M is a tall matrix, and 𝐒 ̃ a wide matrix, which ranks are at most 2kd. Thought the authors do not have a formal proof of being full rank, they never observed a case where either M or 𝐒 ̃ was rank deficient, albeit a bad conditioning number on M was observed. We acknowledge that the case may happen, for instance if fi = fj for i ≠ j. We have now a situation that is equivalent to the instantaneous over-determined mixture that is solved with Blind Source Separation. The algorithm AZCMA [START_REF] Petrochilos | Algebraic algorithms to separate overlapping secondary surveillance radar replies[END_REF] would not be a good candidate since the remaining frequency can be found twice for every initial source. Given that the mode S is pseudo-Gaussian by nature [START_REF] Petrochilos | A zero-cumulant random variable and its applications[END_REF], the usage of ICA or JADE [START_REF] Chaumette | An ICA-based technique for radiating sources estimation: application to airport surveillance[END_REF] seems risky. We did not check if the new sources 𝐒 𝐀 ̃ still comply with the Manchester encoding property [START_REF] Petrochilos | Algebraic algorithms to separate overlapping secondary surveillance radar replies[END_REF], therefore we cannot use MDA. Moreover, with any BSS algorithm, we would have to find in the set of estimated sources which 2k ones are paired to recreate the original d replies. The PA algorithm [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF] overcomes this limitation as it will directly estimate the 2k derived sources together, the pairing is then trivial. In the next paragraph, we remind the PA algorithm, while in the last one, we presents its application to the reshaped equation [START_REF] Petrochilos | Algebraic algorithms to separate overlapping secondary surveillance radar replies[END_REF].

An overview of the Projection Algorithm (PA) [13]

In [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF], two algorithms were proposed to separate from an antenna array overlapping replies: either using the Mode S protocol or the older Mode A/C protocol. For the simplest case of two partially overlapping Mode S replies, the more robust Projection Algorithm (PA) has been demonstrated.

The Extended Projection Algorithm (EPA) considers all other cases of mixture and hence is less robust. One particular case, i.e. a totally overlapping pair of Mode S replies, albeit with a low probability of appearance in the real world, is cumbersome as we have to resort to another algorithm: the Manchester Decoding Algorithm (MDA), [START_REF] Petrochilos | Improved MDA, a case for de-garbling SSR Mode S replies[END_REF].

The projection algorithms need the output of an antenna array. At the time sample n, the i-th element output is xi[n]. After collecting N samples and stacking the samples from each of m array elements into m vectors x[n], the resulting [m×N] data matrix is:

𝐗 = [ x 1 [1] x 1 [2] … x 1 [N] ⋮ ⋮ ⋮ x 𝑚 [1] x 𝑚 [2] … x 𝑚 [N] ] (13)
Let us consider now the reception of signals from d different sources by an m-elements antenna array. The data model is expressed as:

𝐗 = 𝐌 𝐒 + 𝐍 ( 14
)
where M is the [ md  ] mixing matrix that contains the array signatures and the complex gains of the sources, 𝐒 is the [d×N] sources signal matrix (equation ( 4)), and N is the [m×N] noise matrix. In the case considered here, i.e. d ≤ m, the signals separation is a typical blind source separation (BSS)

problem. The noiseless matrix X has a low rank factorization, with a rank of d.

Table 1 shows the steps of the Projection Algorithm [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF]: amplitude with a noise threshold, the presence and the number of sources is estimated as a function of the time. In the following (•) (1) indicates the matrix (or the vector) composed by the columns (or the elements), related to the time interval [t1 : t2] and the same for (•) (2) and the time interval [t3 :

t4].
With this notation and the data model in equation ( 14), the following relations hold:

𝐗 (𝟏) = 𝐌𝐒 (𝟏) + 𝐍 (𝟏) (15.a) 𝐗 (𝟐) = 𝐌𝐒 (𝟐) + 𝐍 (𝟐) (15.b)
By definition, S (1) contains only a part of the first source s1, the row s2 is zero-padded from the beginning till t2, S (2) contains only a part s2 and s1 is zero padded from t3 to t4. Therefore, it is possible to simplify the equations ( 15) in: Note that X (1) and X (2) are rank-one matrices in the noiseless case. The ensuing step of the algorithm consist in the estimation of the [m × 2] mixing matrix: 𝐌 ̂= [𝐦 ̂1 𝐦 ̂2]. The vectors 𝐦 ̂1 and 𝐦 ̂2 are obtained by a SVD of X (1) and X (2) respectively and correspond to the largest singular values. The final step is the derivation of the vectors w1 and w2 by the rows of the Moore-Penrose pseudoinverse of 𝐌 ̂. They works as spatial filters, applied to X to estimate the sources vectors si. If the angle between the two signals impinging on the array tends to zero, 𝐌 ̂ becomes ill conditioned, the noise at the output of the beamformers increase dramatically, and the signal to noise ratio (SNR)

𝐗 (𝟏) =
decreases. The limitation on the angle of arrival is shown in the Appendix I of [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF], to which we kindly refer the reader interested to a more detailed description including the EPA algorithm.

The PASA algorithm

The concept behind the Projection Algorithm Single Antenna (PASA) consists in transforming the received single channel data stream x into a matrix X by reshuffling its elements, i.e. reshaping the vector. The low rank matrix X can then be processed with a modified version of the PA [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF]. The application of PASA is in this paper limited to an overlap of two SSR Mode S signals. Table 2 summarize the PASA: where l is the rounded part of N/m. This matrix contains in its first column the first m samples of the vector x, in its second column the second m samples of x and so forth:

𝐱 = [ x[1] x[2] ⋯ ⋯ ⋯ x[𝑁]] received data stream (17.a) 𝐗 = [ x[1] ⋮ x[𝑚 + 1] ⋮ ⋮ x[𝑚] ⋮ x[2𝑚] ⋯ x[𝑚 + 1] ⋮ ⋱ ⋯ ⋮ x[2𝑚] ] derived data matrix (17.b)
Defining a data model as in equation ( 14), the matrix M is not created by array response but, rather, by source's residual frequency (fr) difference, see equation ( 2). The goal is to estimate M to derive two subspaces on which to project X to recover the signal vectors s1 and s2. The projectors, obtained by the pseudo-inverse of 𝐌 ̂, can be called the signature vectors. The matrices S1 and S2, resulting from the application of the signature vectors onto X, shall contain the two sources, respectively. The estimated signal vectors 𝐬 ̂1 and 𝐬 ̂2 are then obtained by an inverse reshaping of S1 and S2

respectively.

As in the PA, the first step of PASA consists in the estimation of the sources arrival timing, to define the two subset's of the columns of X, in which only one source is present. For this task, the m value, that is the number of rows of X, is a fundamental parameter determining the rank of X and the conditioning number of 𝐌 ̂. Using a too small value (like 2 or 4) the conditioning number of 𝐌 ̂ results to be too high, with the pseudo-inverse being not robust and the noise contribution after the projection too large. Moreover, with m incommensurable with the number of samples contained in a bit time, the rank of 𝐌 ̂ is maximum and the detection of the number of sources is not possible using a whiteness test like for PA. If m is set equal to the number of samples contained in half bit time, (i.e. m= ½ fs Tb, where fs is the sampling frequency and Tb is the bit duration, i.e. 1 µs), each column of X contains either a pulse signal or noise only, or, depending on the synchronization between the observed vector x and the starting sample of the signal, a frame composed by a partial pulse and noise, see figure 4.

So far two signature vectors are required for each source. This means that, adopting m= ½ fs Tb, performing a whiteness test on X to estimate the sources number and timing, two large singular values (above threshold) are expected when one source is present, and four large singular values when two sources overlap, thus making easy the estimation of the ti's. As for PA, the whiteness test is performed on a sliding window of 4 Tb, i.e. 4 µs. Once the timing estimation is done, 𝐗 (1) and 𝐗 (2) , containing only the first and only the second source respectively, are extracted and used to estimate the signature vectors (two vectors for each source). The matrix 𝐌 ̂ is composed with the vectors corresponding to the first and to the second largest singular values of the decomposition of 𝐗 (1) and 𝐗 (2) . The next step consisting in the deriving of the projection matrix: 𝐖 = 𝐌 † = [𝐌 ̂𝐇𝐌 ̂]-1 𝐌 ̂𝐇. The application of W onto matrix X gives the estimation of the sources matrix 𝐒 ̂ from which the separated sources are obtained by an inverse reshaping from the S columns into a row vector. The formulation of the algorithm steps are summarized in the Appendix. As a final remark, we remind that PASA aims at (and is able to) separate exactly two sources. This is a design choice based on the probability to have two Mode S replies overlapping, and considering the much lower probability of overlapping of three or more mode S replies. In fact, considering for example a high Mode S FRUIT rate environment (20k messages/s), the probability of receiving exactly two interfered messages (i.e. 3 %) is ten times the probability to receive three or more interfered messages (i.e. 0.3 %) [START_REF] Piracci | MHz ATC systems capacity improvement[END_REF].

The PASA method: experimental results

In this section the practical implementation of the PASA algorithm is presented. A processing architecture to perform all the algorithm steps is described, and a case of study is shown. Finally the performance of the method is discussed, analysing the success rate with overlapping signals. To evaluate the effectiveness of the algorithm we used Mode S signals recorded by an ad-hoc experimental receiver [START_REF] Galati | Transponder Data Recorder: final implementation and first results[END_REF]. The device is a multichannel receiver (four linear channels and one logarithmic channel), equipped with an array antenna. The logarithmic channel is mostly used as a signal detector, and the four linear channel are sampled at typical rate of 100 Msamples/s directly at the intermediate frequency (21.5 MHz). This receiver was used in single-channel mode to perform the experimental analysis as described in the following paragraphs.

Practical implementation and case of study

It has been previously assessed that, in order to detect the number of sources and their timing based on the reshaped matrix X, it is necessary to set m equal to the number of samples in one half bit time (i.e. 0.5 µs). But, using integer multiples of the value ½ fs Tb, like m=fs Tb or m=2 fs Tb, we have observed better performance for the sources sub-spaces estimation, thanks to a smaller conditioning number for 𝐌 ̂. In order to take advantage of both configurations, a parallel type of processing has been devised and tested. This is shown in Figure 6: the processing performs two data reshaping's in parallel. The matrix obtained by a reshaping using m =fs Tb is suitable for the projection, as it gives better results. The data reshaped with m = ½ fs Tb are the best input for the sources timing estimation. The sources timing is achieved by comparing the singular values of the matrix X with a CFAR threshold [START_REF] Galati | An advanced CFAR technique for SSR reply detection[END_REF]. Once estimated the sources timing, it is possible to extract X (1) and X (2) from X, they contain only the first and the second source respectively. Then, the vectors mi are computed by the SVD of X (1) and X (2) and the matrix X is projected over the pseudoinverse M † of 𝐌 ̂, (𝐌 ̂=[m1 ... mn]). In practice, the matrix X is multiplied by M † . Finally both sources are recovered with an inverse reshaping of the result. confidence bit" being contained in the Appendix I of [START_REF]DO-260A MOPS for 1090 MHz ES ADS-B and TIS-B[END_REF]. However the gain with the respect to a classical decoder is interesting since without PASA it is not possible to detect and decode the trailing reply, and the leading is decoded with several wrong bits. The PASA algorithm has worked thanks to the frequency shift between the sources (0.19 MHz). The resulting conditioning number of the matrix 𝐌 ̂ is equal to 9.5. application of PASA it is possible to correctly detect the preamble of the trailing reply. The preamble detection is performed with the standard routine in compliance to the recommendation of appendix I of [START_REF]ED-117 MOPS for Mode S Multilateration Systems for Use in Advanced Surface Movement Guidance and Control Systems (A-SMGCS)[END_REF]. Since we did not have 'a priori' information about the messages code and the transponders unique address, it was not possible to analyze the decoding errors. Therefore, for each case we have analyzed: the time delay (td), the power ratio between the signals (PR), the frequency shift (fsh) and the conditioning of the matrix 𝐌 ̂ (cond). A post processing analysis has been done to investigate any correlation between the success of the method and these characteristics.

In table 3

Table 3: Characteristics of the superimposed signals involved in PASA evaluation

The application of the PASA method permitted the trailing reply detection in 63 cases, i.e. with a success rate of 60 %, where a failure case is defined when downstream the application of the PASA method it is not yet possible to detect the trailing reply, because of the missed preamble detection.

In order to analyze how the performance are conditioned by the signals characteristics, the parameters of table 1 were evaluated separately on success and on failure cases. It seems that the effectiveness of the PASA is not strongly affected by the relative time delay of the interfering signals, although it was not practically possible to test the method with very small time delays cases. The PASA seems to behave better at low power ratio between replies, order of 0 dB, anyway its behaviour is not significantly correlated with power ratio. The effectiveness of the method seems to be significantly dependent on the frequency relative shift and the conditioning number, as expected. An analysis was done on the mixed signals with an high relative frequency shift (more than 0.15 MHz) but identified as failure cases: they are characterized either by a high power difference between the two sources (the trailing signals are over 10 dB above the leading one), or a poor time delay (less than 10 µs). A typical situation generating failure is the following: if there is a high power difference (more than 3 dB), the projection of the data matrix onto the subspace of the smaller source is not efficient, and the components of the most powerful signal persist on the reconstructed signal that is not detectable. The Mode S receiver in compliance with the requirements of [START_REF]DO-260A MOPS for 1090 MHz ES ADS-B and TIS-B[END_REF] is able to handle this situation: if during the data-block decoding a preamble more than 3 dB greater than the earlier is validated, then the earlier signal is rejected so that the data-block demodulation of the new signal can proceed. By design, the robustness of the PPM permits a free-error decoding (also thanks to the cyclic redundancy check) in the presence of a lowamplitude interference. Another identified cause arises when overlapped signals have a too small time difference, leading to two main complications: i) if the relative time delay is smaller than the window length of the whiteness test (4 µs) it is not possible to recognize the sources timing and the algorithm fails; ii) if the relative time delay is low, but sufficient to estimate the sources timing, there are a few columns of the matrix X on which to estimate the signals signature vectors, then the projection matrix does not permit a good estimate of the two sources sub-spaces and the output SNR is too low to recognize the signals. Summing up: the higher the relative frequency shift between the mixing sources, the smaller the M conditioning and better the PASA performance. A limiting factor is represented by other parameters: when the relative time delay is less than 10 µs, and/or the signal level difference is above 3 dB, the effectiveness of the method is generally insufficient.

The experimental results demonstrate that PASA is useful as signals separation algorithm between replies from similar ranges, which should be a typical scenario of an approach area or of a large airport.

Concluding remarks

In this paper we presented a method called PASA to separate overlapped Mode S signals using a single channel receiver, in order to enhance the Mode S signal-based surveillance of a busy terminal area or a big airport. The method is partly based on the array processing algorithm named PA, presented in [START_REF] Petrochilos | Separation of SSR signals by array processing in multilateration systems[END_REF]. A measurement and test campaign with a dedicated receiver allowed us to test the method with live cases and to analyze its effectiveness and its limits. It has been noticed that PASA permits a significant improvement of the success rate, while there isn't any other reasonably option to solve the problem of signal overlapping. Concerning the experimented failure cases, it was The same for the second source.

Last step is the easy reconstruction of si[n] from 𝐒 ̂𝑖 = [ 𝐬 1 𝐬 2 … 𝐬 𝐿 ], (i = 1, 2).

  situations. This paper is organized as follows. Paragraph 2.1 shows the data model; paragraph 2.2 explains the data adaptation to transform the single channel problem into a multi-channels problem and shows the related data properties. Paragraph 2.3 links with the multi-channel separation methods and explains why the PA principle was chosen to solve the problem. Paragraph 2.4 describes PA in more detail. Paragraph 2.5 shows the practical implementation of PASA. Finally Paragraph 3 presents a case-study and the related performance analysis with recorded live signals. A typical Mode S station for ADS-B and MLAT/WAM is composed by an antenna, an analog front-end and a digital section for signals detection and decoding. In figure 1 a general block diagram of a Mode S receiver is shown; the dotted-line box contains the extra logical functions for the PASA method.

Figure 1 :

 1 Figure 1: 1090 MHz receiver scheme for the PASA (Projection Algorithm Single Antenna)

Figure 2 :

 2 Figure 2: Schematic view of two Mode S replies with a time delay of td = t2 -t1

the

  Appendix). A simple starting point is to reshape this vector s[n] of length N into a matrix SA by stacking m elements of s at a time (being m<<N). The dimensions of the matrix SA are [m × l]

Figure 3 :

 3 Figure 3: Eigenvalues for different reshaping factor, real short mode S reply, N=3400, Fs = 50 MHz, SNR = 24.8 dB According to Figure 3, for most of the reshaping factor values, the matrix SA is not full rank, but

Figure 4 :

 4 Figure 4: How to perform the transformation of x[n] into X

Table 1 :

 1 timing detection, interference-free sub-matrices X1 e X2 extraction 3 Mixing matrix M estimation, (m1 estimated using X1 and m2 estimated using X2) 4 Beamforming matrix computation W = M † 5 Sources estimation by S = W•X Summary of PA Let us consider the case of two (d=2) superimposed signals received with an m-elements array antenna, with m ≥ 2. Let (t1, t2) and (t3, t4) the beginning and ending time of the leading and trailing signal respectively, see fig. 2. Three time slots are defined: [t1 : t2] in which only the leading signal is present, [t2 : t3] with signals superimposing and [t3 : t4] in which only the trailing signal is present.

Figure 5

 5 Figure5shows the amplitude of a signal received in field experiments, being composed by two

  𝐦 𝟏 𝐬 𝟏 (𝟏) + 𝐍 (𝟏) (16.a) 𝐗 (𝟐) = 𝐦 𝟐 𝐬 𝟐 (𝟐) + 𝐍 (𝟐) (16.b) where m1 and m2 are the columns of the [m × 2] mixing matrix M, and si, i=1, 2, are the source vectors, i.e. the rows of the [2× N] matrix S.

Figure 5 :

 5 Figure 5: Amplitude of received signal composed by two Mode S overlapping sources

step 1 Table 2 :

 12 Data vector x acquisition and reshaping to obtain data matrix X (using m = ½ fs Tb) 2 Sources number and timing detection, sub-matrices X1 e X2 extraction 3 Mixing matrix M estimation, (m1 estimated using X1 and m2 estimated using X2) 4 Signature matrix computation W = M † 5 Sources data matrix estimation S = W•X 6 Matrix S inverse-reshaping to recover the separated sources s1 and s2 Summary of PASA Let us consider the superimposition of two signals s1[n] and s2[n] transmitted from two independent sources. The received data stream is collected in a row vector x[n] of length N, and the columns of the matrix X are obtained taking m elements of x at a time (m<<N). The dimensions of the matrix X are [m × l]

Figure 6 :

 6 Figure 6: Proposed architecture for source's timing and extraction

Figure 7 :Figure 8

 78 Figure 7: Singular values of X as function of time, by a sliding window over the columns of X

Figure 8 : 5 3. 2

 852 Figure 8: Results of PASA application, garbled signal shown in figure 5

  the estimated characteristics on the involved signals are shown: the min and max values, the mean and the standard deviation of each parameter have been evaluated. High values of the standard deviation related to the mean values show that the characteristics of the considered signals are rather randomly distributed in their own ranges.

  number ( n ) relative frequency shift (fsh) relative time delay (td) signal power ratio (PR) In Figure 9 the histograms and the mean and standard deviation values of the estimated parameters are shown.

Figure 9 :

 9 Figure 9: Signal characteristics evaluated on the success and failure groups

  verified that they are due to some characteristics of the superimposed signals, like frequency shift, time delay and amplitude levels ratio. A direction for future research on this topic are tests with a large databases to analyze the correlation between the failures and the signals parameters, to finalize future algorithmic improvements.Appendix: formulation of the algorithm stepsLet x[n] the vector of length N of the sampled data stream with two partially overlapped Mode S signals. For the sake of simplicity, let N be large enough to contain both signals and to be an integer multiple of m= ½ fs Tb. We collect from x[n] the following vectors:𝐱 𝑛 = [ x[𝑚 • 𝑛 -𝑚 + 1] ⋮ x[𝑚 • 𝑛] ] of dimension m,with n=1, …, L, and L=N/m. Let 𝐗 ≜ [ 𝐱 1 𝐱 2 … 𝐱 𝐿 ], a [m×L] data matrix composed by vectors xn.The sampled signal x[n] contains the sum of the two sources (s1 and s2) samples and the noise (n):x[n]= s1[n] + s2[n] + n[n].Since the transformation is linear, we have:X = S1 + S2 + N,where S1, S2, and N are [m×L] matrices obtained by s1[n], s2[n] and n[n] respectively, with the same method as X.In the noiseless case S1 and S2 are low rank matrices (i.e. with rank 2 or 3 depending by the initial point of the stacking from s1 and s2).Considering a rank factorization, S1 and S2 may be written as:S1 = M1 V1, S2 = M2 V2,where Mi is a [m×r] matrix and Vi is a [r×L] and r is the rank of Si.Using the introduced factorization we obtain X = M1 V1 + M2 V2 + N, which can be written as:𝐗 = [𝐌 𝟏 |𝐌 𝟐 ] [ 𝐕 𝟏 -𝐕 𝟐 ] + 𝐍 ≜ 𝐌 𝐕 + 𝐍.The idea of PASA is the application of PA to this data model, using the pseudo-inverse of M to extract the sources: 𝐖 = 𝐌 † = [𝐌 𝐇 𝐌] -1 𝐌 𝐇 . The used projectors are: 𝐖 𝟏 = [𝐌 𝟏 𝟎 𝐦×𝐝 𝟏 ] 𝐖, 𝐖 𝟐 = [𝟎 𝐦×𝐝 𝟏 𝐌 𝟐 ] 𝐖, where d1 and d2 are the number of signature vectors used for each source, respectively. Using the projectors on X the sources are extracted. For the first source we have: 𝐒 ̂1 = 𝐖 1 𝐗 = 𝐖 1 𝐌𝐕 + 𝐖 1 𝐍 = [𝐌 𝟏 𝟎 𝐦×𝐝 𝟏 ] 𝐖𝐌𝐕 + 𝐖 1 𝐍 = [𝐌 𝟏 𝟎 𝐦×𝐝 𝟏 ] 𝐌 † 𝐌𝐕 + 𝐖 1 𝐍 = = [𝐌 𝟏 𝟎 𝐦×𝐝 𝟏 ] [ 𝐕 𝟏 -𝐕 𝟐 ] + 𝐖 1 𝐍 = 𝐌 1 𝐕 𝟏 + 𝐖 1 𝐍 = 𝐒 𝟏 + 𝐖 1 𝐍

Acknowledgement

This work has been partially funded by the Company Selex-ES in the frame of Cooperation Agreement with the University of Rome "Tor Vergata", D.R. 1848/2006, Scientific Area ING/INF-03.