Riccardo Spezia
email: spezia@ens.fr

Francois-Xavier Coudert

Anne Boutin

F.-X Coudert

STUDYING NEUTRAL DIPOLAR ATOMS IN WATER VIA ATOMISTIC SIMULATIONS: STRUCTURE, SPECTROSCOPY AND FORMATION OF AN EXCITONIC STATE

Keywords: Excitonic State, Hydrated Neutral Atoms, Solvated Electron, Condensed Phase Ultrafast Reactions, Molecular Dynamics

Introduction

Cation reduction in solution is an everyday process having a fundamental role in a broad range of fields, from industrial to biological sciences. For this reason several efforts were done to understand the elementary processes involved. In particular the electron transfer, the crucial step thermodynamically and kinetically regulating redox reactions, was largely investigated during the last decades 1 . Pulse radiolysis of a metal ion in solution offers a powerful method to reduce solvated metal ions and then study by time-resolved spectroscopy the formation of metal atoms and their coalescence [START_REF] Baxendale | The Study of fast Processes and transient Species by Electron Pulse Radiolysis[END_REF] . Within this technique a solvated electron is formed, able to reduce the metal ions present in solution.

Different ways can be adopted to generate a solvated electron 3 . Hart and Boag pointed it out in water for the first time in 1962 4 , the so-called hydrated electron, its absorption spectrum being a non-structured band with a maximum value at 1.72 eV, confirming previous hypothesis of its existence 5 . Then the hydrated electron was largely characterized experimentally, including studies of its stationary spectrum in different conditions 6 . More recently, its formation dynamics became the subject of several studies 7 . Furthermore, in last ten years it was the subject of several theoretical and computational studies 8,9 .

Since the solvated electron is a very strong reducing agent, when produced in a solution it is able to reduce cations. Of course, the reduction process is determined by the chemical nature of the cation. Considering, for simplicity, only monovalent cations (similar considerations can be made for a generally charged cation), when a solvated electron is produced in a solution containing a soluble salt of this cation, we have two possible processes

M + + e - sol M 0 (1)
and

M + + e - sol (M + , e - sol ) (2)
When the reaction 1 is favorized, then the atoms dimerize or associate with excess ions

M 0 + M + M 2 + (3)
and, by a multistep process, these species progressively coalesce into clusters. The most studied reduction reaction assisted by a solvated electron is probably that of the silver cation. It was in fact the first ion aqueous solution studied by pulse radiolysis 10 and was recently revisited 11,12 . The product of reaction 1 is also interesting, since it is a neutral atom in a polar solvent, water generally. Such a state, before reacting with other cations, can be stabilized by an atomic induced dipole moment due to the hybridization of atomic orbitals, achieving a large dipole moment. This large dipole moment can then polarize its neighbors. The resulting reaction field enhances the polarizability of the atom. If there is a sufficient number of neighbors, the reaction field can overcome the energy cost of hybridizing the atom entirely and an effective dipolar atom would result. This stabilized state is called excitonic state. The transition to an excitonic phase consisting of dipolar atoms was invoked to understand the metal-nonmetal transition in expanded liquid mercury 13 . Fluid lithium can also undergo a similar transition 14 . Another possibility of having a transition to an excitonic state can occur in a variety of impurity-doped matrix systems, like for example low temperature alkali-doped rare gas solids at low impurity concentration. Such a possibility was described by Logan 15 , considering an alkali atom that, when isolated, contains a single ns electron outside a closed shell. He assumed a sufficiently low impurity concentration, so that the impurity-impurity interactions may be neglected. Such conditions are experimentally possible for these systems [START_REF] Schwentner | Electronic Excitations in Condensed Rare Gases[END_REF] . Moreover alkali impurity atoms can be stabilized at low temperature in solid ammonia or in glasses [START_REF] Meier | [END_REF] . Liquid host matrices are also possible: Li atom in liquid ammonia presenting the characteristics of an excitonic state has been observed from path integral quantum Monte Carlo calculations of Sprik et al. 18,19 . A similar situation was found by quantum/classical and Car-Parrinello molecular dynamics in the case of the neutral silver atom generated in aqueous solution as in reaction 1 20 .

In this review we show our recent progresses addressed to study reaction 1 by means of atomistic simulations and mean-field excitonic state theory, examining in particular the case of the silver atom in water. The case of reaction 2 does not lead to the formation of an excitonic state. It can be computationally studied with the same quantum/classical approach, as was done in the case of sodium 21 , but it cannot be interpreted by means of the same mean field theory. Hence, since we are limiting this brief review to the formation of an excitonic state, here we will show only the case of silver, reaction 1. At this end we employ a mixed quantum/classical molecular dynamics approach based on an adiabatic simulation technique 9 where the solvated electron is treated quantum-mechanically and the cation and bulk water classically, as described in the next section. Within this approach one can obtain the absorption spectrum of the system, that can be directly compared with experiments and excitonic state theory. Finally, some insights into the formation process are shown. They can be obtained putting together information coming from quantum umbrella sampling of the electron-cation association with free simulation aimed to see the spontaneous reactive process 21,22 .

This brief review is organized as follows. In section 2 we give an overview of the computational techniques adopted to have microscopical information about the processes. Then in section 3 we give a short description of the mean-field formulation of the excitonic state transition and we apply it to the case of silver atom in water comparing theoretical results with information taken from the atomistic simulation and experiments. Finally, in section 4, the reactive process is investigated. The paper ends with general conclusions of our studies at the present stage and perspectives for future investigations.

Computational approaches

The systems under investigation are composed by a monovalent cation, an excess electron and solvent water molecules. When the cation is reduced by the electron a neutral atom is formed, otherwise a contact pair or, in a limit case, two noninteracting particles are obtained. The formation and the stability of the neutral atom in water can be studied by MD simulations. Different techniques can be employed, generally following the Born-Oppenheimer approximation. A well established way is to use Car-Parrinello molecular dynamics (CPMD) 23 , where all the interactions are obtained from first principles. Within this approach, one has the advantage that the interactions between particles are not dependent on the system, since no ad hoc parametrization of the potential is needed. Of course, there will be limitations on the dimension of both the system and the time length of the simulations. Another possibility is to divide the system in two parts, one described by classical interactions and the other following the Schrödinger equation. Using this quantum/classical molecular dynamics (QCMD) it is possible to avoid some CPMD limitations, although a parametrization of the potential describing the interactions is needed. In this section we describe the basics of QCMD that can be used to study the reduction of a cation by an electron in solution, while CPMD is a well known method largely presented in literature [START_REF] Remler | [END_REF] and hence we refer the reader to that. At this end we first describe the interactions between the particles of the system -electron, cation and solvent -and then we show the basics of the dynamics and finally the quantum umbrella sampling that can be used to study the reaction leading to the cation reduction by the solvated electron.

Definition of the interactions

Before directly addressing the dynamics, a reasonable description of the interactions between the particles is needed. At this end we adopt an additive perspective for the interactions, where the total potential of the system is given by:

V = V el,cat + V el,solv + V cat,solv + V solv,solv (4)
where V el,cat and V el,solv represent the electron/cation and electron/solvent interactions respectively, V cat,solv is the cation/solvent interaction potential and V solv,solv is the solvent/solvent interaction potential. The two former terms are considered as pure classical two-body interactions, composed, as in usual classical force fields, 25 by a Coulombic and a Lennard-Jones (LJ) term

V ij = 1 4π 0 q i q j r ij + 4 ij σ ij r ij 12 - σ ij r ij 6 (5)
where i and j are two interacting atoms holding to different molecules. For water, the solvent, we adopt the SPC model 26 , while for Ag + /water interactions there were no LJ parameters in the literature, and hence we obtained them 21,20 in order to better reproduce structural and energetic properties of silver cation in water [START_REF] Skipper | [END_REF] . The electron/solvent and electron/cation interaction energies are evaluated adding in the time-independent Schrödinger equation of the electron two pseudopotential terms, leading to Te + Vint (r, S) ψ n (r, S) = E n (S)ψ n (r, S)

where Te denotes the electronic kinetic energy operator, r the electronic coordinates, S the solvent/cation configurations, ψ n the electronic wavefunction and Vint the effective electron-solvent/cation interaction potential. The latter is composed by two contributions, the electron/water, V el/wat , and the electron/cation, Vel/cat , pseudopotentials. For V el/wat we use the pseudopotential developed by Turi and Borgis 28 , that is based on a quantum ab initio calculation for one water molecule plus an additional electron confined in a box in the static exchange theory limit. It has the following form

V el/wat = V SE (r, r O , r H1 , r H2 ) + V pol (r, r O ) (7)
where r O , r H1 , r H2 are the coordinates of oxygen and hydrogen atoms of water molecules. The first term, V SE , takes into account the interactions of the excess electron with the frozen molecular orbitals of a given water, having the form:

