Manuel Corral Valero

Bénédicte Prelot

Grégory Lefèvre
email: gregory.lefevre@chimie-paristech.fr

MUSIC Speciation of γ-Al 2 O 3 at the Solid Liquid Interface: How DFT Calculations can Help with Amorphous and Poorly Crystalline Materials

Keywords: Multisite complexation (MUSIC), Density Functional Theory (DFT), alumina, surface speciation, surface reactivity, Solid Liquid Interface, Ab Initio Molecular Dynamics (AIMD)

Transition aluminum oxides, such as γ-Al 2 O 3 or alumina, are widely used in many different technical applications that rely on the surface reactivity of this material at the solid liquid interface. The speciation of surface sites of this material confronts several obstacles. On the one hand, alumina is a poorly crystalline oxide, thus allowing for a limited amount of empirical structural information for an important number of surface sites with different trends in reactivity and, on the other hand, it is a metastable material. In this work, we show several ways in which the Multi Site Complexation Model (MUSIC), combined with atomistic information from Density Functional Theory and Ab Initio Molecular Dynamics, can manage to perform speciation calculations of γ-Al 2 O 3 surface sites at the solid liquid interface. Although the results are in good qualitative agreement with experimental titration curves, and they can serve as a guide for the interpretation of the reactivity of this material at the initial stages of an impregnation experiment and chemical weathering phenomena, this work highlights the need of more complex AIMD simulations to accurately model these phenomena in γ-Al 2 O 3 surface/liquid interfaces.

INTRODUCTION

Among the different transition aluminum oxides, γ-Al 2 O 3 -also referred to as alumina-is one of the most widespread materials in a variety of technical applications and, in particular, for heterogeneous catalysis supports [START_REF] Euzen | Handbook of Porous Solids[END_REF] and water treatment [START_REF] Kasprzyk-Hordern | Chemistry of Alumina, Reactions in Aqueous Solution and its Application in Water Treatment[END_REF] . Hence, predicting the retention of metals and organic species by this compound at the solid liquid interface (SLI) is of paramount importance.

Unfortunately, γ-Al 2 O 3 is also poorly cristallineand extremely difficult to characterize rigorously by most conventional analytical techniques. However, complex Density Functional Theory (DFT) calculations has been successfully performed to obtain a detailed atomic description of physicochemical phenomena (see Raybaud et al. [START_REF] Raybaud | From γ-Alumina to Supported Platinum Nanoclusters in Reforming Conditions: 10 Years of DFT Modeling and Beyond[END_REF][START_REF] Raybaud | First Principles Surface Thermodynamics of Industrial Supported Catalysts in Working Conditions[END_REF] for a review). Moreover, this is a metastable material and its bulk structure changes in presence of water [START_REF] Abi Aad | Chemical Weathering of Alumina in Aqueous Suspension at Ambient Pressure: A Mechanistic Study[END_REF][START_REF] Lefèvre | Hydration of γ-Alumina in Water and Its Effects on Surface Reactivity[END_REF] ; this process takes some time to be detected and the rate and extent to which surface groups are affected by this evolution is not fully understood. Besides, there is evidence that this process may be inhibited by the sorption of inorganic species [START_REF] Abi Aad | Inhibition by Inorganic Dopants of γ-Alumina Chemical Weathering Under Hydrothermal Conditions: Identification of Reactive Sites and their Influence in Fischer-Tropsch Synthesis[END_REF] . Despite these difficulties, there have been attempts to depict the structure and reactivity of γ-Al 2 O 3 surfaces at the SLI [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF][START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF][START_REF] Réocreux | Biomass Derivatives in Heterogeneous Catalysis: Adsorption, Reactivity and Support from First Principles[END_REF] from DFT based Ab Initio Molecular Dynamics (AIMD) simulations. These theoretical results provides insights about the properties of this material in contact with water at short time scales -before the subsequent transformation of alumina by water or its inhibition by the adsorption of a species. Therefore, these studies could provide a starting point for future detailed analysis of the evolution of alumina in contact with water and its inhibition by the adsorption of water at the atomic scale (as in, for example, reference [START_REF] Réocreux | Biomass Derivatives in Heterogeneous Catalysis: Adsorption, Reactivity and Support from First Principles[END_REF] ).

When it comes to the reactivity of amorphous or poorly cristalline materials in water, effective 1-pKa and 2-pKa models are used [START_REF] Bourikas | The Mechanism of the Protonation of Metal (Hydr)Oxides in Aqueous Solutions Studied for Various Interfacial/Surface Ionization Models and Physicochemical Parameters: a Critical Review and a Novel Approach[END_REF][START_REF] Kummert | The Surface Complexation of Organic Acids on Hydrous γ-Al 2 O 3[END_REF] . This allows for a rationalization of sorption experiments, although at the cost of a detailed description of phenomena at the atomic scale. In the case of alumina, works from Raybaud et al. shows that γ-Al 2 O 3 surfaces present about 3 classes of sites with different reactivity [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF][START_REF] Digne | Use of DFT to Achieve a Rational Understanding of Acid-Basic Properties of Gamma-Alumina Surfaces[END_REF] . Hence, atomistic models for this material, at the SLI and based on first principles approaches, could be a starting point for rationalizing their reactivity.

The so-called MUSIC model from Hiemstra and co-workers (see [START_REF] Hiemstra | Multisite Proton Adsorption Modeling at the Solid/Solution Interface of (Hydr)Oxides: A New Approach: II. Application to Various Important (Hydr)Oxides[END_REF][START_REF] Hiemstra | Multisite Proton Adsorption Modeling at the Solid/Solution Interface of (Hydr)Oxides: A New Approach: I. Model Description and Evaluation of Intrinsic Reaction Constants[END_REF] for the original model and [START_REF] Hiemstra | A Surface Structural Approach to Ion Adsorption: The Charge Distribution (CD) Model[END_REF][START_REF] Hiemstra | Intrinsic Proton Affinity of Reactive Surface Groups of Metal (Hydr)oxides: The Bond Valence Principle[END_REF] for refinements of this work) allows for a prediction of the protonation constant (pKa) of surface OH groups on metal (hydr)oxides based on atomic scale information, mainly the number of bonds around surface species and their effective charge. Thus, MUSIC-based pKa constants in conjunction with an electrical double layer model has been successfully used to simulate the surface charge of complex metal oxides in water [START_REF] Tang | Prediction of Isoelectric Point of Manganese and Cobalt Lamellar Oxides: Application to Controlled Synthesis of Mixed Oxides[END_REF] . The bottleneck for the application of this approach to an amorphous or poorly crystalline material remains the lack of information about the structural parameters of the involved surfaces sites. In literature, this problem has been resolved by a so-called pragmatic approach, where the authors described the mean surface reactivity with only one surface site (>SO 1/2-in its deprotonated state) and one equilibrium reaction (>SO 1/2-+ H + = >SOH 1/2+ ) with a pK equal to the Point of Zero Charge (PZC) [START_REF] Mayordomo | Selenium(IV) Sorption Onto γ-Al2O3: A Consistent Description of the Surface Speciation by Spectroscopy and Thermodynamic Modeling[END_REF] . Moreover, Hiemstra et al. [START_REF] Hiemstra | Multisite Proton Adsorption Modeling at the Solid/Solution Interface of (Hydr)Oxides: A New Approach: II. Application to Various Important (Hydr)Oxides[END_REF] have described the surface reactivity with the singly coordinated AlOH found on gibbsite (pK 1,2 = 9), associated to a >SO 1/2-with a pK equal to 7 (chosen from the experimental value of the PZC). In the work by Nagashima and Blum

21

