E Fléchard

Liénard

Mustapha Said

Hamid Berriche

Alaa Abd El-Hady

A L Shahrani

Gilles Ban

Xavier Flechard

Etienne Lienard

Simulation of ion behavior in an open three-dimensional Paul trap using a power series method

Keywords: Paul trap, Higher order potentials, Simulations, Power series method, Buffer gas cooling, Dipolar excitation

Simulation of ion behavior in an open three-dimensional

Paul trap using a power series method

Introduction

Paul traps are widely used in different experiments [1] [START_REF] Knoop | Lifetime, collisional-quenching and jmixing measurements of the metastable 3D levels of Ca+[END_REF][3] [START_REF] Ban | Measurement of the β -ν angular correlation parameter in nuclear β-decay[END_REF]. They are made of few electrodes and the trapping is obtained by the application of a DC and an AC voltage in the radiofrequency domain. Simulations of the trapped ions dynamics are crucial in some experiments such as the experiment of the "Laboratoire de Physique Corpusculaire de Caen (LPC Caen)", where an open three-dimensional Paul trap (the LPCTrap) is used for the determination of the beta-neutrino angular correlation parameter [START_REF] Flechard | Paul trapping of radioactive 6 He + ions and direct observation of their β-decay[END_REF]. The error on this parameter depends on the spatial and velocity distributions of the trapped ions. These last have been determined using commercial softwares.

In real radiofrequency Paul traps such as the LPCTrap, the general expression of the potential in spherical coordinates (ρ, θ, ϕ) is given by [START_REF] Gudjons | Influence of anharmonicities of a Paul trap potential on the motion of stored ions[END_REF] [START_REF] Wang | The non-linear resonance ion trap. Part 2. A general theoretical analysis[END_REF] Φ(ρ, θ) = (U DC + V AC cos(Ωt))

∞ n=0 b n ρ n P n (cosθ) (1)
U DC and V AC cos(Ωt) are the DC and AC voltages applied to the trap. b n are constants and P n are the Legendre polynomials of order n. Using the cylindrical coordinates (r, z, φ), with z along the trap's axis, ρ is replaced by (z 2 + r 2 ) and cosθ = z ρ . That transforms expression (1) into [START_REF] Wang | The non-linear resonance ion trap. Part 2. A general theoretical analysis[END_REF][8] Φ(r, z) = (U DC + V AC cos(Ωt))

∞ n=0 C n H n (r, z) (2)
C n are constants and H n (r, z) are functions of r and z.

For the ideal quadrupole trap all the terms of the sum except those for n = 0 and n = 2 vanish. The equation of motion for a direction u = x, y or z is the Mathieu equation. For every direction the analytical solution exists independently of the 2 other directions [START_REF] Ruby | Application of the Mathieu equation[END_REF]. When trapped, an ion performs oscillations at the frequencies

ν u = (n ± β u 2 )f (3)
f = Ω 2π is the frequency of the trapping field, n an integer and β u is given by a recursion formula [START_REF] March | An Introduction to Quadrupole Ion Trap Mass Spectrometry[END_REF]. ν x and ν y are equal and are noted ν r . The trapping is possible only for some values of β u . This defines the trap's stability diagram. The first region of this diagram corresponds to β z and β x = β y taking the values between 0 and 1.

For real radiofrequency ion traps, terms with n higher than 2 must be included in the expression of the potential [START_REF] Knoop | Lifetime, collisional-quenching and jmixing measurements of the metastable 3D levels of Ca+[END_REF]. The equations of motion for the 3 directions become coupled and non-linear [START_REF] Wang | The non-linear resonance ion trap. Part 2. A general theoretical analysis[END_REF]. When the applied potential has a symmetry along the trap's axis (z axis) and a symmetry about the plane z = 0, only the even terms exist [START_REF] Eades | Nonlinear Resonance Effects During Ion Storage in a Quadrupole Ion Trap[END_REF]. By considering the first 6 terms, the potential (2) becomes [START_REF] Wang | The non-linear resonance ion trap. Part 2. A general theoretical analysis[END_REF][10] [START_REF] Londry | Computer Simulation of Single-Ion Trajectories in Paul-Type Ion Traps[END_REF] 4 -24z 2 r 2 + 3r 4 H 6 = 16z 6 -120z 4 r 2 + 90z 2 r 4 -5r 6 H 8 = 128z 8 -1792z 6 r 2 + 3360z 4 r 4 -1120z 2 r 6 + 35r 8 H 10 = 256z 10 -5760z 8 r 2 + 20160z 6 r 4 -16800z 4 r 6 + 3150z 2 r 8 -63r 10 (5) Where r 2 needs to be replaced by x 2 + y 2 when using the cartesian coordinates. The equations of motion for an ion of mass m and charge Q are

Φ = (U DC + V AC cos(Ωt))(C 0 + C 2 H 2 + C 4 H 4 + C 6 H 6 + C 8 H 8 + C 10 H 10 ) (4) With H 2 = r 2 -2z 2 H 4 = 8z
d 2 x/dξ 2 = T (ξ)xG(x, y, z) d 2 y/dξ 2 = T (ξ)yG(x, y, z) d 2 z/dξ 2 = T (ξ)zH(x, y, z) (6)
Where ξ = Ωt and

T (ξ) = -a z + 2q z cosξ G(x, y, z) = - 1 8 - C 4 16C 2 (-48z 2 + 12r 2 ) - C 6 16C 2 (-240z 4 + 360z 2 r 2 -30r 4 ) - C 8 16C 2 (-3584z 6 + 13440z 4 r 2 -6720z 2 r 4 + 280r 6 ) - C 10 16C 2 (-11520z 8 + 80640z 6 r 2 -100800z 4 r 4 + 25200z 2 r 6 -630r 8 ) H(x, y, z) = 1 4 - C 4 16C 2 (32z 2 -48r 2 ) - C 6 16C 2 (96z 4 -480z 2 r 2 + 180r 4 ) - C 8 16C 2 (1024z 6 -10752z 4 r 2 + 13440z 2 r 4 -2240r 6 ) - C 10 16C 2 (2560z 8 -46080z 6 r 2 + 120960z 4 r 4 -67200z 2 r 6 + 6300r 8 ) (7)
The Mathieu parameters a z and q z are given by

a z = - 16QU DC C 2 mΩ 2 q z = 8QV AC C 2 mΩ 2 (8)
At King Khalid University in Abha, Saudi Arabia, we are using a radiofrequency Paul trap in order to separate between the different calcium isotopes. It is made of six rings with an axis and a mid-plane of symmetry. In order to realize numerical simulations of our work, we use the SIMION software package [12] which gives accurate results but needs long execution time specially for the systematic study of the behavior of tens of ions under the action of different trapping voltages. The trap is similar to the LPCTrap [START_REF] Rodriguez | The LPCTrap Facility : A novel transparent Paul trap for high-precision experiments[END_REF]. This last has been studied experimentally and with simulation at different occasions [START_REF] Duval | Status of the LPCTrap facility at GANIL[END_REF][15] [START_REF] Mery | Measurement of the beta-neutrino angular correlation parameter in the decay of 6He using a Paul trap[END_REF] [START_REF] Rodriguez | The LPCTrap facility for in-trap decay experiments[END_REF]. Mainly the ions spatial and velocity distributions, the mean ion's kinetic energy and the oscillation frequencies have been investigated. In a previous research work, we developed a numerical method to study the behavior of ions in an ideal Paul trap [START_REF] Herbane | Minimized computational time method for the dynamics of ions trapped in an ideal quadrupole ion trap[END_REF]. It is based on the power series solution of the differential equations [START_REF] Boas | Mathematical Methods in the Physical Sciences[END_REF]. It gives accurate results in relatively small calculation time. In this work we present the extension of this method when higher order terms are present in the trapping potential. Our results are compared to the SIMION simulations and to the results of the LPC Caen.