V SE = V O (r, r O ) + V H1 (r, r H1 ) + V H2 (r, r H2 ) (8) V X = Q X erf(A X |r -r X |) |r -r X | + B X erf(C X |r -r X |) -erf(D X |r -r X |) |r -r X | (9)
where the index X represents either O or H, Q X is the partial charge of the molecular sites and A X , B X , C X and D X are adjustable parameters. The second term, V pol , is a polarization contribution, accounting for electronic correlation effects, for which we use the form proposed by Barnett et al. 29

V pol = α W 2(|r -r O | 2 + r 2 p ) 2 (10)
where α W is the water molecular polarizability and r p a typical atomic lenght. For the electron/cation pseudopotential we employed the form proposed by Durand and Barthelat 30 with a set of refined parameters for silver 21 . It has the form:

Vel,cat = - 1 r + Ŵ ps (11)
where -(1/r) is the Coulombic interaction term and W ps is the nonlocal term of the pseudopotential, given by

Ŵ ps = l W l (r) Pl (12)
W l = e -α l r 2 n l i=1 c i,l r n i,l (13)
where Pl is the projector on the spherical harmonics Y lm . The single electron wavefunctions ψ n (r, S) are expanded into a basis of 7×7×7 spherical Gaussian functions centered on the nodes of a cubic boxes. The form of integrals and the parameters used to represent the electron wavefunctions within the above pseudopotentials were presented previously in detail 9,21 .

Quantum/classical molecular dynamics

Having in hand the interaction potentials described in the previous subsection, it is possible to study the electron/cation/water systems using mixed quantum-classical molecular dynamics (QCMD). We consider each system as composed by two subsystems, a classical subsystem, the cation and the solvent (water), following classical Newton's equation of motion, and a quantum subsystem, the excess electron, following the Schrödinger equation. Note that the electron is not, strictly speaking, treated quantum-dynamically, since we obtain the electron wavefunction solving the time-independent Schrödinger equation 6; anyway for each classical configuration of the solvent/cation subsystem, we have a new electron wavefunction, that contributes to the forces acting on the classical subsystem via the Hellman-Feynman theorem [START_REF] Hellmann | Einführung in die Quantenchemie[END_REF] , such that: In the studies we are reviewing here, we always computed these forces for the ground state of the excess electron (n = 0), although an extended usage to excited states is also possible.