, authors assumed that only Al VI are present at the surface and fit the acid-base titration curves to obtain the densities of the surface groups. Thus, in these works, even if crystallography-based surface models have been used, the authors have not followed a prediction approach, whose ability is claimed in the MUSIC model, but have rather used them as a starting point of a fitting approach. In this work we propose to use the DFT model of γ-Al 2 O 3 surfaces developed by Raybaud et al. in [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF][START_REF] Digne | Use of DFT to Achieve a Rational Understanding of Acid-Basic Properties of Gamma-Alumina Surfaces[END_REF] conjunction with the MUSIC approach for the prediction of alumina's pKa values. To the best of our knowledge, there are few examples in the literature combining DFT and MUSIC calculations and amongst them we may cite those of Kubicki and co-workers [START_REF] Machesky | Surface Protonation at the Rutile (110) Interface: Explicit Incorporation of Solvation Structure within the Refined MUSIC Model Framework[END_REF][START_REF] Ridley | Anatase Nanoparticle Surface Reactivity in NaCl Media: A CD-MUSIC Model Interpretation of Combined Experimental and Density Functional Theory Studies[END_REF] . In reference [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF] DFT calculations were successfully used to depict the reactivity of alumina's surface sites at the solid-gas interface. Moreover, AIMD simulations has been performed in reference [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] in order to obtain a description of this material's structure at the solid water interface. In what follows, we are going to use this information [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] , mainly the coordination number for each OH site, the average number of surrounding hydrogen bonds and an estimate of their surface densities. Protonation constants obtained by these means are used with the Basic Stern double layer model and the speciation software PHREEQC [START_REF] Parkhurst | User's guide to PHREEQC (version 2)-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF] in order to compute the surface charge at the SLI as a function of pH; theoretical results are compared with experimental titration curves. [START_REF] Machesky | Surface Protonation at the Rutile (110) Interface: Explicit Incorporation of Solvation Structure within the Refined MUSIC Model Framework[END_REF] Although results from AIMD simulations in references 8-10 could be considered as a snapshot of γ-Al 2 O 3 before the transformations induced by water described [START_REF] Abi Aad | Chemical Weathering of Alumina in Aqueous Suspension at Ambient Pressure: A Mechanistic Study[END_REF][START_REF] Lefèvre | Hydration of γ-Alumina in Water and Its Effects on Surface Reactivity[END_REF] , the results we present here provides insights about the initial stages of this process and opens perspectives for future and more complex works in this area.

MATERIALS AND METHODS

The experimental surface charge used for comparison was obtained from titration methods [START_REF] Ahmad | Sorption Behavior of Co(II)-γ-alumina System in the Presence of Organic Additives[END_REF] .

Measurements were carried out using homemade titration set-up composed of an auto-burette (ABU901) ionometer-pHmeter (PHM250) equipped with GK2401 pH electrode (all from Radiometer), these devices being controlled by the TestPoint procedure. The experimental procedure was adapted from Ali Ahmad et al. [START_REF] Ahmad | Influence of Morphology and Crystallinity on Surface Reactivity of Nanosized Anatase TiO2 Studied by Adsorption Techniques. 2. Solid-Liquid Interface[END_REF] . Alumina (0.1 g of γ-Al 2 O 3 ) was first placed in contact with a 50 mL NaNO 3 background electrolyte. After acidifications, the suspension was equilibrated during 2 h. The titrations were carried out in the acid-to-base direction by adding a 0.02 mol L -1 NaOH titrant. The Background electrolyte was used as reference, with concentration as follows: 10 -1 , 10 -2 , and 10 -3 mol.L -1 NaNO 3 . Experiments were performed in a thermostatic titration cell regulated at 298 K by a thermocryostat and kept under inert gas stream. The surface charge density σ 0 was calculated using the following conventional procedure: (1) where ∆V is the difference between the volume of titrant added to the suspension (V s ) and that added to the reference phase (V e ) (ΔV = V s -V e ), F is the Faraday's constant (96500 C.mol -1 ), C denotes the titrant concentration (mol L -1 ); m and S refer to the mass (g) and the specific area of the solid sample (m 2 g -1 ), respectively. Simulated surface charge for different representations of the particle surfaces have been calculated using the PHREEQC software [START_REF] Parkhurst | User's guide to PHREEQC (version 2)-A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF] following the scheme below (Fig. 1). Fig. 1. Scheme of the method followed to simulate the surface charge curves. The inputs coming from DFT calculations are illustrated in red. See text for details.

Several inputs are needed to make this simulation. First, we considered the two most important orientations of alumina's nanocrystallites, the (100) and (110) faces. The atomistic model of these faces is inherited from the bulk model of γ-alumina developed by Krokidis and co-workers in ref. [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to gamma-Alumina[END_REF] In this work, the authors modeled the topotactic transformation of boehmite (AlOOH) taking place under calcination and leading to γ-alumina. The authors in this paper tested the validity of their model by comparing, on the one hand, experimental and computed XRD patterns and, on the other hand, the Al IV -Al VI distributions in their model with [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to gamma-Alumina[END_REF] Al NMR experiments from previous publications cited therein. The topotactic transformation taking place during boehmite calcination was first reported by Lippens et al. in [START_REF] Lippens | Study of Phase Transformations During Calcination of Aluminum Hydroxides by Selected Area Electron Diffraction[END_REF] based on selected area electron diffraction results. In the latter reference, the authors show that exposed faces of boehmite transform in certain faces of γ-alumina.

The face indexation in [START_REF] Lippens | Study of Phase Transformations During Calcination of Aluminum Hydroxides by Selected Area Electron Diffraction[END_REF] has been used in [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF][START_REF] Digne | Use of DFT to Achieve a Rational Understanding of Acid-Basic Properties of Gamma-Alumina Surfaces[END_REF] by Digne and co-workers, where the bulk model from reference [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to gamma-Alumina[END_REF] is used to build the surface models of the hydroxylated (100) and (110) γ-alumina faces at different water coverages (corresponding to different temperature and water partial pressure exposures); the validity of these surface models was tested by comparing theoretical OH stretch frequencies with experimental FT-IR experiments. In ref. [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] we used the latter results in order to model the reactivity of this material at the SLI. Inspired by the approach of Digne et al., in that paper we also computed the stretching frequencies of alumina's OH groups in presence of water from the velocity autocorrelation function obtained from AIMD simulations.

Both the (100) and (110) orientations account for more than 90 % of the specific surface area and, following the approach in [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF] to describe the reactivity at the solid-gas interface, we considered that each accounts for 20 % and 80 % of the specific area, respectively. Second, we used the detailed structural information of surface groups from DFT calculations in ref. [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] . Here, we briefly summarize the methods used in that work. The average structure of the alumina/water interface was obtained from periodic AIMD simulations in the microcanonical ensemble and a time step of 0.5 fs; kinetic energies were rescaled to obtain average temperatures close to 300 K. These simulations were performed within the framework of the Born-Oppenheimer approximation with the CP2K/Quickstep package based on a hybrid Gaussiand and Plane Wave method described in [START_REF] Vandevondele | Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach[END_REF] , the PBE DFT functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF][START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] with Grimme D2 method for dispersion forces 32 , Goedecker-Teter-Hutter pseudopotentials [START_REF] Goedecker | Separable Dual-Space Gaussian Pseudopotentials[END_REF] and double-ζ plus polarization Gaussian basis sets for all elements. For the plane waves region, an energy cutoff of 400 Ry was used. Considering the results in [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] (in which different approaches to model the SLI interface for thoses system are compared), we chose to use larger simulation boxes - 110) facets, respectively -and shorter simulation times (10 ps after equilibration) in order to obtain structural statistics about the structure of the different surface OH groups (See table 1 for a full account of the different OH groups densities and unit cell vectors in the primitive cell). As for the amount of interfacial water (180 molecules in each cell) it was chosen and kept constant to achieve a density of 1g/L. These water molecules were spread in a volume of 20 Å thickness times the surface area (similar for both cells). It is worth pointing out that our approach gave results in good agreement with ref. [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF] in which a similar analysis was performed about the structure of water at the alumina/water interface.

There is a large variety of hydroxyl groups present in these surfaces; each one is characterized by the coordination of Al atoms (they can be either four coordinated or six coordinated in tetrahedral or octahedral environments, respectively) and the number of Al-O bonds, giving rise to µ 1 , µ 2 and µ 3 groups. The full list is presented in Table 1 and Fig. 2. The charge z of each Al n OH d group (short name for a real structure shown in Fig. 1) is calculated from the number and the type of Al atom around the oxygen (+0.5 for Al VI or +0.75 for Al IV following the bond valence approach [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF][START_REF] Digne | Use of DFT to Achieve a Rational Understanding of Acid-Basic Properties of Gamma-Alumina Surfaces[END_REF] ), and the number of protons (+1 for each proton). The computation of their pKa values using the MUSIC model is based on the empirical formula given in [START_REF] Hiemstra | Intrinsic Proton Affinity of Reactive Surface Groups of Metal (Hydr)oxides: The Bond Valence Principle[END_REF] and retranscripted below in equation 1. Thus, it was found that the protonation constant of a surface site is proportional to the undersaturation of the oxygen valence (the oxygen charge, -2, is not fully neutralized by the electron transfer to the neighbor atoms) and is defined as:

