Noémie Janot

Pascal E Reiller

Marc F Benedetti
email: benedetti@ipgp.fr

Modelling Eu(III) speciation in a Eu(III)/PAHA/α-Al 2 O 3 ternary system

Keywords: Complexation, Adsorption, Fractionation, Modelling, Natural organic matter

Eu(III) behaviour is described under a various range of environmental conditions.

• CD-MUSIC and NICA-Donnan models are used to calculate Eu(III) speciation.

• Modifications of HA reactivity due to its fractionation are taken into account.

• Model shows that Eu(III) is preferentially bound to adsorbed humic acid fraction.

Introduction

Interactions between humic substances (HS) and mineral surfaces have a great influence on contaminant behaviour in the environment. Presence of humic substances greatly modifies the adsorption of metals onto minerals, and may enhance their solubility and potential transport [START_REF] Mccarthy | Mobilization of transuranic radionuclides from disposal trenches by natural organic matter[END_REF][START_REF] Mccarthy | Lanthanide Field Tracers Demonstrate Enhanced Transport of Transuranic Radionuclides by Natural Organic Matter[END_REF]. Several models can describe the different binary systems: Model VI [START_REF] Tipping | Humic Ion-Binding Model VI: An Improved Description of the Interactions of Protons and Metal Ions with Humic Substances[END_REF] or NICA-Donnan [START_REF] Kinniburgh | Ion binding to natural organic matter: competition, heterogeneity, stoichiometry and thermodynamic consistency[END_REF] for metal complexation to humic substances, CD-MUSIC for metal adsorption onto surfaces [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 | A Surface Structural Approach to Ion Adsorption: The Charge Distribution (CD) Model[END_REF], Ligand and Charge Distribution (LCD) model for humic substances retention onto oxides [START_REF] Filius | Modeling the Binding of Benzenecarboxylates by Goethite: The Ligand and Charge Distribution Model[END_REF][START_REF] Filius | Modeling the binding of fulvic acid by goethite: the speciation of adsorbed FA molecules[END_REF]. However, the additivity of the binary systems -i.e., metal/HS, metal/surface, and HS/surface is often not respected for the ternary systems metal/HS/surface [START_REF] Christl | Interaction of copper and fulvic acid at the hematite-water interface[END_REF][START_REF] Vermeer | Metal Ion Adsorption to Complexes of Humic Acid and Metal Oxides: Deviations from the Additivity Rule[END_REF]. Some studies have succeeded in describing ternary systems containing fulvic acids [START_REF] Heidmann | Sorption of Cu and Pb to kaolinite-fulvic acid colloids: Assessment of sorbent interactions[END_REF][START_REF] Weng | Cu2+ and Ca2+adsorption to goethite in the presence of fulvic acids[END_REF], but the larger size and higher heterogeneity of humic acid makes it more difficult to describe [START_REF] Weng | Humic Nanoparticles at the Oxide-Water Interface: Interactions with Phosphate Ion Adsorption[END_REF]. In fact, like real polymers [START_REF] Blaakmeer | Adsorption of Weak Poly Electrolytes on Highly Charged Surfaces. Poly(acrylic acid) on Polystyrene Latex with Strong Cationic Groups[END_REF], humic acids adsorption onto minerals increases with ionic strength [START_REF] Christl | Interaction of copper and fulvic acid at the hematite-water interface[END_REF][START_REF] Fernandes | Spectroscopic Identification of Ternary Cm-Carbonate Surface Complexes[END_REF][START_REF] Weng | Adsorption of Humic Substances on Goethite: Comparison between Humic Acids and Fulvic Acids †[END_REF][START_REF] Weng | Adsorption of humic acids onto goethite: Effects of molar mass, pH and ionic strength[END_REF], which is not observed for fulvic acids [START_REF] Filius | Adsorption of fulvic acid on goethite[END_REF], and is often different from the behaviour of small organic molecules [START_REF] Ali | Competitive sorption of simple organic acids and sulfate on goethite[END_REF]. Recently, Weng et al. [START_REF] Weng | Effects of Fulvic and Humic Acids on Arsenate Adsorption to Goethite: Experiments and Modeling[END_REF] succeeded to describe the behaviour of arsenate in a humic acid-containing ternary system; however, arsenate has a different behaviour towards colloids than metallic ions.

Despite numerous studies of these ternary systems, it is still difficult to predict the speciation of metal ions [START_REF] Reiller | Modelling metal-humic substancessurface systems: reasons for success, failure and possible routes for peace of mind[END_REF][START_REF] Bryan | The effects of humic substances on the transport of radionuclides: Recent improvements in the prediction of behaviour and the understanding of mechanisms[END_REF], a major reason being the fractionation of humic substances constituents when they adsorb onto mineral surfaces [START_REF] Meier | Fractionation of aquatic natural organic matter upon sorption to goethite and kaolinite[END_REF][START_REF] Hur | Molecular weight fractionation of humic substances by adsorption onto minerals[END_REF][START_REF] Reiller | Sorption of Aldrich Humic Acid onto Hematite: Insights into Fractionation Phenomena by Electrospray Ionization with Quadrupole Time-of-Flight Mass Spectrometry[END_REF][START_REF] Claret | Fractionation of Suwannee River Fulvic Acid and Aldrich Humic Acid on α-Al 2 O 3 : Spectroscopic Evidence[END_REF][START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF]. This phenomenon has been characterized by several techniques, including UV-visible spectroscopy, size-exclusion chromatography, time-resolved luminescence spectroscopy (TRLS), and asymmetric flow-field flow fractionation. These techniques have been applied to different ternary systems, with contrasting results depending on the mineral surface, the HS origin, and the solution parameters such as pH and ionic strength. The adsorption of HS onto mineral surfaces induces a mass/size fractionation, as noted by several authors [START_REF] Meier | Fractionation of aquatic natural organic matter upon sorption to goethite and kaolinite[END_REF][START_REF] Hur | Molecular weight fractionation of humic substances by adsorption onto minerals[END_REF][START_REF] Pitois | Humic acid sorption onto a quartz sand surface: A kinetic study and insight into fractionation[END_REF][START_REF] Zhou | Size fractionation upon adsorption of fulvic acid on goethite: equilibrium and kinetic studies[END_REF], but chemical fractionation is also occurring. Spectrophotometric studies show preferential adsorption of the most aromatic moieties [START_REF] Reiller | Sorption of Aldrich Humic Acid onto Hematite: Insights into Fractionation Phenomena by Electrospray Ionization with Quadrupole Time-of-Flight Mass Spectrometry[END_REF][START_REF] Claret | Fractionation of Suwannee River Fulvic Acid and Aldrich Humic Acid on α-Al 2 O 3 : Spectroscopic Evidence[END_REF][START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF][START_REF] Zhou | Size fractionation upon adsorption of fulvic acid on goethite: equilibrium and kinetic studies[END_REF]. The two HS fractions -the one in solution and the one onto the surface -present different composition and affinity towards protons and ions [START_REF] Heidmann | Sorption of Cu and Pb to kaolinite-fulvic acid colloids: Assessment of sorbent interactions[END_REF][START_REF] Claret | Fractionation of Suwannee River Fulvic Acid and Aldrich Humic Acid on α-Al 2 O 3 : Spectroscopic Evidence[END_REF][START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF].

A new major step to successfully describe a metal/HA/surface system is to characterize the reactivity of HA fractions created upon adsorption. However, this is difficult due to the low concentrations of HA in the supernatant from adsorption experiments. Recently, the use of spectrophotometric titrations has facilitated the measurement of environmentally relevant humic acid concentrations, i.e., less than 10 mg/L [START_REF] Dryer | Situ Examination of the Protonation Behavior of Fulvic Acids Using Differential Absorbance Spectroscopy[END_REF][START_REF] Janot | Using Spectrophotometric Titrations To Characterize Humic Acid Reactivity at Environmental Concentrations[END_REF]. In Janot et al. [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF], we applied this method to determine the modification of proton-binding behaviour of HA fractions after adsorption onto α-Al 2 O 3 . Here, Eu(III) speciation is studied in presence of the same purified Aldrich humic acid (PAHA) and α-Al 2 O 3 surface. Eu(III) is considered to be (i) a representative of the lanthanides, which use is increasing in modern industry and which are important for understanding geochemical processes; and (ii) as an analogue of some actinides(III). We have already described the evolution of the Eu(III) chemical environment in the binary and ternary systems using time-resolved luminescence spectroscopy (TRLS) [START_REF] Janot | Influence of solution parameters on europium(III), α-Al2O3, and humic acid interactions: Macroscopic and timeresolved laser-induced luminescence data[END_REF], which guides our modelling hypotheses. Aluminium oxides are not the most common minerals in the environment, but their surface sites are used here as an analogue of the aluminol surface sites of clays [START_REF] Fernandes | Spectroscopic Identification of Ternary Cm-Carbonate Surface Complexes[END_REF], which are ubiquitous in natural systems. Here, we use the CD-MUSIC and NICA-Donnan models to construct a model describing the repartitioning of Eu(III) within the different compartments of the system (solution, surface, adsorbed-PAHA, and dissolved-PAHA).

Materials and methods

Materials

The α-Al 2 O 3 sample was purchased from Interchim (pure 99.99%, size fraction 200-500 nm). As in previous studies [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF][START_REF] Janot | Using Spectrophotometric Titrations To Characterize Humic Acid Reactivity at Environmental Concentrations[END_REF][START_REF] Janot | Influence of solution parameters on europium(III), α-Al2O3, and humic acid interactions: Macroscopic and timeresolved laser-induced luminescence data[END_REF][START_REF] Janot | Colloidal α-Al2O3, Europium(III) and Humic Substances Interactions: A Macroscopic and Spectroscopic Study[END_REF], the solid was washed thrice with carbonate-free NaOH and thrice with milli-Q water before drying and storage at room temperature under N 2 (g) atmosphere. Specific area was measured by N 2 -BET method and found to be 15 m 2 /g. Commercial Aldrich humic acid was purified (PAHA) according to Kim et al. [START_REF] Kim | Characterization of humic and fulvic acids from Gorleben groundwater[END_REF]. Stock suspension at 5 g PAHA /L was prepared by diluting PAHA in NaOH (pH around 10) overnight to completely dissolve the sample.

Potentiometric titrations of a-Al 2 O 3

Potentiometric titrations were performed in a computer-controlled system under N 2 (g) atmosphere in a thermostatic vessel (25 °C) equipped with a stirrer using 50 mL solutions at 100 g/L α-Al 2 O 3 . Ionic strength I was fixed using KNO 3 (for 0.01 M and 1 M titrations) or NaClO 4 (0.1 M titration). The quoted ionic strengths are the initial values before the addition of any titrant. Experimental values of I were calculated for every data point, accounting for both background electrolyte ions and free H + and OH -. The pH values of the solutions were controlled during titrations by addition of 0.1 mol/L acid and base solutions. Base titrant (titrisol for 0.1 mol/L solution) was prepared with degassed Millipore water. The pH values were recorded with two pH Metrohm 6.0133.100 glass electrodes and a single Metrohm 6.0733.100 reference electrode. The pH electrodes were calibrated by performing a blank titration of the background electrolyte prior to the sample titration. The aqueous suspension was purged with N 2 (g) at pH 4 for 90 min. The suspension was then titrated by adding small volumes of titrant, and pH was recorded as a function of the titrant volume added to the suspension. After each addition, a drift criterion for pH was used (ΔmV/min < 0.1) and a maximum waiting time of 60 min was set for acquiring each data point. Usually stabilization of pH took less than 2 min. A similar procedure was followed for the blank solution titration. Two forward and backward titrations were performed to eliminate the hysteresis effect.

Adsorption experiments

The methodology of these experiments have been described in a previous study [START_REF] Janot | Colloidal α-Al2O3, Europium(III) and Humic Substances Interactions: A Macroscopic and Spectroscopic Study[END_REF]. Briefly, batch experiments were carried out at ambient temperature with 10 -6 mol/L Eu(III) concentration at different pH, ionic strengths (0.01 M and 0.1 M NaClO 4 ) and various concentrations of organic and/or inorganic compounds. To minimize carbonation of the systems, and formation of Eu(CO 3 ) n 3-n complexes, headspace was purged with nitrogen prior to closing the tubes. The pH values of the solutions were adjusted by addition of fresh 0.01 or 0.1 mol/L HClO 4 or NaOH solutions. They were measured using one combined glass electrode (Radiometer Analytical) connected to a Radiometer Analytical pH Meter. To prevent KClO 4 precipitation in the electrode frit, and to minimize junction potential, the electrode filling solution was modified to use 0.1 M NaCl. The pH electrode was calibrated using buffer solutions (pH 4.01, 7.01, and 10.00), yielding σ pH of 0.05-0.08. After 3 days of equilibration, the final pH was checked and samples were ultracentrifuged at 60,000 rpm during 2 h. Eu(III) concentration in the supernatant was measured by ICP-AES or by standard addition method in TRLS, diluting an aliquot of the supernatant in 3 mol/L K 2 CO 3 [START_REF] Reiller | Dual Use of Micellar-Enhanced Ultrafiltration and Time-Resolved Laser-Induced Spectrofluorimetry for the Study of Uranyl Exchange at the Surface of Alkylsulfate Micelles[END_REF].

Spectroscopic measurements

We studied Eu(III) speciation in binary and ternary systems by TRLS. A detailed description of the installation can be found elsewhere [START_REF] Janot | Colloidal α-Al2O3, Europium(III) and Humic Substances Interactions: A Macroscopic and Spectroscopic Study[END_REF]. Briefly, a luminescence signal was collected after a gate delay D = 10 μs after the excitation by a laser flash, during a gate width W = 300 μs. To increase the signal to noise ratio, 1000 accumulations were performed for each spectrum. To avoid aggregation and settling of samples in the α-Al 2 O 3 -containing systems, samples were manually shaken between two measurements. All luminescence measurements were performed at ambient temperature (21 ± 2 °C). The excitation wavelength was set at λ exc = 393.8 nm, which corresponds to the 7 F 0 → 5 L 6 transition of Eu(III). The observed luminescence corresponds to the transitions from the 5 D 0 excited state to the ground 7 F j manifold [START_REF] Bünzli | Luminescent Probes[END_REF]. They are the 5 D 0 → 7 F 0 transition, forbidden for magnetic and electric reasons 3/12 10.1016/j.colsurfa.2013.02.052 (λ max ≈ 580 nm), the 5 D 0 → 7 F 1 magnetic dipole transition (λ max ≈ 593 nm), and the 5 D 0 → 7 F 2 hypersensitive electric dipole transition (λ max ≈ 615 nm) [START_REF] Judd | Optical Absorption Intensities of Rare-Earth Ions[END_REF][START_REF] Ofelt | Intensities of Crystal Spectra of Rare-Earth Ions[END_REF][START_REF] Jørgensen | Hypersensitive pseudoquadrupole transitions in lanthanides[END_REF]. Photo-degradation [START_REF] Monsallier | Influence of photochemical reactions on the complexation of humic acid with europium(III)[END_REF] of the humic acid can be neglected due to the relatively low amount of energy provided to the system, as described previously [START_REF] Janot | Colloidal α-Al2O3, Europium(III) and Humic Substances Interactions: A Macroscopic and Spectroscopic Study[END_REF][START_REF] Brevet | Spectral and temporal luminescent properties of Eu(III) in humic substance solutions from different origins[END_REF][START_REF] Reiller | Bi-exponential decay of Eu(III) complexed by Suwannee River humic substances: Spectroscopic evidence of two different excited species[END_REF][START_REF] Reiller | Europium(III) complexed by HPSEC size-fractions of a vertisol humic acid: Small differences evidenced by time-resolved luminescence spectroscopy[END_REF].

Results and Discussion

Surface charge of a-Al 2 O 3

The experimental charging curves are presented in Fig. 1, together with the results of the fits obtained using CD-MUSIC model [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 | A Surface Structural Approach to Ion Adsorption: The Charge Distribution (CD) Model[END_REF]. Calculations were done using the ECOSAT software [START_REF] Keizer | ECOSAT, a computer program for the calculation of chemical speciation and transport in soil-water systems[END_REF]. The pH PZC was found to be 9.6, meaning that the storage procedure managed to avoid carbonation of the surface. The oxide surface is positively charged below pH PZC , with charge increasing with ionic strength.

The application of Pauling valence rules to α-Al 2 O 3 structure, in which the Al ions distribute their charge to the six surrounding O atoms, results in a charge of -0.5 per ≡AlOH site. Those singly coordinated surface groups can be protonated following the reaction

≡ AlOH 1/ 2 + H + ⇔ ≡ AlOH 2 +1 /2 ¿ K H in which K H is the protonation constant.
In the α-Al 2 O 3 structure, doubly coordinated surface sites ≡Al 2 OH also exist, but they are uncharged in the pH range studied and are not taken into account for the description of surface charging behaviour [START_REF] Hiemstra | Interfacial charging phenomena of aluminum (Hydr)oxides[END_REF]. The fitted parameters are summarized in Table 1. The electrostatic double layer of the water-oxide interface was described using the Basic Stern model [START_REF] Westall | A comparison of electrostatic models for the oxide/solution interface[END_REF]. Ion pair formation of surface groups with electrolyte ions is taken into account, and described with the following equations

≡ AlOH -1/ 2 + C + ⇔ ≡AlOH -1 /2 -C + ¿ K C ≡ AlOH 2 +1/ 2 + A -⇔ ≡ AlOH 2 + 1/ 2 -A -¿ K A
where C + and A -are the cation and the anion from the electrolyte, respectively.

The electrolyte ions were all positioned in the outer plane of the Stern layer. Evidence of asymmetric binding of electrolyte ions on gibbsite has been previously reported [START_REF] Hiemstra | Interfacial charging phenomena of aluminum (Hydr)oxides[END_REF], so different affinity constants have been determined for the different electrolyte used. Their values were derived from those found for the charging of goethite [START_REF] Rietra | Electrolyte Anion Affinity and Its Effect on Oxyanion Adsorption on Goethite[END_REF]. The values obtained for cation pair formation are in accord with the ones of Rietra et al. [START_REF] Rietra | Electrolyte Anion Affinity and Its Effect on Oxyanion Adsorption on Goethite[END_REF]. However, NO 3 -ions were found to have very weak interaction with the alumina surface (logK NO3 -= -1.6), which is much lower than the observations made by Rietra et al. [START_REF] Rietra | Electrolyte Anion Affinity and Its Effect on Oxyanion Adsorption on Goethite[END_REF] who found (logK NO3 -= -1) with goethite. The affinity constants for ClO 4 -are 4/12 coherent in the two studies: -1.4 in our case, -1.7 for Rietra et al. [START_REF] Rietra | Electrolyte Anion Affinity and Its Effect on Oxyanion Adsorption on Goethite[END_REF].

The obtained capacitance (C = 2.2 F/m²) is quite high compared to the values proposed by Hiemstra and van Riemsdijk [START_REF] Hiemstra | On the relationship between charge distribution, surface hydration, and the structure of the interface of metal hydroxides[END_REF] for metal (hydr)oxides surfaces (C ≈ 1 F/m²). It is however similar to the value found by Rabung et al. [START_REF] Rabung | Sorption of Eu(III) on a Natural Hematite: Application of a Surface Complexation Model[END_REF] in the case of hematite (C = 2.24 F/m²). A previously observed surface transformation of α-Al 2 O 3 occurred after hydration, which leads to a surface structure closer to gibbsite -Al(OH) 3 -surface [START_REF] Eng | Structure of the hydrated α-Al2O3 (0001) surface[END_REF]. After 10 h titrations in solution, our sample of α-Al 2 O 3 may develop gibbsite-like, or hydroxide-like, surface properties.

Eu(III) adsorption onto α-Al 2 O 3

Eu(III) adsorption onto α-Al 2 O 3 has been measured in Janot et al. [START_REF] Janot | Influence of solution parameters on europium(III), α-Al2O3, and humic acid interactions: Macroscopic and timeresolved laser-induced luminescence data[END_REF] at different pH and two ionic strengths, 0.01 M and 0.1 M NaClO 4 (Fig. 2). As expected, it shows a great dependency to pH and a low dependency to ionic strength. Together with TRLS data showing a loss of water molecules in the first coordination sphere of the ion when surface complexes are formed [START_REF] Janot | Influence of solution parameters on europium(III), α-Al2O3, and humic acid interactions: Macroscopic and timeresolved laser-induced luminescence data[END_REF], this indicates inner-sphere complexation [START_REF] Rabung | Sorption of Eu(III) on a Natural Hematite: Application of a Surface Complexation Model[END_REF][START_REF] Wang | Sorption and complexation of Eu(III) on alumina: Effects of pH, ionic strength, humic acid and chelating resin on kinetic dissociation study[END_REF][START_REF] Tan | Adsorption of Eu(III) onto TiO2: Effect of pH, concentration, ionic strength and soil fulvic acid[END_REF]. The influence of metal/surface ratio has also been studied at pH = 6.2 ± 0.2 by changing α-Al 2 O 3 concentration in the system (see Figure S1 To confirm the accuracy of this model, the parameters were then used to describe independent data sets from literature of Eu(III) binding to hydrous alumina [START_REF] Reiller | Modelling metal-humic substancessurface systems: reasons for success, failure and possible routes for peace of mind[END_REF] and γ-Al 2 O 3 [START_REF] Bryan | The effects of humic substances on the transport of radionuclides: Recent improvements in the prediction of behaviour and the understanding of mechanisms[END_REF]. Fig. 3 shows that the results of modelling are in good agreement for Eu(III) concentrations close to our conditions, i.e., between 6 × -6 and 6 × 10 -7 mol/L, even if the mineral surfaces used in these studies were different from our surface, especially with higher specific area (above 100 m 2 /g compared to our 15 m 2 /g surface). For higher concentrations of Eu(III) (above 10 -5 mol/L), the model is slightly overestimating Eu(III) adsorption, with pH-edges found at lower pH values than in experiments (see blue line on Fig. 3a and blue, red and green lines on Fig. 3b). The ability of our model to describe literature data then validates a range of concentrations in which this model can be used.

Eu(III) binding to humic acid

Eu(III) binding to PAHA has been studied by TRLS at pH 4 and different PAHA concentrations. Eu(III) luminescence spectra present an hypersensitive transition, 5 D 0 → 7 F 2 , the intensity of which is enhanced in the metal-ligand complex relative to its intensity for the aquo ion [START_REF] Jørgensen | Hypersensitive pseudoquadrupole transitions in lanthanides[END_REF]. Consequently, the change of relative intensities of 5 D 0 → 7 F 2 and 5 D 0 → 7 F 1 transitions indicates a change in Eu(III) environment symmetry. This asymmetry ratio 5 D 0 → 7 F 2 / 5 D 0 → 7 F 1 (noted 7 F 2 / 7 F 1 ) evolves with the relative concentration of metal and humic acid [START_REF] Janot | Influence of solution parameters on europium(III), α-Al2O3, and humic acid interactions: Macroscopic and timeresolved laser-induced luminescence data[END_REF]. Above 20 mg/L of PAHA in the system, the asymmetry ratio is constant, suggesting that all the Eu(III) is complexed to humic acid. When PAHA is absent from the system, the asymmetry ratio value is the value corresponding to free Eu(III). It is then possible to use this evolution as a spectral titration curve [START_REF] Dobbs | A novel approach to metal-humic complexation studies by lanthanide ion probe spectroscopy[END_REF] (see Figure S2 of SI). [START_REF] Bryan | The effects of humic substances on the transport of radionuclides: Recent improvements in the prediction of behaviour and the understanding of mechanisms[END_REF]. Symbols are the data from literature, and lines are the results of our modelling.

6/12

10.1016/j.colsurfa.2013.02.052

We used these asymmetry ratios as indicators of the proportion of Eu(III) complexed, to adjust the NICA-Donnan parameters for Eu(III)/humic acid binding [START_REF] Milne | Generic NICA-Donnan Model Parameters for Metal-Ion Binding by Humic Substances[END_REF][START_REF] Marang | Combining Spectroscopic and Potentiometric Approaches to Characterize Competitive Binding to Humic Substances[END_REF], as shown in Fig. 4. Modelling has been done using ECOSAT software. The NICA-Donnan equation for metal binding to humic substances can be found elsewhere, together with the definition of the different parameters [START_REF] Milne | Generic NICA-Donnan Model Parameters for Metal-Ion Binding by Humic Substances[END_REF].

The adjustment of Eu(III)-binding parameters was done at pH 4 only, because 100% complexation was immediately reached when PAHA was added in the system at pH 6. At higher pH, due the presence of hydrolyzed species of Eu(III) in solution, it is impossible to have the asymmetry ratio value corresponding to the absence of complexation. Thus, we only adjusted the parameters for the low affinity sites (S1). Results are given in Table 3.

logK

Eu(III) speciation in the ternary system

A ternary system is difficult to describe due to HA fractionation and modification of reactivity after adsorption onto the mineral surface [START_REF] Meier | Fractionation of aquatic natural organic matter upon sorption to goethite and kaolinite[END_REF][START_REF] Hur | Molecular weight fractionation of humic substances by adsorption onto minerals[END_REF][START_REF] Reiller | Sorption of Aldrich Humic Acid onto Hematite: Insights into Fractionation Phenomena by Electrospray Ionization with Quadrupole Time-of-Flight Mass Spectrometry[END_REF][START_REF] Claret | Fractionation of Suwannee River Fulvic Acid and Aldrich Humic Acid on α-Al 2 O 3 : Spectroscopic Evidence[END_REF]. In this work, spectrophotometric titrations of PAHA supernatants presented in previous studies [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF], [START_REF] Janot | Using Spectrophotometric Titrations To Characterize Humic Acid Reactivity at Environmental Concentrations[END_REF] are used to quantify proton-binding behaviour of humic fractions. Two pools of PAHA are defined: one staying in solution after centrifugation of the suspension (PAHA sol ) and one adsorbed onto the mineral (PAHA ads ), each fraction showing different NICA-Donnan parameters. Titrations of original PAHA material and modelling of its charge before fractionation using the NICA-Donnan model have been published elsewhere [START_REF] Janot | Using Spectrophotometric Titrations To Characterize Humic Acid Reactivity at Environmental Concentrations[END_REF]. In a previous study on the binary PAHA/α-Al 2 O 3 system, we used spectrophotometric titrations of supernatant from adsorption experiments to determine modifications of the proton binding parameters of the two PAHA fractions [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF]. Here, we take these modifications into account for modelling the behaviour of Eu(III) in presence of both PAHA and α-Al 2 O 3 .

In our conditions -at 0.1 M NaClO 4 and an initial PAHA concentration of 28 mg/L -between 90 and 30% of PAHA is adsorbed onto α-Al 2 O 3 in the ternary system, depending on pH [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF]. Except at pH < 4, fractionation rate of PAHA is therefore moderate. In order to test the validity of our approach, we then chose to modify the proton-binding parameters of the two fractions to values corresponding to a moderate fractionation rate.

Spectrophotometric titrations of the dissolved fraction after adsorption did not show any change of the proton-binding site density of high-affinity sites Q 2 , neither of the apparent heterogeneity parameters m 1 and m 2 compared to the initial material [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF]. Therefore, we did not change these parameters. However, the proton-binding site density of the low affinity sites Q 1 and the median values of the affinity distribution for proton for both types of sites logK H,1 and logK H,2 were modified: Q 1 and logK H,1 were lowered, and logK H,2 increased compared to the original compound, as seen by spectrophotometric titrations of nonsorbed fractions [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF]. The same calculation was performed for PAHA ads : from spectrophotometric titrations and mass balance calculations, we ob-7/12 The logK H,1 values for protonation varied for the two fractions compared to the original material. Consequently, the logK Eu,1 values for Eu(III)/ PAHA binding determined previously (see Table 3) were modified in the same proportions. Following the same reasoning, the logK Eu,2 value for Eu(III) binding to PAHA sol was updated as well. The NICA-Donnan parameters taken into account for modelling Eu(III) binding to dissolved and adsorbed PAHA fractions are given in Table 3. For metal binding to humic acid, generic values of intrinsic heterogeneity parameters p 1 (0.62) and p 2 (0.41) were used [START_REF] Milne | Generic NICA-Donnan Model Parameters for Proton Binding by Humic Substances[END_REF]. The m i values for proton binding were the same in original PAHA and fractions, so values of the non-ideality ion-binding parameters n Eu,i for Eu(III) binding were kept unchanged.

Using these parameters, the surface concentration of Eu(III) in each batch experiment was calculated, using as inputs of the model the values of pH, Eu(III) and α-Al 2 O 3 concentrations, as well as amounts of both PAHA fractions. Eu(III) and α-Al 2 O 3 concentrations may slightly vary between two experiments due to modifications of solution volume when adjusting the pH value. These differences within a series explain the otherwise apparently peculiar variations of model for some data. Model values were calculated for the experiments performed at different pH and I (see circles in Fig. 2), and for experiments performed at constant pH and various initial PAHA concentrations (see Fig. 5). They show theoretical and actual surface con-centration of Eu(III), i.e. Eu(III) bound to the surface and to adsorbed PAHA.

Under these modelling hypotheses, the behaviour of Eu(III) is well described in a ternary Eu(III)/PAHA/α-Al 2 O 3 system depending on pH, especially at 0.1 M. Modification of PAHA binding parameters with fractionation are determined from adsorption experiments performed at pH ≈ 7 and I = 0.1 M [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF], thus results of the modelling provide better fits to experimental data at this ionic strength value and neutral pH. The trends of Eu(III) behaviour depending on initial PAHA concentration are captured, but the values are not always accurate. The isotherm at pH ≈ 6.2 is well described, whereas the Eu(III) adsorption at pH ≈ 7.4 seems underestimated. The behaviour of Eu(III) at pH 4 is less well described, even if the general trend is present. These conditions corre-8/12 Table 4. NICA-Donnan proton-binding intrinsic parameters for PAHA before contact with the oxide (PAHA i ), PAHA fraction in solution (PAHA sol ), and adsorbed PAHA fraction (PAHA ads ), determined from spectrophotometric titrations [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF], [START_REF] Janot | Using Spectrophotometric Titrations To Characterize Humic Acid Reactivity at Environmental Concentrations[END_REF].

In bold are the modified parameters used to model the PAHA fractions compared to PAHA i . spond to a high fractionation rate, with a large percentage of PAHA being adsorbed and PAHA reactivity being greatly modified [START_REF] Janot | Characterization of humic acid reactivity modifications due to adsorption onto α-Al2O3[END_REF], [START_REF] Janot | Using Spectrophotometric Titrations To Characterize Humic Acid Reactivity at Environmental Concentrations[END_REF]. The disagreement between experimental and calculated concentrations under these conditions may be due to the values chosen for adapted protonation parameters of both PAHA fractions, which correspond to a median fractionation rate (around 70%). Fig. 6 shows the model results of pH-dependent Eu(III) speciation. It is observed that between pH 4 and 8, Eu(III) is mainly bound to PAHA ads . The influence of the mineral surface becomes significant only above pH 9, which is in agreement with spectroscopic observations [START_REF] Janot | Colloidal α-Al2O3, Europium(III) and Humic Substances Interactions: A Macroscopic and Spectroscopic Study[END_REF]. TRLS spectra of Eu(III) in the ternary system showed an apparition of surface-like features for pH values higher than 8: decrease of asymmetry ratio, broadening of 5 D 0 → 7 F 1 transition, shifting of 5 D 0 → 7 F 2 apparent maximum towards higher wavelengths [START_REF] Janot | Influence of solution parameters on europium(III), α-Al2O3, and humic acid interactions: Macroscopic and timeresolved laser-induced luminescence data[END_REF].

At 0.01 M NaClO 4 , macroscopic observations showed less adsorbed PAHA, so the surface concentration of Eu(III) is lower, but Eu(III) repartitioning is similar (data not shown). Eu(III) repartitioning at different surface coverage ratios is shown in Figure S3 of SI. At pH 6.2 and 7.4, Eu(III) is partially complexed to the surface for surface coverage ratios below 10 mg PAHA /g α-Al2O3 . When PAHA concentrations increase, Eu(III) is totally bound to PAHA, mostly to the adsorbed fraction. At pH 4, the surface has no influence, in agreement with binary Eu(III)/α-Al 2 O 3 observations. Eu(III) surface concentration is only due to its complexation with adsorbed PAHA. Below 20 mg PAHA /g α-Al2O3 , a significant proportion of Eu(III) is free in solution, which is not seen at higher pH values.

Conclusions

In this model of contaminant behaviour in a ternary metal/HA/mineral system, we took into account the modification of humic acid reactivity due its fractionation after adsorption onto a mineral surface. We used results from spectrophotometric titrations of PAHA supernatants from adsorption experiments performed at pH ≈ 7 and I = 0.1 M. The results of our approach show a good description of experimental observations at different pH, ionic strength and humic acid concentration values, thus providing model validation. However, the description is less accurate when fractionation rate of the humic acid is high, such as at low pH. Performing spectrophotometric titrations of supernatants from adsorption experiments at other pH values would allow a better characterization of PAHA reactivity modifications. Moreover, the protonation parameters used for the modelling were calculated from titrations made on supernatant from binary PAHA/α-Al 2 O 3 system. The presence of Eu(III) may have an influence on the fractionation of humic moieties as complexation presumably primarily involves the carboxylic 9/12 10.1016/j.colsurfa.2013.02.052 groups at low pH values [START_REF] Naber | Scanning transmission X-ray and laser scanning luminescence microscopy of the carboxyl group and Eu(III) distribution in humic acid aggregates[END_REF]. Still, the results of this study show the importance of the reactivity modifications of the humic substances, and provide a step towards a better understanding and description of these complex ternary systems.

Figure 1 .

 1 Figure 1. Surface charge of α-Al 2 O 3 at different pH and ionic strengths. The dotted lines are the results of modelling.

  of Supplementary Information). At high Eu(III)/α-Al 2 O 3 ratio, with 0.5 g/L α-Al 2 O 3 , 58% of initial Eu(III) is retained onto the surface. The proportion of adsorbed Eu(III) increases with α-Al 2 O 3 concentration, and reaches 100% at approximately 5 g/L α-Al 2 O 3 . Modelling in the framework of the CD-MUSIC model has been done using ECOSAT software. Data for Eu(III) speciation in solution has been taken from Hummel et al. [53]. Spectroscopic re-5/12

Figure 2 .

 2 Figure 2. Surface concentration of Eu(III) in the binary Eu(III)/α-Al 2 O 3 (squares) and ternary (circles) systems depending on pH and I [adapted from 32]. Lines are the results of modelling. [Eu(III)] = 10 -6 mol/L; C(α-Al 2 O 3 ) = 1 g/L; Error bars correspond to 2σ.

Fig. 3 .

 3 Fig. 3. Eu(III) surface concentration on (a) 4.4 g/L hydrous alumina (105 m²/g), for I = 0.1 M KNO 3 and [Eu(III)] = 4.3 10 -5 mol/L (squares) and 0.1 M NaClO 4 and [Eu(III)] = 6.0 10 -7 mol/L (diamonds) [21]; (b) 3.6 g/L γ-Al 2 O 3 (119 m 2 /g) for I = 0.1 M NaClO 4 and various Eu(III) concentrations[START_REF] Bryan | The effects of humic substances on the transport of radionuclides: Recent improvements in the prediction of behaviour and the understanding of mechanisms[END_REF]. Symbols are the data from literature, and lines are the results of our modelling.

Fig. 4 .

 4 Fig. 4. Results of adjustment of NICA-Donnan parameters for Eu(III)/PAHA binding at pH 4;[Eu] total = 10 -6 mol/L; I = 0.1 M NaClO4 . The generic parameters come from Milne et al.[START_REF] Milne | Generic NICA-Donnan Model Parameters for Metal-Ion Binding by Humic Substances[END_REF].

Fig. 5 .

 5 Fig. 5. Eu(III) surface concentration in the ternary system depending on pH and coverage ratio. Lines are the results of modelling. [Eu(III)] = 10 -6 mol/L; C(α-Al 2 O 3 ) = 1 g/L; I = 0.1 M NaClO 4 , Error bars correspond to 2σ.

Fig. 6 .

 6 Fig. 6. Eu(III) speciation in the ternary system depending on pH. I = 0.1 M NaClO 4 , [PAHA] = 28 mg/L, C(α-Al 2 O 3 ) = 1 g/L, [Eu(III)] = 10 -6 mol/L.

  10.1016/j.colsurfa.2013.02.052

	Specific area (m²/g)	15
	Site density (sites/nm²)	8
	Stern layer capacitance (F/m²)	2.2
	logKH+	9.6
	logKNa+	-0.8
	logKK+	-1.5
	logKNO3-	-1.6
	logKClO4-	-1.4

Table 1 .

 1 Fitted parameters used for describing the charging behaviour of α-Al 2 O 3 surface.

Table 2 .

 2 Fitted CD-MUSIC parameters for Eu(III) complexation to α-Al 2 O 3 . 0 and Dz 1 represent the charge distribution on the inner and outer planes of the Stern layer, respectively.

	Dz

Table 3 .

 3 Adjusted NICA-Donnan parameters used for Eu(III) binding to PAHA and PAHA fractions after contact with mineral surface (dissolved and adsorbed). a Milne et al. [55]. 10.1016/j.colsurfa.2013.02.052 served an increase in Q 1 and logK H,1 values regarding the values fitted for non-fractionated PAHA [27]. The calculated parameters for both fractions PAHA sol and PAHA ads are reported in Table 4.

Acknowledgment. We thank Benjamin Kocar (SLAC) for editing the language and two anonymous reviewers for their comments which improved the manuscript. This work was financed through the MRISQ (DEN/DISN/AvC) project of CEA.

Supporting Information Available.


B Lopez

P Antonelli

S Wolf

S Lagarde

W Jaffe

R Navarro

U Graser

R Petrov

G Weigelt

Y Bresson

K H Hofmann

U Beckman

T Henning

W Laun

Ch Leinert

S Kraus

S Robbe-Dubois

F Vakili

A Richichi

P Abraham

J.-C Augereau

J Behrend

Ph Berio

N Berruyer

O Chesneau

J M Clausse

C Connot

K Demyk

W C Danchi

M Dugué

G Finger

S Flament

A Glazenborg

H Hannenburg

M Heininger

Y Hugues

J Hron

S Jankov

F Kerschbaum

G Kroes

H Linz

J.-L Lizon

Ph Mathias

R Mathar

A Matter

J L Menut

K Meisenheimer

F Millour

N Nardetto

U Neumann

E Nussbaum

A Niedzielski

L Mosoni

J Olofsson

Y Rabbia

T Ratzka

F Rigal

A Roussel

D Schertl

F.-X Schmider

B Stecklum

E Thiebaut

M Vannier

B Valat

K Wagner

L B F M Waters

MATISSE, perspective of imaging in the mid-infrared at the VLTI

Keywords: Interferometry, Mid-Infrared Spectral Domain, VLTI

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The Very Large Telescope of the European Southern Observatory (ESO) recently opened a new era of astronomy for Europe. In particular, mid-infrared interferometry with the Mid-Infrared Interferometric Instrument for the VLTI, MIDI, operating since the end of the year 2002. It allows spatially resolved observations of reemission regions of hot dust in circumstellar disks, AGB stars, winds of hot stars, and the tori of AGNs with a resolution of 10-20 mas.

MIDI was proven to be very successful in interferometric and spectroscopic observations since it allows the comparison of the chemical composition of dust on very different spatial scales. However, the investigation of small-scale spatial structures in general, and the quantitative analysis of spectroscopic observations in particular, are strongly limited due to the small number of visibility points measured in a reasonable amount of time (over one or a few nights) and due to the lack of phase information.

Furthermore, the comparison between models and visibility points in the Fourier plane is currently the main strategy for the interpretation of MIDI observations. The selection of more realistic models and the removal of possible ambiguities in the model fitting and resulting interpretations require a sufficient sampling of the uv-plane and the reconstruction of model-independent images. Moreover, the classical analysis usually relies on two-dimensional models with rotation symmetry. This approach, which is only justified by large-scale --if at all existent --symmetries of the considered objects, is expected to be strongly misleading or simply wrong on the size scales investigated with MIDI.

In contrast, MATISSE will for the very first time allow image reconstruction of the small-scale regions traced with MIDI and thus finally allow an investigation of these structures based on an unprecedented level of constraints.

Our study underlines the interest for a higher spectral resolution and a large spectral coverage : L, M, and N bands. MATISSE will allow to perform interferometric spectroscopy with three different spectroscopic resolutions in the range of R~30-1000.

In particular, the following spectral resolutions are considered :

• R=30 (amorphous dust features),

• R=100--300 (crystalline dust features, emission features of Polycylic Aromatic Hydrocarbons --PAHs),

• and R=500--1000 (e.g., molecular bands).

These spectral resolutions will allow a fundamental analysis of the composition of gases and dust grains in various environments. For example, MATISSE will be the first instrument ever able to probe the spatial distribution of CO ice in protoplanetary disks and thus provide a key in the understanding of the complex disk chemistry and planet formation.

In most astrophysical domains which require a multi-wavelength approach, MATISSE will be a perfect complement of forthcoming international high angular resolution facilities such as ALMA (the Atacama Large Millimeter Array). MATISSE covers the mid-infrared spectral domain, between the near-infrared domain, for which many interferometric facilities are developed, and ALMA which will operate at (sub) millimeter wavelengths. With the extended wavelength coverage from the M to the N band, MATISSE will not only allow to trace different spatial regions of the sources, but also different physical processes and thus provide insights into previously unexplored areas (such as the investigation of the distribution of volatiles in addition to that of the dust).

The mid-infrared spectral domain is very relevant for the study of the environment of various astrophysical sources. The MATISSE characteristics and expected sensitivity are given in the first part of the presentation. Our science case studies show the wide field of applications of MATISSE. They will be illustrated in the second part of this presentation through the perspective of imaging the circumstellar environments/discs of young stellar objects. In a third part, a general view of the MATISSE modules is presented (the detailed design study being in progress).

INSTRUMENT CHARACTERISTICS

The instrument characteristics derive from the goals defined in the 'Statement of Work and Technical Specifications for the Phase A Study of the MATISSE Instrument' (VLT-SPE-ESO-15860-0001), from the 'MATISSE Phase~A Science Case' study (VLT-TRE-MAT-15860-4325), and, from the feasibility study and choices developed by the 'MATISSE Phase A System Design' (VLT-TRE-MAT-15860-4327).

MATISSE has the following characteristics :

The performance goal in term of sensitivity is the following :

THE ADVANTAGE OF IMPROVING THE U-V PLANE COVERAGE FOR THE PERSPECTIVE OF IMAGING IN THE MID-INFRARED

Observations in the different mid-infrared bands (L, M, N) will allow to trace regions with different characteristic temperatures, allowing to study the physical conditions and chemical environment in different spatial regions, such as in proto-planetary disks or stellar winds. Furthermore, these observations will provide images with different spatial resolution, providing the basis for a comparison with lower-resolution images obtained at large telescopes with adaptive optics --tracing the large-scale structure of the scientific targets --in different wavelength regions (e.g., L/M band: NACO; N/Q band: VISIR).

Beside an individual analysis of measurements in different bands, a combination is thus expected to provide much stronger constraints on models of any of the potential science targets than N band observations alone. Depending on the individual band, unique spectral features (dust/gas) are accessible.

THE SCIENTIFIC POTENTIAL OF THE SPECTROSCOPIC MODES

MATISSE will allow to perform interferometric spectroscopy with three different spectroscopic resolutions. In particular, the following spectral resolutions are considered:

• Spectroscopic resolution R=30 (amorphous dust features),

• R=100--300 (crystalline dust features, emission features of PAHs),

• and R=500--1000 (molecular bands, selected atomic lines). The potential of the different modes, which are focussed on different scientific aspects, are outlined in brief below.

Gas: The 3 to 13 m spectral region is extremely rich in spectral diagnostics of gas covering a huge range in physical and chemical conditions. Strong permitted transitions of HI, HeI, and HeII probe ionized gas covering a wide range in density. The MATISSE spectral region contains a large number of strong ro-vibrational bands of abundant molecules such as CO, OH, H 2 O, SiO, C 2 H 2 , HCN, CS, SiC, that probe molecular gas near young and old stars, as well as in molecular clouds. The gas-phase diagnostics require spectral resolutions of the order of 10 3 for optimal detection sensitivity.

Dust: The MATISSE spectral region contains the strongest resonances of abundant carbon-rich and oxygen-rich dust species. Oxygen-rich dust shows prominent bands at 9.7 µm of amorphous silicates, and several bands of crystalline silicates, such as the forsterite Mg 2 SiO 4 bands at 11.3 and 12.0 µm. In addition, several simple oxides have resonances, such as SiO 2 at 8.6 µm, amorphous Al 2 O 3 at 11.5 µm, and the 13 µm band which is possibly due to Spinel (MgAl 2 O 4 ). Carbon-rich dust has resonances of SiC at 11.3 µm, sometimes attributed to TiC. Furthermore, the famous UIR bands, now commonly attributed to PAHs have their strongest vibrational bands in the MATISSE spectral range (e.g. 3.3, 3.4, 7.7, 8.6, 11.3, 12.7 µm), due to C-H and C-C stretching and bending modes of these aromatic molecules. The 3.52 µm band is believed to be due to nano-diamonds. The spectral resolution required to detect these dust species is of the order of 50 to 100, while studies of the mineralogy of the dust require a spectral resolution of 10 3 . Volatile dust or ice species have resonances somewhat displaced in wavelength compared to the gasphase molecules they form. The most prominent ice band is due to H 2 O and is near 3.1 µm. The shape of the band is sensitive to the ice lattice structure, i.e. whether the ice is crystalline or amorphous. A second water ice band can be found near 12 µm. A second ice band accessible for MATISSE is CO at 4.7 µm. MATISSE will be the first instrument ever to probe the spatial distribution of these ices, providing a key for understanding the complex chemistry of these disks.

SCIENCE CASES FOR MATISSE

The following areas were identified as those where MATISSE will have a major impact by addressing qualitatively new questions (the scientific potential of MATISSE is outlined in the corresponding document, 'MATISSE Phase A Science Case' study (VLT-TRE-MAT-15860-4325, 1 st of June 2007) and in the 'Complement to the Science Case Document' (VLT-TRE-MAT-15860-4336, 15 th of September 2007) :

• Planet formation,

• Star formation (low and massive stars, multiple stars),

• Circumstellar disk evolution, ,

• Active Galactic Nuclei,

• Extrasolar planets,

• Solar System Minor Bodies,

• Dust and winds from evolved stars,

To illustrate the potential of MATISSE, we outline the example of planet formation studies in the following discussion.

Example: Planet Formation

Current star formation models predict the existence of circumstellar disks around young stellar objects as a natural consequence of the star formation process.

Given the importance and the strong interest in the astronomical community to study circumstellar disks in order to reveal the origin of planetary systems such as our own, circumstellar disks are key subject of observations by all available means. Unfortunately, the innermost region of disks with a radius of a few AU, where planet formation is expected to take place, can only be marginally investigated so far by the study of the spectral energy distribution, single visibility measurements (e.g. with MIDI or AMBER), high-resolution images obtained with the Hubble Space Telescope or by using adaptive optics. This will remain the case until ALMA will be in operation with its longest baselines (~ 2012), unless one can provide high angular resolution images in the mid-infrared wavelength range.

In contrast to the typically ambiguous interpretations of visibilities as obtained with MIDI, MATISSE is aimed at answering pivotal questions about the inner regions of young circumstellar disks, such as:

[1] What is the surface brightness profile in circumstellar disks around T Tauri HAeBe stars ?

MATISSE will allow to consider that circumstellar disks which are not seen exactly face-on do show a brightness profile which cannot be described by iso-brightness contours, which are centered on the location of the central star, as assumed in the previous data analysis. It will permit to derive the radial temperature profile of the hot dust on the disk surface and the inner disk rim. This profile will provide information about the radial and vertical structure of the disk, and thus on the physical conditions in the potentially planet-forming regions of circumstellar disks.

[2] Does the brightness profile show evidence for dust grain growth and sedimentation ?

MATISSE will allow to provide strong observational constraints for the region, where dust grain growth is presumably fastest.

As simulations of dust settling in circumstellar disks show, the disk flaring and thus the ability to absorb stellar radiation even at large distances from the star depends on the grain size distribution remaining in the upper disk layers --an effect that can also be investigated with MATISSE (Dullemond and Dominik 2004).

[3] Is the inner disk structure modified by early stages of planet formation ?

The inner region of circumstellar disks is expected to show large-scale (sub-AU -AU sized) density fluctuations / inhomogenities. Locally increased densities are expected to show up as local brightness variation (due to increased absorption/ shadowing effects) in the mid-infrared images (see Figure 1).

Fig. 1: Simulated 10 µm images of the inner region of a circumstellar T Tauri disk, with an embedded Jupiter-mass planet at a distance of 5.2 AU from the central star (Wolf and Klahr, in prep.). The hot region around the planet above the center of the disk is indicated as bright areas in these reemission images. Assuming a distance of 140 pc, a 20 mas scale is indicated.

[4] Is there indirect or even direct evidence for the presence of planets ?

Once planets have been formed, they may significantly alter the surface density profile of the disk and thus cause signatures in the disk that are much easier to find than the planets themselves.

The appearance and type of these signatures depend on the mass and orbit of the planet, but even more on the evolutionary stage of the circumstellar disk. MATISSE would for example allow to observe a gap, which is expected to be created along the orbit of a young giant planet within the inner few AU of a young gas-dominated disks [START_REF] Wolf | MATISSE Phase A Science Case' study[END_REF]).

A BRIEF VIEW OF THE DESIGN

MATISSE is composed of several sub-systems. The first one is called "Warm optics". It represents the instrumental part between the VLTI feeding optics and the cryostat windows. Located on an optical table of 1.5x2.1 m size, it is composed by 5 modules : the field rotation compensator, the beam commuting device, the spectral separator, the anamorphic optics and the delay lines.

The field rotation compensator optimizes the field orientation in regards to the slits of the spectrograph and compensates the field and the polarization rotation existing in the focal laboratory.

The beam commutating device is used as a calibration tool for instrumental effects (including chromatic and detector effect) by inverting them thanks to beam commutations. A spectral separation between the L (and M) band and the N band is necessary because the beam combination modes, the materials of optics and the detectors are different. This separation is made by a set of 4 dichroïcs.

The anamorphic optics is used to make an image magnification in the spatial direction in order to optimize the fringe sampling.

The delay lines (one set per spectral band) transfer the pupil at the cold pupil mask level of the spectrograph and are used for fine OPD adjustment.

The cold optics represents the MATISSE part between the cryostat windows and the detector. It is made up of 2 separate cryostats (one for each spectral band). The functionalities of each cryostat are similar. They are composed by several elements:

The shutters select the chosen beams (permitting observations with 2, 3 or 4 telescopes).

The pupil masks reduce the thermal background (in particular the thermal background emitted by the optical mounts outside the cryostat) and avoid parasitic light.

Focusing mirrors produce an intermediate image to insert pinholes or slits used as spatial filter. After the spatial filter, collimating mirrors collimates the beams. The pair separator composed by beam splitters and mirrors module separates the 4 beams in 6 pairs of beam. By combining different positions of shutters and beam splitter, several configurations of observation can be obtained: 4 telescopes without specific photometric images, 3 telescopes with specific photometric images and 2 telescopes with or without photometric images. The π phase shifter separates in two parts each pair of beam and produces a π phase between these two outputs.

In case of observations with 4 telescopes, 12 images are produced (2 images per pair and 6 pairs of beam). 12 tilted flat mirrors introduce a slight angle between the beams in order to position all these images at different locations on the detector.

Prism and grisms can be chosen to obtain the desired spectral resolution. The camera optics focalize the beam on the detectors and is the beam combiner (multi-axial combination) of the instrument.

The Figure 2 gives the general implementation of the MATISSE instrument. The warm optical table in the centre of the figure has a size of about 1.2 x 2.4 m (TBC). The two cryostats are located symmetrically on the two sides of the warm optical bench. Cryostat L&M Warm Optics

VLTI

The Figure 3 is a 3D view of the cryostat. On the left side behind the "grey wall" we can identify the four input windows. Then the beams are deviated by four mirrors towards four couples of parabolic mirrors (image plane, spatial filtering). After the spatial filtering module, in the centre of the figure we can identify the pair creator module followed by the OPD module (centre -right). At the bottom right we can see the 0 -π  module followed by the tip/tilt adjustment mirrors for re-positioning the pairs onto the detector. Several flat mirrors are placed before the camera optics and the detector.

CONCLUSION

There exist only two other 10 µm ground based long-baseline interferometers (ISI 3 and Keck I) and two projects presently under study for mid-infrared space interferometry (Darwin and TPF I).

The Infrared Stellar Interferometer ISI 3, equipped with three 1.6 m telescopes, allows imaging and operates in heterodyne mode at a moderate sensitivity level, related to the very narrow band used.

The Keck Interferometer, on the other hand, has two large apertures on fixed baselines. Its sensitivity is equivalent to that of the VLTI. However, its imaging capability is only based on the supersynthesis effect which allows the exploration of a limited range of spatial frequencies in the Fourier plane.

Darwin and TPF I (Terrestrial Planet Finder Interferometer) are space projects, particularly focused on detecting earthlike extrasolar planets and the search for bio-signatures in their atmospheres.

MATISSE is a unique instrument (see Figure 4). It offers new scientific breakthroughs and will contribute to the evolution of mid-infrared interferometry.

Fig. 2 :

 2 Fig. 2: General top view of the optical design of MATISSE

Fig. 4 :

 4 Fig. 4: Uniqueness of MATISSE at theVLTI represented in a Spatial Resolution versus Wavelength diagram.


Laurent Marichal

Fabrice Cousin

Gaël Giraudon-Colas

Serge Pin

Jean Philippe Renault

Jean Philippe

Assessing the in situ protein corona structure

Assessing the in situ protein corona structure Laurent Marichal 1 , Fabrice Cousin 2 , Gaël Giraudon--Colas 1 , Serge Pin 1 , and Jean Philippe Renault 1 1 NIMBE, CEA Saclay, CNRS, Université Paris-Saclay, Gif-sur-Yvette, France 2 Laboratoire León-Brillouin, CEA Saclay, CNRS, UniversitéParis-Saclay, Gif-sur-Yvette, France

As soon as nanoparticles are placed in a biological medium, dynamic interactions happen between their surfaces and neighbouring biomolecules. Proteins are particularly prone to interactions and can form the so-called protein corona around these particles. This gives a biological identity to nanoparticles that will condition their biodistribution and potential toxicity. We thoroughly studied the interactions happening between model proteins (myoglobin and hemoglobin) and nanoparticles (silica). We combined a series of physico-chemical techniques such as oxygenation studies, isothermal titration calorimetry, circular dichroism, and small-angle neutron scattering (SANS). First, we found out that hemoglobin have a stronger affinity to silica surface compared to myoglobin, partly because of its larger size [1]. Secondly, adsorbed and free hemoglobin have different oxygenation properties, especially in terms of oxygen affinity (increase of oxygen affinity for adsorbed hemoglobin) [2]. The structure of adsorbed proteins was then assessed by circular dichroism and SANS. Very subtle secondary and tertiary structure changes could be found while the quaternary structure was totally preserved. This indicates that even without altering their structure, proteins adsorbed on nanoparticles can have their function drastically altered. Finally, SANS studies allowed us to probe the complex structures of the protein corona (Fig. 1) [2]. Myoglobin and hemoglobin form a monolayer of proteins around nanoparticle. We could also determine that, in every system studied, protein/nanoparticle interactions lead to a reaction-limited aggregation process. Besides, it was found that even incomplete coronas were organized by short-distance repulsive interactions between adsorbed proteins. For standard proteins such as the ones studied here the term "protein corona" is particularly well-suited even though, this is not a universal phenomenon since very large proteins tend to form open lattices when interacting with nanoparticles [1].

Figure 1 .

 1 Figure 1. Schematic representation of the protein/nanoparticle system studied by small-angle neutron scattering. Solutions can be in H 2 O (left) or D 2 O (right) in order to contrast match the nanoparticles to the solvent.


Victor Jung
email: victor.jung@univ-cotedazur.fr

Jean-Charles Régin
email: jean-charles.regin@univ-cotedazur.fr

Checking Constraint Satisfaction

Keywords: Multi-valued Decision Diagram, Inclusion, Constraint Learning

We address the problem of verifying a constraint by a set of solutions S. This problem is present in almost all systems aiming at learning or acquiring constraints or constraint parameters. We propose an original approach based on MDDs. Indeed, the set of solutions can be represented by the MDD denoted by M DDS . Checking whether S satisfies a given constraint C can be done using M DD(C), the MDD that contains the set of solutions of C, and by searching if the intersection between M DD(S) and M DD(C) is equal to M DD(S). This step is equivalent to searching whether M DD(S) is included in M DD(C). Thus, we give an inclusion algorithm to speed up these calculations. Then, we generalize this approach for the computation of global constraint parameters satisfying C. Next, we introduce the notion of properties on the MDD nodes and define a new algorithm allowing to compute in only one step the set of parameters we are looking for. Finally, we present experimental results showing the interest of our approach.

Introduction

Many works in Constraint Programming try to improve the quality of a model by adding new implicit constraints [START_REF] Picard-Cantin | Learning parameters for the sequence constraint from solutions[END_REF], redundant constraints [START_REF] Bessière | Learning implied global constraints[END_REF] or global constraints [START_REF] Leo | Globalizing constraint models[END_REF]. All these works face a common problem: the verification of the satisfaction of constraints by a given set of solutions. Some choose a brute force approach [START_REF] Leo | Globalizing constraint models[END_REF][START_REF] Picard-Cantin | Learning parameters for the sequence constraint from solutions[END_REF], others prefer a more specific but ad-hoc approach [START_REF] Beldiceanu | A constraint seeker: Finding and ranking global constraints from examples[END_REF]. In all cases, these methods go through the solutions to test if they satisfy constraints. Constraints are not necessarily tested individually, but the solutions can be considered one after the other.

In this paper, we propose a more global and efficient method to test whether a set of solutions verifies one or a set of constraints. Multi-valued decision diagrams (MDDs) are a very efficient data structure to represent a set of solutions in a compressed way and for which many operators are available to combine MDDs without decompressing them. We therefore propose to use M DD(S) the MDD which corresponds to the set of solutions S. We show that we can simply test if S satisfies a constraint C, by using M DD(C), the MDD that represents the solutions of C, and then by performing the intersection between M DD(S) and M DD(C). This method is simple to implement since it only requires constructing the two MDDs and performing their intersection, for which efficient algorithms are available, and then testing whether the resulting MDD is similar to M DD(S). However it has an important flaw: it will create and calculate a MDD even if the intersection will not be equal to M DD(S). It therefore risks doing many operations unnecessarily. To avoid this we introduce an inclusion operator between MDDs since this is what we want to test: is M DD(S) included in M DD(C)?

This operator is efficient when it is a question of verifying a precise and unique constraint, but not very efficient when it is a question of searching for the parameters of a constraint such that the resulting constraint is satisfied by a set of solutions. Finding the parameters of a global constraint so that it is satisfied by a set of solutions is a recurrent problem at present ( [START_REF] Picard-Cantin | Learning parameters for the sequence constraint from solutions[END_REF], [START_REF] Picard-Cantin | Learning the parameters of global constraints using branch-and-bound[END_REF], [START_REF] Leo | Globalizing constraint models[END_REF], [START_REF] Beldiceanu | A constraint seeker: Finding and ranking global constraints from examples[END_REF], [START_REF] Bessière | Learning implied global constraints[END_REF]). To solve this problem we propose to work with M DD(S) which we enrich by introducing the notion of node properties. Then, a process called "the parent-child propagation of the parameters" is performed through the MDD. More precisely, the global constraint is expressed by properties including the parameters and these properties are propagated in M DD(S) in order to compute for each sub-tree of the MDD those which are compatible with the constraint. Thus, we determine the most restrictive parameters of the constraint that are satisfied by the MDD.

This article is organized as follows. First, we give some basic definitions. Then, we present a general scheme to check the satisfaction of a constraint and improve it by defining a new operation between MDDs. Next, we address the problem of finding parameters of global constraints by introducing the notion of node properties. Finally, we provide benchmarks and results testing the different approaches described in this article, and we conclude.

Preliminaries

Constraint Programming

A finite constraint network N . is defined as a set of n variables X = {x 1 , . . . , x n }, a set of current domains D = {D(x 1 ), . . . , D(x n )} where D(x i ) is the finite set of possible values for variable x i , and a set C of constraints between variables. We introduce the particular notation D 0 = {D 0 (x 1 ), . . . , D 0 (x n )} to represent the set of initial domains of N on which constraint definitions were stated. A constraint C on the ordered set of variables

X(C) = (x i1 , . . . , x ir ) is a subset T (C) of the Cartesian product D 0 (x i1 ) × • • • × D 0 (x ir ) that specifies the allowed combinations of values for the variables x i1 , . . . , x ir . An element of D 0 (x i1 )ו • •× D 0 (x ir ) is called a tuple on X(C). A value a for a variable x is often denoted by (x, a). Let C be a constraint. A tuple τ on X(C) is valid if ∀(x, a) ∈ τ, a ∈ D(x). C is consistent iff there exists a tuple τ of T (C) which is valid. A value a ∈ D(x) is consistent with C iff x ∈ X(C)
or there exists a valid tuple τ of T (C) with (x, a) ∈ τ . We denote by #(a, τ ) the number of occurences of the value a in a tuple τ .

We present some constraints that we will use in the rest of this paper. A global cardinality constraint (gcc) constrains the number of times every value can be taken by a set of variables. This is certainly one of the most useful constraints in practice. Note that the alldiff constraint corresponds to a gcc in which every value can be taken at most once.

Definition 1 A global cardinality constraint is a constraint C in which each value a i ∈ D(X(C)) is associated with two positive integers l i and u i with l i ≤ u i defined by gcc(X, l, u) = {τ |τ is a tuple on X(C) and ∀a i ∈ D(X(C)) : l i ≤ #(a i , τ ) ≤ u i } Definition 2
Given X a set of variables, l and u two integers with l ≤ u and V a set of values. The among constraint ensures that at least l variables of X and at most u will take a value in V , that is among(X, V, l, u) = {τ | τ is a tuple on X(C) and l ≤ a∈V #(a, τ ) ≤ u} This constraint has been introduced in CHIP [START_REF] Beldiceanu | Introducing global constraints in CHIP[END_REF]. The sequence constraint [START_REF] Beldiceanu | Introducing global constraints in CHIP[END_REF] is a conjunction of sliding among constraints.

Definition 3 Given X a set of variables, q, l and u three integers with l ≤ u and V a set of values. The sequence constraint holds if and only if for 1 ≤ i ≤ n q +1 among({x i , ..., x i+q-1 }, V, l, u) holds. More precisely Sequence(X, V, q, l, u) = { τ | τ is a tuple on X(C) and for each sequence S of q consecutive variables: l ≤ v∈V #(v, τ, S) ≤ u}

Multi-valued Decision Diagram

The decision diagrams considered in this paper are reduced, ordered multi-valued decision diagrams (MDD) [START_REF] Kam | Multi-valued decision diagrams[END_REF][START_REF] Srinivasan | Algorithms for discrete function manipulation[END_REF][START_REF] Bergman | Decision Diagrams for Optimization[END_REF], which are a generalization of binary decision diagrams [START_REF] Bryant | Graph-based algorithms for boolean function manipulation[END_REF]. They use a fixed variable ordering for canonical representation and shared sub-graphs for compression obtained by means of a reduction operation. An MDD is a rooted directed acyclic graph (DAG) used to represent some multivalued functions f : {0...d -1} n → true, f alse. Given the n input variables, the DAG contains n + 1 layers of nodes, such that each variable is represented at a specific layer of the graph. Each node on a given layer has at most d outgoing arcs to nodes in the next layer of the graph. Each arc is labeled by its corresponding integer. The arc (u, a, v) is from node u to node v and labeled by a. Sometimes it is convenient to say that v is a child of u. All outgoing arcs of the layer n reach tt, the true terminal node (the false terminal node is typically omitted).

There is an equivalence between f (a 1 , ..., a n ) = true and the existence of a path from the root node to the tt whose arcs are labeled a 1 , ..., a n . The reduction of an MDD is an important operation that may reduce the MDD size by an exponential factor. It consists in removing nodes that have no successor and merging equivalent nodes, i.e. nodes having the same set of neighbors associated with the same labels. This means that only nodes of the same layer can be merged.

MDD of a constraint. Let C be a constraint defined on X(C). The MDD associated with C, denoted by M DD(C), is an MDD modeling the set of tuples of C. More precisely, M DD(C) is defined on X(C), such that the arc labels of the layer of the variable x correspond to values of x, and a path of M DD(C) (that is a path from the root node to the tt node) where a i is the label of layer i corresponds to a tuple (a 1 , ..., a n ) on X(C).

Operators We reproduce here the description of the generic Function apply [START_REF] Perez | Efficient operations on MDDs for building constraint programming models[END_REF][START_REF] Perez | Decision Diagrams : Constraints and Algorithms[END_REF] because we will see that the inclusion can be easily modelled thanks to it1 . From the MDDs mdd 1 and mdd 2 it computes mdd r = mdd 1 ⊕ mdd 2 , where ⊕ is union, intersection, difference, symmetric difference, complementary of union and complementary of intersection. Function apply is mainly based on the possible combinations of labeled arcs. It proceeds by associating nodes of the two MDDs operands. Each node x of the resulting MDD is associated with a node x 1 of the first MDD and a node x 2 of the second MDD represented by a pair (x 1 , x 2 ). First, the root is created from the two roots. Then, the layers are successively built. From the nodes of layer i -1 the nodes of layer i are built as follows.

For each node x = (x 1 , x 2 ) of layer i -1, the arcs outgoing from nodes x 1 and x 2 and labeled by the same value v are considered. We recall that there is only one arc leaving a node x with a given label. Thus, there are four possibilities depending on whether there are y 1 and y 2 such that (x 1 , v, y 1 ) and (x 2 , v, y 2 ) exist or not. The action that is performed for each of these possibilities will define the operation performed for the given layer. For instance, a union is defined by creating a node y = (y 1 , y 2 ) and an arc (x, v, y) each time one of the arcs (x 1 , v, y 1 ) or (x 2 , v, y 2 ) exists. An intersection is defined by creating a node y = (y 1 , y 2 ) and an arc (x, v, y) when both arcs (x 1 , v, y 1 ) and (x 2 , v, y 2 ) exist. Thus, these operations can be simply defined by expressing the condition for creating a node and an arc.

Function apply, given in Algorithm 1 takes as parameters the two MDDs, two arrays op, V having as many elements as layers, and typeOp the operation type (i.e. intersection, union...). For each layer i, op[i] contains 4 entries, each one representing the fact that we create an arc or not for a combination of arc existence in the two MDDs and V [i] represents the set of values needed by the complementary set. If it is equal to nil then V [i] will be equal to the union of the values of the neighbors of the considered nodes. At the end the resulting MDD is reduced by calling pReduce algorithm [START_REF] Perez | Efficient operations on MDDs for building constraint programming models[END_REF]. The values of op[i] defining the binary operations are defined as follows for the different combinations: If this intersection is equal to M DD(S) then this means that all solutions of S satisfy the constraint C. The proof of the soundness of this approach is quite immediate : an MDD is a set of solutions, so if the intersection does not modify M DD(S) then it means that any solution of M DD(S) is also a solution of M DD(C) and therefore this solution satisfies the constraint. This approach is not particularly efficient, because it systematically requires the intermediate calculation of an intersection between MDDs. However, we are not interested in this intersection 2 . What matters is to know whether this intersection is similar to the initial MDD. We can reduce what we are trying to do in a single step: we check a relation of inclusion. Indeed, answering the question M DD(S) ∩ M DD(C) = M DD(S)? is equivalent to answering the following question: M DD(S) ⊆ M DD(C)? As this operator does not exist in the literature, we propose to create it.

op[0] op[1] op[2] op[3] ¬a1 ∧ ¬a2 ¬a1 ∧ a2 a1 ∧ ¬a2 a1 ∧ a2 layer [1..r-1] r [1..r-1] r [1..r-1] r [1..r-1] r A ∩ B F F F F F F T T A ∪ B F F T T T T T T A -B F F F F T T T F Algorithm 1 Generic Apply Function. apply(mdd1, mdd2, op, V, typeOp): (MDD, bool) // L[i] is the set of nodes in layer i. root ← createNode(root(mdd1), root(mdd2)) L[0] ← {root} for each i ∈ 1..r do L[i] ← ∅ for each node x ∈ L[i -1] do get x1 and x2 from x = (x1, x2) if V [i] = nil then V [i] ← values(ω + (x1) ∪ ω + (x2)) for each v ∈ V [i] do if ∃(x1, v, y1) ∈ ω + (x1) then if ∃(x2, v, y2) ∈ ω + (x2) ∧ op[0] then createArc(L, i, x, v, w[i]) if ∃(x2, v, y2) ∈ ω + (x2) ∧ op[1] then addArcAndNode(L, i, x, v, nil, y2) else if ∃(x2, v, y2) ∈ ω + (x2) ∧ op[2] then addArcAndNode(L, i, x, v, y1, nil) if ∃(x2, v, y2) ∈ ω + (x2) ∧ op[3] then addArcAndNode(L, i, x, v, y1, y2) if typeOp = Inclusion then if ∃(x1, v, y1) ∈ ω + (x1)∧ ∃(x2, v, y2) ∈ ω + (x2) then return (nil, f alse) if typeOp = Inclusion then return (nil, true) merge all nodes of L[r] into t pReduce(L) return (root, true) addArcAndNode(L, i, x, y1, v, y2) if ∃y ∈ L[i] s.t. y = (y1, y2) then y ← createNode(y1, y2) add y to L[i] createArc(L, i, x, v, y) 3

Operator of Inclusion

The inclusion operator between MDDs is easily done using the generic function apply. Let's consider that we want to know if M DD 1 is included in M DD 2 . We use the same rules as the intersection operator with a notable exception: if an edge a is in M DD 1 but not in M DD 2 , then we end the algorithm by returning false. In this case there is at least one solution in M DD 1 which is not in M DD 2 so M DD 1 cannot be included in M DD 2 . For the three other cases we can easily find those of the intersection. Using the terminology of the preliminaries, we have clearly: ¬a 1 implies that no arc is created and a 1 ∧ a 2 implies an arc creation, since the solutions are common.

We notice that it is not necessary to keep in memory the MDD that is built. Indeed, we just need to know if we can create each new level. To do this, only the last level that has just been built is useful and must be kept in memory, the others being no longer useful can be destroyed. The reduction of the built MDD is no longer useful either since we are only interested in the ability to build an MDD from the root to tt. Function apply must therefore return true instead of performing the reduction at the end. This allows us to save time compared to the previous method.

Inferring parameters of global constraints

The inclusion operator allows to answer in an efficient way to the question of the satisfaction of a constraint by a set of solutions S. However, in practice, the question that is often asked is more general: given a global constraint C involving a set of parameters, for which parameters S satisfies this constraint?

Let us consider P a set of parameters and C(P ) a constraint defined using these parameters. Formally, we can present the problem in the following way: What are the sets P such that ∀s ∈ S, s satisfies C(P )?

A first way to proceed is to check for each set of parameters P if we have M DD S ⊆ M DD C(P ) .

We propose to add additional information, called properties, to each node of an MDD. This idea has similarities with the scheme introduced by J. Hooker et al. [START_REF] Hoda | A systematic approach to MDD-based constraint programming[END_REF]. This information is used to memorize the valid parameters from the root to the node in relation to the constraint under consideration. When we will reach tt, we will know the parameters that are checked by all the solutions. In addition, retaining all the parameter sets is superfluous and we can be satisfied with retaining the more restrictive parameters. Other acceptable parameters may be derived from these restrictive parameters. For example, if the constraint sequence(X, V, q = 3, l = 1, u = 2) is satisfied then the constraints sequence(X, V, q = 3, l = 0, u = 2), sequence(X, V, q = 3, l = 1, u = 3, ) and sequence(X, V, q = 3, l = 0, u = 3) are also satisfied. So, the more restrictive parameters are (q = 3, l = 1, u = 2). We present on an example the ideas of our algorithm. We will use the binary representation of a sequence constraint. Indeed, for a sequence constraint, we can abstract any X and V into a binary problem with V = {1}. If x i = 1 then it means that x i takes its value in V (x i ∈ V ) otherwise we have x i =0. So we are in the presence of only binary variable and we are looking for the parameter values (q, u, l) which are satisfied by S. Consider the red node f in M DD(S) (Fig. 1). There are two paths to access this node from the root: take 0, 0 then 1 (which corresponds to the blue path -left-on the diagram) or take 1, 1 then 0 (which corresponds to the path in green -right-). In total, there are 5 sub-MDDs (i.e. smaller MDDs contained in the main MDD) having this red node f as the terminal node: 1 having as starting point the general root ({0-0-1, 1-1-0}), then 2 having as root the two nodes of the first layer ({0-1}, {1-0}) and finally 2 having as root the two nodes of the second layer ({1}, {0}). Strictly speaking, there are actually 6 sub-MDDs, since the MDD consisting only of the red node f exists.

Properties (i.e. satisfied sequences) are added to nodes:

-Node R. This node contains only the basic information, i.e. (q = 0, l = 0, u = 0), since the only way to reach this node is to start from it and take no edge. This is the basic case, viable for all nodes. -Node a (blue). In addition to the basic case, it is possible to reach this node starting from the root and taking the value 0. As we take 0 times the value 1, the satisfied Sequence is (q = 1, l = 0, u = 0).

-Node b (green). Same as for the previous node, except that we take once the value 1. The satisfied Sequence is therefore (q = 1, l = 1, u = 1). -Node c (blue). We retrieve the information from the parent node (there is only one here). So we have (q = 0, l = 0, u = 0) and (q = 1, l = 0, u = 0). There is only one edge that can be traversed, with a value of 0. If we add 0 to the preceding satisfied sequences, we obtain (q = 1, l = 0, u = 0) and (q = 2, l = 0, u = 0). We retain these sequences. -Node d (green). By the same reasoning, we obtain (q = 1, l = 1, u = 1) and (q = 2, l = 2, u = 2). -Node f (red). We start by looking at parent Node c (blue). By adding the fact that we can reach the red node f by taking the value 1, we have (q = 1, l = 0, u = 1), (q = 2, l = 0, u = 1) and (q = 3, l = 1, u = 1).

In the same way, looking at the side of parent Node d (green), we obtain (q = 1, l = 0, u = 1), (q = 2, l = 1, u = 2) and (q = 3, l = 2, u = 2). We notice that Sequence constraints of size 2 and size 3 are not compatible. In this case, the union of the two Sequence constraints is performed (since both are satisfied). We thus obtain (q = 1, l = 0, u = 1), (q = 2, l = 0, u = 2) and (q = 3, l = 1, u = 2). We can then check that for each path leading to the red node f , we take between 0 and 2 times the value 1 for a path of size 2, and between 1 and 2 times the value 1 for a path of size 3.

Fig. 2 shows a slightly more complete example.

Fig. 2: Satisfied sequence constraints for each node. Value q corresponds to the index in the array associated with a node.

However, one thing is noticeable: it is possible to lose information (Fig. 3). For example, the node b contains the information (q = 0, l = 0, u = 0), (q = 1, l = 0, u = 0) and (q = 2, l = 1, u = 1), but we can see that for paths of size 1 belonging to the MDD it is possible to take between 0 and 1 times the value 1. To solve this problem, we make the union of all the nodes of a layer (for each layer) that we store in an accumulator. This operation can be performed at the same time as the information is constructed. Fig. 3: Simple example to show information loss. Value q correspond to the index in the array associated with a node.

We remind that it is not necessary to retain all the information as represented on the diagram: we can simply keep the last two layers since the construction is done in a sliding way.

Implementation

The information is represented in the form of a property associated with each node. Each constraint will define its own property. The propagation of information from the MDD nodes is performed using a breadth-first approach, because we need all parent nodes to be correctly defined before propagating to the children. A depth-first approach would be strictly speaking impossible because all children would have to be re-explored each time the parents are updated (which is highly inefficient). During propagation, one looks to see if the child is already carrying a property or not. Two cases are possible.

-The child's property is already defined. In this case, the information already present must be merged with the new information provided by the parent. Function mergeWithProperty of a property is in charge of this.

-The child property is not defined. In this case we simply create new information based only on the parents' information. Function createProperty of the property performs this operation.

Algorithm 2 is a possible implementation of this mechanism. The important Functions are mergeWithProperty and createProperty. They depend on the type of constraint that the property represents, so it is difficult to define a general way to represent the information. Technically, Function createProperty is quite simple. We believe that the real difficulty lies in the definition of Function mergeWithProperty, because the information must be complete and valid, i.e. it must represent the state of the constraint in a correct way for the node that contains it at any moment of the propagation. Space complexity. We retain (i+1) information for each node of layer i. The space complexity thus depends both on the layer where we are located (the last two layers in reality) and on the number of nodes present in the layer. As it is not possible to predict the number of nodes in a layer for any MDD, the complexity remains rather vague. Let L i be the number of nodes in the layer i and L the number of layers. The space complexity is in: O(Max(i

× |L i | + (i + 1) × |L i+1 |)), ∀i s.t. 0 ≤ i < L

Properties definitions

We provide the definition of properties for sum ( x∈X x = [min, max]), sequence, and cardinality constraints as examples. These constraints have been chosen in relation to the problems we are working on, and not in relation to the difficulty of implementation. These are also quite common constraints.

Sum Constraint. The sum property is quite simple to model (Algorithm 3). As the result is an interval, we just need to use a pair (min, max) to represent the data. Switching from a parent to a child is simply the addition of an interval with an integer (the value of the arc), and the merge operation is a simple union between two intervals. Each node contains the minimum and maximum value that can be obtained by a path from the root to that node. The end node tt therefore contains the minimum and maximum value that can be obtained by taking any path through the MDD.

Algorithm 3 Sum Property createProperty(label)

(min, max) ← (thismin + label, thismax + label) return (min, max) mergeWithProperty(property, label)
thismin ← min(thismin, propertymin+ label) thismax ← max(thismax, propertymax+ label)

getResult() return this

Cardinality Constraint We represent the property of a global cardinality constraint as a values matrix of size |V | × 2. Each value considered in the constraint is associated with a pair (min, max) representing the minimum and maximum number of times the value is taken. Moving from a parent to a child through an arc labeled by a, amounts to incrementing by 1 the number of times the value a is taken, which is a simple addition operation on the intervals. On the other hand, performing the merge amounts, as for the sum, to performing a union operation. Algorithm 4 is a possible implementation.

Sequence Constraint. The sequence property is the most complex of the three (Algorithm 5). In itself, Function createProperty and mergeWithProperty correspond, as for the cardinality constraint, to perform respectively an addition on intervals and a union operation. The difference comes from the fact that, for the sequence, it is necessary to have separate information, represented here by accumulator. This peculiarity comes from the sliding and global aspect of the sequence constraint: each node contains information about a local sequence, i.e. sequences that contain it as the final value. However, it is possible in this case to lose information -as shown in Fig. 3. This accumulator can be implemented in different ways: either by building it on the fly (memory gain but time loss), in which case it is not necessary to retain the information on more than two layers (the last two in a sliding manner), or by building it once the propagation of the property is completed, but all the information of all the nodes must be in memory (time gain but memory loss). Here, the accumulator is built on the fly.

Experiments

Testing environment

The algorithms have been implemented in Java 12. The experiments were performed on a Windows 10 machine using a Ryzen 2600 AMD CPU and 32 GB of RAM for the car sequencing problem and on a machine having four E7-4870 Intel processors, each having 10 cores with 256 GB of memory and running under Scientific Linux for the nurse rostering problem. The different tests comparing the methods presented in this paper were performed using solutions from instances of Car Sequencing and Nurse Rostering (represented as a MDD). The details of this instances can be found in appendix. However, we give some important information:

Car Sequencing The Car Sequencing MDD contains 25942 nodes, 53985 arcs and represents 2.6 × 10 14 solutions. There are two options with capacity 1/2 and 2/3, 4 car classes (configurations) and a total of 100 cars.

Nurse Rostering

The Nurse Rostering MDD contains 128325 nodes, 220600 arcs and represents 1.2 × 10 28 solutions. There are 6 nurses, 28 days and 3 shifts. The scheduling has some predefined data : 1 means that the nurse is working that day, -means that the nurse is not working that day and 0 means that it is yet to be determined. Each day, each shift has a minimum number required of nurse working that shift (we can have more, but never less). We have various constraints, such that a nurse cannot work more than 7 days straight and must have at least 2 free days in a row in a 2-weeks window.

Testing each raw solution individually requires testing 2.6×10 14 and 1.2×10 28 solutions (respectively for Car Sequencing and Nurse Rostering problems). Doing such benchmark in a reasonable amount of time is out of the question.

Comparison between inclusion and intersection based inclusion

We will compare the two methods to compute the inclusion presented in this paper. We do not take into account the time and space needed to create and store the MDDs representing the constraints -we only focus on time and space needed to compute the inclusion between two MDDs.

We can see from results in Table 1 that the inclusion method is better in every way than the intersection based inclusion. We observe improvements by a compute all sequences with the other methods. Furthermore, a very big part of the process is to build all the MDDs of the constraints, resulting in an increase in both time and memory consumption, as expected. Once again we find the factor 2 that we had in our previous comparison between the two methods of inclusion.

Conclusion

This article sheds light on a new aspect of the interest of using MDDs in the context of constraint programming. We have introduced a new inclusion operator that allows to answer the question of the satisfaction of a constraint more efficiently than by using the classical sequence of operations on MDDs (intersection, reduction, comparison). In addition, we have shown that adding properties to the nodes of a MDD allowing to represent locally the state of a constraint is a very efficient way to obtain the parameters of a global constraint (we presented GCC, sum and sequence), provided that this constraint can be formalized as a node property. This method is much more advantageous than a succession of inclusion operations, both from a temporal and spatial point of view, because it does not require any constraint construction in the form of an MDD. The use of inclusion is nevertheless of interest when the constraint is particularly complex, unique, or very difficult to formalize in the form of a property. The question of the use of properties for other constraints (other than global) seems to be the next step in order to answer in more detail the problem of extracting the parameters of a constraint from a set of solutions.

Acknowledgments

This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

Appendix

Car Sequencing The Car Sequencing MDD contains 25942 nodes, 53985 arcs and represents 2.6 × 10 14 solutions. There are two options with capacity 1/2 and 2/3, 4 car classes (configurations) and a total of 100 cars, defined as follow : Nurse Rostering The Nurse Rostering MDD contains 128325 nodes, 220600 arcs and represents 1.2 × 10 28 solutions. There are 6 nurses, 28 days and 3 shifts.

The scheduling has some predefined data : 1 means that the nurse is working that day, -means that the nurse is not working that day and 0 means that it is yet to be determined. Each day, each shift has a minimum number required of nurse working that shift (we can have more, but never less). We have various constraints, such that a nurse can't work more than 7 days straight and must have at least 2 free days in a row in a 2-weeks window.

  Checking constraint satisfactionA first solution to test whether a set of solutions S satisfies a constraint C is to represent S by an MDD, denoted M DD(S), then use M DD(C) the MDD of the constraint C and calculate the intersection between M DD(S) and M DD(C).

Fig. 1 :

 1 Fig. 1: Sequences for M DD(S)

  As these functions depend on the constraint and the implementation, we can't give any further details.

	Algorithm 2 Propagation
	propagatePropertyMDD, property

MDD.root.addP roperty(property) for each layer L in MDD do for each node in L do transferProperty(node) return MDD.tt.getP roperty().getResult() transferProperty(node) for each (label, child) in node.children do if child.hasP roperty() then child.getP roperty().mergeWithProperty(node.getP roperty(), label) else child.addP roperty(node.getP roperty().createProperty(label)) Time complexity: in O((|A| -|N |)× O(mergeWithProperty) + |N |× O(createProperty)), where |A| is the number of edges and |N | the number of nodes. Function createProperty is called only once for each node, that is |N | times globally. Function mergeWithProperty is globally called |A| -|N | times. O(mergeWithProperty) and O(createProperty) are the time complexity of the functions mergeWithProperty and createProperty.

Table 2 :

 2 Car Sequencing Problem (Time in ms, memory in MB).

			GCC	Sum		Sequence
	Methods Time Memory Time Memory Time Memory
	Intersection 6 080 3 240 2 137 1 328 11 306 6 493
	Inclusion 3 371 2 204 1 391 1 021 6 114 3 867
	Properties	48	10	22	4	82	37
			GCC		Sequence
	Methods		Time Memory Time Memory
	Intersection 663 647 170 281 704 405 137 625
	Inclusion		283 994 137 625 235 996 65 768
	Properties		1 952	348	5 677	1 436

Table 3 :

 3 Nurse Rostering Problem (Average time in ms, memory in MB).

Table 4 :

 4 Input of the Car Sequencing Problem

	Option Capacity	Car Class
	Name	m/N	1 2 3 4
	Option 1	1/2	0 0 1 1
	Option 2	2/3	0 1 0 1
	Number of cars	15 37 21 27

Unlike Perez and Régin[START_REF] Perez | Efficient operations on MDDs for building constraint programming models[END_REF], the complementary of an MDD M is computed by making the difference between the universal MDD and M . This avoids the need of a dedicated algorithm.

We could also perform the intersection between M DD(S) and the negation of M DD(C) and check whether it is empty or not. However the computation of the negation is required so it does not improve the classical intersection.

factor between 2 and 3 in time and between 1.5 and 3 in memory for GCC and Sum constraints, and almost 5 in time for the sequence constraint. This result was clearly expected, as the inclusion operation is at worst an intersection, without the reduce and compare part of the intersection based inclusion.

In the Car Sequencing problem, the GCC constraint expresses the number of time a car has to be produced and the Sequence constraint represents the maximum capacity of each option (that is, for any subsequence of q consecutive cars, the maximum number of cars that can have this option). In the Nurse Rostering problem, the GCC constraint expresses the demands for a shift (that is the minimum number of nurses required for a given shift) and the Sequence constraint expresses the fact that a nurse must have at least a certain number of breaks in a sliding time window. The Sum constraint does not have any powerful meaning in these problems, but we decided to test it on at least one problem. We did not test the sum constraint for the Nurse Rostering problem, hence the dash symbol.

Learning parameters of a global constraint

We compare the different methods (inclusion, intersection based inclusion and properties) in order to determine the parameters of a global constraint. The tested constraints are the sequence, sum and GCC constraints. To determine the parameters with the inclusion and intersection based methods, a dichotomous search on the parameters is performed. For example, if the constraint SUM(a, b) is satisfied, we test if the constraint SUM(a, b/2) is satisfied: if it is, we test SUM(a, b/4), otherwise we test SUM(a, b * 3/4). We modify the parameters one by one until they are fixed. However, for both methods, we do not compute the sequence parameters for all possible values of q because it would take too much time (we stop at q = 11). The time and memory used when building the MDDs are integrated into the results.

As we can see from these results, using properties to compute the parameters of a global constraint is better than performing successive inclusion checks by at least a factor 65 in time (1 391ms vs 22ms for the sum constraint for the Car Sequencing, in Table 2) and at most a factor 145 in time (283 994ms vs 1 952ms for the GCC constraint for the Nurse Rostering, in Table 3). For the sequence constraint, we reach a factor of 75 with the properties even if we do not