F q S = -∇ S E 0 (S) = - +∞ -∞ dr ψ 0 (r, S) 2 ∇ S Vint (r, S) (14
For the classical subsystem we immersed the cation in 300 or 800 SPC water molecules, with the excess electron treated with the Born-Oppenheimer dynamics described above. Standard periodic boundary conditions are employed to simulated bulk behaviour with the Ewald summation technique to calculate electrostatic interactions [START_REF] Allen | Computer Simulation of Liquids[END_REF] . The simulations are performed in the NVT ensemble using the Nosé-Hoover thermostat to ensure isothermal conditions [START_REF] Nosé | [END_REF] . The integration of translation and rotation equation of motion is performed using the Gear predictor-corrector algorithm [START_REF] Gear | Numerical Initial Value Problems[END_REF] , fixing the time-step at 0.5 fs, and the simulations runs have a typical time length of a few tens of picoseconds after a stabilization time.

Reaction free energy profile

The cation reduction reaction path can be studied obtaining, from QCMD simulations, free energy curves as a function of the electron/cation distance, r n (S), defined as

r n (S) = | Ψ n |r|Ψ n -R cat | = ξ (15)
The distance ξ is the easiest microscopic property given by simulations we can relate to the formation of an excitonic state, coupled with a wavefunction analysis. The limiting situation ξ = 0 corresponds to the electron located on the atom. Increasing ξ, we have first a distortion of the atom leading to a dipolar atom (the excitonic state, largely described in the next section), then an electron/cation contact pair is formed, while for larger values of ξ we have two distinct solvated species, the hydrated electron and the cation (see fig 1).

Along this coordinate, the free energy function will be given by

F n (ξ) = -k B T ln P n (ξ) ( 16
)
where P n (ξ) is the probability to find the electron in its electronic state n for the electron/cation distance ξ,

P n (ξ) = δ(r n (S) -ξ) n (17)
The solvent average • • • n and the data points r n (S) for the reaction coordinate are calculated by QCMD with the total Hamiltonian

H n (S) = E n (S) + H(S) (18)
where H(S) is the classical kinetic and potential energy due to the solvent and the cation. To sample the needed electron/cation distances we use the Umbrella Sampling (US) approach, adding the following quadratic potential to the Hamiltonian of equation 18:

U ξ0 = 1 2 k(r n (S) -ξ 0 ) 2 (19)
The US technique was adapted for the case of a system including a quantum particle by Borgis and Staib [START_REF] Borgis | [END_REF] and this methodology was succesfully used in the case of cations reduction process 21,22 .

Dipolar excitonic state

In this section we briefly review the basic concepts of the mean field dipolar excitonic state theory based on works of Logan and coworkers 15,36 with a particular focus on the application to a neutral atom in a liquid host matrix. Notably the case of silver in water is analyzed looking for its properties. In particular from excitonic state theory it is possible to obtain the hybridized atomic orbitals and hence the atomic wavefunctions of the states considered. In this way the induced silver dipole and electronic transitions between s and p states can be evaluated and directly compared with QCMD results and experimental data.

General theory

Considering an atom containing a single ns electron outside a closed shell, like elements of groups 1 (alkali atoms) and 11 (for example silver), it is possible to write the isolated atom Hamiltonian considering simply a single ns orbital with energy 0 = -ω 0 and three np orbitals, with energy 1 = ω 0 . Neglecting higher excited states, the atomic Hamiltonian is therefore reduced to a simple four level system (FLS) that can be represented by

Ĥ0 = ω 0 -|s s| + α |p α p α | (20)
where α = x, y or z. An external electric field will introduce the perturbation ĤE = -μ • E 0 , where μ is the dipole moment operator. In the above basis the perturbation Hamiltonian is given by

ĤE = α (|s p α | + |p α s|) ξ α (21)
where ξ α = M E 0 α represents the energy contribution due to the external field and M = s|μ α |p α is the transition dipole moment. The new Hamiltonian, Ĥ0 + ĤE , can be diagonalised leading to four eigenvalues, ω 0 -, ω 0 (doubly degenerate) and ω 0 + , where

ω 0 ± = ± ω 2 0 + ξ 2 (22)
and the corresponding normalized, orthogonal eigenfunctions are given by

|Ψ 0 -= |s + λp 0 [1 + λ 2 ] -1/2 ( 23
)
|ψ 0 0 = |p 1 and |p -1 ( 24
)
|ψ 0 + = |p 0 -λs [1 + λ 2 ] -1/2 (25)
where

λ = ω 2 0 + ξ 2 1/2 -ω 0 |ξ| (26)
Since the new (perturbed) wavefunctions are obtained, one can express the new (perturbed) transition dipole moments in this basis, giving

Ψ 0 -| μ|ψ 0 + = 1 -λ 2 1 + λ 2 M (27) Ψ 0 -| μ|ψ 0 0 = 1 (1 + λ 2 ) 1/2 M (28)
We can now consider that the atom (the impurity) is immersed in an inert (non reactive) condensed phase matrix, solid or liquid. Assuming a sufficiently low concentration of the impurity, we can neglect impurity-impurity interactions, thus considering a single matrix-bound impurity. The mean field theory we are going to shortly review here is based on the following assumptions: (i) the impurity atom valence electron remains strongly localized; (ii) the ns → np excitations in the matrix-bound impurity are of Frenkel type; (iii) there is no overlap with the host matrix absorption bands.

Following the continuum dielectric theory [START_REF] Fröhlich | Theory of Dielectrics[END_REF] , we can express the local field acting on the matrix-bound impurity, f (E), as given by two contributions:

f (E) = G(E) + R(E; µ ) ( 29
)
where G(E) is the electric field inside the cavity due to external sources and R(E; µ ) is the reaction field due to the total impurity dipole moment, µ, that polarizes its sourronding. Hence, the mean field interaction Hamiltonian, Ĥint , can be expressed as

Ĥint = -μ • f (E) (30) = -μ • 3 m 2 m + 1 E + g(ρ m ) µ (31)
where g(ρ m ) is the reaction field factor

g(ρ m ) = 8πρ m 3 [ m -1] [ m + 1] (32)
with ρ m and m being the matrix density and dielectric constant respectively. Hence, the effective Hamiltonian, Ĥi , will be given by Ĥi = Ĥ0 + Ĥint (33) To obtain Ĥint from equation 31 we should know µ , the impurity dipole moment, that, as the expectation value of any quantum operator, is given by

µ = Q -1 i Tr μ exp -β Ĥi ( µ ) ( 34
)
where Tr is over the eigenstates of Ĥi = Ĥ0 + Ĥint .

We have hence a self-consistent equation for the matrix-bound impurity dipole moment. The solution for µ will be a function of both the atomic properties of the free impurity and the thermodynamic properties of the host matrix, i.e. the environment. In the absence of an applied external electric field, (E = 0), the self-consistent equation for µ is

µ 0 y( µ 0 ) = µ 0 tanh [βs(ρ m )y( µ 0 )] 1 + exp [-β ω 0 ] sech [βs(ρ m )y( µ 0 )] (35)
with

y( µ 0 ) = [α 0 g(ρ m )] -1 1 + α 0 ω 0 [g(ρ m ) µ 0 ] 2 (36)
where α 0 is the polarizability of the FLS and

s(ρ m ) = ω 0 α 0 g(ρ m ) > 0 (37)
Note that to solve equation 35 one has only to know properties that characterize the isolated impurity atom, α 0 and ω 0 , and the host matrix, ρ m and m . It is evident that one possible solution of equation 35 is the trivial one, µ 0 = 0, corresponding to the absence of an atomic dipole. Moreover, given the atomic properties α 0 and ω 0 , at a fixed temperature, exists a critical density, ρ m,c such that for ρ m < ρ m,c only the trivial solution µ 0 = 0 can fulfil the self-consistent equation 35. On the other hand, for ρ m > ρ m,c two solutions are possible, µ 0 = 0 and µ 0 > 0, where the latter is favoured thermodynamically. Hence, for ρ m = ρ m,c the system has a transition to a dipolar excitonic state, µ 0 > 0.

A simplified expression of the self-consitent equation 35 can be obtained for the limit situation y( µ 0 ) → 1, leading to

µ 2 0 M 2 = [α 0 g(ρ m )] 2 -1 [α 0 g(ρ m )] 2 (38)
that preserves all the features of a transition to a dipolar excitonic state, and that was conveniently used to analyze the presence of a dipolar excitonic state in liquid systems 19,20 .

Dipolar atom properties

Once the excitonic state is obtained, the original atomic properties will of course be changed. This changes can be thought in terms of hybridization of original s and p orbitals, leading to the wave-function

|Ψ( λ) = |S + λP z 1 + λ2 -1/2 (39)
and the sp-hybrid state will have the following dipole moment

µ α 0 = Ψ( λ)|μ α |Ψ( λ) (40) = 2 λM [1 + λ2 ] δ αz ( 41
)
where λ is a mixing coefficient that here can be regarded as a parameter. It will be useful to express the total free energy change of the system associated to the formation of the dipolar atomic state as a function of λ, taking as reference the normal state, µ 0 = 0, where λ = 0. This total free energy change, ∆A( λ), is composed by two contributions: (i) the hybridization free energy cost, A h ( λ), and (ii) the free energy arising from the interaction with the environment, A s ( λ). Then the total free energy change, ∆A( λ) = A h ( λ) + A s ( λ), will be given by

∆A( λ) = 2 ω 0 λ2 [1 + λ2 ] 1 - α 0 g(ρ m ) [1 + λ2 ] (42)
Note that for ρ m < ρ m,c ∆A( λ) has a minimum at λ = 0, i.e. the normal state is the stable state; on the other hand for ρ m ρ m,c , ∆A( λ) has a minimum for

λ 2 = [α 0 g(ρ m ) -1] [α 0 g(ρ m ) + 1] ( 43
)
Once the dipolar excitonic state is formed (ρ m > ρ m,c ) the energies of the perturbed state described by |Ψ -and |ψ 0 + in equations 23 and 25 will be splitted by

∆ EI = ω 0 1 + α 0 ω 0 [g(ρ m ) µ 0 ] 2 1/2 -1 (44) such that ω ± = ± (ω 0 + ∆ EI ) ( 45
)
Atomistic simulations can be seen in the spirit of this theory looking, e.g., for properties related to the impurity atom that can be directly addressed by them. In this way one can interpret microscopic observations taken from simulations in terms of a mean-field theory based on experimental observables (i.e. an essentially macroscopic theory). In particular the average dipole moment of the impurity can be obtained from MD simulations. This study was done originally by Logan 19 analyzing a simulation study of Li in ammonia of Sprik et al. 18 , where he pointed out that the simple equation 38 can predict the formation of an atomic dipole of Li obtaining a good quantitative agreement with results found in the atomistic simulation. More recently, the behaviour of neutral silver atom in water was studied by two indipendent and different atomistic simulations, CPMD and QCMD 20 . Also in this case a neutral dipolar atom is formed, i.e. an excitonic state, providing a good agreement between theory and simulations (notably 1.96 D from the theory versus 1.9 D and 2.4 D obtained from CPMD and QCMD respectively). This state is particulary interesting also from the point of view of the solvent structures. In fact it mainly presents the typical features of an hydrophobic species, like similar oxygen and hydrogen peaks in the electron/water radial distribution function, but it is in part organized to stabilize the atomic induced dipole, that adiabatically follows the electric field applied by the solvent 20 . As already pointed the silver cation is reduced by the electron in water solution leading to a dipolar atom, that can be well interpreted as an excitonic state . From the same kind of atomistic simulations we also noticed that replacing silver cation with sodium cation a different situation is obtained 21 . Changing the chemical nature of the cation, sodium now, it is no more reduced, as expected, and only a metastable cation/electron contact pair is formed with an average distance of ∼2 Å. This metastable contact pair was also found by Laria in micelles [START_REF] Laria | [END_REF] and its presence was suspected by experiments 39 -of course being a metastable state, it is experimentally difficult to be detected. In this case the basic assumptions of the excitonic theory are no more valid, for example simply because the s electron is no more localized on the atom.

Electronic transitions

Systems where an excess electron is formed in solution in presence of cations are generally investigated experimentally via electronic spectroscopy 40,12 . It is evident, from equations 23-25 and 44-45, that excitonic theory is in principle able to provide transition energies between ns and np states of the impurity. To fully investigate the spectroscopic consequences of a transition to an excitonic state, at least two effects must be incorporated in the theory: the spin orbit coupling interaction and the spatially local, primarly repulsive, interactions (corresponding to the usual van der Waals classical repulsive short range interactions). Hence, the Hamiltonian Ĥi of equation 33 is replaced by

Ĥ int = Ĥi + Ĥso + V0 ( 46
)
where Ĥso is the spin-orbit coupling Hamiltonian and V0 is a spherically simmetric interaction term taking into account the local interactions that do not mix the ns and np orbitals. Then, it is possible to write simple expressions for the total resulting transition energies

∆E 1 = E 0 + δ + 3 2 ∆ EI - 1 4 ζ - 1 2 ∆ 2 EI + ζ∆ EI + 9 4 ζ 2 1/2 (47)
∆E 2 = E 0 + δ + ∆ EI + 1 2 ζ (48)
∆E 3 = E 0 + δ + 3 2 ∆ EI - 1 4 ζ + 1 2 ∆ 2 EI + ζ∆ EI + 9 4 ζ 2 1/2 (49)
where E 0 = 2 ω 0 , ζ is an effective np spin-orbit coupling constant and δ is the energy shift between ns and np orbitals due to V0 (i.e. δ = p| V0 |p -s| V0 |s ). Note that equations 47-49 contain the normal state (∆ EI = 0, ∆E 2 = ∆E 3 ) and the excitonic state (∆ EI = 0). To directly address the transition energies from equations 47-49, δ is probably the most problematic term to estimate. It is also a crucial value, since, disregarding it, the theory can provide only blue-shifted spectra. If a red-shifted spectrum is observed, it can be justified by a negative δ that bilances a positive ∆ EI . Unfortunately, δ cannot be obtained analitically but only via semiempirical considerations. In the case of silver, the explicit case under our interest, we can make a first drastic, but reasonable, assumption considering δ p << δ s and then δ = -δ s . This term, δ s , can be further estimated from the energy of |S state in vacuo, E |S vac , and in solution, E |S sol , as obtained from simulations reported in ref. 21 , and using ∆ EI via equation 44 as the electrostatic contribution of the solvent. Thus,

δ s = E |S sol -E |S vac + ∆ EI (50)
The spin orbit contribution can be taken from recent experimental data 41 . Using the values listed in tab 1 for Ag 0 in water, we can obtain the three transition energies: ∆E 1 = 3.33 eV (≡ 372 nm), ∆E 2 = 3.51 eV (≡ 353 nm) and ∆E 3 = Table 1. Parameters adopted and values obtained to study the excitonic state of silver atom in water at ambient conditions. (a) experimental value from ref. 3548.9 cm -1(c)

3.75 eV (≡ 331 nm). All these values are in good agreement with experimental observations 12 . Moreover from QCMD simulations we can obtain as well the overall spectrum of Ag 0 in water that is in good agreement with experiments and also with the previous theoretical considerations (see figure 2).

Formation reactions

In this section we show how it is possible to obtain via atomistic simulations some information on the reactive process between the solvated electron and a cation in water solution. Generally speaking, one can be interested in equilibrium thermodynamics or kinetics properties. In the former case the quantum Umbrella Sampling described in section 2.3 can be used in order to have the free energy profile of the reaction. This is useful mainly to locate the stable state along an electron/cation distance coordinate, as defined in equation 15. As expected, this free energy curve presents a clear minimum for electron/silver distance ∼ 0.5 Å, as shown in figure 3 (for comparison, the free energy curve for a sodium cation is also plotted). This corresponds to the average distance obtained from unconstrained simulations of the excitonic state described previously. It is also possible to show, from a simple thermodynamic cycle, that the free enthalpy of reduction, ∆G * , is +0.96 eV and -0.74 eV for sodium and silver respectively 22 . In addition, another important thermodynamic information comes out from umbrella sampling calculations, namely the free energy barrier to escape from the stable state. In the case of silver it was interesting to note that the atomic excitonic state presents a high barrier (over 40 k B T ), confirming that this state is both thermodynamically and kinetically stable.

Unfortunately, a full quantitative kinetic information cannot be directly addressed from these curves, since they correspond to an equilibrium solvation regime 42 that assumes the solvent instantaneously relaxing into an equilibrium state. As a consequence, they lack a solvent coordinate that, consistently with Mar- cus theory for electron transfer 1 , should be employed to fully describe the reaction process. Moreover, as we noted recently, the barrier high changes as a function of the number of water molecules considered 22 , preventing us from extracting quantitative kinetic information.

Despite of the previously mentioned limitations, the choise of the electron/cation distance to describe the free energy profile of the reactive process has some advantages. As already noticed, we can determine which is the thermodynamically stable state from the same intuitive definition done in terms of electron/cation distance used in the previous sections. It is also useful since an excitonic stable state can come out directly and it can be related to equilibrium properties of other intermediate states. Moreover, it can be related to experiments done with different cation concentrations leading to different stationary absorption spectra. 39,43,12 As already remarked, from these curves a full kinetic information cannot be directly addressed. Whereas a large barrier can tell us, qualitatively, that the thermal reaction leading to a separated cation/electron state is kinetically unfavoured, a quantitative calculation of kinetic constants needs more effort. Anyway, in the present case we can be sufficiently convinced that this state is kinetically stable other than thermodynamically. Doing independent free QCMD simulations at different starting distances d 0 4 Å, we can obtain some insights into the reduction process (see figure 4). In fact, it is possible to decompose the silver reduction process in the following steps: (i) a diffusive (brownian) step where the electron/cation distance reaches a value of ∼4 Å, (ii) a solvent/solute organization step that allows the decrease of distance down to ∼1 Å, occurring in ∼0.15 ps and (iii) a relatively slow decrease of this distance until the formation of the excitonic state, occurring in ∼0.8 ps. Note that the time scale of the non diffusive processes we identified above (i.e. the last two steps) are in agreement with the experimental hypothesis that the reaction should occurs in less than 1 ps 12 .

To identify the nature of the system between steps (ii) and (iii), that we will call the reduction intermediate state, we obtained its absorption spectrum (figure 5) and the radial density distributions around silver (figure 6) by an average over the first 200 fs of step (iii). Comparing the spectrum with those of a free hydrated electron and of the silver excitonic state, we see that its maximum is close to the silver excitonic state, while its width is closer to the free hydrated electron. Moreover, the study of radial density distributions reveals that this intermediate state has a weak solvation structure, with oxygen atoms still ordered around the silver but hydrogen atoms with less defined position. We interpret these results as a rotational reorganization of the solvation structure, with water molecules still presenting a configuration reminiscent of the cation solvation structure. Thus, we can think that the electron rapid jump on the silver cation (step (ii)) is initiated by a favorable configuration of the fluctuating orientation of water molecules around the cation.

Conclusions and perspectives

Here, we have reviewed our recent progresses in understanding the formation and stability of neutral atoms in water. In particular, we have shown as from atomistic simulations and theoretical analysis it is possible to better understand the physico-chemical basis of such phenomena. At this end we have studied the simplest reduction process, i.e. the direct interaction between a cation, the reducible species, and the electron, the reducent, in solution. We have found, in the case of silver, that the stationary state we obtained from simulations, the excitonic state, was in good agreement with experimental studies, done in the same thermodynamic conditions, and theoretical analysis. In particular, atomistic simulations are able to well reproduce experimental UV-Vis spectra of such systems. This is particulary encouraging in order to extend the research to other similar systems, like other cations in solution and in confined systems. At this end, the study of cations reduction by pulse radiolysis in zeolites is an intriguing perspective, that could help to better understand recent advanced experiments (REF). Moreover, the identification of this excitonic state is particularly interesting since it allows us to analyze it also via theoretical formulation. This analysis is particulary important because we can have another confirm of the good quality of the data obtained from simulations and we can be able to rationalize them in terms of physical macroscopic conditions. Furthermore, we have shown that studying this kind of systems by means of mean field excitonic theory it is also useful to understand some spectroscopic data, obtaining a nice agreement between experiments, theory and simulations. This result can give us the basis to investigate with a sufficient degree of confidence the formation of such neutral atoms in different thermodynamic states (of particular interest can be the supercritical conditions) and in perspective it can be also used to predict the spectroscopic behaviour of other atoms.

Finally, we have shown as we can have a general picture of the overall reactive processes, from both constrained -the umbrella sampling studies -and unconstrained simulations. This gives us a first picture of the elementary reduction reaction. To better understand this point, we are going to complete the studies we have done with other techniques aimed to simulate complex reactions (REF) and also, from the theoretical point of view, with a more detailed study of non-equilibrium solvation effects.

Fig. 1 .

 1 Fig. 1. Schematic representation of different states as a function of the electron/cation distance.

10 R

 10 . Spezia, F.-X. Coudert & A. Boutin

  12 R. Spezia, F.-X. Coudert & A. Boutin

Fig. 2 .

 2 Fig. 2. Absorption spectrum of silver atom in bulk water: experimental (dashed line), calculated from QCMD simulations (solid line) and obtained from excitonic state theory (vertical black lines). The three components obtained from QCMD simulations are separately shown (dotted, dasheddotted and dashed-dotted-dotted) corresponding to s→p transitions. The in vacuum transition is also shown (vertical open line).

Fig. 3 .

 3 Fig.3. Free energy profiles as a function of electron-cation distance obtained for silver (solid) and sodium (dashed) cations.

Fig. 4 .

 4 Fig. 4. Electron/cation distance as a function of time during unconstrained simulations of the reduction process. On the right panel, trajectories have been shifted in time; the bold curve is the mean of all trajectories.

Fig. 5 .

 5 Fig. 5. Absorption spectrum of the excess electron in the reduction intermediate state (dashed line), compared to the Ag excitonic state (solid line) and the free hydrated electron (dotted line).

Fig. 6 .

 6 Fig.6. Radial density functions for oxygen (upper panel) and hydrogen (lower panel) around the silver atom in the reduction intermediate state (thick line), compared to the silver excitonic state (thin line) and the Ag + cation (dashed line).

  41 ; (b) calculated from equation 44;(c) calculated from equation 50.

	Property	Value
	α Ag	7.9	Å3
	M Ag	4.8 D
	ω 0	15082.9 cm -1(a)
	ζ	613.8 cm -1(a)
	m	78.3
	ρ m	1.000 g/cm 3
	∆ EI	1380.3 cm -1(b)
	δ s

R. Spezia, F.-X. Coudert & A. Boutin

Acknowledgements

We wish to acknowledge C. Nicolas, R. Vuilleumier and P. Archirel for fruitful collaborations in this project and D. Borgis, M. Mostafavi, B. Lévy and A. H. Fuchs for several interesting discussions. We thank the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) of the CNRS for generous allocations of computer time and the CNRS/DFG and European Union through Marie Curie project for fundings.