Log K = -19.8 ( s Al + d.s H + a.(1-s H ) -2) (1)
Where s Al is +0.5 or +0.75 for VI and IV coordinated Al atoms, respectively, d is the number of protons (0, 1 or 2), s H is +0.8, a is the number of accepting hydrogen bonds with water molecule (see Fig. 1); a and d depends on the geometric constraints around the oxygen atom. In the original work, it is assumed that a+d = 2 for simply coordinated groups (µ 1 -OH), 1 for triply coordinated groups (µ 3 -OH), and 1 or 2 for doubly coordinated groups (µ 2 -OH), and that each proton forms one (value of p in the Fig. 1) hydrogen bond with a water molecule.

The interfacial structure with the evolution of charges and potential near the surface should be taken into account to properly reflect surface reactions. Several models have been proposed such as, for example, Triple Layer Model, Triple Plane Model or Basic Stern Model [START_REF] Lützenkirchen | Surface Complexation Models of Adsorption[END_REF] . We chose to use the latter, as proposed by Boily [START_REF] Boily | Modeling Proton Binding at the Goethite (α-FeOOH)-Water Interface[END_REF] , due to its having fewer adjustable parameters. The value of the capacitance has been chosen at 0.92 F/m 2 , in the range of the typical values found for modelling metal (hydr)oxides.

RESULTS AND DISCUSSION

3.1 Description of surface from DFT calculations

Calculation of protonation constants

Protonation constants were calculated for each site using the MUSIC model from Hiemstra et al. [START_REF] Hiemstra | Multisite Proton Adsorption Modeling at the Solid/Solution Interface of (Hydr)Oxides: A New Approach: II. Application to Various Important (Hydr)Oxides[END_REF][START_REF] Hiemstra | Multisite Proton Adsorption Modeling at the Solid/Solution Interface of (Hydr)Oxides: A New Approach: I. Model Description and Evaluation of Intrinsic Reaction Constants[END_REF] with corrections accounting for the influence of hydrogen bonds at the solid-liquid interface based on the bond-valence principle in reference [START_REF] Hiemstra | Intrinsic Proton Affinity of Reactive Surface Groups of Metal (Hydr)oxides: The Bond Valence Principle[END_REF] . For the latter, we considered two approaches: firstly, by choosing to set the number of hydrogen bonds around the oxygen atom of the hydroxyl group arbitrarily to the one in ice (as in reference [START_REF] Hiemstra | Intrinsic Proton Affinity of Reactive Surface Groups of Metal (Hydr)oxides: The Bond Valence Principle[END_REF] ); secondly, by using for Al n OH surface species at the Point of Zero Charge the average number of hydrogen bonds obtained from AIMD simulations in reference [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] . The results of both approaches are presented in tables 2 and 3, respectively. Table 2. Calculation of protonation constants using the number of H bonds around oxygen atoms in hydroxyl groups (a+d) as defined in [START_REF] Hiemstra | Multisite Proton Adsorption Modeling at the Solid/Solution Interface of (Hydr)Oxides: A New Approach: I. Model Description and Evaluation of Intrinsic Reaction Constants[END_REF] . The deprotonated surface group with a (negative) charge z, Al n O z , can react twice with one proton, allowing the calculation of two successive protonation constants. The difference between some of these constants are more than 10 log units, and only those in the pH range 0-14 are listed in the Table 2.

OH

The protonation constants of surface groups involving Al atoms of octahedral geometry are similar to those calculated on gibbsite surfaces in [START_REF] Yong | Interfacial Charging Phenomena of Aluminum (Hydr)oxides[END_REF] ; surface sites involving Al atoms in tetrahedral positions are more acidic, since there is a stronger Al to O charge transfer in O bonded to Al IV sites (charge 3/4) than to Al VI sites (charge 3/6).

Using the number of water molecules obtained by AIMD simulation in ref. [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] , the undersaturation of the oxygen atoms can be refined. For instance, sites with the same local structure µ n become different for each face, since DFT simulations have shown that the structure of the water molecule in interaction with the surface depends on the orientation of the face (see [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] ).

The values of protonation constants illustrate the strong effect of the number of hydrogen bonds.

Thus, for µ 1 (Al IV ) on the (110) surface, compared to the ice structure, one more hydrogen bond has been found. This leads to a decrease of the second protonation constant of about 4 log unit, and allows to the first protonation to occur below pH 14. An opposite effect is seen for µ 2 (Al VI -Al VI ), with a lower mean number of hydrogen bonds (between 1 and 1.5 instead of 2), and an increase of log K for the two protonations.

In the calculation of pKas, another assumption is made concerning the charge of the proton of the hydroxyl group, in interaction with surrounding water molecules by hydrogen bonds. In the initial model, one Al n OH…OH 2 bond is formed for each proton of surface OH groups. The charge +1 of the proton decreased to +0.8 due to the formation of one H donating bond. However, this choice for the number of hydrogen bonds comes from an analogy with ideal bulk water, and the real number of H donating bonds could be different. AIMD calculations allow to determine the mean number of water molecules involved in H donating bonds and can be used instead of the arbitrarily set value of 1 (p in Figure 1). This value is important because it has an impact on the second protonation constant since the calculation of this value is based on the geometry of the surface group after the first protonation step.

In Table 4 it can be seen that the arbitrarily and generally used set value of 1 hydrogen bond is slightly overestimated. This can be explained by steric effects at the surface. Thus, instead of +0.8 units, the charge of H will increase following the relation 1-0.2p. However, these results show that a 25% decrease in the H charge leads to a small decrease of about 1 log unit for K. Al VI µ 3 O -0.5 + H +  µ 3 OH 0.5 4.75 n.a. * for H bonds around oxygen, see Table 3. n.a.: non applicable since the site is deprotonated.

One of the surface sites obtained by DFT simulations (Table 1, first line), consists on a surface aluminum atom with a chemisorbed water molecule. We consider µ 1 -H 2 O groups as interfacial water molecules and thus neglect them in the chemical equations considered in the interfacial model. Indeed, if we assume that the interaction of H 2 O molecules with the support is weak and therefore the pKa of these molecules is the same as those in bulk water. The model described above (Table 4) has been implemented in PHREEQC and the curves of surface charge have been simulated for three different ionic strengths (Fig. 3). The result follows the general trends found in aluminum (hydr)oxides titration curves with a PZC around 8-9 [START_REF] Kosmulski | Chemical Properties of Material Surfaces[END_REF] . The surface speciation is then detailed for a given ionic strength (1 mM). The protonation fraction of each site for both faces as a function of pH is shown in Fig. 4. Several groups neither release nor accept H from interfacial water. Therefore, for both faces, fluctuations in surface charge can be accounted for by the interaction of µ 1,VI -OH with interfacial water. Besides, and in accordance with the analysis in reference [START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] , the interaction of these sites with the solvent is stronger in the (110) orientation. Al VI-Al IV (mono-/de-), (4) µ 2 -Al VI-Al VI (mono-/de-) and ( 5) µ 3 -Al VI (mono-/de-). The indicated acid/base couple (mono-or di-/ de-or mono-) correspond to the equilibria detailed in Table 3.

Moreover, at all pH values considered in this work, µ 3 -OH groups lose their proton readily. This is an interesting result since there is evidence that µ 3 -OH are present at the solid-gas interface [START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF] . Hence, µ 3 may tend to capture protons from liquid water during the drying process of alumina. This is a demonstration that the reactivity of a surface sites depends on the surrounding media and therefore, in order to rationalize adsorption processes at the solid liquid interface, materials characterizations need to be done in the presence of the encompassing liquid phase.