Series solution of the equations of motion

We look for the solutions u(ξ) (u = x, y or z) of equation ( 6) as power series [START_REF] Boas | Mathematical Methods in the Physical Sciences[END_REF]

u = ∞ n=0 A u n (ξ -ξ 0 ) n (9)
Where ξ 0 is an arbitrary constant. The second derivative of u(ξ) is then

d 2 u dξ 2 = ∞ n=0 (n + 1)(n + 2)A u n+2 (ξ -ξ 0 ) n (10)
To apply the power series method, we replace the function T (ξ) by its Taylor expansion around ξ 0 .

T (ξ) = ∞ n=0 A T n (ξ -ξ 0 ) n (11)
The coefficients A T n are calculated as follows

A T 0 = 2q z cosξ 0 A T 1 = -2q z sinξ 0 A T n = - A T n-2 n(n -1) f or n ≥ 2 (12) A T 0 is then replaced by A T 0 -a z to get A T 0 = 2q z cosξ 0 -a z (13)
The functions G and H of equation ( 6) need also to be replaced by power series. For that we use the coefficients of x(ξ), y(ξ) and z(ξ) to calculate the coefficients of x 2 (ξ), y 2 (ξ), z 2 (ξ). These last are used to calculate the coefficients of r 2 (ξ) and those of all the powers and products appearing in G and H. We have the property that the product of the series

F 1 = ∞ n=0 A (1) n (ξ -ξ 0 ) n and F 2 = ∞ n=0 A (2) n (ξ -ξ 0 ) n is the series F 3 = ∞ n=0 A (3) n (ξ -ξ 0 ) n such that A (3) n = n i=0 A (1) i A (2) n-i (14)
and

F 1 + F 2 gives the series F 4 = ∞ n=0 A (4)
n (ξ -ξ 0 ) n with

A (4) n = A (1) n + A (2) n ( 15
)
The third property we need to calculate the series of G and H is that if we add a constant term C to a series, F 1 for example, we get the series

F 5 = ∞ n=0 A (5)
n (ξ -ξ 0 ) n such that

A (5) 0 = C + A (1) 0 A (5) n = A (1) n f or n > 0 (16)
So, if we have the series of of x(ξ), y(ξ) and z(ξ) using equations ( 14), ( 15) and ( 16), we find the series of H and G then we calculate the series of xG, yG and zH. Finally by multiplying these last by the series of T (ξ) given by equation [START_REF] Londry | Computer Simulation of Single-Ion Trajectories in Paul-Type Ion Traps[END_REF], we get the three second terms of equation [START_REF] Gudjons | Influence of anharmonicities of a Paul trap potential on the motion of stored ions[END_REF]. We write them

∞ n=0 S x n (ξ -ξ 0 ) n , ∞ n=0 S y n (ξ -ξ 0 ) n and ∞ n=0 S z n (ξ -ξ 0 )
n for the equations of x, y and z respectively.

Using equation [START_REF] Eades | Nonlinear Resonance Effects During Ion Storage in a Quadrupole Ion Trap[END_REF] and equating both terms of equation ( 6) we find the relations

A x n+2 = S x n (n + 1)(n + 2) A y n+2 = S y n (n + 1)(n + 2) A z n+2 = S z n (n + 1)(n + 2) ( 17
)
These are recursion relations. In fact, knowing A x i , A y i and A z i for the order i = 0, 1, ..., n, one calculates S x n , S y n and S z n then uses equation [START_REF] Rodriguez | The LPCTrap facility for in-trap decay experiments[END_REF] to calculate A x n+2 , A y n+2 and A z n+2 . The coefficients of order 0 (proportional to the initial position) and order 1 (proportional to the initial velocity) are needed to get the others at any order.

Practically, we terminate the sum in equation ( 9) at n = n max . That means that the coordinates of an ion are given by polynomials of degree n max and not by infinite series. If we want to calculate the trajectory of an ion from t = 0 to a relatively large time t max , a large value of n max is necessary. However, this will not work in practice because the processor will truncate the numbers when high powers of time or phase are calculated and gives wrong results. The solution is to divide the large time interval into small intervals of width ∆t or in an equivalent way, divide the corresponding phase interval into small intervals of width ∆ξ. For every small interval, ξ 0 of equation ( 9) is taken equal to l∆ξ with l an integer having the value 0 for the first interval, 1 for the second one and so on. The n max coefficients of the polynomial for x(ξ), y(ξ) and z(ξ) are calculated for every interval. The first 2 of them needed to apply the recursion relations [START_REF] Rodriguez | The LPCTrap facility for in-trap decay experiments[END_REF] correspond to

A u 0 = u(ξ 0 ) A u 1 = du dξ | ξ=ξ0 = v u (t = ξ0 Ω ) Ω ( 18
)
Where v u is the component of the velocity in the u direction.

For the first interval, A u 0 is the u coordinate of the initial position and A u 1 is the component of the initial velocity divided by Ω. For the following intervals, A u 0 and A u 0 are obtained by imposing the continuity of u(ξ) and its derivative du dξ .

When we applied this method and solved numerically the equation of motion of an ion in an ideal Paul trap [START_REF] Herbane | Minimized computational time method for the dynamics of ions trapped in an ideal quadrupole ion trap[END_REF], we found that n max = 15 and ∆ξ = 0.38π, give accurate results with a relatively small calculation time. The position and the velocity of the ion are calculated for ξ = 0, 0.38π, 0.76π, ... means with a time step equal to 19% of the RF period. In the following we use the same values for n max and ∆ξ.

Application to the LPCTrap

The LPCTrap is made of six rings having the same axis of symmetry (Z axis) and it has a median plane of symmetry (z = 0 plane). Its general scheme is shown in figure 1. When using the trap, the voltages are applied symmetrically to the rings R1 and R2 while the others are grounded.

Z X z max r max R1 R2 R3 R4 R5 R6
We used SIMION7 and calculate the potential for z between -z max +0.01 mm and z max -0.01 mm and r between 0 and r max -0.01 mm when 1000 V is applied to R1 and R2 and the others grounded. The potential was calculated in points spaced by 0.01 mm in z and r. Fitting the potential by equation (1) for U DC = 1000 V , V AC = 0 and up to the 100 th term where we consider only the even terms, gives the constants

C 0 = 0.378473 C 2 = -3111.505 m -2 C 4 = -9.973175 × 10 4 m -4 C 6 = -1.51631875 × 10 10 m -6 C 8 = -6.1841484375 × 10 10 m -8 C 10 = 3.17524609375 × 10 16 m -10 (19)
The fitting is based on the χ 2 minimization. 1000C 0 = 378.473 V is the potential at the trap's center when 1000 V is applied to R1 and R2 and the others grounded. SIMION gives 378.579377 V. When fitting with a degree between 60 and 130, the maximum variations of these constants relative to their values given in ( 19) are 2.6 × 10 -4 % for C 0 , 4.8 × 10 -4 % for C 2 , 3.5 × 10 -2 % for C 4 , 9.5 × 10 -3 % for C 6 , 9 % for C 8 and 0.1 % for C 10 . Because C 2 is negative , the parameter a z is positive if a positive DC voltage is applied to the inner rings and vice versa (equation ( 8)). q z is always negative. In the experimental work of the LPC Caen, q z was given positive [START_REF] Rodriguez | The LPCTrap Facility : A novel transparent Paul trap for high-precision experiments[END_REF]. When using our convention, that means that it is in fact the absolute value of q z which was given.

