Ionela Rusu

Ioan Cosmin Manolache-Rusu

Andrei Diaconu

Oleg Palamarciuc

A Il'ya

Gural'skiy

Gabor Molnar

Rotaru Aurelian
email: aurelian.rotaru@usm.ro

Pressure Gradient Effect on Spin-Crossover Materials: Experiment vs. Theory

We studied the effect of a non-hydrostatic pressure on the hysteretic spin crossover in coordination complexes.

By introducing into an Ising-like model a double distribution of the interactions and gap energy, respectively, we were able to generate the major hysteresis loop and the FORC diagram for spin-crossover systems of 10 6 hysterons (like-spin domains). We show that, for high pressure gradients around the spin-crossover system, the thermal hysteresis loop takes an asymmetric shape, in good agreement with the experimental data on pressure effect recorded at low temperatures, below the solidification of the pressure transmitting medium. Interestingly, the FORC diagram method seems to be much more sensitive to local changes than the "bulk" parameters, which characterize the major hysteresis loop.

Introduction

Spin-crossover (SCO) complexes form an interesting class of switchable materials that can exist in two different spin states, known as the low-spin (LS) and high-spin (HS) states. These complexes may be switched between This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582

the two spin states under various stimuli including temperature, pressure, light, magnetic and electrical fields [START_REF] Bousseksou | Molecular spin crossover phenomenon: recent achievements and prospects[END_REF][START_REF] Halcrow | Spin-Crossover Materials. Properties and Applications[END_REF]. Among these stimuli, pressure is one of the most appealing thermodynamic parameters used to probe the spin transition properties of spin-crossover materials since it allows for tuning both the crystal field as well as the interactions between SCO centers. Since the pioneering works of Drickamer [START_REF] Slichter | Pressure-Induced Electronic Changes in Compounds of Iron[END_REF][START_REF] Drickamer | Electronic Transitions and the High Pressure Chemistry and Physics of Solids[END_REF], Kӧnig [START_REF] Konig | The Effect of Pressure on the Thermal hysteresis of the 1 st -Order Spin Transition in Bis (1,10-Phenanthroline-2-Carbaldhyde Phenylhydrazone) Iron (II) Complexes[END_REF] and Gütlich, [START_REF] Ksenofontov | Spin crossover behavior under pressure of Fe(PM-L)2(NCS)2 compounds with substituted 2'-pyridylmethylene 4-anilino ligands[END_REF][START_REF] Garcia | Influences of temperature, pressure, and lattice solvents on the spin transition regime of the polymeric compound [Fe(hyetrz)3]A2•3H2O (hyetrz = 4-(2 '-hydroxyethyl)-1,2,4-triazole and A-= 3nitrophenylsulfonate)[END_REF],

the pressure effect has continuously attracted the interest of the scientific community, being intensively studied until nowadays. The pressure effect on spin-crossover materials has been investigated using various techniques, such as Mössbauer spectroscopy [START_REF] Konig | The Effect of Pressure on the Thermal hysteresis of the 1 st -Order Spin Transition in Bis (1,10-Phenanthroline-2-Carbaldhyde Phenylhydrazone) Iron (II) Complexes[END_REF], X-Ray and neutron diffraction [START_REF] Granier | High-pressure single-crystal x-ray diffraction study of two spin-crossover iron(II) complexes: Fe(Phen)2(NCS)2 and Fe(Btz)2(NCS)2[END_REF][START_REF] Guionneau | High pressure and very low temperature effects on the crystal structures of some iron(II) complexes[END_REF][START_REF] Shepherd | Pressure-induced two-step spin transition with structural symmetry breaking: X-ray diffraction, magnetic, and Raman studies[END_REF][START_REF] Shepherd | High-pressure spin-crossover in a dinuclear Fe(II) complex[END_REF][START_REF] Legrand | Neutron powder diffraction studies of two spin transition Fe-II complexes under pressure[END_REF][START_REF] Tissot | Temperature-and Pressure-Induced Switching of the Molecular Spin State of an Orthorhombic Iron(III) Spin-Crossover Salt[END_REF], X-Ray spectroscopies [START_REF] Roux | Pressure-induced spin-state crossovers in six-coordinate Fe(II)LnL'm(NCS)2 complexes with L=L' and L not equal L': A XANES investigation[END_REF][START_REF] Boillot | Pressure-induced spin-state crossovers at room temperature in iron(II) complexes: comparative analysis; a XANES investigation of some new transitions[END_REF], Raman spectroscopy [START_REF] Molnar | Pressure tuning Raman spectroscopy of the spin crossover coordination polymer Fe(C5H5N)2[Ni(CN)4][END_REF][START_REF] Molnar | Raman spectroscopic study of pressure effects on the spin-crossover coordination polymers Fe(pyrazine)[M(CN)4]•2H2O (M = Ni, Pd, Pt). First observation of a piezo-hysteresis loop at room temperature[END_REF], optical reflectivity [START_REF] Jeftic | Helium gas pressure cell for pressures up to 1 kbar (0.1 GPa) in conjunction with the cold head of a closed-cycle He refrigerator[END_REF][START_REF] Codjovi | Pressure and temperature hysteresis in the spin-transition solid Fe(btr)2(NCS)2•H2O, pure and diluted in Ni matrix[END_REF][START_REF] Tanasa | Piezo-and thermo-switch investigation of the spin-crossover compound [Fe(PM-BiA)2(NCS)2[END_REF][21][START_REF] Jureschi | Thermo-and piezochromic properties of [Fe(hyptrz)]A2•H2O spin crossover 1D coordination polymer: Towards spin crossover based temperature and pressure sensors[END_REF][START_REF] Rotaru | Pressure effect investigated with first-order reversal-curve method on the spin-transition compounds FexZn1x(btr)2(NCS)2•H2O (x=0.6,1)[END_REF][START_REF] Levchenko | The pressure-induced spin transition in the Fe(phen)2(NCS)2 model compound[END_REF][START_REF] Levchenko | Electrical Voltage Control of the Pressure-Induced Spin Transition at Room Temperature in the Microporous 3D Polymer Fe(pz)Pt(CN)4[END_REF], magnetometry [START_REF] Ksenofontov | Spin crossover behavior under pressure of Fe(PM-L)2(NCS)2 compounds with substituted 2′-pyridylmethylene 4-anilino ligands[END_REF][START_REF] Ksenofontov | Pressure-induced spin state conversion in antiferromagnetically coupled Fe(II) dinuclear complexes[END_REF][START_REF] Gutlich | Pressure effect studies in molecular magnetism[END_REF][START_REF] Gutlich | Pressure effect studies on spin crossover systems[END_REF][START_REF] Galet | Pressure effect and crystal structure reinvestigations on the spin crossover system: [Fe(bt)2(NCS)2] (bt=2,2'-bithiazoline) polymorphs A and B[END_REF][START_REF] Agusti | Thermal, pressure and light induced spin transition in the two-dimensional coordination polymer {Fe(pmd)2[Cu(CN)2]2}[END_REF] and dielectric spectroscopy [START_REF] Diaconu | Piezoresistive Effect in the [Fe(Htrz)2(trz)](BF4) Spin Crossover Complex[END_REF]. In general, the application of an external pressure induces a transformation of the SCO complex to the low-spin state, i.e. to the state of lower molecular volume, and has the effect of shifting the switching temperatures towards higher values. However, the interpretation of the pressure effect on molecular SCO phenomena has been, in some cases, compromised by the relatively poor or incomplete information resulted from experimental difficulties related to the need for hydrostatic conditions at low temperatures. The practical implementation consists of immersing the sample in a pressure transmitting fluid meant to support no shear. However, the melting point of fluids is usually increased under pressure and solidification inevitably occurs at some pressure. Beyond this point, the pressure across the experimental volume is generally inhomogeneous and differential as well as shear stresses appear [START_REF] Murata | Pressure transmitting medium Daphne 7474 solidifying at 3.7 GPa at room temperature[END_REF]. Depending on the type of measurement, this can lead to a more or less dramatic decrease in the quality and accuracy of the data and often to the presence of 'anomalies', which might be erroneously ascribed to new physical phenomena. For example, the hysteretic properties of spin-crossover microparticles, embedded in glass-forming or semicrystalline matrices, reveal a larger hysteresis loop as compared to the bulk [START_REF] Tanasa | Hysteretic behavior of Fe(phen)2(NCS)2 spin-transition microparticles vs. the environment: A huge reversible component resolved by first order reversal curves[END_REF][START_REF] Tanasa | Unraveling the Environment Influence in Bistable Spin-Crossover Particles Using Magnetometric and Calorimetric First-Order Reverse Curves[END_REF]. This effect has been attributed to an increase of the cooperativity of the SCO system due to particle-matrix interactions.