The simulated curves have been compared to experimental titrations carried out on -alumina (Fig. 7). The detailed outcomes for experimental surface charge are described in reference [START_REF] Machesky | Surface Protonation at the Rutile (110) Interface: Explicit Incorporation of Solvation Structure within the Refined MUSIC Model Framework[END_REF] . Keeping in mind that the simulated curves are not based on any fitting parameters, they can be considered as in good qualitative agreement with experiment at pH below 9. Moreover, effects induced by the ionic strength are well reproduced and PZC are in reasonable good agreement with experiments. Beyond the pH value of 9, the match between the experimental and theoretical curves worsens, due to the asymmetry of the experimental titration curves, whose slopes are higher above pzc, as this has been previously observed by Hiemstra et al. This trend cannot be captured by a model based only on protonation surface reactions, which give intrinsically concentration evolutions symmetrical on both sides of pKa. To solve this problem, authors have used the possibility to form surface ion pairs between charged surface hydroxyls and electrolyte ions according to the following general equations, where A -and C + are respectively inert anion (Cl -, NO 3 -) or cation (Na + , K + ), respectively.

Al

n OH d z+d + A - Al n OH d A z+d-1 K c,A with z+d > 0 (2) Al n OH d z+d + C +  Al n OH d C z+d+1 K c,C with z+d < 0 (3)
The formation of surface ion pairs decreased the surface potential, in absolute value, and so the electrostatic repulsion between surface and H + /OH -ions. It allows to increase the surface density of (de)protonated groups.

The direct measurement of the constants remains difficult, and they are generally obtained from the fitting of acid-base curves and a difference in the behavior of cations and anions can be at the origin of the asymmetry. Thus, Hiemstra et al. have suggested that Cl -are able to reach subsurface hydroxyl sites, modifying the impact of counter-ions on the surface concentration of positive or negative groups. Mayordomo et al. [START_REF] Mayordomo | Selenium(IV) Sorption Onto γ-Al2O3: A Consistent Description of the Surface Speciation by Spectroscopy and Thermodynamic Modeling[END_REF] have used a more common approach, with a higher value for K c,Na (0.8 vs. 0.67 for K c,Cl ). In a first attempt to take into account surface ion pairs in the modelling, the same simulation has been done including the equilibria ( 2) and ( 3) for all the charged surface species.

Values from Mayordomo et al. [START_REF] Mayordomo | Selenium(IV) Sorption Onto γ-Al2O3: A Consistent Description of the Surface Speciation by Spectroscopy and Thermodynamic Modeling[END_REF] This calculation illustrates the importance of surface ion pair formation in the prediction of surface charging. The current approach is rather rough, due to the lack of independent data on adsorption of counter-ions. The modelling of these reactions would be an important step in the understanding of the interface structure. According to our simulations, the pH of interfacial water seems to affect mostly the speciation of alumina µ 1,VI -OH groups. These results could be used for the interpretation of complexation data or a refinement of a heterogeneous catalyst preparation by an impregnation technique. Moreover, future and more complex AIMD simulations for the description of the chemical weathering phenomena of this material (as in [START_REF] Réocreux | Biomass Derivatives in Heterogeneous Catalysis: Adsorption, Reactivity and Support from First Principles[END_REF] ) could benefit from this result allowing researchers to focus on the most reactive sites of this complex system.