Anne Boutin

Riccardo Spezia

François-Xavier Coudert

Mehran Mostafavi

Molecular Dynamics Simulations of the Temperature and Density Dependence of the Absorption Spectra of Hydrated Electron and Solvated Silver Atom in Water

We report mixed quantum-classical molecular dynamics simulations of the optical absorption spectrum of the solvated silver atom and electron in liquid water. The simple one electron model is shown to be able to reproduce the strong temperature dependence of the absorption spectra of hydrated electron as well as the much weaker dependence for the silver atom. A qualitative explanation is provided for this experimental fact. When extending these simulations to very low densities corresponding to supercritical conditions the results display a progressive "desolvation" of the hydrated electron. Two other distinct theoretical models lead to results similar to those of the present QCMD simulations.

Introduction

Understanding the mechanism of elementary chemical events such as charge transfer in solution is one of the major challenges in chemical physics today. Owing to its importance in many domains of chemistry and biology, the hydrated electron has been the subject of numerous studies. This is also true for the radiation-induced reduction of metal ions in solution. One of the most studied elementary reactions is the reduction of a silver ion in water, following a radiolysis pulse: A full understanding of the mechanism of elementary reactions requires a precise knowledge of the nature of the reactive species. In order to better characterize short-lived species in solution, their optical spectra have been thoroughly investigated. In the case of the hydrated electron, the observation of a broad structureless absorption band centered at 715 nm dates back to the early 1960s [1][2][3][START_REF] Hart | The Hydrated Electron[END_REF]. Soon after, this species was detected in other solvents and the optical spectrum of the solvated electron was found to depend strongly on the nature of the solvent [START_REF] Jou | [END_REF]. The maximum of the absorption band can be shifted to as high as 2200 nm in tetrahydrofuran. The temperature effect on the optical spectrum has also been studied [6][7][8][9][10], and it has been recently investigated in supercritical water [11,12].