When simulating the ions behavior in the trap, we considered that an ion is trapped as long as it remains inside the effective trapping volume we define as |z| < z max and r < r max . If at some time step, an ion does not fulfill these conditions, it is considered as lost and its trajectory is no longer calculated.

For some of our simulations, we took for the ions uniform random initial positions inside the effective trapping volume. For the velocities we used what is known about the LPCTrap to be as close as possible to the experiment. In fact, the ions kinetic energy distribution before the injection in the trap has been measured [START_REF] Darius | A linear radiofrequency quadripole for the cooling and bunching of radioactive ion beams[END_REF]. It is a gaussian with the standard deviation ∆E = 3 eV F W HM . We then start by considering ions with random velocities following a Maxwellian distribution with the temperature T = ∆E kB . k B is the Boltzmann constant. However, most of the ions generated with these initial conditions are lost during the first microseconds of the trapping. We then consider only those which remain trapped after 100 microseconds.

Comparison with SIMION7

In order to compare our method to SIMION, we studied the evolution of the position of an ion of mass 40 amu and charge 1 e up to 50 ms. We considered separately the cases where the ion was moving in the axial z direction and in the radial x direction. For all the cases we used the frequency of the AC field equal to 1 M Hz.

Figure 2 shows the positions for V AC = 329 V and U DC = 0 which corresponds to (q z = -0.5, a z = 0) calculated using equation [START_REF] March | An Introduction to Quadrupole Ion Trap Mass Spectrometry[END_REF].

There is a good overlap between the positions calculated by SIMION and those with our method for the first periods of the ion's motion. Figure 2 (a) corresponds to the an ion starting from rest at t = 0 and z = 5 mm, x = y = 0. The overlap between the 2 curves ceases for long times. There is a phase shift between the 2 curves which increases with time and is equal to about π around t = 3 ms. Figure 2 (b) shows the curve in the same conditions but with the ion starting from z = 1mm. The agreement between the 2 curves is good for longer times and the phase shift is less than π for t = 30 ms. is observed between our results and those of SIMION. Similar to the motion in the z-direction, the rate of the increase of this phase shift is the smaller for the smallest oscillation amplitude. This can be explained by the fact that when the excursion of the ion is large, it enters in regions of space where terms in the potential of order higher than 10, and which we do not include in our simulations, become important. These terms exist in SIMION and they give the difference between the results of the 2 methods (SIMION and our method).

Figure 3 shows other examples corresponding to V AC = 870 V and U DC = -178 V that is (q z = -1.32, a z = -0.54). This working point lies in the limits of the stability diagram with the β x = β x and β z close to unity. The ion starts from the origin at t = 0 and with initial kinetic energy of 0.1 and 0.01 eV for the oscillation in the z-direction and 0.1 eV for the x-direction. The difference between our results and those of SIMION are seen faster than for the cases of Figure 2. We also observe that the agreement between the curves of SIMION and ours depend on the amplitudes of oscillation. The smaller they are the better it is. In all cases, the amplitudes and the frequencies of oscillation we get, are close to those obtained by SIMION. The main advantage of our method is that it allows the calculation of an ion's trajectory in relatively short time. Using a DELL computer of 3.1 GHz frequency, we get the positions and velocities up to 50 ms in less than 4 s while several minutes are needed by SIMION7. Our codes are written in the C language under Windows and SIMION was used in its default computational quality.

Since we have similar oscillation amplitudes and frequencies when using our method and the SIMION software, we expect to get similar results with the 2 methods when considering the maximum excursion of the ions and their frequencies. The observed phase shift does not have an important effect.

The stability Diagram

In order to get the first stability domain of the trap, we start by generating 500 ions of mass 40 amu and charge 1 e. These are Ca + ions. We take for them uniform random positions with -z max < z < z max and r < r max . The velocities are those of a maxwellian distribution with a temperature of 1000 K. These can be ions created by the ionization of hot atoms. For every working point, the trajectories are calculated up to a trapping time of 2 ms. At the end, the number of ions N (q z , a z ) which remain in the trap is counted. These are those which satisfy the trapping condition given in section 3. The RF frequency is 1 M Hz.

Figure 4 gives the result when V AC is scanned between 0 and 800 V with a step ∆V AC of 40 V . U DC is taken 0. That is a scan of q z between 0 and -1.2167 with a step q z = -0.06 and a z = 0. Our result are compared to those we get when using SIMION. We see that the 2 curves have the same general shape. Specially there is a relative minimum for both of them at the same value of V AC = 440 V (q z = -0.67). In figure 5, we show a more detailed investigation of the trap's stability diagram. For that we scanned V AC between 0 and 1000 V with a step ∆V AC = 5 V . That is q z between 0 and -1.52 with ∆q z = -0.0076. For every value of V AC , U DC is scanned between -230 V and 80 V with 2 V step which corresponds to a z between -0.7 and 0.21 with ∆a z = 0.006. Figure 5 shows a 2D plot of N (q z , a z ) as a gray scale map. One recognizes the classical shape of the first stability domain. In the domain we see lines where the number of ions is reduced. They are the non linear resonances [START_REF] Wang | The non-linear resonance ion trap. Part 2. A general theoretical analysis[END_REF] which occur when

n z ν z + n r ν r = nf (20)
n z , n r and n are integers.

A non linear resonance is induced by the terms of the trapping potential having an order larger than 2. When it occurs, the ions oscillation amplitudes become larger and some of them can be lost.

We identify the resonances 4ν z = f and ν z -ν r = 0. This is confirmed by the analytical calculation of these resonances obtained using the recursion relation to get ν z and ν r [START_REF] March | An Introduction to Quadrupole Ion Trap Mass Spectrometry[END_REF]. The analytical results are shown by symbols in the figure. The line 4ν z = f is induced by the absorption by the ions of energy from the RF field and has been observed in some experiments with real traps [START_REF] Alheit | Observation of instabilities in a Paul trap with higher-order anharmonocities[END_REF] [START_REF] Champenois | Characterisation of a miniature Paul-Straubel trap[END_REF]. The line ν z -ν r = 0 corresponds to the coupling between the ions different degrees of freedom and has also been observed experimentally [START_REF] Vedel | Evidence of radial-axial motion couplings in an rf stored ion cloud[END_REF].

The buffer gas cooling

The collision of the trapped ions with the atoms or the molecules of a buffer gas induces their cooling or heating [START_REF] Major | Exchange-Collision Technique for the RF Spectroscopy of Stored Ions[END_REF]. When the mass of the ion is much bigger than the mass of the buffer gas atoms, its energy is reduced (cooling).

The Monte Carlo simulation of this phenomenon has been implemented since decades [START_REF] Londry | Computer Simulation of Single-Ion Trajectories in Paul-Type Ion Traps[END_REF]. The most simple method is based on the Langevin collision theory where the collision probability is independent on the ion's velocity. The other method is the hard sphere model (HS1) which considers the ions and the gas particles as spheres of radii r i and r g respectively. The collision cross section is given by