(

  2x3) and (2x2) 4-slab thick supercells containing 1164 (96 Al 2 O 3 units, 48 adsorbed water molecules and 180 water molecules in the liquid phase) and 1004 (64 Al 2 O 3 units, 48 adsorbed water molecules and 180 water molecules) atoms for the (100) and (

Fig. 3 .

 3 Fig.3. Simulated surface charge curves at different ionic strengths from the DFT-MUSIC model

Fig. 4 .

 4 Fig. 4. Protonated fraction of each sites as a function of pH from calculation by the DFT-MUSIC model

Fig. 7 .

 7 Fig. 7. Surface charge curves at different ionic strengths (1) 1 mM, (2) 10 mM and (3) 100 mM from

  have been used for adsorption of cations or anions, along with the electrostatic model they used (TLM). The result is shown in Fig. 8. As expected, two effects can be seen: a global increase of the surface charge (in absolute value), stronger in the pH range where surface is negative.

Fig. 8 .

 8 Fig.8. Simulated surface charge curves at an ionic strength of 100 mM using DFT-MUSIC model(1)

  presented in this work show to what extent atomic scale structural dataobtained from Density Functional Theory calculations-can be used with the MUSIC model to simulate the speciation in solution of surface sites on the poorly crystalline material, γ-Al 2 O 3 . The results presented are in qualitative agreement with experimental data, despite the absence of any fitting parameter. Hence, although more elaborate theoretical calculations need to be performed in order to get a better match between simulation and experiment, the atomistic model for γ-Al 2 O 3surfaces from Digne and co-workers[START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF] can be used with the methodology presented herein to obtain the speciation of this material in suspension at the initial stages of the impregnation or chemical weathering of γ-Al 2 O 3 and seems an appropriate starting point to improve MUSIC based results with more demanding theoretical calculations in the future as the ones done in[START_REF] Sulpizi | Acidity Constants from Vertical Energy Gaps: Density Functional Theory Based Molecular Dynamics Implementation[END_REF] at the DFT level.

Table 1

 1 and figure2present, for the main orientations of γ-Al 2 O 3 , the different sort of OH groups present in the surface models of Digne and co-workers at the highest water coverage considered in reference[START_REF] Digne | Hydroxyl Groups on Gamma-Alumina Surfaces: A DFT Study[END_REF] (0.17 and 0.18 OH/nm 2 for the (100) and (110) orientations, respectively) and used in the Ab Initio Molecular Dynamics study of Ngouana-Wakou et al.Fig. 2: From left to right, (1x1) representation of the (100) and (110) γ-Al 2 O 3 surface orientations considered in this work (only the uppermost atoms are shown). Aluminum atoms are shown in green, hydrogen atoms in white and oxygen in surface and OH groups are shown in blue and red, respectively. Al atoms are in VI coordination unless otherwise stated.

[START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF]

. Chemisorbed water molecules are referred to as H 2 O sites; OH groups are referred to as µ n sites, where n stands for the number of aluminum atoms bonded to the central oxygen atom. Transition aluminas contain two different types of aluminum atoms: six-coordinated and four-coordinated atoms in octahedral (Al VI ) and tetrahedral (Al IV ) environments, respectively. The nature of the Al involved in the different OH groups are listed in column 2 of Table

1

and illustrated in Fig.

2

.

Table 1 .

 1 Characteristics of γ-Al 2 O 3 surface OH groups from references[START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF][START_REF] Mayordomo | Selenium(IV) Sorption Onto γ-Al2O3: A Consistent Description of the Surface Speciation by Spectroscopy and Thermodynamic Modeling[END_REF] at water coverages of 0.17

	and 0.18 OH/nm 2 for the (100) and (110) orientations, respectively (see figure 2)
	OH site	Al involved	Amount in simulation box *	Surface Density (nm -2 )
			(100)	(110)	(100)	(110)
	H 2 O	Al VI	2	3	4.31	4.48
	µ 1	Al VI	2	1	4.31	1.49
	µ 1	Al IV	0	1	0	1.49
	µ 2	Al VI and Al IV	1	2	2.16	2.98
	µ 2	Only Al VI	0	1	0	1.49
	µ 3	Only Al VI	1	1	2.16	1.49
	* unit cell vectors (u,v) for the (100) surface: 8.36 and 5.55 Å; for the (110) surface: 8.03 and
	8.36 Å.

Table 3 .

 3 Calculation of protonation constants using a number of H bonds around oxygen atoms in hydroxyl groups (a+d) obtained from AIMD simulations in[START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] .

	Al involved	Protonation reaction	Log K Number of H bonds

Table 4 .

 4 Calculation of protonation constants using AIMD simulations to determine the number of H bonds around oxygen atoms in hydroxyl groups (a+d), and the number of H bonds around their respective protons (p).

	OH site Al involved	Protonation reaction	Log K	Number of H
					bonds with
					proton*
		Surface (100)
	µ 1	Al VI	µ 1 OH -0.5 + H +  µ 1 OH 2 0.5	8.16	0.86
	µ 2	Al VI-Al IV	µ 2 O -0.75 + H +  µ 2 OH 0.25	10.32	n.a.
	µ 3	Al VI	µ 3 O -0.5 + H +  µ 3 OH 0.5	5.94	n.a.
		Surface (110)
	µ 1	Al IV	µ 1 OH -0.25 + H +  µ 1 OH 2 0.75	0.12	0.75
	µ 1	Al IV	µ 1 O -1.25 + H +  µ 1 OH -0.25	12.97	n.a.
	µ 1	Al VI	µ 1 OH -0.5 + H +  µ 1 OH 2 0.5	9.16	0.95
	µ 2	Al VI-Al IV	µ 2 O -0.75 + H +  µ 2 OH 0.25	10.45	n.a.
	µ 2	Al VI -Al VI	µ 2 OH 0.25 + H +  µ 2 OH 2 1.25	0.29	0.57
	µ 2	Al VI -Al VI	µ 2 O -0.75 + H +  µ 2 OH 0.25	13.86	n.a.
	µ 3

ACKNOWLEDGMENTS

We would like to acknowledge Dr. P. Raybaud (IFP Energies nouvelles) is also greatly acknowledged for fruitful discussions and a review of an initial version of this manuscript. This work was supported by the French National Research Agency within the framework of the ANR-14-CE08-0019 SLIMCAT project. Calculations were performed using HPC resources from GENCI-CINES (Grant A0020807386).


Matteo Conforti

Fabio Baronio

Costantino De Angelis

Guido Sanna

Davide Pierleoni

Paolo Bassi

Spectral shaping of femtosecond pulses in aperiodic quasi-phase-matched gratings

Keywords: Quasi phase matching, Second harmonic generation, Pulse shaping, Genetic algorithms PACS: 42.65.Ky, 42.65.Re, 42.65 02.60.Pn

Introduction

Important practical applications of high-energy femtosecond laser sources are emerging in a variety of fields, including material processing for optical [1] or mechanical [2] applications, biotechnology and medicine [3], environmental monitoring [4], and ultra-fast chemistry [5]. The development of compact and reliable sources with the desired temporal and spectral characteristics is then of current high technological interest.

In this framework, self-induced effects occurring during propagation of light in nonlinear media can be exploited to modify light propeties in the spatial, temporal and spectral domains [START_REF] Kivshar | Optical solitons[END_REF]. Specifically, nonlinear spectral manipulation could be used to design devices enabling to modulate the amplitude, phase profile, and carrier frequency of femtosecond pulses over a broad range. These devices, potentially much simpler than optical parametric amplifiers [START_REF] Cerullo | [END_REF], could enhance the performances of sources such as fibre lasers [8], which are very compact and efficient but do not offer frequency tunability. In this context a key role is played by the choice of the nonlinear medium where self-induced effects take place: periodically poled optical crystals offer unique possibilities since they provide high effective nonlinearity and the flexibility of engineering the nonlinear poled structure with a photolithographic mask. Transverse quasi-phase-matched gratings have been made both for shaping of second harmonic beams [9] and to extend the spectral coverage of optical parametric oscillators [10]. Pulse compression, gain conversion [11,12] and spectral shift [13] enhancement in frequency-conversion schemes in longitudinally chirped quasi-phase-matched gratings have been demonstrated. More recently the issue of arbitrary second-harmonic pulse shaping in aperiodically poled crystals was solved by means of both optimal control and simulated annealing techniques [14][15][16].

In this Paper, we analyze the use of aperiodically poled nonlinear crystals to tune the carrier frequency and spectral profile of fundamental frequency (FF) broadband near-infrared femtosecond pulses, by means of cascading second harmonic (SH) interactions. To this goal, we use a genetic algorithm for properly designing piecewise uniform periodically poled ferroelectric crystals capable of moulding femtosecond infrared pulsed radiation. After briefly recalling the basic equations describing the phenomenon we are investigating and the features of genetic algorithms, we will report and comment results concerning some possible applications of the proposed structures.

Aperiodic quasi-phase-matched gratings for spectral shaping

We consider the process of parametric quadratic wave mixing in a non uniform quasi-phase-matching (QPM) structure with variable QPM period. Averaging over the QPM periods, at the lowest order the coupled nonlinear partial differential equations governing the propagation of the complex electric field amplitudes E m (z, t) (m = 1 fundamental field FF, m = 2 second harmonic field SH) of two plane waves of central frequencies ω i (ω 2 = 2ω 1 ) and wave numbers k m travelling along the z axis read:

i ∂E 1 ∂z = β 1 2 ∂ 2 E 1 ∂t 2 -χ 1 exp[-i∆k(z)z]E * 1 E 2 , i ∂E 2 ∂z = iδ ∂E 2 ∂t + β 2 2 ∂ 2 E 2 ∂t 2 -χ 2 exp[i∆k(z)z]E 2 1 , (1)
where ∆k(z) = 2k 1 -k 2 + 2π/Λ(z) is the effective wave-number mismatch, Λ(z) is the variable poling period, δ is the difference between the inverse group velocity at FF and SH, β 1,2 is the group velocity dispersion at FF and SH, and χ m = 2ω 1 χ (2) /πcn m , being χ (2) the second-order susceptibility. In the following, since the refractive indices at FF and SH are similar (n 1 ≈ n 2 ), we assume

χ 1 = χ 2 = χ.
The problem of the spectral control of the FF field at output by means of mismatch engineering can be stated as follows: having fixed the FF wave at the input of the nonlinear crystal, determine the piecewise constant effective mismatch ∆k(z) that minimizes a cost (or fitness) function.

For example we may want to minimize the distance between the FF spectrum at output [ Ê1 (z = L, ω)] and the target [ Ê1,target (ω)], where Ê1 (ω) stands for the Fourier transform of E 1 (t). As a measure of the distance between the target and the FF wave at output, we can use the following fitness function:

J = +∞ -∞ Ê1 (z = L, ω) -Ê1,target (ω) 2 dω. ( 2
)
Another case of interest is the maximization of the shift of the center of mass of the FF spectrum toward red or blue. In this case the fitness function reads

J = ± +∞ -∞ ω | Ê1 (z = L, ω)| 2 dω +∞ -∞ | Ê1 (z = L, ω)| 2 dω , ( 3
)
where the + (-) sign imposes redshift (blueshift).

Once the mismatch function ∆k(z) is obtained, it is straightforward to determine the nonuniform aperiodic QPM structure that implements the desired ∆k(z), by using the formula Λ

(z) = 2π/[∆k(z) -2k 1 + k 2 ].
The goal here is thus to find the optimum ∆k(z) function. Searching among a collection of candidate solutions for a desired one is a very common problem and many different algorithms can be used for this kind of optimization (as the Dijkstra or the Floyd algorithms). In this framework genetic algorithms are a class of optimization procedures that apply principles of natural evolution to optimize a given objective function [START_REF] Davis | Genetic Algorithms and Simulated Annealing[END_REF][START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF], and are able to solve a variety of problems that are not well suited for standard optimization algorithms, including problems in which the objective function is discontinuous, nondifferentiable, stochastic, or highly nonlinear. Genetic algorithms were recently successfully used even for solving very complex optical design problems [START_REF] Goh | [END_REF].

The genetic procedure that we used to find the optimal ∆k(z) can be summarized in the following iterative algorithm:

(1) Given a crystal of length L, choose a suitable number of uniform domains characterized by ∆k n , with n = 1, . . . , N . The variables to be optimized (∆k n ) are stored in a vector which is usually called individual. Each individual represents an aperiodic structure. (2) Randomly generate an initial population of M individuals. In order to have a sufficiently broad population it is advisable to take M ≥ N . (3) Compute the fitness function J(m) for each individual m in the current population, and select parents based on their fitness. This steps requires to solve M times Eqs. (1) (we employed a standard finite-difference fourth order Runge-Kutta method). (4) Generate the new population from the parents. Children are produced either by making random changes to a single parent (mutation) or by combining the vector entries of a pair of parents (crossover ). ( 5) If a generated individual fits the target, or more generally it meets a predefined convergence criterion, the iterative procedure stops, otherwise it goes back to point 3.

Results and discussion

In order to show the capabilities of our method to design devices able to perform spectral shaping of femtosecond infrared pulses, we consider Periodically Poled Stoichiometric Lithium Tantalate (PPSLT) crystals. For FF input pulses with center wavelength of 1500 nm, the parameters to be used in Eqs.

(1) are:

χ 0 = 1.87 • 10 -5 V -1 , δ = 4 • 10 -10 s/m, β 1 = 0.94 • 10 -25 s 2 /m and β 2 = 3.7 • 10 -25 s 2 /m.
As a first example we apply our design method to blue shifting spectral shaping. We consider a L = 2.5 cm long crystal, a gaussian transform limited FF input with peak intensity I = 20 GW/cm 2 and spectral full width at half maximum intensity (FWHMI) of 55nm (corresponding to a temporal FWHMI of T F F = 60 f s) centered at λ = 1500 nm. As a target we choose a blue shifted FF gaussian wave with the same spectral bandwidth centered at λ = 1460 nm (blue shifting) and a FF maximum depletion of 50%. We divide the crystal in N = 10 uniform domains, which turned out to be a sufficiently large number to obtain reasonable agreement with the target. Using a higher number of domains does not considerably improve the performances of the method. In all the examples we obtained convergence to the optimal solution after around 500 iterations (or generations) of the genetic algorithm. These runs requires a CPU time of a few hours (1-4) on a standard personal computer.

The results of the genetic algorithm based design are shown in Fig. 1. In Fig. 1(a) we show the piecewise constant effective mismatch ∆k(z) of the PPLST aperiodic grating poling period. Figure 1(b) shows the FF spectral evolution in the (λ, z) plane: it is clear that the frequency shifting mainly takes place in the first half of the crystal, whereas in the last part the spectral shape is modified to fit the gaussian target. In the second example, using the same input signal, we apply our design procedure to red shifting spectral shaping. As a target we thus choose a FF gaussian wave of spectral full width at half maximum intensity (FWHMI) of 55nm centered at λ = 1540 nm and a FF maximum depletion of 50%.

The results are shown in Fig. 2. In Fig. 2(a) we show the piecewise constant effective mismatch ∆k(z) of the PPLST grating. Figure 2(b) shows the FF spectral evolution in the (λ, z) plane. Once again we can note that the spectral shift takes place mostly in the first half of the device, while the last part of the crystal engineers the spectral shape of the output. Figure 2(c) reports the spectral intensity of input FF field (thin blue line) and output FF field (thick black line) compared to the FF target field (red circles): the agreement is very good. Note that in this example the effective cascading nonlinearity is of self-focusing type [21], owing to the positive values of ∆k(z). This regime usually prevents to obtain clean spectral and temporal profiles, because of the combined effect of positive dispersion and effective self-focusing cascading nonlinearity, which eventually leads to pulse breakup. However the engineering of the grating prevents these undesirable phenomena and gives substantially good results. Finally Fig. 2(d) shows the FF (blue line) and output (black line) temporal intensity pulse profiles. The temporal profile is smooth and does not exhibit breakup phenomena. As demonstrated by the above examples, the spectral shift caused by a phasemismatched second harmonic generation process under conditions of strong group-velocity mismatch, can be efficiently controlled by acting on mismatch engineering. However spectral shifting can also be affected by other param-eters, mainly by the initial FF intensity and chirp [21]. Our genetic device design is well suited also to take into account these parameters. In the following we will present an example where the FF input chirp is among the accessible parameters to optimize in order to obtain the maximum spectral shift of the FF pulse. Thus the genetic variables to be optimized will include the mismatch ∆k(z) and the chirp C of the input pulse, i.e. each individual in the genetic algorithm now represents an aperiodic structure and a particular chirp of the input pulse.

Let us now consider a L = 2.5 cm long crystal and a FF wave with a temporal gaussian profile with T F F = 60 f s, input peak intensity I = 20 GW/cm 2 , centered at λ = 1500 nm. As a target, we choose to maximize the spectral FF blue-shift by optimizing the mismatch grating and the initial FF pulse chirp using the fitness function (3). A similar result was obtained in Ref. [13], where a linear variation of ∆k(z) was assumed. Figure 3(a) shows that the linear dependence is indeed the best choice and also identifies the best values of the initial and final mismatch. The optimized initial chirp is C = -β 1 z c with z c = -0.013 m which means that the waist plane of the input FF gaussian pulse is roughly located in the middle of the crystal. Figure 3(b) shows the FF spectral evolution in the (λ, z) plane. Figure 3(c) reports the spectral intensity of chirped input FF field (thin blue line) and normalized output FF field (thick black line). Even if the target pulse shape is not imposed by the fitness function, the genetic algorithm gives a gaussian-shaped profile. The maximum attainable shift is 65 nm, with a strong pump depletion (about 88%). This strong pump depletion occurs since the values of the shift are related to energy exchange between FF and SH: as a rule of thumb the more energy is converted from FF to SH, the higher is the frequency shift. Finally Fig. 3(d) shows the FF chirped input (solid blue line) and normalized FF output temporal pulse profile (thick black line). As clearly shown in the picture, by properly acting on the chirp of the input field it is possible to simultaneously control the output fundamental field spectral shape and to compensate dispersive effects to get an output temporal FWHMI smaller than the input one.

Conclusions

In conclusion a genetic technique has been used to mould the spectral properties of femtosecond infrared optical pulses, propagating in longitudinally aperiodic poled crystals, by mismatched second order parametric processes.

In our work we account for the engineering of the non uniform QPM structure and the optimization of the input pulse characteristics (i.e., the initial chirp).

Our results suggest a valuable means for shaping infrared femtosecond pulses in a way that is simpler and cheaper than the state of art. In fact the proposed devices are potentially much simpler than optical parametric amplifiers and could enhance the performances of sources such as fibre lasers, which are very compact and efficient but emit at a single fixed frequency. Moreover the method described here is quite general and can be exploited in other areas of optics whenever the spectral control of the optical pulse is required.

Figure 1 (

 1 c) shows the spectral intensity of input FF field (blue line) and output FF field (thick black line) compared to the FF target field (red circles): the agreement is very good. Finally Fig.1(d)shows the FF input (blue line) and output (black line) temporal pulse profile. The temporal profile is very clean, almost symmetric and delayed of about 100f s. This delay is caused partially by the spectral blueshift, which moves the FF pulse in a spectral region where the group velocity is decreased owing to dispersion, an partially induced by an effective steepening due to the dragging by the slower SH pulse[20].

Fig. 1 .

 1 Fig. 1. (Color online) a) Effective mismatch profile ∆k(z). b) FF spectral evolution in the (λ, z) plane. c) Spectral intensity profile of FF input field (thin blue line) and FF output field (thick black line) compared to the FF target (red circles). d) Input (thin blue line) and output (thick black line) FF temporal intensity profiles.