An indispensable condition for generating a high-quality hydrostatic pressure is that the pressure transmitting medium should be a fluid in the whole experimental P-T range. When the pressure homogeneity is poor, it is sometimes difficult to distinguish real coexistence of two phases from a simple mixture of two phases due to pressure distribution around the sample. This adversity is likely to occur when the pressure region under discussion is beyond the solidification of the pressure medium.

The purpose of the present study is to investigate the effect of a non-hydrostatic pressure on the thermally induced spin transition in spin-crossover complexes.

Theory: the Ising-like model

A useful microscopic model of SCO, the so-called Ising-like model, has been proposed by Bousseksou et al. [START_REF] Bousseksou | Ising-like model for the 2-step spin crossover[END_REF] nearly 3 decades ago. This approach takes into account the degenerate energy levels and it was successfully used to describe both static and dynamic properties of SCO compounds under different external stimuli (temperature, pressure, magnetic and electric fields) [START_REF] Boukheddaden | One-dimensional Ising-like systems: an analytical investigation of the static and dynamic properties, applied to spin-crossover relaxation[END_REF][START_REF] Chiruta | Lattice architecture effect on the cooperativity of spin transition coordination polymers[END_REF]. The Hamiltonian in the Ising-like model can be written as follows:

( )

( ) , ( ) ˆˆˆˆ1 ln 2 B i ij i j i i i j i H k T g J G       =  - - -    ; (1)
where 𝜎 ̂𝑖 is the fictitious spin operator associated to molecule i having the eigenvalues +1 (when the molecule is in the HS state) and -1 (when the molecule is in the LS state), Δ is the energy gap between the HS and LS states, G stands for a long range interaction parameter that is expressed here in a mean-field approach with   being the average fictitious spin.

The pressure effect has been introduced considering that the energy gap is changing due to the elastic energy change pV   when increasing the applied pressure [START_REF] Slichter | Pressure-Induced Electronic Changes in Compounds of Iron[END_REF]:

( ) ( ) Δ Δ 0 p p p V  = = + ; (2)
where V  is the volume variation of the SCO molecules during the spin transition phenomenon and p is the externally applied pressure. In the mean-field approach, the Hamiltonian HMF can be written as follows:

( )

(n ˆ2 l ) p k T g H qJ G MF B     = +  - - ; ( 3
)
where q is the number of nearest neighbours. The average fictitious spin  of a hysteron will be further denoted by m that can be expressed as:

( ) ln () tanh 2 B BB k T g p qJ L m m k T k T -  + =+   . ( 8
)
For the sake of simplicity, in the following calculations we use several normalized parameters: the energy gap

( ) / B pk  = 
, the interaction parameter In order to closer represent a real system, we considered a spin crossover system of 10 6 hysterons (i.e. spin-like domains), characterized by a double Gaussian distribution of energy gaps and interaction parameters of the form:

( ) ( ) ( ) 22 00 1 , exp exp 2 2 2 J J J JJ PJ              -- - - =             (9)
Thus, the high spin fraction can simply be calculated by using the normalization condition:

( , ) ( , )

N HS i HS P J n J n N    =  (10)
Due to the high number of hysterons for our system, we employed high performance calculation (HPC) by using a Nvidia GeForce GTX760 GPU (GPU = Graphics Processing Unit) with a parallel architecture. Thus, the 1152 hysterons have been simultaneously simulated by using the 1152 CUDA cores of the GPU, allowing us to reasonably reduce the calculation time (from 14 h using sequential compilation to 20 msec. using HPC) of such systems with a macroscopic number of hysterons, close to that of a real system.

Results and discussions

The pressure gradient has been naturally introduced into the Ising-like model through the standard deviation of the reduced energy gap distribution, ()   , expressed in Kelvin. The following computational parameters have been used: 0 1200 K = , ln 7 g = .

We can see that for low pressure gradients, i.e. under hydrostatic pressure, the thermal hysteresis loop has a symmetric shape and becomes strongly asymmetric as the pressure gradient is increased. The same effect, i.e. asymmetric thermal hysteresis loop for a high pressure gradient, has been observed for various values of the interaction parameters. Such effects have been experimentally observed mostly in magnetic measurements under pressure performed on complexes characterized by low switching temperatures [START_REF] Ksenofontov | Spin crossover behavior under pressure of Fe(PM-L)2(NCS)2 compounds with substituted 2'-pyridylmethylene 4-anilino ligands[END_REF] or in nanoparticles embedded in various matrixes [START_REF] Tanasa | Hysteretic behavior of Fe(phen)2(NCS)2 spin-transition microparticles vs. the environment: A huge reversible component resolved by first order reversal curves[END_REF][START_REF] Tanasa | Unraveling the Environment Influence in Bistable Spin-Crossover Particles Using Magnetometric and Calorimetric First-Order Reverse Curves[END_REF]. Indeed, in magnetic measurements, a clamp-type pressure cell is typically used [START_REF] Gaspar | Pressure effect investigations on spin-crossover coordination compounds[END_REF].

In such pressure cells, the pressure is transmitted to the sample through silicon oil such as Daphne 7373 or Flourinert 77/70.

Experiment vs Theory

In most of the previous works published on this topic, the hydrostatic character of the pressure cell is evaluated from the abruptness of the magnetic transition curve (from diamagnetic to superconductive state) of a superconducting manometer (usually Pb or Sn thin wires). In this context, we measured the superconducting transition of a Pb wire by using two pressure media: distilled water and Daphne 7373 (see Figure S1 in ESI). Surprisingly, the transition curves recorded on Pb in both water and Daphne 7373 have the same shape. However, a shift to lower values of the transition temperature of Pb in water, compared with that one of Pb in Daphne 7373 has been observed. This denotes that water applies a higher pressure around the Pb wire due to its negative coefficient of compressibility.