σ hs = πr 2 hs = π(r i + r g ) 2 (21)
Its implementation is more complicated than the Langevin theory but it has been shown that the convenience of a model or the other depends on the parameter [START_REF] Smith | Unimolecular decomposition of a polyatomic ionin a variable-temperature selected-ion-flow-drift tube : Experiment and theoretical interpretation[END_REF] ǫ = e 2 α e 2(4πε 0 ) 2 k B T g r 4 hs [START_REF] Champenois | Characterisation of a miniature Paul-Straubel trap[END_REF] α e is the gas particle's electric polarizability, ε 0 the electric permittivity of vacuum and T g is the gas temperature. If ǫ ≪ 1, the hard sphere collisions dominate. In our case we want to study the cooling of 6 Li + ions by molecular Hydrogen at room temperature (300 K). Using α e = 0.8 Å3 [START_REF] Miller | Atomic and Molecular Polarizabilities[END_REF] and r hs = 5.3 Å (this value will be justified bellow), the parameter ǫ is about 0.27. We then choosed the model HS1.

In the HS1 model, the collision probability per unit time can be calculated by [START_REF] Zipkes | Kinetics of a single trapped ion in an ultracold buffer gas[END_REF] dP dt = nσ hs v (23) n = p kB Tg is the buffer gas density, v is the mean ion's velocity relative to the buffer gas particles and p the buffer gas pressure. When considering a Maxwellian distribution for the gas velocities, v has the approximate expression [START_REF] He | Simulation of external ion injection, cooling and extraction processes with SIMION 6.0 for the ion trap/reflectron time-offlight mass spectrometer[END_REF]

v = v 2 i + v g 2 ( 24
)
Where v i is the ion's speed, v g = 8kB Tg πmg is the mean speed of the gas particles and m g their mass.

The collisions between an ion and a buffer gas particle are supposed to be elastic. To get the ion's velocity after the collision, a reference frame other than the laboratory frame is more convenient. SIMION8 uses a frame where the gas particle is stationary [29]. Parks and Szoke use the center of mass (CM) frame [START_REF] Parks | Simulation of collisional relaxation of trapped ion clouds in the presence of space charge fields[END_REF]. We follow their procedure. To get the velocity of the ion after the collision, its velocity and the velocity of the buffer gas particle in the CM frame before the collision are calculated. The z-axis is rotated and made parallel to the ion's velocity. The collision then occurs in the z direction. The angle α between the line connecting the centers of the 2 colliding spheres and the z-axis defines the impact parameter. It varies between 0 and π 2 . The orientation of the 2 spheres about the rotated z-axis is defined by an angle ϕ r taking the values between 0 and 2π. Knowing these angles, the velocity of the ion after the collision is calculated in the rotated CM frame then transformed back to the laboratory frame.

For every ion, the collision probability is calculated every time step, and a random number between 0 and 1 is generated. If this number is smaller or equal to the calculated probability, the collision is supposed to occur. A random velocity with components following maxwellian distributions with the temperature T g are affected to the buffer gas particle. The new velocity of the ion is calculated according to the procedure explained above.Two additional random numbers are needed. One between 0 and 1 whose value is affected to sin 2 α, and a second between 0 and 2π affected to ϕ r [START_REF] Parks | Simulation of collisional relaxation of trapped ion clouds in the presence of space charge fields[END_REF].

In the LPC Caen experiment, the main buffer gas is molecular hydrogen. The cooling of singly charged ions of mass 6 amu has been simulated for a trapping frequency of 1.3 M Hz and V AC = 80 V [START_REF] Duval | Status of the LPCTrap facility at GANIL[END_REF] with a buffer gas pressure of 5 × 10 -4 mbar. The same reference gives the experimental cooling times τ of 6 Li + for the pressures of 6 × 10 -6 , 10 -5 and 4.3 × 10 -5 mbar. They are 12.9, 7.2 and 2.5 ms respectively. The final temperature of the ions has also been measured by applying electric pulses to the rings R1, R2 and R3 and measuring the time of flight (TOF) to a microchannel plate detector [START_REF] Flechard | Temperature measurement of 6He+ ions confined in a transparent Paul trap[END_REF]. The comparison of the TOF to SIMION8 simulations give an ion's mean kinetic energy of 0.11 eV .

In our simulations, the test of the occurrence of a collision is realized every time step. When a collision occurs, the velocity given by equation ( 18) is changed according to the described procedure. For every working point, the initial positions and velocities of the ions were generated according to the procedure described in section 3. We first studied the evolution of 1000 6 Li + up to 90 ms. For that, we start by fixing a time step for the calculation of some physical quantities of the ion cloud. The latter are the standard deviations σ x , σ y , σ z , σ vx , σ vy and σ vz for the positions and the velocities and the ion's mean kinetic energy E . They are calculated every simulation time step during the first 5 periods of the trapping field following every calculation step. We then calculate the mean values of the 7 quantities for the 5 periods. We get σ x , σ y , σ z , σ vx , σ vy , σ vz and E for every calculation step.

We varied the value of the collision cross section σ hs until getting the closest cooling time to the experimental value for the pressure of 10 -5 mbar. τ is obtained by fitting the ion's mean kinetic energy by a function exponentially decaying in time. We fixed for that U DC to 0, V AC to 65 V and f to 1.15 M Hz. These are the usual working parameters in the LPC Caen [START_REF] Rodriguez | The LPCTrap Facility : A novel transparent Paul trap for high-precision experiments[END_REF]. We get σ hs = 90 Å2 , corresponding to a cooling time of 7.37 ms. This is shown in figure 6. When used with the pressures of 6 × 10 -6 and 4.3 × 10 -5 mbar, we find cooling times of 11.8 and 2 ms respectively, which are also close to the experimental values. The figure shows that after a trapping time of about 55 ms, the ions are thermalized. The average value of the mean kinetic energy considered for the points after thermalization is 0.128 eV for 6 × 10 -6 mbar, 0.129 eV for 10 -5 mbar and 0.128 eV for 4.3 × 10 -5 mbar. This is in good agreement with the experimental results [START_REF] Flechard | Temperature measurement of 6He+ ions confined in a transparent Paul trap[END_REF].

From σ hs = 90 Å2 we calculate r hs = 5.3 Å. this last gives ǫ = 0.27 (equation ( 22)) which justifies the use of the hard sphere model.

In reference [START_REF] Duval | Status of the LPCTrap facility at GANIL[END_REF] simulated results for the behavior of 6 Li + ions in residual H 2 buffer gas are presented as well. Realistic 6 Li + -H 2 interaction potential was used. The simulations focus on the spatial and velocity distributions of the trapped 6 Li + . The trapping voltages were U DC = 0 and V AC = 80 V with a frequency of 1.3 M Hz and a residual pressure of 5 × 10 -4 mbar. In order to test our method, we used it to reproduce the results of these simulations. We have a calculation time of 100 µs and the total time 3000 µs. In figure 7 we show the evolution of σ x , σ y and σ z for the five periods following a trapping time of 1000 µs. The mean values of these standard deviations as function of the trapping time are given in figure 8. We observe a thermalization after 1000 µs. The average value of the standard deviations for the thermalized ions are σ x = 1.16 mm, σ y = 1.14 mm, σ z = 0.6 mm, σ vx = 1121 m/s, σ vy = 1094 m/s and σ vz = 1217 m/s. They are close to the results of reference [START_REF] Duval | Status of the LPCTrap facility at GANIL[END_REF] which are σ x = 1.2 mm, σ z = 0.63 mm, σ vx = 1115 m/s and σ vz = 1200 m/s.

The dipolar excitation