Fig. 2 .

 2 Fig. 2. (Color online) a) Effective mismatch profile ∆k(z). b) FF spectral evolution in the (λ, z) plane. c) Spectral intensity profile of FF input field (blue line) and FF output field (black line) compared to the FF target (red circles). d) Input (blue line) and output (black line) FF temporal intensity profiles.

Fig. 3 .

 3 Fig. 3. (Color online) a) Effective mismatch profile ∆k(z). b) FF spectral evolution in the (λ, z) plane. c) Spectral intensity profile of FF chirped input field (thin blue line) and normalized FF output field (thick black line). d) Input chirped (solid blue line) and normalized output (thick black line) FF temporal pulse profile.The results of the iterative algorithm are shown in Fig.3. In Fig.3(a) we show the piecewise constant effective mismatch ∆k(z) of the PPLST grating. The values of ∆k(z) evolve in a linear fashion from -20000m -1 to 20000 m -1 . A similar result was obtained in Ref.[13], where a linear variation of ∆k(z) was assumed. Figure3(a) shows that the linear dependence is indeed the best

Acknowledments

Part of this work has been developed in the framework of the project PRIN 2005 "Spectral and temporal control of femtosecond pulses by second order non linear processes".


est destinée au dépôt

Evidence of a link between resting energy expenditure and bone remodelling, glucose homeostasis and adipokine variations in adolescent girls with anorexia nervosa