We have investigated the pressure effect on the magnetic properties of [Fe(PM-BiA)2(NCS)2] spin-crossover complex by using a clamp pressure cell. The choice of this complex has been motivated by the low values of the switching temperatures (around 170 K), bellow the melting point of the pressure transmitting medium. The sample was sealed in a Teflon sample holder surrounded by Daphne 7373 as the pressure transmitting medium. The magnetic measurements were performed using a MPMS3 SQUID magnetometer (Quantum Design Inc., USA), in DC mode, under a static magnetic field of 1000 Oe with a temperature sweep rate of 2 K/min. The pressure value inside the pressure cell has been determined by measuring the superconducting transition temperature of a Pb wire. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. crossover complex (polymorph I) shows a wider, more asymmetric thermal hysteresis loop (Fig. 2a). A similar behavior was noticed in our simulations for higher values of the pressure gradient (Fig. 2b). Thus, the asymmetric shape of the thermal hysteresis loop recorded on the [Fe(PM-BiA)2(NCS)2] complex might be attributed to the presence of pressure gradient around the sample. Indeed, Daphne 7373 has a freezing/melting point, at ambient pressure, of ca. 180 K (see Figure S2 in the ESI) which increases with increasing applied pressure. It is worth noting herein that the [Fe(PM-BiA)2(NCS)2] spin crossover complex (polymorph I) under external pressure undergoes a crystallographic phase transition [START_REF] Ksenofontov | Spin crossover behavior under pressure of Fe(PM-L)2(NCS)2 compounds with substituted 2'-pyridylmethylene 4-anilino ligands[END_REF]21], which can amplify the deformation of the hysteresis loop.

PLEASE CITE THIS ARTICLE AS

For a more in depth analysis of this effect, we employed the first-order reversal curve (FORC) diagram method.

FORC diagram method

The FORC diagram method has been proposed as an alternative characterization tool to the classical major hysteresis loop (MHL) analysis [START_REF] Roberts | First-order reversal curve diagrams: A new tool for characterizing the magnetic properties of natural samples[END_REF]. It is already known that the FORC diagram method provides detailed information from within the major hysteresis loop and enables determination of the distribution of switching temperatures and interaction fields for all "particles" (such as the like-spin mono-domains) that contribute to the hysteresis loop [START_REF] Rotaru | Pressure effect investigated with first-order reversal-curve method on the spin-transition compounds FexZn1x(btr)2(NCS)2•H2O (x=0.6,1)[END_REF][START_REF] Pike | The effect of magnetic interactions on low temperature saturation remanence in fine magnetic particle systems[END_REF][START_REF] Enachescu | FORC method applied to the thermal hysteresis of spin transition solids: first approach of static and kinetic properties[END_REF][START_REF] Tanasa | Physical parameter distribution in spin transition systems derived from FORC data[END_REF][START_REF] Enachescu | Kinetic hysteresis in spin crossover solids analyzed using FORC diagrams[END_REF]. In the case of a temperature induced spin transition, the FORCs are a specific class of minor hysteresis loops, for which the sweep of the temperature is reversed at a specific temperature, Ta, on one of the branches of the major thermal hysteresis loop, back to a specific temperature at which the domain structure is saturated. The temperature steps are chosen such that Ta and Tb are regularly spaced, which means that This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. similar with that one experimentally recorded in the 1D spin crossover compounds, at ambient pressure [START_REF] Dirtu | Insights into the Origin of Cooperative Effects in the Spin Transition of [Fe(NH2trz)3](NO3)2: the Role of Supramolecular Interactions Evidenced in the Crystal Structure of [Cu(NH2trz)3](NO3)2•H2O[END_REF][START_REF] Weber | Influence of Hydrogen Bonding on the Hysteresis Width in Iron(II) Spin-Crossover Complexes[END_REF], to +0.9 for high pressure gradients.

PLEASE CITE THIS ARTICLE AS

MHL vs FORC

A remarkable finding that we would like to highlight here is the sensitivity of the FORC diagram method regarding small local changes compared to the classical major hysteresis loop.

Conclusions

In this work we analyzed the effect of a non-hydrostatic pressure on the spin transition properties of molecular spin crossover complexes in terms of major thermal hysteresis loop and FORC diagram methods. We have shown that for high pressure gradients, the thermal hysteresis loop becomes asymmetric. Remarkably, the FORC diagram method has been proved to be much more sensitive to local changes than the "bulk" parameters, which characterize the major hysteresis loop. This makes the FORC diagram method a suitable technique to characterize the hydrostaticity degree of a pressure cell. These results aim to shed light on a series of reported results on the pressure effect on spin crossover materials and on hysteretic properties of spin-crossover microparticles embedded in glass-forming matrices, analyzed at cryogenic temperatures, below the solidification of the medium surrounding the SCO particles.

SUPPLEMENTARY MATERIAL

See supplementary material for the superconducting transition of a bulk Pb wire recorded at ambient pressure and under a pressure of ca. 7.7 kbar, in water and Daphne oil and for Heat capacity profiles for Daphne 7373.

Bk

  is Boltzmann constant (=1.38064852(79)×10 -23 J⋅K -1 ), g is the degeneracy ratio between HS and LS energy levels (>1), ij J stands for the short-range interaction parameter representing the cooperative interaction that only exists between the nearest neighbouring pairs ( ij J includes the nearest neighbors number), , ij  denotes the fact that the sum is only made over the nearest neighbors pairs,

10 -

 10 stands for the conversion of the pressure units from Pascal to bar), and the high spin fraction ( 1) / 2 HS nm =+ .

Figure 1

 1 shows the temperature dependence of the high spin fraction, nHS, simulated for various values of the pressure gradient and different values of the effective interaction parameter.

Figure 1 .

 1 Figure 1. Thermal hysteresis simulated for various values of the pressure gradient and different values of the effective interaction parameter: (a) 0 300 K = ; (b) 0 400 K = and (c) 0 500 K = respectively.

Figure 2 .

 2 Figure 2. (a) Temperature dependence of the magnetic susceptibility recorded on [Fe(PM-BiA)2(NCS)2] polymorph I, at various pressures. (b) Simulated thermal hysteresis curves for different pressure gradients.

  can be plotted on a regular grid. The FORC distribution ( ) derivative of the fitted surface and it can be assigned to the center of grid as a representation of the density of the FORC distribution ( ),ab TT  at that point. A FORC diagram is a contour plot of a FORC distribution with Ta=Tup, and Tb=Tdown on the vertical and horizontal axes, respectively.

Figure 3

 3 Figure3shows the FORC diagram simulated for various pressure gradients. As it can be seen, the FORC diagrams

Figure 3 .

 3 Figure 3. FORC diagrams simulated for various values of the pressure gradient. The following computational

Figure 4 .

 4 Figure 4. Pressure gradient influence on the statistical parameters obtained from the statistical analysis of the simulated FORC diagrams, for various values of the effective interaction parameter.

Figure 5

 5 shows the effect of small pressure gradients on the thermal hysteresis loop and on the FORC diagram, respectively. If the major hysteresis loop is only slightly affected by a pressure gradient of ( ) 30 K  = , the behavior inside the thermal hysteresis loops, expressed by the FORC diagram, suffers strong modifications even in the presence of small pressure gradients.

Figure 5 .

 5 Figure 5. Comparative analysis of the pressure gradient effect on the major thermal hysteresis loop (a) and on the FORC diagram (b-c) simulated for various values of the pressure gradient. The following computational parameters have been used: 0 1200 K = , 0 300 K = and ln 7 g = . 4. Conclusions

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582

Acknowledgments

This work was partially founded by the European Commission through the SPINSWITCH project (H2020-MSCA-RISE-2016, Grant Agreement No. 734322). The financial support of EXCALIBUR project Contract

No. 18 PFE/16.10.2018 is also acknowledged.

Chemistry, 2011 (2011) 3193-3206. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0042582

Conflicts of Interest:

The authors declare no conflict of interest.


André Quinquis
email: quinquis@ensieta.fr

Cornel Ioana
email: ioanaco@ensieta.fr

Emanuel Radoi

POLYNOMIAL PHASE SIGNAL MODELING USING WARPING-BASED ORDER REDUCTION