The dipolar excitation is a common experimental method used for the determination of the ions oscillation frequencies (see for example [START_REF] Vedel | New Schemes for Resonant Ejection in R.F. Quadrupole Ion Traps[END_REF]). It consists on the superposition of a dipolar alternating voltage of small amplitude and variable frequency (tickle) to the trapping field. When the tickle's frequency is equal to the ions oscillation frequencies, there is a resonance. Ions absorb energy and the amplitude of their motion increases. They may be lost from the trap.

The additional voltage is given by C 1 cos(ωt)z where C 1 is a constant. It induces the additional force -C 1 Qcos(ωt) in the z direction. In the system of equations ( 6), the equations for x and y remain unchanged, while in the

with

e 0 = cos(ω rel ξ 0 ) e 1 = -ω rel sin(ω rel ξ 0 ) e n = -ω 2 rel e n-2 n(n -1) f or n ≥ 2 (26)
and ξ 0 being equal to N ∆ξ for the N th step. This affects the recursion relation [START_REF] Rodriguez | The LPCTrap facility for in-trap decay experiments[END_REF] where the term -F d en (n+1)(n+2) has to be added to the expression of A z n+2 . In the case of the LPCTrap, the dipolar excitation has been realized by applying an additional voltage V t cos(ωt) to one of the inner rings, R1, and -V t cos(ωt) to the second one R2 [START_REF] Mery | Measurement of the beta-neutrino angular correlation parameter in the decay of 6He using a Paul trap[END_REF] [START_REF] Rodriguez | The LPCTrap facility for in-trap decay experiments[END_REF]. The tap is loaded with 6 Li + ions, the tickle is applied for 10 ms with some frequency ω and then the number of the ions remaining in the trap is counted with an MCP detector. The operation is repeated for different ω. When there is no resonance, the number of the ions remaining in the trap is almost constant. At the resonance, this numbers clearly decreases as is shown in figure 9(a).

To simulate the effect of such a field, we used SIMION7 to calculate the electric potential inside the trap when 1 V is applied to R1 and -1 V to R2. These additional voltages, produce in the central region of the trap, an electric potential which can be approximated by a linear function of z. The excitation is then mainly axial. The region where the potential has been calculated has -5 mm ≤ z ≤ 5 mm,0 ≤ r ≤ 5 mm and we used 0.01 mm step for z and r. This potential has been fitted by the function c 1 z. We find

c 1 = 75.39 m -1 (27) C 1 is equal to c 1 V t .
In order to approach the experimental conditions as well as possible we considered, besides the dipole excitation, also the cooling by buffer gas. Every time step then, the collisions of every ion with the molecules of the buffer gas were implemented by Monte Carlo simulations according to the method described in section 6. The buffer gas was hydrogen at 300 K and 6 × 10 -6 mbar.

For the LPCTrap, the tickle was applied for 10 ms after the cooling of the ions by the buffer gas. The trapped ions are 6 Li + . The RF has 60 V amplitude and 1.15 M Hz frequency [START_REF] Mery | Measurement of the beta-neutrino angular correlation parameter in the decay of 6He using a Paul trap[END_REF]. To simulate this, we start for this working point by considering the evolution of 500 ions up to 55 ms under a hydrogen pressure of 6 × 10 -6 mbar. The positions and the velocities of 100 of the remaining ions were stored and used as initial conditions for the simulation of the dipole excitation with the different tickle frequencies which has a constant amplitude of 0.6 V . We varied the tickle's frequency with a 1 kHz step and we counted the number of ions remaining in the trap after 10 ms excitation. The result is shown in figure 9 (b).

Fitting the resonance peaks with lorentzians, gives for the fundamental frequency 186.7 ± 0.3 kHz. The experiment gives 188 kHz. In figure 9 (b), the second peak has the frequency 961.6 ± 0.05 kHz. The sum of the frequencies of the two peaks, i.e. 1148.3 ± 0.3 kHz, is almost equal to the frequency of the trapping field.

Conclusion

We developed a method for the simulation of the dynamics of ions trapped in a Paul trap with the presence of terms in the potential up to the order 10. The equations of motions of the ions are solved numerically. Every time step, the coordinates of the ion are represented by a 15 degree polynomials of time. The coefficients of the polynomials are determined by imposing the continuity of the position and the velocity and by using the power series method. When applied to particular cases, the method gives oscillation amplitudes and frequencies close to those obtained by SIMION7, however our calculations are much faster.

The method was then applied to an existing trap. When knowing the coefficients of the trapping potential, we could consider the evolution of the position of hundreds of ions up to several tens of milliseconds. This is the time scale for the trap. By scanning the Mathieu parameters and counting the number of ions remaining in the trap after 2 ms for every working point, we get the first stability domain which shows the nonlinear resonances.

The buffer gas cooling is introduced by a Monte Carlo method where The hard sphere model is used, here again the results are close to the experimental ones and those obtained by simulations.

At the end we introduced a dipole excitation to ions cooled by the buffer gas. The results show good agreement with the experiment.

Aknowledgements

The present work has been accomplished under the project No 09-ADV826-07 funded by KACST (King Abdul Aziz City for Science and Technology) through the Long Term Comprehensive National Plan for Science, Technology and Innovation program in Saudi Arabia. The authors thank Yvan Merrer from LPC Caen for the technical informations about their trap.

Figure 1 :

 1 Figure 1: Scheme of the LPCTrap. zmax = 8.8 mm, rmax = 10.5 mm.

Figure 2 :

 2 Figure 2: Evolution of the position of an ion starting from rest at different initial positions (x 0 , y 0 , z 0 ). Working point (qz = -0.5, az = 0) (a) x 0 = y 0 = 0, z 0 = 5 mm. (b) x 0 = y 0 = 0, z 0 = 1 mm. (c) x 0 = 5 mm, y 0 = z 0 = 0. (c) x 0 = 1 mm, y 0 = z 0 = 0.

Figure 3 :

 3 Figure 3: Evolution of the position of an ion starting from the origin with different initial kinetic energies E 0 . Working point (qz = -1.32, az = -0.54) (a) E 0 = 0.1 eV with the motion in the z-direction. (b) E 0 = 0.01 eV with the motion in the z-direction. (c) E 0 = 0.1 eV with the motion in the x-direction.

Figure 4 :

 4 Figure 4: Number of trapped ions as function of V AC for U DC = 0

Figure 5 :

 5 Figure 5: The first stability domain of the LPCTrap as obtained by the simulations for Ca + ions. The symbols give the theoretical limits of the domain and the positions of the non linear resonances.

Figure 6 :

 6 Figure 6: Cooling times for Ca ions in a H 2 buffer gas at three different pressures.

Figure 7 :

 7 Figure 7: Standard deviation of the positions around a trapping time of 1000 µs.

Figure 8 :

 8 Figure 8: Evolution of the mean value of the standard deviations of the positions and velocities

Figure 9 :

 9 Figure 9: (a) Experimental dipole excitation spectrum of 6 Li + for V AC = 60 V and f = 1.15 M Hz obtained from Ref. [16]. (b) The simulated spectrum for the same conditions

[ 1 ]

 1 C. Hao, R.E. March. Ion Trap Mass Spectromety, A survey of recent research activity in quadrupole ion trap mass spectrometry. Int. J. Mass Spectrom. 212 (2001) 337-357


M A Veganzones

M Simões

G Licciardi

J Bioucas

J Chanussot

HYPERSPECTRAL SUPER-RESOLUTION OF LOCALLY LOW RANK IMAGES FROM COMPLEMENTARY MULTISOURCE DATA