Laurent Maïmoun, S. Guillaume, P. Lefebvre, P.

  Philibert, H. Bertet, M.C. Picot, L. Gaspari, F. Paris, M. Seneque, A.-M. Dupuys, et al. To cite this version: Laurent Maïmoun, S. Guillaume, P. Lefebvre, P. Philibert, H. Bertet, et al.. Evidence of a link between resting energy expenditure and bone remodelling, glucose homeostasis and adipokine variations in adolescent girls with anorexia nervosa. Osteoporosis International, 2016, 27 (1), pp.135 -146. 10.1007/s00198-015-3223-x. hal-01800536


J Castillo-Navarro
email: javiera.castillonavarro@onera.fr

N Audebert

A Boulch

B Le Saux

S Lefèvre

What Data are needed for Semantic Segmentation in Earth Observation?

Keywords: Deep Learning, Supervised Learning, Semantic Segmentation, Land Use/Land Cover Mapping

This paper explores different aspects of semantic segmentation of remote sensing data using deep neural networks. Learning with deep neural networks was revolutionized by the creation of ImageNet. Remote sensing benefited of these new techniques, however Earth Observation (EO) datasets remain small in comparison. In this work, we investigate how we can progress towards the ImageNet of remote sensing. In particular, two questions are addressed in this paper. First, how robust are existing supervised learning strategies with respect to data volume? Second, which properties are expected from a large-scale EO dataset? The main contributions of this work are: (i) a strong robustness analysis of existing supervised learning strategies with respect to remote sensing data, (ii) the introduction of a new, large-scale dataset named MiniFrance.

I. INTRODUCTION

The ever-growing amount of remote sensing imagery data in the last two decades has allowed new developments in the fields of ecology, urban planning or natural disaster response. Those data are also more easily available, even openly as in the Copernicus or Landsat programs. However, data exploitation still remains a bottleneck. It requires human interprets, for example to identify tree species and study deforestation in a local ecosystem, or to find new buildings and measure growth of urban areas.

Thanks to the new deep learning methods developed for processing multimedia images in recent years, it is now possible to automate most of these tasks for Earth Observation (EO) data. Indeed, many state-of-the-art algorithms for object detection and image segmentation or classification [START_REF] Audebert | Beyond RGB: Very High Resolution Urban Remote Sensing With Multimodal Deep Networks[END_REF], [START_REF] Rey | Detecting animals in African Savanna with UAVs and the crowds[END_REF] have been successfully transferred to aerial and satellite images. It allows to produce quickly and without human intervention precise semantic maps, in both urban and rural contexts.