The high-order ambiguity function (HAF) was introduced for the estimation of polynomial-phase signals (PPS). Currently the HAF suffers from noise-masking effects and from the appearance of undesired cross terms in the presence of multi-components PPS. The multi-lag product HAF concept was then proposed as a way to improve the performances of the HAF. Nevertheless, performances of the new methods are affected by the error propagation. This effect is due to the technique used for polynomial order reduction, common for current approaches : signal multiplication with the complex exponentials formed with the estimated coefficients.

In this paper, we introduce an alternative method to reduce the polynomial order, based on the successive unitary signal transformation, according to each polynomial order. We will prove that this method considerably reduces the effect of error propagation.

INTRODUCTION

It is well known that there is no transformation from the Cohen's class which can produce the complete concentration along the instantaneous frequency law (IFL) when this one is a nonlinear function of time. Therefore, different high order distributions have been developed in order to better match the non-linear timefrequency behavior of the analyzed signal [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF], [START_REF] Porat | Digital Processing of Random Signals[END_REF], [START_REF] Boashash | Polynomial Wigner-Ville distributions and their relationship to time-varying higher order spectra[END_REF], [START_REF] Stankovic | L-class of time-frequency distributions[END_REF]. For example, the polynomial phase signal constitutes a good model in a variety of applications (e.g. Radar, Communication, etc ) [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF].

One of the first approaches to estimate the parameters of the PPSs [START_REF] Porat | Digital Processing of Random Signals[END_REF] provides good results for high signal-tonoise ratio (SNR). Nevertheless, since the HAF is a nonlinear method, it suffers from three basic problems : 1) noise-masking effects for low SNR, 2) cross terms in the presence of multi-component PPSs (mc-PPSs) and 3) the propagation of the approximation error from an order to other.

Recently, different methods have been proposed in order to eliminate the first two limitations. The key point is to use the multi-lag concept in the HAF computing procedure [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF]. Moreover, multiplying the HAFs obtained for some lag sets (product HAF -PHAF), the performances related to 1) and 2) are considerably improved with respect to multi-lags HAF [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF].

Nevertheless, the effect of propagation error remains a serious limitation when we try to estimate a deeply nonlinear IFL (underwater transitory signals, digital modulation, etc). Therefore, we propose a new procedure for polynomial order compensation, based on the recursive signal warping. In fact, using the order reduction property of the warping technique, the idea is to iteratively apply this method to reduce the phase orders.

This paper is organized as follows. In section 2, we present the Product HAF method. The major limitation of the PHAF, related to the error propagation effect, is described in section 3. Afterwards, using the warping technique [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF], a new method for order compensation is proposed in section 4. Some examples will be presented in section 5. We will finally present some remarks in "Conclusion and Perspective" (section 6).

PRODUCT HIGH-ORDER AMBIGUITY FUNCTION

As it was illustrated in [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF], [START_REF] Porat | Digital Processing of Random Signals[END_REF], the classical HAF algorithm presents some limitations, related to the noise robustness and the cross-terms presence. In order to solve these aspects, the multi-lag HAF (mlHAF) concept has been initially proposed in [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF]. In fact, the mlHAF is based on the generalization of the high order instantaneous moment HIM [START_REF] Porat | Digital Processing of Random Signals[END_REF] :

( ) ( ) ( ) * 1 1 1 2 1 1 2
; ; ;

K K - K K K - K K K - HIM s t HIM s t HIM s t τ τ - - - - = + -             τ τ τ (1)
where

( ) i τ τ τ ,..., , 2 1 = i τ
is the lag set. Applying the Fourier transform to (1), we obtain the ml-HAF of the signal s(t) : [ ] ( )

; , ;

j t K K m lH A F s H IM s t e dt α α ∞ - -∞ =     ∫ τ τ (2)
The main property of HIM is that, assuming a PPS model for the analyzed signal, i.e.

( ) ( )

0 exp exp K k k k s t A j t A j a t φ =   = =     ∑ (3)
the K th order HIM is reduced to a harmonic with amplitude

( ) ( ) 2 1
;

e x p ω φ - = ⋅+     % % k k k k HIM s t A j t τ (4)
where

k K k a k 1 ! ~-= τ ω (5).
Based on these results, Porat [START_REF] Porat | Digital Processing of Random Signals[END_REF] has proposed an algorithm which estimates sequentially the coefficients {a k }. At each step, using a spectral analysis method, we estimate the spectral peak and, using the HAF, we compute an estimation value ( k a ˆ) of a k . Using this value, the effect of the phase term of the higher order is removed:

( ) ( ) ( ) ( ) 1 ˆk k k k ja t s t s t e - - = ⋅ (6)
Using the ml-HIM concept (relation ( 1)), Barbarossa and al [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF] introduced the Product HAF: the mlHAFs computed, via relation [START_REF] Porat | Digital Processing of Random Signals[END_REF], for different lag sets

( ) ( ) { } 1, ; 1, 1 l L l l K 1 K 1 i i K τ = - -     =       = = - T τ τ (7)
are multiplied, obtaining also a more robust method and a cross-term free representation :