Remote sensing hyperspectral images (HSI) are quite often locally low rank, in the sense that the spectral vectors acquired from a given spatial neighborhood belong to a low dimensional subspace/manifold. This has been recently exploited for the fusion of low spatial resolution HSI with high spatial resolution multispectral images (MSI) in order to obtain super-resolution HSI. Most approaches adopt an unmixing or a matrix factorization perspective. The derived methods have led to state-ofthe-art results when the spectral information lies in a low dimensional subspace/manifold. However, if the subspace/manifold dimensionality spanned by the complete data set is large, the performance of these methods decrease mainly because the underlying sparse regression is severely ill-posed. In this paper, we propose a local approach to cope with this difficulty. Fundamentally, we exploit the fact that real world HSI are locally low rank, to partition the image into patches and solve the data fusion problem independently for each patch. This way, in each patch the subspace/manifold dimensionality is low enough to obtain useful superresolution. We explore two alternatives to define the local regions, using sliding windows and binary partition trees. The effectiveness of the proposed approach is illustrated with synthetic and semi-real data.

Introduction

In recent years, there has been a huge improvement of spectral and spatial resolutions in the design of remote sensing sensors. However, it seems that it is not possible to achieve both of them simultaneously. This is due, on the one hand, to the system trade-off related to data volume and signal-to-noise ratio (SNR) limitations and, on the other hand, to the specific requirements of different applications [START_REF] Songa | Spatio-spectral fusion of satellite images based on dictionary-pair learning[END_REF]. Then, there is a need for super-resolution techniques that fuse high spectral resolution images, such as hyperspectral images (HSI), with high spatial images such as multispectral images (MSI) or panchromatic images, in order to obtain high spectral and spatial (super-resolution) images.

Recently, some techniques dedicated to the fusion of HSIs and MSIs have been proposed. A general trend is to associate this problem with either linear spectral unmixing, which assumes that the underlying data can be described by a mixture of a relatively small number of "pure" spectral signatures, corresponding to the materials present in the scene or, with the learning of a spectral dictionary that codifies the information present on the images. Since both HSIs and MSIs capture the same scene, the underlying materials (the so-called endmembers) or the dictionaries should be the same. Therefore, the spectral information extracted from one of the images should also be able to explain the other one. Due to the high spectral resolution of the HSIs, the endmembers or the dictionary are extracted from these data, and are then used to reconstruct the MSI. Since MSIs have high spatial resolution, the final reconstructed image will be at super-resolution.

Related work

Zurita et al. [START_REF] Zurita-Milla | Unmixing-based landsat tm and meris fr data fusion[END_REF], introduced one of the first unmixing-based approaches to the fusion of remote sensing multiband images. A related approach is proposed in [START_REF] Kawakami | High-resolution hyperspectral imaging via matrix factorization[END_REF], where a very high-resolution hyperspectral image is estimated from a lowerresolution hyperspectral image and a high-resolution RGB image. The method starts by identifying an unmixing matrix used to represent the hyperspectral spectra and then uses this matrix in conjunction with the RGB input to compute, via sparse regression, representation coefficients for the high-resolution hyperspectral image. This methodology can be viewed as a factorization of the input into a mixing matrix and a set of maximally sparse coefficients. An approach with similar flavour is proposed in [START_REF] Charles | Learning sparse codes for hyperspectral imagery[END_REF]. The main difference is that the mixing matrix is replaced by a dictionary learnt using a non-negative matrix factorization with sparsity regularization on the code. In [START_REF] Huang | Spatial and spectral image fusion using sparse matrix factorization[END_REF], the hyperspectral data is unmixed via the K-SVD algorithm, and the multispsectral data is reconstructed using orthogonal matching pursuit to induce sparsity. Authors in [START_REF] Songa | Spatio-spectral fusion of satellite images based on dictionary-pair learning[END_REF] proposed a method where two dictionaries were learned from the two different data, and then used a dictionary-pair learning method to establish the correspondence between them. A similar and older technique introduced in [START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF] alternately unmixes both sources of data to find the signatures and the abundances of the endmembers.

Contribution

In real world HSI, it is very likely that in a small spatial neighbourhood the number of different materials is small, i.e., these images are locally low rank. We take advantage of this property to propose a local super-resolution methodology. Firstly, we propose two approaches, one using a square sliding window of fixed size and a second using a binary partition tree representation, to obtain a set of patches from the image, such that the set of pixels with indices in each patch span a subspace of lower dimensionality than that of the whole image. Secondly, we exploit the locally low rank property of the pixels of each patch in two different ways: by using an endmembers induction algorithm (EIA) to find the low rank subspace spanned by the spectral signatures of the materials present in each patch; and, by a singular value decomposition (SVD)-based approach to find the local low-rank subspace of the data of each patch. The proposed methodology does not need to estimate the spatial blur, as far as it is constant across bands.

Outline

The remainder of the paper is organized as follows: in Sec. 2 the super-resolution problem is formulated, in Sec. 3 the proposed local super-resolution methodology is introduced and finally, in Sec. 4 experimental results and some conclusion remarks are provided.

Problem formulation

Let X ∈ R n h ×n denote a HSI with n h spectral bands (rows of X) and n = n x ×n y pixels (columns of X). We may interpret X either as a collection of n h images (or bands) of size n x ×n y , each one associated to a given wavelength interval, or as a collection of n spectral vectors of size n h , each one associated with a given pixel. In this work, we are concerned with the estimation of X, which we term the original HSI, from two degraded observations of X: a) a low spatial resolution HSI, Y h ∈ R n h ×(nx/d)×(ny/d) , where d > 1 denotes a spatial downsampling factor, and b) a MSI, Y m ∈ R nm×nx×ny , where n m ≪ n h . We assume that Y h is generated as

Y h = XBM + N h , (1)
where B ∈ R n×n is a matrix modeling band independent sensor blur, M ∈ R n×(n/d 2 ) is a masking matrix accounting for spatial downsampling of size d on both spatial dimensions, and N h is an additive perturbation. Concerning the MSI, Y m , we assume the generation model

Y m = RX + N m , (2)
where the matrix R ∈ R nm×n h holds in its columns the n h spectral responses of the multispectral sensor, and N m is an additive perturbation.

Let us suppose that it is possible to learn a dictionary D ∈ R n h ×n d from the hyperspectral image Y h , and that the columns X, denoted by x i , for i = 1, . . . , n, may be sparsely represented as linear combinations of the columns of D. That is, given x i for i ∈ S, there is a sparse vector α i ∈ R n d (i.e., only a few components of α i are non-zero) such that:

x i = Dα i . (3)
By replacing (3) in [START_REF] Zurita-Milla | Unmixing-based landsat tm and meris fr data fusion[END_REF], we obtain:

Y m = RDA + N m , (4)
where A ≡ [α 1 , . . . , α n ], is often termed as code in dictionary learning and sparse regression applications. If equation ( 4) can be solved with respect to A, then we may plug its solution into (3) and thereby obtain an estimate of X.