However, these learning algorithms rely heavily on the availability of annotated image databases. Even if collaborative cartographic resources such as OpenStreetMap can be used as annotations [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the INRIA aerial image labeling benchmark[END_REF], these are restricted in terms of semantics (only roads, buildings, etc.). Thus, the question of quantifying the influence of existing datasets on the models we learn arises. Moreover, we need to define what it requires to make a good dataset for training EO data classification algorithms, and thus J. Castillo-Navarro's work is partially funded by a CNES grant. N. Audebert's work is partially funded by ONERA-Total Naomi project. make a first step towards the ImageNet (a massive, annotated image dataset) of remote sensing.

Thus, the contributions of this paper are twofold: (i) An experimental analysis of the amount of data necessary to successfully achieve supervised learning, and moreover, of the required data variability; (ii) The constitution of a large-scale dataset of various remote sensing images and annotations from different sources, which overcomes limitations of standard reference datasets.

The rest of the paper is organized as follows. First we explain the problem statement and explore some related work in Sec. II. We then discuss the robustness of supervised learning in two cases: on current, small-scale datasets in Sec. III and at large-scale in Sec. IV. Finally, we conclude and draw perspectives in Sec. V.

II. PROBLEM STATEMENT AND RELATED WORK

Convolutional Neural Networks (CNNs) have greatly benefited to computer vision and remote sensing for different tasks including image segmentation such as land use and land cover mapping. Furthermore, the pixel-wise predictions of Fully Convolutional Networks (FCN) and derivatives offer an appealing approach for Earth Observation data analysis. These models have reached the state-of-the-art for remote sensing data analysis including: land cover mapping in urban areas [START_REF] Marmanis | Semantic segmentation of aerial images with an ensemble of CNNs[END_REF], [START_REF] Paisitkriangkrai | Effective semantic pixel labelling with convolutional networks and conditional random fields[END_REF], object segmentation [START_REF] Audebert | Segment-before-Detect: Vehicle Detection and Classification through Semantic Segmentation of Aerial Images[END_REF] or semantic segmentation for multi-modal and multi-scale remote sensing data [START_REF] Audebert | Beyond RGB: Very High Resolution Urban Remote Sensing With Multimodal Deep Networks[END_REF]. Deep networks achieve top results on most of the public benchmark EO datasets [START_REF] Campos-Taberner | Processing of extremely high-resolution LiDAR and RGB data: Outcome of the 2015 IEEE GRSS Data Fusion Contest-Part A: 2-D Contest[END_REF], [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the INRIA aerial image labeling benchmark[END_REF], [START_REF] Mou | Multitemporal Very High Resolution from Space: Outcome of the 2016 IEEE GRSS Data Fusion Contest[END_REF], [START_REF] Rottensteiner | The ISPRS benchmark on urban object classification and 3D building reconstruction[END_REF], with excellent performances over 80 to 90% accuracy. However, these methods are supervised, and still require a large amount of dense pixel-wise annotated data during the training phase. Most of the time, there exists only a small number of annotations on which to train models. In this context, we aim to answer the following question: how much data are needed to train a supervised neural network for semantic segmentation in Earth Observation?

This problem can be tackled in three ways: a first category tries to reduce the number of parameters of a model, while keeping data unchanged [START_REF] Marcos | Land cover mapping at very high resolution with rotation equivariant cnns: Towards small yet accurate models[END_REF], a second one trains models with large, incomplete ground truth [START_REF] Maggiolo | Improving maps from CNNs trained with sparse, scribbled ground truths using fully connected CRFs[END_REF] and a third one is to focus on a small amount of well annotated data. In this paper, we consider the third option. In addition, applying a model trained on one dataset to another one yields in a serious drop of accuracy. This is known as the domain adaptation problem [START_REF] Yokoya | Open data for global multimodal land use classification: Outcome of the 2017 IEEE GRSS Data Fusion Contest[END_REF], and leads us to wonder what are the desired properties of a good training dataset.

To this aim, we set up several experiments to investigate the behaviour of network training with respect to the data. First, we aim at quantifying the amount of data required to obtain good performances in a standard set-up. Second, we renew the experiment with a more large-scale dataset to understand the relationship between training and genericity of the resulting models. In all experiments, we use as a baseline an efficient and versatile FCN: SegNet [START_REF] Audebert | Beyond RGB: Very High Resolution Urban Remote Sensing With Multimodal Deep Networks[END_REF], [START_REF] Badrinarayanan | SegNet: A Deep Convolutional Encoder-Decoder Architecture for Image Segmentation[END_REF]. It presents an encoderdecoder architecture, its encoder part is similar to VGG16 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] and thus can be initialized with pre-trained weights. Training is achieved with Stochastic Gradient Descent and the standard cross-entropy loss function weighted over pixels.

III. ANALYSIS OF SUPERVISED LEARNING ON SMALL-SCALE DATASETS

The ISPRS 2D Semantic Labeling [START_REF] Rottensteiner | The ISPRS benchmark on urban object classification and 3D building reconstruction[END_REF] Vaihingen dataset consists of 33 infrared-red-green tiles with a spatial resolution of 9cm/px and an average size of 2000 × 1500px. Dense annotations are available on 16 tiles for 6 classes of interest: impervious surfaces, buildings, low vegetation, trees, cars and clutter. The associated benchmark being now closed, we perform experiments using 12 annotated tiles for training (train) and 4 tiles for validation (denoted by val).

In our first experiment, we aim to test the sensitivity of supervised learning to training. So we reduce the amount of annotated images used for training, from 12 tiles to only one, while val remains unchanged. Additionally, two different training settings are compared on the course of this experiment. First, we train our network by initializing the encoder with the pre-trained weights from VGG-16 on ImageNet. Secondly, we repeat the experiment, but weights are randomly initialized using the policy from [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on ImageNet classification[END_REF]. This allows us to study if we can benefit from transfer learning. Results are presented in Fig. 1. The outcomes of this experiment are somehow surprising. However, going from 90% to 80% of overall accuracy and 65% to 50% of mIoU when reducing the number of training tiles from 12 to 1 is much better than one would expect. Indeed, we supposed that reducing the number of training tiles would seriously impact the performance of the network. One possible reason is that all the images in the Vaihingen dataset are alike, thus, to generalize on them is a relatively easy task. On the other hand, the mean IoU decreases faster than the accuracy showing that one must be careful when interpreting these scores.

To better understand the quantitative scores from Fig. 1 in terms of segmentation quality, Fig. 2 shows the different predictions obtained for tile 30. We can observe that the quality of the segmentation map decreases notably when less annotated tiles are used during the training phase. It is interesting to note that there is not a considerable difference between training with 10 tiles and with 6 tiles, however there is a greater difference when training with 1 tile: borders are less smooth and little objects (such as cars) are not well learned, which explains why the mean IoU decreases faster on Fig. 1.

To assess the idea that the Vaihingen dataset has much redundancy, we observed its statistical distribution. In Fig. 4(a) and (b), we compare the color histograms over the 3 channels for train and val. Indeed, they are almost identical, which indicates that learning on a single location might not be so challenging. Actually, this is even promising in terms of practical business applications, since mapping an area can be achieved after labeling only a few images.

IV. SUPERVISED LEARNING AT LARGE-SCALE

The previous section stressed out a limitation of standard datasets for semantic mapping. If some already take into account multiple locations, they are devoted to a single class (such as buildings [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the INRIA aerial image labeling benchmark[END_REF], [START_REF] Mou | Multitemporal Very High Resolution from Space: Outcome of the 2016 IEEE GRSS Data Fusion Contest[END_REF]) or to land cover classes [START_REF] Yokoya | Open data for global multimodal land use classification: Outcome of the 2017 IEEE GRSS Data Fusion Contest[END_REF], but do not offer generic land use classes. Consequently, we introduce the MiniFrance dataset with the purpose of adding variety to the existing datasets for semantic segmentation.

MiniFrance consists of aerial images of 16 cities or conurbations from different regions in France (see Fig. 3). It is built using the 2 following data sources:

• Orthophotos from the BD ORTHO® free access database from the IGN (French National Mapping Agency) 1 , acquired between 2012 and 2014. Each tile size is 10, 000×10, 000 pixels and the ground sampling distance of these images is 50cm/px. • Labelled ground truth obtained from the Copernicus Urban Atlas 2012 database 2 . In this paper we focus on the second level of the semantic hierarchy, and 15 land use classes are considered, including urban fabric, industrial and transport units, arable lands, pastures, forests, open spaces, waterlands, and water. We also propose a fixed partition for coherency of comparisons: 8 cities are used for training and the remaining 8 ones for test, keeping diversity in terms of architecture and urban design in both subsets. All in all, it is constituted of 2121 images, each of them of size 10, 000 × 10, 000 pixels. Therefore, MiniFrance dataset is 2719 times larger than Vaihingen.

Similarly to Sec. III, we first test the influence of the amount of training data over the classification. However, due to computational times 3 , we conduct more focused experiments. We train with the whole dataset, then only consider 10% of images on the dataset (we make sure to pick 10% of images from each conurbation to conserve the diversity of the dataset), and finally use only one city for training (the seaside town of Caen, which represents a similar amount of data: 12.5%). Test set remains the same. Following results of previous experiments (see Fig. 1), weights are initialized as the pre-trained weights for the VGG-16 network. Hence, results are shown in Table I 4 . Performances are not reaching the same level than on Vaihingen, which could be expected since working with this dataset is still at early stages and the land use classes are more abstract and difficult than land cover ones. However, considering our current issue, it is worth noting that training with all and 10% data leads to similar scores, both in accuracy and IoU. By picking our 10% sample images all over the dataset, we preserved the diversity of the training set and did not degrade the results too much (even if more data is better). On the contrary, training with a single location implies a 10% loss in accuracy and 5% less of mIoU. Clearly, the training dataset does not then offer enough variety to encompass all the potential images of the test set.

In a second experiment, we apply the model trained on the whole dataset to each city or conurbation of the test set. Results are shown in Table II. It is interesting to observe that the performance of the network varies notably between some conurbations, revealing differences between cities and a lack of generalization capacity from the model. Indeed, we observe in Fig. 4 that the statistical distribution of the pixel colors differs from train to test. Thus MiniFrance is a much more diverse dataset than many others, and offers exciting challenges to overcome. To better understand these numbers, Fig. 5 presents semantic maps obtained during testing. The first row shows a quite accurate prediction. Indeed, the model correctly identifies most classes present on the image. However, the predicted map appears more fragmented than the ground-truth, which shows that the network is sensitive to color variations of the image, and sometimes misses some abstract semantic classes. In the second row, we can observe an image where groundtruth is not reliable and for which the prediction obtained by our method seems actually more accurate.

V. ANALYSIS AND CONCLUSION

Experiments in sections III and IV aimed at shedding some light on the relationship between data and semantic segmentation networks. On small-scale, single-location datasets such as Vaihingen, these networks are robust to a drastic decrease of the amount of training data when testing on similar data. This is promising as practical applications can be solved with less labeled data than what is usually thought, but meanwhile is limited because models are likely to overfit. Going large-scale raises new issues. Indeed, images are more diverse, corresponding to various geographical idiosyncrasies and various acquisition conditions and times. Also, large-scale annotations are often approximate, due to crowd-sourcing or semi-automated processes. Additionally, a new challenge is the higher semantic level of the classes with respect to preexisting datasets which focus on structural classes (building, roads, etc.). Indeed, in crowdsourced maps such as OSM, classes are more symbolic and the classes are more land use oriented, related to the parcel and based on the interpretation of structural classes (e.g. building density in an area).

To promote research on large-scale EO analysis, we constituted MiniFrance, a new dataset for urban semantic segmentation. By using different cities for training and testing, we augment the variety of the dataset, which will require learning algorithms with good generalization capacities and the right bias/variance compromise [START_REF] Geman | Neural networks and the bias/variance dilemma[END_REF]. MiniFrance raises new challenges including predicting high level semantics classes from approximate ground truth or structured learning. In addition, in this paper we also proposed a first FCN baseline for further experimental comparison.

Fig. 1 :

 1 Fig. 1: Influence of the training set size (number of tiles) on the network performances, with and without weight initialization.

Fig. 2 :

 2 Fig. 2: Visualization of results of the different training phases over the Vaihingen dataset.

Fig. 3 :

 3 Fig. 3: MiniFrance dataset overview.

Fig. 4 :

 4 Fig. 4: Per channel color histograms over the Vaihingen and the MiniFrance datasets, comparing train and validation/test subsets.

Fig. 5 :

 5 Fig. 5: MiniFrance semantic segmentation results. From left to right: RGB image, ground truth and prediction. Legend of main classes: Urban fabric , Industrial, commercial, public, military, private and transport units , Complex and mixed cultivation patterns , Herbaceous vegetation associations .

TABLE I :

 I Classification performances with respect to amount of data.

	Train set	OA	mIoU
	100 %	52.40	15.79
	10 %	50.14	15.25
	Caen only (∼ 12.5%)	42.09	10.05

TABLE II :

 II Results by conurbation, over the entire train set. Marseille Rennes Angers Quimper Vannes Clermont Lille Cherbourg OA 46.13 51.56 44.85 50.82 49.51 46.51 61.35 67.54 mIoU 12.77 15.05 13.15 13.93 12.66 11.40 16.93 15.82

	Score

http://professionnels.ign.fr/bdortho

https://land.copernicus.eu/local/urban-atlas/urban-atlas-2012/view

Using a Titan X GPU, training over MiniFrance takes

hours, while testing takes 25 hours

On this table and for the rest of the document, OA stands for Overall Accuracy and mIoU is Mean IoU.