(

; ,

K l i L i K K l i i l K PHAF f mlHAF s f τ τ - = - = = -       =       ∏ ∏ ∏ T τ (8)
The simulation results given in [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF] proved that the PHAF solves both noise robustness and ambiguity problems, providing also a correct IFL estimation. Nevertheless, the analyzed signals were characterized by smooth time-frequency behavior (usually, 3 th order PPs). If this condition is not verified, one of the major limitation of the PHAF based approaches, related to the error propagation phenomenon, acts. This phenomenon is studied in the next section.

ERROR PROPAGATION IN POLYNOMIAL PHASE MODELING

Let consider the signal given in (3) and we denote with K a ˆ the estimation of the k th order polynomial coefficient.

In real applications [START_REF] Porat | Digital Processing of Random Signals[END_REF], since a spectral estimation of a discrete sequence is involved, this value differs from the theoretic one by

K K K a a - = ε
-the approximation error. This error is directly related to the number of points in Fourier transform and to the SNR ( [START_REF] Barbarossa | Product High-Order Ambiguity Function for Multicomonent Polynomial-Phase Signal Modeling[END_REF]). Using this estimate, we remove, via (6), the k order phase component

( ) ( ) (
)

1 1 0 0 1 K K k k k k K k k K j at a a t j a t K j t K K K s t Ae Ae e ε - - = =     +   -   - ∑ ∑ = = ⋅ (9)
The (k-1) th order HIM of s (k) (t) is expressed as

( ) ( ) ( ) 2 1 1 1 ! 1 1 ; K K K j K a t K K HIM s t e τ φ τ - - -   - + -   -  =   (10)
This expression is a consequence of the main HAF property previously presented : the (k-1) th order HIM of a PPS is a sinusoid with an angular frequency related, via [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF], to the (k-1) th order polynomial coefficient. Nevertheless, due to the measurement error existing for the k th order, 1 ˆK a -is not the correct value of the (k-1) th order polynomial coefficient of signal s (k-1) . In order to find the relation between the errors at order k and (k-1), we can evaluate ( )

1 1 ;τ - -    K K HIM s
using the recurrent definition given in (1) and the relation ( 4) and [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF]. We obtain the following expression :

( ) ( )( ) 2 1 1 1 ! 1 1 ; τ ετ τ - - -   - + + -   -  =   K K K K j K a K t R K K HIM s e (11)
where R is a residual term which does not depend on t.

Since the mlHAF evaluation supposes the computation of Fourier transform of HIM, we focus only on the coefficient of t. Therefore, if we compare the corresponding terms which appear in (10) and ( 11) we obtain the following relation between the errors for two consecutive orders :

2 1 2 1 1 τ ε τ τ - - - - - + = K K K K K K a K a (12)
and, with the notation

1 1 1 ε - - - = - K K K
a a , we obtain :

1 K K K ε ε τ -= (13)
Furthermore, using the optimal value of the lag proposed in [START_REF] Porat | Digital Processing of Random Signals[END_REF] i.e. K N / = τ (N-sample number), we get:

1 K K N ε ε -= (14)
This relation shows that the error existing at a given order is transmitted at the inferior order by multiplication of N. The next figure illustrates this dependence for N=10. From this figure, it can be observed that even if the measurement error for the highest order is insignificant, its effect through the lower orders becomes deeply disturbing. It explains why the error propagation effect does not affect the polynomial estimation when a small approximation order is required (3 or 4). Nevertheless, there are many situations which impose a high Error K=8, N=10 Orders approximation order : digital modulations, transitory signals, etc. One example is given in figure 2 where we process, via PHAF-based phase modeling method, a frequency shift keying (FSK) modulation. The theoretical IFL is plotted in figure 2.b. Note that the SNR is about 30 dB. The PHAF-based estimation procedure is applied, starting with the order 9. At the superior orders (9, 8, 7) PHAF performs quite well : the propagation error is insignificant, but its effect is accumulated and it becomes disturbing for inferior orders (down to 6). This can be observed in the figure 2.a. , where the peak locations give information about the estimation quality. Consequently, the estimation of the polynomial coefficients is not correct (figure 2.b); the evaluated IFL does not match the correct time-frequency behavior of the FSK modulation.

A potential solution to reduce the error propagation is proposed in next section.

WARPING-BASED PHASE ORDER REDUCTION

In time-frequency analysis, the warping operator principle has recently become a very useful concept for generation of non-linear time-frequency representation (TFR) [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF]. Mathematically, a warping operator is defined as a unitary transform U on l 2 (ℜ), whose effect is x-axis deforming :

( )( ) ( ) ( ) 1/ 2 ' s x w x s w x =     U (15)
where w is a smooth, one-to-one function, comprising a large subclass of unitary transformations ( [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF]). Generally, these functions, called warping functions, are chosen to ensure the "linearization" of signal time-frequency behavior. The design of a warping function associated to a signal is possible if the signal modulation nature is known. [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF]. This is not always possible, especially in the case of passive systems. For this case, the polynomial phase modeling seems to be a potential solution for signal IFL describing.

To reduce the error propagation effect, we propose the replacement of this technique with a warping-based method. More precisely, we define a particular warping operator as :

( ) ( ) K K K w a t t t t / 1 ˆ :         = = →  K w K w w K (16)
where ˆk a is the k th order estimated polynomial coefficient.

The following example illustrates the effect of this warping operator for a 3 th order PPS given by : ( ) The Wigner-Ville Distribution (WVD) of this signal is depicted in figure 3. Using the warping operator obtained by the particularization of the relation (16) (for k=3 and Therefore, in order to reduce the polynomial order we apply, in an iterative manner, a warping operator designed as shown in (16). The effect of this operator on the PPS defined in (3) is depicted by : Therefore, the result of the warping transform of a K th order PPS consist in a (K-1) th order PPS for new temporal variable ( ) t w K . The (K-1) th order PHAF of this signal, with respect to variable ( ) t w K , peaks to a frequency location related, via relation [START_REF] Baraniuk | Unitary Equivalence : a new twist on signal processing[END_REF], to the a K-1 coefficient. Once a K-1 estimated, we construct the (K-1) th order unitary operator U K-1 as indicated in (16). The process is iterated until all polynomial coefficients are estimated.

( ) ( ) ( ) ( ) 1/ 1 0 exp exp ˆ-=               = ⋅                       ∑ K K K m K K w K m w m K t y t A ja j a t a K U ( ) ( ) th

A

  , frequency ω ~ and phase φ ~:

Figure 1 .

 1 Figure 1. Error propagation effect.

Figure 2 .

 2 Figure 2. PHAF-based phase modeling of a FSK signal

  a linear time-frequency structure (i.e. a 2 th order PPS) as illustrated in figure3.

Figure 3 .

 3 Figure 3. Order reduction using the warping operator (16)

As the relation (18) shows, the (K-1) th coefficient of the signal (U K y) depends only on the a K-1 coefficient, avoiding also the process of error propagation to the lower orders. This will be practically proved in the next section.

SIMULATION RESULTS

In this section, we firstly give the performances in terms of estimated variances as a function of the SNR. We assumed a 3 th order PPS embedded in white gaussian noise. Two methods were compared -classical PHAF estimation method and PHAF-based estimation method with warping-based phase compensation (denoted "WarpComp" method). Each variance was computed for 500 trials and, for each order, it was compared with the Cramer-Rao bound (CRB), theoretically evaluated in [START_REF] Porat | Digital Processing of Random Signals[END_REF]. The first plot proves that, for the highest order, the performances of both methods are similar : the estimation of the highest order coefficient depends on the noise only. The next two pictures show that, using the warping-based phase compensation, the estimation performances remains close to the CRB as in the case of the highest-order coefficient. Consequently, the performances of this method depend only on the noise, whereas in the PHAF case they are affected also by the error propagation phenomenon.

The error propagation effect reduction is also illustrated in the next figure, using the signal presented in the figure 2. This figure shows, the proposed method provides a much more accurate estimation of the IFL of a FSK signal. It is illustrated by an existence of a single peak at each polynomial order. Consequently, the polynomial order reduction through the procedure based on the warping operators (section 4) considerably improves the performances of the PHAF-based estimation procedure.

The statistical analysis and the result for a FSK signal prove that this method can successfully be used to estimate the polynomial model of a general class of signals.

CONCLUSION AND PERSPECTIVES

In this paper we have proposed a method, for polynomial order reduction, based on the warping principle, applied to each polynomial order. This method constitutes an attractive way to attenuate the effect of error propagation which inherently appears in any estimation algorithm of the PPSs. As shown by the results, this method, associated with a modern procedure for polynomial phase modeling, provides accurate timefrequency information about the analyzed signal.

In further works, we intend to apply this procedure in real contexts, such underwater signal processing or digital modulation recognition.


Ahed Alboody

Matthieu Puigt

Gilles Roussel

Vincent Vantrepotte

Cédric Jamet

Trung Kien

EXPERIMENTAL COMPARISON OF MULTI-SHARPENING METHODS APPLIED TO SENTINEL-2 MSI AND SENTINEL-3 OLCI IMAGES

Keywords: Image fusion, Remote sensing, Sentinel-2 MSI, Sentinel-3 OLCI, Simulations, Real data

Multi-spectral images are crucial to detect and to understand phenomena in marine observation. However, in coastal areas, these phenomena are complex and their analyze requires multi-spectral images with both a high spatial and spectral resolution. Unfortunately, no satellite is able to provide both at the same time. As a consequence, multi-sharpening techniques-a.k.a. fusion or superresolution of multi-spectral and/or hyper-spectral images-were proposed and consist of combining information from at least two multi-spectral images with different spatial and spectral resolutions. The fused image then combines their best characteristics. Various methods-based on different strategies and tools-have been proposed to solve this problem. This article presents a comparative review of fusion methods applied to Sentinel-2 MSI (13 spectral bands with a spatial resolution ranging from 10 to 60 m) and Sentinel-3 OLCI (21 spectral bands with a spatial resolution of 300 m) images. Indeed, both satellites are extensively used in marine observation and, to the best of the authors' knowledge, the fusion of their data was partially investigated (and not in the way we aim to do in this paper). To that end, we provide both a quantitative analysis of the performance of some state-of-the-art methods on simulated images, and a qualitative analysis on real images.

INTRODUCTION

The satellite observation of our planet knew significant instrumental advances for several decades, with consequent developments in terms of spatial resolution-e.g., in water remote sensing with high spatial resolution (10-60 m)-and in terms of spectral resolution (hyperspectral imagery). However, the Signal-to-Noise Ratio (SNR) of a Multi-Spectral or Hyper-Spectral Imaging (MSI/HSI) sensor is proportional to the ratio between the sensor area and the number of observed spectral bands. Therefore, to maintain a constant SNR, increasing the number of spectral bands in an hyper-spectral image implies a decrease in spatial resolution. As a consequence, our planet is currently observed by MSI systems having a very good spatial resolution but a low spectral resolution and by HSI systems having a very good spectral resolution but a low spatial resolution.

In remote sensing applied to water color, Sentinel-2 MSI and Sentinel-3 Ocean and Land Color Instrument (OLCI) are extensively used. As shown in Fig. 1-which exhibits the spectral band positions and widths of both instruments-for the considered application, they have 13-with 10, 20, or 60 m spatial resolution-and 21

This work is partially funded by CNES within the TOSCA "OSYNICO" project. Experiments presented in this paper were carried out using the CAL-CULCO computing platform, supported by SCoSI/ULCO. usable spectral bands-with 300 m spatial resolution-in the visible and near-infrared range, respectively. Both instruments allow to map several variables describing the biogeochemical dynamics of the marine environment, e.g., Chlorophyll A, Suspended Matter, or Organic Carbon. However, the study of specific environments-such as coastal or estuarine areas-requires observations at a sufficiently fine spatial scale and, at the same time, at a fine spectral resolution, which is not possible using Sentinel-2 or Sentinel-3 alone.

Multi-sharpening [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF] or multi-sensor image fusion is the process of combining relevant information from two or more MSI/HSI images into a single image with complementary spatial and spectral resolution characteristics [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF][START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF]. To that end, several methods have been proposed in the literature and are based on, e.g., component substitution-using, e.g. Gram-Schmidt [START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF], the Sylvester equation [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF], or Principal Component Analysis [START_REF] Kwarteng | Extracting spectral contrast in landsat thematic mapper image data using selective principal component analysis[END_REF]-coupled Nonnegative Matrix/Tensor Factorization (NMF/NTF)-e.g., [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF][START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF][START_REF] Li | Fusing hyperspectral and multispectral images via coupled sparse tensor factorization[END_REF]and more recently on deep learning, e.g., in [START_REF] Dong | Image superresolution using deep convolutional networks[END_REF]. However, to the best of the authors knowledge, fusing Sentinel-2 MSI and Sentinel-3 OLCI images was not yet considered in the literature. Indeed, even if in [START_REF] Korosov | Fusion of data from Sentinel-2/MSI and Sentinel-3/OLCI[END_REF], the authors suggested a method for the fusion of data from Sentinel-2 and Sentinel-3, they proposed in practice an approach to combine Sentinel-2 and Moderate Resolution Imaging Spectroradiometer (MODIS) images. The fusion of Sentinel-2 and Sentinel-3 data is let for future work in the conclusion of [START_REF] Korosov | Fusion of data from Sentinel-2/MSI and Sentinel-3/OLCI[END_REF].

Moreover, the authors of [START_REF] Wang | Spatio-temporal fusion for daily sentinel-2 images[END_REF] considered time series of Sentinel-2 and Sentinel-3 images to create daily Sentinel-2 data, i.e., they do not aim to create new multi-spectral images combining the best characteristics of Sentinel-2 and Sentinel-3. On the contrary, they use the same spectral bands in Sentinel-2 and Sentinel-3 in order to generate a better spatial resolution of Sentinel-3 images. As a consequence, the objective of these authors slightly differs from ours.

As Sentinel-2 has several bands with coarse to fine spatial resolutions, the super-resolution of Sentinel-2 images has also been considered, e.g., in [START_REF] Palsson | Sentinel-2 image fusion using a deep residual network[END_REF]. The proposed approach is based on a deep convolutional neural network including a residual neural network for the fusion problem modeling. However, taking into account Sentinel-3 data was not considered by the authors of [START_REF] Palsson | Sentinel-2 image fusion using a deep residual network[END_REF].

As a consequence, the goal of this paper is to provide a first comparison of the performance of state-of-the-art multi-sharpening techniques applied to Sentinel-2 and Sentinel-3 data. While an objective comparison of such methods on real data is challengingexcept if we aim to cross-validate their performance, as proposed in, e.g., [START_REF] Puigt | Cross-validation of blindly separated interstellar dust spectra[END_REF] for HSI unmixing-it remains possible on semi-realistic simulations. We thus aim to propose (i) a first quantitative analysis of the performance of some state-of-the-art techniques on simulated data and (ii) a qualitative analysis of their performance on real data.

The remainder of the paper reads as follows. In Section 2, we briefly introduce the tested methods in this paper. Section 3 (respectively, Section 4) then presents the simulated dataset (respectively, some real dataset) and the reached multi-sharpening performance. Lastly, we conclude and discuss about future work in Section 5.

TESTED MULTI-SHARPENING METHODS

In this paper, we evaluate the performance of image fusion methods of hyperspectral and multispectral images. To that end, we consider classical and modern multi-sharpening techniques, i.e., (i) based on component substitution-and particularly Gram-Schmidt (GS) [START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF] and Gram-Schmidt Adaptative (GSA) [START_REF] Aiazzi | Improving component substitution pansharpening through multivariate regression of ms + pan data[END_REF]-(ii) coupled Non-negative Matrix Factorization (NMF) [START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF], (iii) Fast fusion based on Sylvester equation (FUSE) [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF], (iv) Generalized Laplacian pyramid with hypersharpening (GLPHS) [START_REF] Aiazzi | MTF-tailored multiscale fusion of high-resolution MS and pan imagery[END_REF], (v) Smoothing filter-based intensity modulation with hypersharpening (SFIMHS) [START_REF] Liu | Smoothing filter-based intensity modulation: A spectral preserve image fusion technique for improving spatial details[END_REF], (vi) Lanara's work presented in ICCV'15 (ICCV'15) [START_REF] Lanaras | Hyperspectral super-resolution by coupled spectral unmixing[END_REF], and (vii) Maximum a Posteriori estimation with a Stochastic Mixing Model (MAP-SMM) [START_REF] Eismann | Resolution enhancement of hyperspectral imagery using maximum a posteriori estimation with a stochastic mixing model[END_REF]. All these algorithms are provided in the HSMSFu-sionToolbox1 [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF]. We decided not to take into account both the "ECCV'14" and "HySure" which were also provided in this toolbox. Indeed, we faced some parameter tuning and execution time issues with the former while the later could not provide any enhancement in our experiments that we do not show because of space constraints. All the tests presented in the remainder of the paper were performed using MATLAB R2019a on Windows 10 x64 machines on a laptop equipped with an Intel i9-9880H CPU @ 2.30GHz with 32.0GB RAM on Windows 10 x64.

The GS transformation is fast and easy to implement and generates fused images with high integration quality color and spatial detail. GS is a commonly used method [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF][START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF] when a low-spatialresolution image is sharpened by adding spatial details obtained by multiplying the difference between a high-spatial-resolution image and a synthetic intensity component by a band-wise modulation coefficient. The improvement lies in computing the synthetic intensity component by performing a linear regression between a highresolution image and lower-resolution bands to mitigate spectral distortion. GSA integrates this technique into the GS algorithm [START_REF] Aiazzi | Improving component substitution pansharpening through multivariate regression of ms + pan data[END_REF]. CNMF [START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF] unmixes the low-spatial-resolution HSI data and the highspatial-resolution MSI data alternately using the relation between sensor properties for the NMF initializations. CNMF alternately unmixes these two inputs images by NMF to estimate the spectral signatures of endmembers and the high-resolution abundance maps, respectively. CNMF starts by unmixing the low-resolution image using vertex component analysis (VCA) to initialize the endmember signatures. The final high-resolution HSI data is obtained as the product of the spectral signatures and the high-resolution abundance maps. FUSE combines a high-spatial low-spectral resolution image and a low-spatial high-spectral resolution image [START_REF] Wei | Fast fusion of multi-band images based on solving a sylvester equation[END_REF]. FUSE utilizes a Sylvester equation to solve the maximization problem of the likelihoods obtained from the forward observation models. A closed-form solution for the Sylvester equation improved computational performance. In GLPHS fusion, spatial details of each low-resolution band are obtained as the difference between a high-resolution image and its low-pass version multiplied by a global gain factor, which can be computed globally [START_REF] Aiazzi | MTF-tailored multiscale fusion of high-resolution MS and pan imagery[END_REF]. A Gaussian filter, matching the modulation transfer function (MTF) of a lower-resolution sensor, is used for low pass filtering. SFIMHS [START_REF] Liu | Smoothing filter-based intensity modulation: A spectral preserve image fusion technique for improving spatial details[END_REF] is based upon a simplified model for solar radiation and land surface reflection, SFIMHS sharpens the low-resolution image by multiplying an upscaled lower resolution image by a ratio between a higher resolution image and its low-pass filtered image (with a smoothing filter) on a pixel-by-pixel basis [START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF]. Spatial details can be modulated to a co-registered lower resolution multispectral image without altering its spectral properties and contrast. SFIMHS can be performed on individual HS bands. Similar to CNMF [START_REF] Yokoya | Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion[END_REF], Lanaras's algorithm jointly unmixes the two input images into the spectral signatures of endmembers and (pure reflectance spectra of the observed materials) the associated fractional abundances (mixing coefficients). A projected gradient method was proposed to alternately update the endmember signatures and the high-resolution abundances by solving the two unmixing problems of the input HS-MS images, respectively. Simplex identification via split augmented Lagrangian (SISAL) is used to initialize the endmembers and sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL) is adopted to obtain initial abundances. The formulation leads to a coupled matrix factorization problem, with a number of useful constraints imposed by elementary physical properties of spectral mixing. The MAP-SMM algorithm [START_REF] Eismann | Resolution enhancement of hyperspectral imagery using maximum a posteriori estimation with a stochastic mixing model[END_REF] adopted an SMM to estimate the underlying spectral scene statistics or, more specifically, the conditional mean vector and covariance matrix of the high-resolution HS image with respect to the MS image. The average spectrum, covariance matrix, and abundance map of each endmember are derived from the low-resolution HS image. A MAP objective function is formulated to optimize the high-resolution HS data relative to the input images based on the SMM statistics. The MAP-SMM algorithm is performed in the principal component subspace of the low-resolution HS image.

It is worth mentioning that all the above state-of-the-art methods were proposed for the fusion of two multi-or hyperspectral images with different spatial and spectral resolutions. However, Sentinel-2 MSI data already consist of several spatial resolutions (10, 20, and 60 m). As a consequence, applying the above multi-sharpening techniques to Sentinel-2 and Sentinel-3 can only be performed for one given Sentinel-2 resolution. As a consequence, we will show below the fusion performance for the three spatial resolutions.

PERFORMANCE OF THE TESTED METHODS ON SYNTHETIC DATA

Many studies have been published during the last years focusing on single-date, intra-sensor data fusion. As presented in [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF][START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF], the researchers evaluate the performance of fusion algorithms mainly on images of the following data sets: AVIRIS Indian Pines, AVIRIS Cuprite, and ROSIS-3 University of Pavia [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF][START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF]. In our study, we aim to evaluate the performance of fusion algorithms on Sentinel-2 MSI and Sentinel-3 OLCI remote sensing images. This section thus provides a brief description of the synthetic dataset used in the experiments. To validate the feasibility of Sentinel-3 and Sentinel-2 image fusion, we propose to consider a purely linear model to generate synthetic Sentinel-3 and Sentinel-2 images using AVIRIS2 endmembers for Ocean water and Seawater from USGS spectral library3 [START_REF] Kokaly | USGS spectral library version 7[END_REF]. We then choose 5 endmembers from this library, each of them being related to an ocean or a seawater spectrum. From these endmembers, we extract 21 spectral bands (respectively, 13) which are the closest to Sentinel-3 (respectively, Sentinel-2) ones.

Then, for each endmember, we generate a synthetic abundance map using the Gaussian Fields method [START_REF]Hyperspectral imagery synthesis (EIAs) toolbox[END_REF]. Each of theses maps is of dimension 3000 × 3000 and is assumed to have a 10 m spatial resolution. Combining them with the above endmembers allows generate three high-resolution Sentinel-3-like images, i.e., a 3000×3000×21 datacube. We then apply a Gaussian filter to downsample the spatial content and we derive both a 1500 × 1500 × 21 and a 500 × 500 × 21 datacubes that we aim to estimate from their downsampled versions. We then generate the observed data by (i) downsampling-by a factor 30-the 3000 × 3000 maps by using a Gaussian filter, in order to derive the observed 100 × 100 × 21 synthetic Sentinel-3, with 300 m spatial resolution. We then derive the observed synthetic Sentinel-2 datacube by applying to the 3000 × 3000 abundance maps a downscale factor of two (with 20 m spatial resolution) and six (with 60 m spatial resolution) and by considering the corresponding spectral bands used in the real Sentinel-2 imagers, i.e., four bands with 10 m, six bands with 20 m, and three bands with 60 m spatial resolution of Sentinel-2.

In order to assess the performance of the tested methods, we use some classical quantitative performance measures [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF][START_REF] Jagalingam | A review of quality metrics for fused image[END_REF], i.e., (i) the Peak Signal-to-Noise Ratio (PSNR)-which is the ratio between the highest possible signal energy and the noise energy-(ii) the Spectral Angle Mapper (SAM)-which is a pixelwise measure of the angle between the reference spectrum and the fused one. SAM values near zero indicate local high spectral quality and we use the average SAM value with respect to pixels for the quality index of the entire data set-(iii) the ERGAS measure-i.e., a normalized average error of each band of processed image-and (iv) Q2n which is a generalization of the universal image quality index (UIQI) and an extension of the Q4 index to HS images based on hypercomplex numbers [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF].

As already explained, the tested state-of-the-art methods are designed to fuse two datacubes with two distinct resolutions. However, Sentinel-2 already has 3 different spatial resolution at different wavelengths. One must thus choose a target spatial resolution among the three available, i.e., 10, 20, or 60 m. Table 1 shows the obtained results for each of these resolutions. They show that FUSE is almost always outperforming the other techniques. Then-except for 60 m resolution where GLPHS outperforms FUSE-GLPHS, MAP-SMM, and SFIMHS provide some performance slightly lower than FUSE. Lastly, the other tested methods provide some performance significantly lower than the former. Let us emphasize again that these results were obtained with simulations in which the purely linear mixing model is assumed. However, such an assumption might not be satisfied in practice.

Method

PERFORMANCE OF THE TESTED METHODS ON REAL DATA

We now investigate the performance reached by the tested methods on real Sentinel-2 and Sentinel-3 datasets4 . The selected study area is located in the islands of The Bahamas in the Atlantic Ocean. The Tongue of the Ocean is a deep-water basin in the Bahamas that is surrounded to the East, West, and South by a carbonate bank known as the Great Bahama Bank. The deep blue water of the Tongue is a stark contrast to the shallow turquoise waters of the surrounding Bank. Generally, waters that are optically shallow (e.g., Grand Bahama Bank) appear blue-green due to high bottom reflectance contributions while optically deep waters appear dark blue. The centre of the Bahamas image is located at the following coordinates: Lat: 26 • 35'58.50"N, and Lon: 77 • 28'29.93"W (DMS), Projection UTM, Zone 18 N, and World Geodetic System 1984.

To assess the performance of the tested methods, as we do not have any ground truth information, we propose to compare their multi-sharpening outputs with those obtained using the ENVI implementation of GS [START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF]. To do so, we use the same performance criteria as in the previous section. However, these results must be carefully interpreted. Indeed, a good performance index-i.e., a high PSNR, a low SAM, a low ERGAS or a high Q2n-only indicates that the fusion result obtained with a tested method is close to the one reached by GS. On the contrary, the fact that one performance index is "poor" only states that the reached performance is different from GS. At this stage, we are not able to decide whether the performance is "better" or "worse" than with GS. tion, all the tested methods (except MAP-SMM for 10 m of spatial resolution) provide a quite similar performance with respect to these indices. We thus can state that they globally provide a quite similar performance. However, when we visually look at the fused images, one might see some specific behaviors, as we now discuss. Due to space constraints, we only show the performance reached for 60 m of spatial resolution. As the tested methods generate some datacubes which are not easy to draw in a compact way, we propose the following procedure. For each pixel, we average the amplitudes along the spectrum axis. We then derive an image that we draw according to the same scale. All these images are shown in Fig. 2. In particular, Fig. 2(a) represents the original 1830×1830×3 Sentinel-2 image that we aim to fuse with the 366 × 366 × 16 Sentinel-3 one plot in Fig. 2(b). Figures 2(c) to 2(i) show the obtained fused images with the different tested methods. One may notice that most of them provide some artefacts. In particular, FUSE seems to provide the worst visual performance. Indeed, one notice some visible artefacts providing a kind of texture. Moreover, it also provides some negative amplitudes which have no physical meaning. Then, please notice that the Sentinel-3 image seem to provide more spatial information that the Sentinel-2 one. This is probably due to the fact that, even if both images were taken the same day, they were not taken at the same time and that the considered Sentinel-2 data only contain 3 spectral samples. However, CNMF and GSA seem to provide a higher visual quality than the other methods. Both methods can also partially remove clouds from Sentinel-3 images. Then, both SFIMHS and MAP-SMM provide some similar images, with a visually good spatial information. Lastly, ICCV'15 seems to stick to Sentinel-2 spatial information. However, this is not the case when the spatial resolution is 10 m-not shown for space considerationwhere it provides the best visual information. Unfortunately, in the absence of ground truth, it is hard to provide more information about these results, which is a perspective discussed in the next section.

CONCLUSION AND PERSPECTIVES

In this paper, we presented a first comparative study of Sentinel-2 and Sentinel-3 image fusion methods. As the ratio between high and low spatial resolutions increases, and for higher spatial resolution (10m), the fusion performance decreases. According to the considered experiments, CNMF and ICCV'15 seem to be the more robust tested methods. In future work, we aim to propose a new method able to take into account all the bands of Sentinel-2 in order to provide a new multi-spectral image with 10-m spatial resolution and Sentinel-3 spectral bands. Moreover, we aim to compare the above methods and our future work with in situ measurements. To that end, atmospheric correction [START_REF] Steinmetz | Sentinel-2 MSI and sentinel-3 OLCI consistent ocean colour products using POLYMER[END_REF] will be applied to Sentinel-2 and 3 data.

Fig. 1 .

 1 Fig. 1. Examples of Sentinel-2 and Sentinel-3 reflectance spectra.

Fig. 2 .

 2 Fig. 2. Fused Sentinel-2 (with 60 m spatial resolution) and Sentinel-3 images, with several methods: (a) original Sentinel-2 image, (b) original Sentinel-3 image, (c) GSA, (d) FUSE, (e) SFIMHS, (f) CNMF, (g) GLPHS, (h) ICCV'15, (i) MAP-SMM.

Table 1 .

 1 Performance of the fusion methods on simulations.

		PSNR	SAM	ERGAS	Q2n
	Perf. obtained with 60 m spatial resolution
	GSA	11.3	11.8811	8.969	1
	FUSE	23.5	0.26654	1.6697	1
	CNMF	15.3	0.49605	4.3611	1
	GLPHS	23.7	0.26208	1.6329	1
	SFIMHS	22.8	0.28683	1.8134	1
	ICCV15	18.8	0.39035	2.8472	1
	MAP-SMM	23.2	0.28777	1.7412	1
	Perf. obtained with 20 m spatial resolution
	GSA	9.8	9.4837	2.9385	1
	FUSE	21.6	0.34359	0.70224	1
	CNMF	8.4	0.47334	3.4464	0.67949
	GLPHS	21.4	0.36691 0.71696	1
	SFIMHS	21.1	0.35972 0.74190	1
	ICCV15	18.2	0.45740	1.0388	1
	MAP-SMM	21.0	0.36554 0.75128	1
	Perf. obtained with 10 m spatial resolution
	GSA	17.4	0.8112	0.56697 0.99774
	FUSE	19.9	0.4321	0.42869	1
	CNMF	15.2	0.60021 0.74173	1
	GLPHS	19.3	0.48046 0.46072	1
	SFIMHS	19.6	0.44427 0.44506	1
	ICCV15	17.3	0.55533 0.57638	1
	MAP-SMM	19.2	0.43245 0.44757	1

Table 2

 2 summarizes the obtained results. One may see that, for each given spatial resolu-

	Method	PSNR	SAM	ERGAS	Q2n
	Perf. obtained with 60 m spatial resolution
	GSA	26.7	0.2518 20.9157 0.88583
	FUSE	25.8	0.2154 19.8935 0.81807
	CNMF	27.5	0.2046 16.7684 0.91749
	GLPHS	25.6	0.2166 21.3205 0.93843
	SFIMHS	25.6	0.2197 21.7276 0.87596
	ICCV15	28.5	0.1887	14.054	0.95898
	MAP-SMM	25.6	0.2196	21.48	0.91038
	Perf. obtained with 20 m spatial resolution
	GSA	23.9	0.2427	5.2917	0.86474
	FUSE	18.8	0.2500	5.6916	0.80446
	CNMF	21.7	0.2570	6.9271	0.75592
	GLPHS	23.2	0.2406	5.8331	0.90253
	SFIMHS	23.0	0.2420	5.9424	0.86716
	ICCV15	23.3	0.2278	5.6937	0.93979
	MAP-SMM	23.1	0.2465	5.8853	0.85204
	Perf. obtained with 10 m spatial resolution
	GSA	24.3	0.2550	2.6741	0.89361
	FUSE	21.4	0.2596	2.7337	0.86143
	CNMF	23.9	0.2492	2.7735	0.92034
	GLPHS	23.3	0.2495	2.9307	0.86243
	SFIMHS	23.2	0.2521	2.9953	0.86188
	ICCV15	23.9	0.2289	2.6807	0.92013
	MAP-SMM	14.7	0.2008	1.3722	0.99306

Table 2 .

 2 Reached performance wrt GS output on real data.

The toolbox may be downloaded at https:// openremotesensing.net/wp-content/uploads/2017/11/ HSMSFusionToolbox.zip.

See https://aviris.jpl.nasa.gov/.

This library is accessible at https://crustal.usgs.gov/ speclab/AV14.php.

They are publicly available at https://scihub.copernicus. eu/dhus/#/home.