The success of a dictionary-based approach depends fundamentally on the ability to solve (4) with respect to A. The difficulty in solving this system comes from the fact that the system matrix RD ∈ R nm×n d is often fat, i.e., n m < n d , yielding an undetermined system of equations. A typical multispectral sensor has less than 10 bands, quite often 4 in the wavelength interval [0.4, 2.5] microns where most hyperspectral sensors operate, whilst n d is often of the order of a few tens. The ill-posedness of (2) when n m < n d may be cured by exploiting the sparsity of the codes α i , for i ∈ S, which opens the door to all sort of sparse regression techniques, many of them recently introduced in compressive sensing [START_REF] Baraniuk | Compressive sensing [lecture notes[END_REF] applications. Nevertheless, we still face a difficulty: in hyperspectral applications, the columns of D tend to be highly correlated, implying that the mutual coherence of the columns of RD is close to 1. This makes (2) ill-posed even when codes α i , for i ∈ S, are sparse [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF][START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF].

In the next section, we introduce two local dictionary-based techniques conceived to cope with the ill-posedness with origin in the matrix system RD. The main idea, in the vein of the local approaches to image restoration, is to decompose the HSI into patches and build patch-dependent dictionaries such that n m ≥ n d in each patch. Thereby, the inverse problem of solving (2) in each patch is well-posed.

3 Local dictionary based HSI super-resolution

Local spatial patches definition

Let P j denote the set of patches obtained from an image, Y h , such that the set of pixels with indices in each patch, Y j ≡ [y i , i ∈ P j ], span a subspace of lower dimensionality than that of Y h . Next, we introduce two approaches to obtain such patches. The first approach uses a sliding window and the second approach relies on a binary partition tree (BPT) representation [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF].

Using a sliding window of fixed size

The sliding window methodology has been broadly used in image processing. It consists in defining a square window of fixed size so the pixels lying inside the window form a patch. The window then slides over the whole image in a standard zig-zag way, be it with some overlapping or not, eventually covering the whole image and defining the image patches, P j .

Using a BPT representation

The BPT is a hierarchical region-based representation of an image in a tree structure [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF]. In the BPT representation, the leaf nodes correspond to an initial partition of the image, which can be the individual pixels, or a coarser segmentation map. From this initial partition, an iterative bottom-up region merging algorithm is applied until only one region remains. This last region represents the whole image and corresponds to the root node. All the nodes between the leaves and the root result of the merging of two adjacent children regions.

Two notions are of prime importance when defining a BPT, the region model, M R , which specifies how a region R is modelled, and the merging criterion, O(M Rα , M R β ), which is a distance measure between the region models of any two regions R α and R β . Each merging iteration involves the search of the two neighbouring regions which achieve the lowest pair-wise similarity among all the pairs of neighbouring regions in the current segmentation map. Those two regions are consequently merged. To build the BPT representation from the hyperspectral image [START_REF] Valero | Comparison of merging orders and pruning strategies for binary partition tree in hyperspectral data[END_REF][START_REF]Hyperspectral image representation and processing with binary partition trees[END_REF][START_REF] Veganzones | Hyperspectral image segmentation using a new spectral mixture-based binary partition tree-representation[END_REF], Y h , we use the first-order parametric model M R :

M R d = x = 1 N R N R i=1 x i ,
where N R is the number of pixels on the region; and, in order to merge regions, the spectral angle distance:

O M Rα , M R β d = d SAM (x α , xβ ) = arccos xα xβ xα xβ .
Once the BPT representation is built, the BPT is pruned to achieve a partition of the image such that the regions of the partition define the image patches, P j . We use the vote-based pruning approach from [START_REF] Tochon | Binary partition tree as a hyperspectral segmentation tool for tropical rainforests[END_REF], where each region of the BPT representation is voted by its children. A region receives a number of votes equal to the leaves of a children node if the distance between the region model and the children model is above some threshold. Once all the regions have collected their number of votes, the BPT is pruned so the partition cardinality is maximized constrained to the regions of the partition having received at least half the total number of votes.

Local dictionary-based super-resolution

EIA-based approach

Given a set of patches P j , j = 1, . . . , P ; for each P j , we identify a mixing matrix D j by applying an EIA to extract a set of endmembers from the set of pixels

Aj ≥0 Y m,j -RD j A j 2 F j = 1, . . . , P, (5)
where A j ∈ R n d ×|Pj | and Y m,j ∈ R nm×|Pj | gathers, respectively, the columns of A and Y m with indices corresponding to the pixels X j in P j . The inequality A j ≥ 0 is to be understood in the component-wise sense, and • F denotes the Frobenius norm. The constraint A j ≥ 0 in ( 5) is used because, in the linear mixing model, the codes A j represent abundances of materials which are necessarily non-negative [START_REF] Ma | A signal processing perspective on hyperspectral unmixing: Insights from remote sensing[END_REF]. It is possible to add the abundances sumto-one constraint, A T j I = 1 |Pj | , but this constraint is usually dropped due to possible scale model mismatches [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. We remark that, since our observations are spectral vectors and thus non-negative (apart from the noise contribution), the non-negativity constraint in ( 5) is equivalent to a form of constrained ℓ 1 regularization and therefore, to some kind of sparsity enforcement (see [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equations[END_REF] for the details).

In order to solve (5), we use the SUnSAL algorithm [START_REF] Bioucas-Dias | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF], which is an instance of the C-SALSA methodology introduced in [START_REF] Afonso | An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF] to effectively solve a large number of CLS problems sharing the same matrix system. Fig. 1 shows the pseudocode of the proposed HSI Super Resolution via Local Dictionary Learning using EIAs (HSR-LDL-EIA). The algorithm operates on each patch independently. Each patch indexes the set of pixels in the hyperspectral image, Y h,j , and in the multispectral image, Y m,j . From the hyperspectral set of pixels, the subspace matrix, D j , is obtained using any EIA. Then, the SUnSAL algorithm is used to estimate the fractional abundances, A j , from the multispectral set of pixels. The fractional abundances are linearly combined with the estimated subspace to estimate the super-resolution pixels of the patch, X j . The estimated super-resolution pixels of all the patches are combined to form the estimated super-resolution image, X.

Algorithm HSR-LDL-EIA Input:

1. Y h , Y m , R, P j , n d 2. for j = 1 to P 3. Y m,j := [Y m,i , i ∈ P j ] 4.
Y h,j := [Y h,i , i ∈ P j ] 5.

D j := EIA(Y h,j , n d ) 6.

A

j := arg min Aj ≥0 Y m,j -RD j A j 2 F
7.

X j := D j A j 8. X := X j , j = 1, . . . , P

SVD-based approach

In this method, we carry out the dictionary learning by solving the following optimization problem:

min rank(D j ) | min β j ,Dj Y h,j -D j β j 2 F ≤ δ, (6)
where Y h,j are the observed hyperspectral vectors with indexes in P j , and β j are weighting coefficients (the code). The optimization problem (6) aims at finding the dictionary D j of lower rank, desirably lesser than n m , yielding a representation error for Y h,j smaller than δ. Using the local low rank argument, the range of the subspace spanned by the super-resolution pixels X j in patch P j should be close to that of D j and thus, we may obtain the code for X j by solving:

min Aj Y m,j -RD j A j 2 F + τ A j 2 F j = 1, . . . , P. (7)
Comparing ( 7) with (5), we have discarded the non-negativity constraint and have introduced the quadratic term τ A j 2 F , where τ ≥ 0 is a regularization parameter. We have discarded the non-negativity constraint because the columns

Algorithm HSR-LDL-SVD Input: 1. Y h , Y m , R, P j 2. for j = 1 to P 3. Y m,j := [Y m,i , i ∈ P j ] 4. Y h,j := [Y h,i , i ∈ P j ] 5.
U kmax := SVD(Y h,j , k max ) 6.

D j := U k , k := min k : U T k Y h,j 2
F ≤ δ 7. A j := D T j R T RD j + τ I -1 D T j R T Y m,j 8.
X j := D j A j 9. X := X j , j = 1, . . . , P of D j are no more interpretable as spectra of pure materials and thus A j are no more interpretable as abundances. The introduction of a quadratic regularizer aims at improving the noise-amplification we may have owing to the very low singular vectors of RD j .

Although problem ( 6) is nonconvex, it can easily be solved based on the SVD of Y h,j . Let U k ∈ R n h ×k be a matrix holding the k left singular vectors of Y h,j ordered by decreasing value of the corresponding singular values. Then, for rank(D i ) = k, it holds that the solution to: min

β j ,Dj Y h,j -D j β j 2 F = U T k Y h,j 2 F (8)
is achieved for D j = U k [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF]. Therefore, the solution of ( 6) is given by:

D j := U k | k = arg k min U T k Y h,j 2 F ≤ δ
Fig. 2 shows the pseudocode of the proposed HSI Super Resolution via Local Dictionary Learning using Single Value Decomposition (HSR-LDL-SVD). The algorithm works in a similar way to the HSR-LDL-EIA algorithm. The main difference is that the subspace matrix, D j , is obtained from the k columns of the singular value decomposition of the hyperspectral pixels. Then, the fractional abundances, A j , are calculated from the multispectral pixels using the analytic solution to the regularized inversion problem [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF]. Finally, the super-resolution pixels of the patch, X j , are calculated; and, the estimated super-resolution pixels of all the patches are combined to form the estimated super-resolution image, X.

Experimental Results and conclusions

In order to assess the performance of the proposed local super-resolution methods, we used the following experimental methodology: 1) Given a super-resolution image, X, we simulate a low-spatial resolution hyperspectral image, Y h , by applying a Gaussian blurring (B) and downsampling the blurred image by a factor of four (M). 2) We also simulate a high-spatial resolution multispectral image, Y m , by applying a spectral response matrix (R) to the original super-resolution image. 3) White noise was added to both images, Y h and Y m , with 30db and 40db respectively. 4) Then, we estimate the super-resolution image, X, by applying a super-resolution methodology to both images, Y h and Y m . 5) Finally, we compare the original and the estimated super-resolution images, using the RMSE, the SAM, the ERGAS and the Q quality measures [START_REF] Songa | Spatio-spectral fusion of satellite images based on dictionary-pair learning[END_REF].