The metallic silver ion/solvated atom (eq.1) is being considered as a model system for the radiolysis-induced charge transfer in solution [13][14][15][16][17][18][19][20][21]. The optical absorption spectra of the transient silver atom, as well as that of the charged dimer (produced by an aggregation process with an excess neighboring silver ion), have been observed by pulse radiolysis in various solvents [START_REF] Belloni | Radiation Chemistry[END_REF]. An indication of the influence of the atom-solvent interaction is given by the solvent dependence of the optical absorption spectra in solution. In ammonia, the absorption spectra of silver atoms and dimers were observed at two different temperatures: 296 K [START_REF] Belloni | Radiation Chemistry[END_REF] and 223 K [START_REF] Farhataziz | [END_REF]. The maximum of absorption band was found at 435 nm and 450 nm, respectively. More recently [21], the rate constant for the reduction of silver ions by hydrated electrons was measured at different temperature and is found to fit with an Arrhenius relationship (from 398 K to 573 K at 200 atm) with an apparent activation energy of 17.8 kJ.mol -1 . The measurement of transient absorption spectra showed that the maximum of the absorption spectrum of silver atom in water is shifted from 355 nm to 370 nm when the temperature is increased from 298 K to 473 K.

In this paper, we concentrate on the (P, ρ, T ) dependence of the optical spectra of the hydrated electron and silver atom, rather than the solvent dependence. Until very recently [24], the red shift of the maximum of the absorption band of the hydrated electron with increasing temperature had not been investigated theoretically. Based on a mixed quantum-classical molecular dynamics simulation, Nicolas et.al. [24] reported that a "cavity picture" of the hydrated electron could qualitatively capture most of the known experimental data. It was also suggested that the red shift of the absorption spectrum could be a density, rather than a temperature effect. No such investigation exists, to our knowledge, in the case of the hydrated silver atom.

In this work, we have performed adiabatic molecular dynamics simulations of both the hydrated electron (extending our previous simulations to supercritical conditions) and hydrated silver atom species. In the ordinary liquid conditions, we show that simple one electron, mixed quantum-classical simula-tions are able to reproduce the strong thermodynamic state point dependence of the optical spectra of the hydrated electron, as well as the weak shift of the absorption spectrum in the case of silver atom. However, this is not the case under supercritical conditions where the simulation results are found to be at variance with recent pulse radiolysis experiments [12].

Methods

We have performed mixed quantum classical molecular dynamics (QCMD) simulations of an excess electron in bulk water and in presence of a silver cation. The excess electron only is treated quantum mechanically, using the Born-Oppenheimer approximation. The method and the simulation details are fully described elsewhere [24,25].

For the excess electron-water interaction we used the pseudo-potential developed by Turi and Borgis [26] while for electron-silver interaction we used the form initially proposed by Durand and Barthelat [27], which we have adapted to a Gaussian grid representation of the electron [25]. Water and cations interact via classical interaction potentials made by coulombic and 6-12 Lennard Jones (LJ) terms. Water-water interactions are described by the SPC model [START_REF] Berendsen | Interaction models for water in relation to protein hydratation[END_REF].

The simulations were carried out in a box of 300 SPC water molecules, an excess electron and/or a silver cation. Standard periodic boundary conditions were used together with the Ewald summation technique for computing the electrostatic interactions. The MD time step was fixed at 0.5 fs and simulation runs lasted a few tens of picoseconds after the stabilization period.