We compared the two proposed local super-resolution approaches, the HSR-LDL-EIA and the HSR-LDL-SVD, to the global approach based on EIAs, HSR-GDL-EIA. The SVD-based approach has been defined explicitly for the local approach. Its global version performs poorly and thus, there is no point in showing the comparison. For the proposed local approaches, we also compared the sliding window and BPT-based approaches to obtain the image patches. In the case of the EIA-based approach, we used the Vertex Component Analysis (VCA) EIA [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF]. We made the comparison over two datasets:

1. A synthetic image composed of multiple geometric shapes (ellipses and rectangles) of different sizes and orientations, where each geometrical element and the background are formed using a different linear mixture of 5 endmembers from the USGS spectral library, for a total of 45 different endmembers in the image. This image simulates a high spectral variability scenario. In this case, we fixed the number of induced endmembers and selected SVD components to p = 2 for the local approaches, and p = 45 for the global one.

2. The well-known Pavia University dataset. This dataset is a real image representing an urban area mainly composed of buildings, urban vegetation and parking lots. The multi-spectral image was simulated using the spectral response of the first three bands of the Ikonos sensor. For the local approaches, we fixed the number of induced endmembers and selected SVD components to p = 3. For the global approach, this value was set to p = 10.

Tables 1-2 show the results obtained for the synthetic and Pavia University datasets respectively. For each quality measure, it is shown the mean and standard deviation of Monte Carlo runs. Results show that the local approaches outperform the global one in the synthetic scene, and that the EIA-based local approach does the same in the real scene. From the point of view of the patches definition, the BPT representation shows a slightly higher performance, hinting at the relevance of the patching definition process. However, the windowing approach shows slightly lower standard deviation probably due to the effects of noise in the BPT partitioning. The lack of space prevents us from showing the obtained super-resolution images for visual comparison. We encourage the interested reader to access the results extended report at first author's website1 . We'll further work on the validation of the proposed methodologies and on the study of the patches definition by means of other BPT pruning strategies. We'll also face the problem of estimate the spectral response matrix, R, from the data.

Figs. 3 and4 show false color representations of the actual Pavia University dataset and the super-resolution estimated images obtained by the compared approaches.

Figure 1 :

 1 Figure 1: HSI Super Resolution via Local Dictionary Learning using Endmember Induction Algorithms (HSR-LDL-EIA).

Figure 2 :

 2 Figure 2: HSI Super Resolution via Local Dictionary Learning using the Singular Value Decomposition (HSR-LDL-SVD).

Figure 3 :Figure 4 :

 34 Figure 3: False color representations of the Pavia University dataset.

Table 1 :

 1 Mean and standard deviation values of the average RMSE, average SAM, ERGAS and average Q quality measures for super-resolution methodologies comparison using the synthetic dataset.

				Synthetic
			RMSE	SAM	ERGAS	Q
	HSR-GDL-EIA	0.0552 ± 0.0230	4.7777 ± 3.5114	4.7583 ± 1.6699	0.9629 ± 0.0252
	HSR-LDL-EIA	Windows BPT	0.0294 ± 0.0007 0.0290 ± 0.0008	2.3944 ± 0.0610 2.2067 ± 0.0714	2.1254 ± 0.0450 2.0574 ± 0.0447	0.9877 ± 0.0007 0.9878 ± 0.0008
	HSR-LDL-SVD	Windows BPT	0.0295 ± 0.0009 0.0288 ± 0.0008	2.3877 ± 0.0751 2.1810 ± 0.0691	2.1093 ± 0.0468 2.0429 ± 0.0535	0.9875 ± 0.0009 0.9881 ± 0.0008
	X j , identified by the corresponding patch. Having identified the dictionaries,
	D j , the code A is estimated by solving the following constrained least squares
	(CLS) problems:
		min

Table 2 :

 2 Mean and standard deviation values of the average RMSE, average SAM, ERGAS and average Q quality measures for super-resolution methodologies comparison using the real Pavia University dataset.

				Pavia University
			RMSE	SAM	ERGAS	Q
	HSR-GDL-EIA	0.0136 ± 0.0017	2.9654 ± 0.3068	2.2056 ± 0.2994	0.9897 ± 0.0024
	HSR-LDL-EIA	Windows BPT	0.0169 ± 0.0001 0.0138 ± 0.0030	2.9804 ± 0.0077 2.6766 ± 0.4147	2.6302 ± 0.0082 2.1924 ± 0.3118	0.9832 ± 0.0001 0.9887 ± 0.0051
	HSR-LDL-SVD	Windows BPT	0.0901 ± 0.1101 0.0153 ± 0.0015	2.7874 ± 0.1449 3.1503 ± 0.1907	15.042 ± 19.988 2.9803 ± 0.3274	0.8445 ± 0.1345 0.9844 ± 0.0037

http://www.gipsa-lab.grenoble-inp.fr/page pro.php?vid=1728

Acknowledges Thanks for funding to Fundação para a Ciência e Tecnologia (FCT), Portuguese Ministry of Science and Higher Education, projects PEst-OE/EEI/0008/2013 and PTDC/EEI-PRO/1470/2012, and grant SFRH/BD/87693/2012.

* This work was supported by the Delegation Générale de l'Armement (DGA) under contract PDE2012600079.