The method used for computing the electron absorption spectrum is fully described elsewhere [24,25,[START_REF] Spezia | [END_REF]. Different thermodynamic conditions were simulated, from ambient to supercritical conditions. Since the SPC model fails at perfectly reproducing the equation of state of liquid water in the whole thermodynamic state space, we decided to perform (N,V,T) ensemble simulations, and to fix the density of the system to the value found in the NIST WebBook [START_REF] Linsstrom | NIST Chemistry WebBook[END_REF], for each given temperature and pressure. For instance, in order to compare our hydrated silver atom simulations to experiments, we used the density values corresponding to the thermodynamic state points investigated by Mostafavi et al. [21], namely 1.00, 0.95 and 0.86 g.cm -3 which correspond to 298, 376 and 476 K respectively, at a pressure of 20 MPa.

Results and discussion

The absorption spectra of the hydrated silver atom at 20 MPa and three different temperatures are shown in Figure 1. The maximum absorption wavelengths as well as the bandwidths are fairly well reproduced by the QCMD simulations. It is interesting to mention that the absolute maximum intensity of the simulated spectra is only weakly dependent on temperature (Table 1). Experimentally, the absorption intensity is often presented as the product of G. (where G is the radiolytic yield and is the absorption extinction coefficient). As shown in Table 1, the experimental maximum intensity does not depend markedly on temperature, at constant pressure. The comparison between experiments and simulations seems to provide an indication that the change in intensity of the experimental spectra could be mainly attributed to a decrease of the radiolytic yield (since the density of the system decreases with increasing temperature), rather than a decrease of . This is also true for the hydrated electron simulations where the maximum intensity of the absorption spectrum is found to vary within less than 20% whatever the density or temperature is.

In Figure 2, we show the evolution with temperature (at constant pressure) of the maximum of the absorption spectrum of the hydrated silver atom species. As already noted above, the simulation results are in good agreement with experiments. In the same figure, we report the previously published shift in the case of the hydrated electron alone [24]. It is clear from these data that the temperature shift is much weaker in the case of silver atom that it is for the hydrated electron. In the latter case, it was suggested that the red shift of the absorption spectrum could be a density, rather than a temperature effect [24].

Here, we have carried out QCMD simulations of the hydrated silver atom at several thermodynamic state points in order to shed some light on the weak shift observed in this case. The results are given in table 2.

At a constant density of 1.0 g.cm -3 , a temperature effect (red shift) is obtained. On the other hand, at a constant temperature of 573 K, a density effect (blue shift) is observed. At constant pressure, these two effects compensate each other and this leads to the weak temperature effect seen in Figure 2.

In a previous QCMD work, we have shown that the hydrated silver atom behaves as a dipolar excitonic state [START_REF] Spezia | [END_REF]. This is to say that the active electron and the silver cation act as a strong, fluctuating, dipole. The average distance between the electron and the cation can be computed from the QCMD simulations. It is shown, in Table 2, that a strong correlation exists between the change in the maximum of the absorption spectrum and the electron-cation distance.

From these findings, we suggest that the behavior of the active electron in the hydrated silver atom species is mainly governed by the presence of the neighboring silver cation. The thermodynamics state point slightly affects the electron-cation distance, and this has an effect on the absorption spectrum. Two distinct temperature and density effects are observed, but these effects are weak. This is in clear contrast with the behavior of the hydrated electron, for which the thermodynamic state point effect was shown to be much stronger. It was suggested to be more or less related to the size of the cavity in which the electron was sitting, and in turn to the density of the system.

We have reinvestigated the density effect on the optical spectra of the hydrated electron, in order to compare the QCMD data to other simulation results of the literature. Our aim was also to extend our simulations to the supercritical conditions, for which some new experimental data became available recently [12,[START_REF] Katsumura | [END_REF]. In Figure 3, we report the QCMD computed evolution of the maximum of the absorption spectrum of the hydrated electron as a function of density, and compare our data with the results obtained using other methods, namely Path Integral MD (PIMD [START_REF] Laria | [END_REF]) and Car-Parrinello MD (CPMD [33]). Interestingly enough, the PIMD and CPMD simulations are in rather good agreement with the present QCMD data and were performed at different temperatures. This seems to confirm the existence of a strong density effect, the maximum of the absorption spectrum decreasing linearly with the density of the system, to as low as 0.1 g.cm -3 . It is worth mentioning that PIMD, CPMD and Adiabatic (this work) simulations rely on rather different methods and approximations.

In a recent experimental work, however, Bartels and coworkers [12] clearly observe a "levelling off" of the density effect. At densities lower than 0.5 g.cm -3 (supercritical conditions), the maximum of the absorption spectra is near 1 eV, and does not change anymore with decreasing densities. An analysis of their data leads Bartels and coworkers to suggest that the inner solvation shell of the electron remains constant in this density region [12]. This is in contrast with the QCMD results which display a progressive and continuous "desolvation" of the electron as the density decreases. This is shown in the radial distribution functions reported in Figure 4. The structure of the solvation shell is progressively lost when the density decreases and the first peak of the radial distribution functions is shifted to larger distances. According to the unpublished results of the Katsumura group [START_REF] Katsumura | [END_REF], even in the supercritical domain a slight decrease of the maximum of absorption is observed from 1.04 eV to 0.97 eV, corresponding to densities from 0.5 to 0.2 g.cm -3 . Such a slight shift has also been observed by Bartels and coworkers [12]. However, our calculations show a marked decrease of E max from 1.15 to 0.87 eV for the same range of density.

Conclusions

We have shown that QCMD simulations are able to capture most of the important features of the optical spectrum of the silver atom and electron species solvated in ordinary liquid water. The strong thermodynamic state point dependence of the absorption spectra of hydrated electron as well as the much weaker dependence for the silver atom are well reproduced, and a qualitative explanation is provided for this experimental evidence. When extending the QCMD calculations to very low densities corresponding to supercritical conditions the simulation results display a progressive "desolvation" of the hydrated electron. This is at variance with recent pulse radiolysis experiments [12] in which the solvation shell is suggested to "level off". The reason for this discrepancy is not clear, since two other distinct theoretical models lead to results similar to those of the QCMD simulations.

aq

  Article published in Chem. Phys.Lett. 409, 219 (2005)

Fig. 1 .

 1 Fig.1. Experimental (dashed lines[21]) and simulated (full lines) normalized absorption spectrum of hydrated silver atom at 20 MPa and 298 K (a), 373 K (b) and 473 K (c).

Fig. 2 .Fig. 3 .

 23 Fig.2. Temperature dependence (at constant pressure) of the maximum of the absorption spectrum of hydrated silver atom. Full squares: experiments[21], open squares: simulations (this work). Also shown are the data for the hydrated electron. Full circles: experiments[11], open circles: simulations[24].

Fig. 4 .

 4 Fig. 4. Radial distribution functions for the electron-oxygen (upper panel), and electron-hydrogen (lower panel) distances at different densities. Solid lines: 1.0 g.cm -3 . Dashed lines: 0.48 g.cm -3 . Dotted lines: 0.20 g.cm -3 .

Table 1

 1 Relative intensity of the maximum of the absorption spectrum of the hydrated silver atom at three temperatures (at 298 K, 360 nm =16000 mol -1 cm -1 ). The experimental and simulation data are normalized with respect to the ambient thermodynamic conditions data.

	Temperature / K	298		373	473
	Density / g.cm -3	1.00	0.95	0.86
		Exp [21]	1.0		0.97	0.87
		QCMD	1.0		0.97	0.92
	Table 2
	Maximum of the absorption spectrum of hydrated silver atom (in eV) as a function
	of temperature and density, from QCMD simulations. The corresponding average
	electron-silver cation distance (in Å) is given in parenthesis.
	Temperature		Density / g.cm -3
	/ K	1.00	0.95	0.86		0.73	0.50
	298	3.59 (0.52)	-	-		-	-
	373	3.52 (0.61) 3.49 (0.59)	-		-	-
	473	3.49 (0.63)	-	3.47 (0.59)	-	-
	573	3.47 (0.66)	-	-	3.47 (0.52) 3.57 (0.47)

Acknowledgements

We are most grateful to David Bartels and Yosuke Katsumura for fruitful discussions, and for providing us with recent supercritical water results (ref. [12,31]) prior to publication.


Yuxiang Zhou

Jean-Baptiste Caillau

Marc Antonini

Michel Barlaud
email: barlaud@i3s.unice.fr

Robust classification with feature selection using alternating minimization and Douglas-Rachford splitting method

This paper deals with supervised classification and feature selection. A classical approach is to project data on a low dimensional space with a strict control on sparsity. This results in an optimization problem minimizing the within sum of squares in the clusters (Frobenius norm) with an 1 penalty in order to promote sparsity. It is well known though that the Frobenius norm is not robust to outliers. In this paper, we propose an alternative approach with an 1 norm minimization both for the constraint and the loss function. Since the 1 criterion is only convex and not gradient Lipschitz, we advocate the use a Douglas-Rachford approach. We take advantage of the particular form of the cost and, using a change of variable, we provide a new efficient tailored primal Douglas-Rachford splitting algorithm. We also provide an efficient classifier in the projected space based on medoid modeling. The resulting algorithm, based on alternating minimization and primal Douglas-Rachford splitting, is coined ADRS. Experiments on biological data sets and computer vision dataset show that our method significantly improves the results obtained with a quadratic loss function.

Introduction

This paper deals with supervised classification with feature selection in high dimensional spaces. In this paper we consider methods where feature selection is embedded in a classification process [START_REF] Saeys1 | A review of feature selection techniques in bioinformatics[END_REF]. However, classification in high dimension suffers from the curse of dimensionality: As dimensions increase, vectors become indiscernible and the predictive power of the aforementioned methods is drastically reduced [START_REF] Barlaud | On k-anonymity and the curse of dimensionality[END_REF][START_REF] Radovanovic | Hubs in space : Popular nearest neighbors in highdimensional data[END_REF]. In order to overcome this issue, the main idea of the following methods is to project data in a low dimensional space. A popular approach for high-dimensional data is to perform Principal Component Analysis (PCA) prior to classification. This approach is however not relevant in general [START_REF] Wei-Chien | On using principal components before separating a mixture of two multivariate normal distributions[END_REF]. Partial least squares method closely related to principal component regression is designed to cope with this issue with high dimensional correlated features [START_REF] Wold | Pls-regression: a basic tool of chemometrics[END_REF][START_REF] Rosipal | Overview and recent advances in partial least squares[END_REF]. An alternative approach is to perform dimension reduction by means of Linear Discriminant Analysis (LDA) [START_REF] De La Torre | Discriminative cluster analysis[END_REF][START_REF] Ding | Adaptive dimension reduction using discriminant analysis and k-means clustering[END_REF]. The authors of [START_REF] Bach | Diffrac: a discriminative and flexible framework for clustering[END_REF] and [START_REF] Flammarion | Robust discriminative clustering with sparse regularizers[END_REF] propose a convex relaxation of this approach in terms of a suitable semi-definite program (SDP) at the cost of an increased computational complexity. Nie et al proposed a feature selection based on 2,1 norm minimization [START_REF] Nie | Efficient and robust feature selection via joint l2,1-norms minimization[END_REF]. A popular approach for selecting sparse features in supervised classification or regression is the Least Absolute Shrinkage and Selection Operator (LASSO) formulation [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Friedman | Regularization path for generalized linear models via coordinate descent[END_REF][START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF][START_REF] Ng | Feature selection, l 1 vs. l 2 regularization, and rotational invariance[END_REF][START_REF] Zou | Sparse principal component analysis[END_REF]. The LASSO formulation uses the 1 norm [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF][START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF][START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via 1 minimization[END_REF][START_REF] Donoho | Signal recovery and the large sieve[END_REF] as an added penalty term instead of an 0 term. However, an issue is that using an 2 norm for the data is not robust to outliers. In this paper, we cope with this problem by minimizing an 1 norm both on the penalty term and the loss function. In this case, the criterion is convex but not gradient Lipschitz. We propose so to use a Douglas-Rachford splitting method. This splitting was successfully used in signal processing [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF][START_REF] Combettes | A douglas-rachford splitting approach to nonsmooth convex variational signal recovery[END_REF][START_REF] Douglas | On the numerical solution of heat conduction problems in two or three space variables[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems withÂapplications to imaging[END_REF][START_REF] Setzer | Split bregman algorithm, douglas-rachford splitting and frame shrinkage[END_REF][START_REF] Steidl | Removing multiplicative noise by douglas-rachford splitting methods[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. However, for classification, we cannot apply straightforwardly Douglas-Rachford since the proximal operator for the affine transform Y -XW involved in the criterion is not available. We take advantage of the particular form of the original cost: The sum of two 1 norms is equal to a single 1 norm after a change of variables, and the linear constraint can be integrated into the cost. We also provide an efficient classifier in the projected space based on medoid modelling.

2 Robust supervised classification

A robust framework

Let X be the nonzero m × d matrix made of m line samples x 1 , . . . , x m belonging to the ddimensional space of features. Let Y ∈ {0, 1} m×k be the label matrix where k 2 is the number of clusters. Each line of Y has exactly one nonzero element equal to one, y ij = 1 indicating that the sample x i belongs to the j-th cluster, see [START_REF] Bach | Diffrac: a discriminative and flexible framework for clustering[END_REF][START_REF] Flammarion | Robust discriminative clustering with sparse regularizers[END_REF] . Let W ∈ R d×k be the projection matrix, k d. The classical quadratic loss criterion in the projected space is:

min W 1 2 Y -XW 2 F + λ W 1 (1)
where . F stands for the Frobenius norm. Projecting the data in lower dimension is crucial to be able to separate them. Besides, it is well known that the quadratic Frobenius loss criterion is not robust to outliers, so we propose to minimize instead the 1 loss cost, with an 1 penalty regularization to promote sparsity and induce feature selection. So, given the matrix of labels Y , we consider the following convex supervised classification problem:

Problem 1 min W Y -XW 1 + λ W 1 .
Here, the 1 norm of an m by n matrix A denotes the 1 norm of its vectorization:

A 1 := A(:) 1 = m i=1 n j=1 |a ij |.

An equivalent formulation

The cost in Problem 1 is the sum of two 1 norms, which suggests to use a splitting algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Such methods are very efficient to minimize the sum of convex functions and do not require differentiability properties. Note indeed that, having replaced the squared Frobenius norm by the 1 norm of Y -XW , it is not possible to use forward-backward splitting [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] as none of the functions is differentiable. In order to use the more general Douglas-Rachford scheme (see next subsection) that is able to cope with mere convex functions, one still has to be able to compute their proximity operators. Now, while the prox of the 1 norm is well known and expressed in terms of soft thresholding, there is no explicit expression for the prox of the 1 norm of the affine transform Y -XW . We propose to introduce the auxiliary variable ζ := (Y -XW )/λ in R m×k and to minimize

min W W 1 + ζ 1 under the affine constraint XW + λζ = Y.
The sum of the two 1 norms is equal to the single 1 norm of the augmented variable W := (W, ζ) ∈ R (d+m)×k . Let C be the affine subset of (d + m) × k matrices such that X W = Y where

X := [X λI m ] ∈ R m×(d+m) .
The problem can be recast as min Then problem 1 is equivalent to

Problem 2 min W W 1 + i C ( W ).
This new problem is readily equivalent to the previous one as

W 1 = W 1 + ζ 1 .
(Note that the cost has been rescaled by a factor λ.) Problem 2 involves two convex functions whose proximity operators can be efficiently computed (see, e.g., Lemma 1), which allows to use a Douglas-Rachford scheme.

Remark 1 Our approach can be seen as a special case of the alternating-direction method of multipliers (ADMM, see [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]). Indeed, Problem 1 falls into the following class:

min (x,y) f (x) + g(y)
under the affine constraint y = Ax, where x and y are vectors of finite dimensional spaces. Then ADMM considers the augmented Lagrangian (γ > 0 is a scalar parameter),

L γ (x, y, z) := f (x) + g(y) + 1 γ (z|Ax -y) + 1 2γ Ax -y 2 .
The Lagrangian is minimized over x, then over y, and a third step of the method updates the Lagrange multiplier z using a proximal maximization step. See also [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems withÂapplications to imaging[END_REF] for a similar approach. In our case, we take advantage of the fact that the two functions f and g are simple convex functions ( 1 norms) whose sum is again a known convex function (yet another 1 norm). This eliminates one function and allows to penalize the affine constraint to integrate it into the cost. The resulting proximal operator is a simple projection operator on a linear subspace, which is easily computed (see, e.g., Lemma 1).

Douglas-Rachford splitting

The method was initially proposed by [START_REF] Douglas | On the numerical solution of heat conduction problems in two or three space variables[END_REF] for solving matrix equations and used in signal and image processing [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF][START_REF] Combettes | A douglas-rachford splitting approach to nonsmooth convex variational signal recovery[END_REF][START_REF] Setzer | Split bregman algorithm, douglas-rachford splitting and frame shrinkage[END_REF][START_REF] Steidl | Removing multiplicative noise by douglas-rachford splitting methods[END_REF]. Problem 2 amounts to minimizing the sum of two convex functions F ( W ) + G( W ), and the (constant-step) Douglas-Rachford scheme is the following [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]:

Fix ε > 0, τ > 0, γ ∈ (ε, 2 -ε), V 0 ∈ R (d+m)×k
, and define

W n := prox τ F (V n ), (2)
V n+1 := V n + γ(prox τ G (2 W n -V n ) -W n ). (3)
Theorem 1 [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] For ε ∈ (0, 1), τ > 0 and γ ∈ [ε, 2 -ε], every sequence generated according to (2-3) converges to a solution to Problem 2.

We recall that, given a lower semicontinuous convex function f from a vector space to R ∪ {+∞}, the proximity operator prox f of f maps a given x to the unique minimizer of f (y) + (1/2) x -y 2 . In our case, F is the 1 norm, so prox τ F is given by soft thresholding (parameterized by τ ). Indeed, one can separate the variables and use the prox of the scaled absolute value dimension one:

soft(x, τ ) = x + τ if x < -τ, = 0 if x ∈ [-τ, τ ], = x -τ otherwise.
Since G is an indicator function, whatever τ the associated prox simply is the projection operator on the affine subspace C of (d + m) × k matrices. For the sake of completeness, we recall the associated expression of this projection.

Lemma 1 Let A be an m × n matrix of rank m < n, and let b be a vector in R m . The orthogonal projection of z ∈ R n on the affine subspace {x

∈ R n | Ax = b} is proj(z, A, b) = z -A T (AA T ) -1 (Az -b). (4)
Proof. The projection z is characterized by the fact that z -z belongs to (ker A) ⊥ = Im A T , plus the fact that Az = b. So there is some ȳ in R m such that z = z -A T ȳ, and

y = Az -AA T ȳ.
Since AA T is invertible, ȳ = (AA T ) -1 (Az -y) which gives the desired expression.

The lemma can readily be applied in R (d+m)×k to obtain the prox of G, independently of τ , and

(2-3) can be translated into Algorithm 1.

Algorithm 1 Supervised classification Douglas-Rachford algorithm. The operators soft and proj denote the soft thresholding and projection operators, respectively. 1:

Input: X, λ, Y , γ, τ , V , N 2: X ← [X λI m ] 3: for n = 0, . . . , N do 4: W ← soft(V, τ ) 5: V ← V + γ(proj(2 W -V, X, Y ) -W ) 6: end for 7: (W, ζ) ← W 8: Output: W
Feature selection and scalability. In applications, particularly in biological ones, the issue of feature selection may be even more important than classification itself. This selection is achieved by means of the 1 sparsity inducing penalty in Problem 1 cost. A feature i ∈ {1, . . . , m} is then selected if the corresponding line in the matrix of weights W is not zero ( W (i, :) = 0). The set of nonzero columns is interpreted as the signature of the corresponding cluster. In order to evaluate this signature, we compute the projected matrix XW . Each block shows the correlation between the data and the signature. We precompute X T (XX T ) -1 ) only once. Note that matrix (XX T ) is a small (m × m) matrix, so we can use fast adapted algorithm to compute the inverse. The complexity of the resulting iterations is O(d × k × d) for the proximal part, plus O(d × k) for the projection.

Classification using medoid

In order to build a robust classifier, we compute a medoid for each cluster. For the j-th cluster, a medoid µ j is any member of the class minimizing the average dissimilarity inside the class in the projected space: µ j := (XW )( ī, :) where

ī ∈ arg min i s.t. yij =1 s.t. y j =1, =i (XW )(i, :) -(XW )( , :) 1 .
Let us define µ ∈ R m×k , the medoid matrix of the clusters. Computing medoids minimization can be reformulated as minimizing the following 1 norm:

min µ Y µ -XW 1 .
Then, a new query x (a dimension d row vector) is classified according to the following rule: It belongs to the (supposedly unique) class j such that

j ∈ arg min j=1,k xW -µ j 1 .
The cost of medoid computation is expected to be O(m × k) in average [START_REF] Bagaria | Medoids in almost linear time via multi-armed bandits[END_REF].

Alternating Minimization

In this section, we iterate the previous procedure. Let define now the within cluster sum of absolute difference criterion:

Problem 3 min W,µ Y µ -XW 1 + λ W 1 .
We propose an alternating (or Gauss-Seidel) scheme. Given the medoid matrix µ, the first subproblem in W is solved by using our primal Douglas-Rachford algorithm. Conversely, given the matrix of weights W , the second subproblem in µ is solved by medoid computation on the projected data as explained in the previous section. Algorithm 2 summarizes the resulting alternating minimization.

Algorithm 2 Supervised classification with alternating minimization and Douglas-Rachford splitting algorithm.

1: Input:

X, µ = I m , λ, Y , γ = 1, τ , V , N , L 2: X ← [X λI m ] 3: for l = 0, . . . , L do 4:
for n = 0, . . . , N do 5:

W ← soft(V, τ ) 6: V ← V + γ(proj(2 W -V, X, Y µ) -W ) 7:
end for 8:

(W, ζ) ← W 9:
µ ← medoids(Y, XW ) 10: end for 11: Output: W, µ Proposition 1 As each step of the alternating minimization scheme decreases the 1 norm Y µ -XW 1 , which is nonnegative, the following readily holds. The 1 norm Y µ -XW 1 converges as the number of iterates L in Algorithm 2 goes to infinity. This property is illustrated in the next section on real data.

Application to real datasets 3.1 Experimental settings

We compare the labels obtained from our classification with the true labels to compute the MCE (misclassification error, i.e. the number of misclassified observations divided by the number of observations) on the test set. We provide TSNE results for visual evaluation [START_REF] Van Der Maaten | Visualizing high-dimensional data using t-sne[END_REF]. In order to evaluate the signature, we plot the heat map of matrix XW in order to visualize the correlation between the data and the signature. We compare our ADRS method with four state of the art methods: Filter method using Ttest (we rank the features according to their p-values using a PLS classifier (Partial least squares) [START_REF] Wold | Pls-regression: a basic tool of chemometrics[END_REF]), PLS with a sequential selection of features [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF], and the Frobenius norm criterion (1) using the classical forward-backward algorithm. We compare algorithms on three public datasets which are available at https://archive.ics.uci.edu/ml/datasets. In our experiments, we set γ = 1 and N = 10 (see Figure 1 for preliminary observations on convergence). Results are reported after a 4-fold cross-validation.

Dataset: Ovarian [21]

The data available on UCI data base were obtained from two sources: The National Cancer Institute (NCI) and the Eastern Virginia Medical School (EVMS). All the data consist of mass-spectra obtained with the SELDI technique. The samples include patients with cancer (ovarian or prostate cancer). Healthy or control patients form a set of 216 samples with 15000 features. Applying classical Ttest, there are 40% p-values smaller than 0.05, thus 6000 features that have strong potential discrimination.

Conclusion

In this paper, we propose to minimize the 1 norm related to the data term with an 1 regularization. This problem can be solved by classical ADMM algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] O'connor | Primal-dual decomposition by operator splitting and applications to image deblurring[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems withÂapplications to imaging[END_REF]. Our contribution is twofold. First, we provide a new primal Douglas-Rachford splitting algorithm adapted to the criterion. Second, we propose a new efficient classification algorithm (ADRS) using medoid and alternating minimization. Experiments on biological datasets and a computer vision dataset show that our method significantly improves the results of methods using a quadratic loss function.

5 Appendix: Minimize Frobenius norm using a gradient-projection splitting method

The criterion is given by:

min W 1 2 Y -XW 2 F + λ W 1
To solve this problem, we use a gradient-projection method. It belongs to the class of splitting methods ( [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Mosci | Solving structured sparsity regularization with proximal methods[END_REF][START_REF] Sra | Optimization for Machine Learning[END_REF]) using separately the convexity properties of the Frobenius cost on one hand, and of the convexity of the set C on the other. We use the following forward-backward scheme to generate a sequence of iterates. For any fixed step γ ∈ (0, 2/σ 2 max (X)), the forward-backward scheme applied to the above criterion converges towards a solution.

Algorithm 3 Forward-Backward splitting

Input: X, Y, W 0 , N, γ, τ for n = 0, . . . , N do V ← W -γX T (XW -Y ) W ← prox(V, τ ) end for Output: W where prox(V, τ ) denotes the proximity operator of the 1 norm (soft thresholding) already define in section 2.

W W 1

 1 under the affine constraint X W = Y Let moreover i C be the indicator function of the set C, vanishing on C and equal to +∞ outside C.

Figure 1 :

 1 Figure 1: Convergence of our primal Douglas-Rachford algorithm: The decay of Y µ -X W 1 versus the number of iterations N emphasizes the fast convergence.

Figure 2 :

 2 Figure 2: Ovarian: The data in the projected space XW show that the centroids and the medoids are very close.

Figure 3 :

 3 Figure3: The MCE versus the number of selected features selected on Ovarian shows that MCE decreases with the number of genes for ADRS, Frobenius and PLS. Note that Ttest assume that data are independent, which is not the case for proteomics dataset.

Figure 4 :

 4 Figure 4: The heatmap of correlation between data and signature, left with ADRS, right with a quadratic (Frobenius) cost.

Figure 5 :

 5 Figure 5: The Tsne of the the data in the projected space Xw shows that the centroids and the medoids are very close in the minority class while they are far in the majority class.

Figure 6 :

 6 Figure 6: The MCE and the number of selected features selected on Facedrive shows that 1500 features are necessary to obtain the best MCE.

Figure 7 :

 7 Figure 7: heatmap of correlation, Left using ADRS , Right using Frobenius

Figure 8 :

 8 Figure 8: The Tsne of the data in the projected space Xw shows that the distance between centroids and medoids depends on the clusters.

Figure 9 :

 9 Figure 9: The MCE versus the number of selected features selected on Cancer dataset shows that only 20 features are necessary to ADRS to obtain an accurate MCE.

Table 1 :

 1 MCE on Ovarian dataset (216 samples, 15000 features, k = 2 clusters which have 121 and 95 samples, respectively): For 100 selected features, ADRS outperforms Ttest by 11.1% , Frobenius by 7% and PLS by 7.9%).

	Ovarian	Ttest Frobenius PLS ADRS
	MCE % 100 features 14.82	10.65	11.57	3.70
	MCE % 200 features 12.03	6.48	10.65	2.78
	MCE % 300 features 13.43	6.48	9.72	1.85

Table 2 :

 2 This table shows the mean of the correlation on the diagonal of the previous heatmap. The signature obtains with ADRS readily improves the one obtained using a Frobenius norm.The FaceDrive dataset contains image sequences of subjects while driving in real scenarios. It is composed of 606 samples with 4800 features each, acquired over different days from 4 drivers with several facial features like glasses and beard. A set of labels assigning each image into 3 possible clusters are given. Note that this dataset is highly unbalanced.

	Ovarian	Cancer Normal
	ADRS	0.43	0.55
	Frobenius	0.12	0.16
	3.3 Dataset:FaceDrive [13].

Table 3 :

 3 MCE on Facedrive (606 samples, 6400 features, k = 3 clusters which have 27, 546 and 33 samples, respectively): For 1500 selected features, ADRS outperform Frobenius by 2.3%, and PLS by 3.6%).

	Facedrive	Ttest+svm Frobenius PLS ADRS
	MCE % 1000f eatures	7.75	8.89	8.41	5.87
	MCE % 1500f eatures	8.25	6.6	7.92	4.29
	MCE % 2000 features	7.43	6.08	7.76	4.01

Table 4 :

 4 This table shows the mean of the correlation on the diagonal of the previous heatmap. Again, ADRS signature improves the signature obtained with a Frobenius cost. This data is a subset of the RNA-Seq (HiSeq) PANCAN dataset, of patients having different types of tumor: BRCA, KIRC, COAD, LUAD and PRAD. It is available at https://archive.ics.uci.edu/ml/datasets and https://www.synapse.org. This dataset is composed of 801 cells with 20531 genes. There are k = 5 clusters which have 300, 78, 146, 141, and 136 samples, respectively.

	FaceDrive Cluster1 Cluster2 Cluster3
	ADRS	0.42	0.04	0.32
	Frobenius	0.17	0.02	0.10
	3.4 Dataset: Gene expression cancer RNA-Seq Data Set.

Table 5 :

 5 MCE on Cancer dataset ( 801 cells and 20531 genes with k = 5 clusters ): For 20 genes, ADRS outperforms Frobenius by 0.8% , and PLS by 6.4%.Figure 10: heatmap of correlation, Left using ADRS , Right using Frobenius

	Cancer	Ttest+svm Frobenius PLS ADRS
	MCE % 10 genes	10.74	5.12	14.11	4.12
	MCE % 20 genes	7.99	2.11	7.74	1.34
	MCE % 40 genes	10.49	1.5	7.49	1.3
	MCE % 80 genes	5.99	1.1	4.74	1.0

Table 6 :

 6 This table shows the mean of the correlation on the diagonal of the previous heatmap, it shows that ADRS signature outperforms Frobenius signature.

	Pancancer BRCA COAD KIRC LUAD PRAD
	ADRS	0.63	0.90	0.82	0.82	0.83
	Frobenius	0.55	0.75	0.73	0.67	0.74