Alessandro Genoni
email: alessandro.genoni@univ-lorraine.fr

Molecular Orbitals Strictly Localized on Small Molecular Fragments from X-Ray Diffraction Data

Keywords: Extremely Localized Molecular Orbitals, experimental wavefunction, electron density, transferability, quantum chemistry, X-Ray diffraction

Nowadays, the electron density is recognized as a fundamental property that contains most of the information concerning the electronic structure of molecules and, therefore, its determination from high-resolution X-Ray diffraction data is becoming more and more important.

In this context we propose a new strategy for the charge density analysis, strategy in which the chemical interpretability of the multipole model is combined with the quantum mechanical rigor of the wavefunction-based approaches. In particular, the novel technique aims at extracting Molecular Orbitals strictly localized on small molecular fragments (e.g., atoms, bonds or functional groups) from a set of measured structure factors amplitudes. Preliminary tests have shown that their determination is really straightforward and, given their reliable transferability, we envisage the possibility of constructing new Extremely Localized Molecular Orbitals databases as alternative to the existing pseudoatoms libraries.

It is well known that the electron density is of paramount importance to fully define the ground state properties of many-electron systems. Its prominence derives from the Hohenberg & Kohn theorem [START_REF] Hohenberg | Inomogeneous Electron Gas[END_REF] which has shown the existence of a functional relation between the ground state electron distribution and the ground state wave function of electronic systems. This is one of the reasons why the determination of the electron density from high-resolution X-ray diffraction measurements has become a mature field of science over the years. [START_REF] Coppens | X-Ray Charge Densities and Chemical Bonding[END_REF][START_REF] Koritsanszky | Chemical Applications of X-ray Charge-Density Analysis[END_REF] In this context, the multipole model [START_REF] Hirshfeld | Difference Densities by Least-Squares Refinement: Fumaramic Acid[END_REF][START_REF] Stewart | Electron Population Analysis with Rigid Pseudoatoms[END_REF][START_REF] Hansen | Testing Aspherical Atom Refinements on Small-Molecule Data Sets[END_REF] is by far the most widely used technique and it consists in fitting the X-Ray diffraction data to a model electron distribution written as a sum of aspherical atomic density functions (also known as pseudoatoms) depending on proper refinable parameters.

Different features have made this strategy really appealing. For instance, we can mention its linear scaling behavior and the generally high observations/parameters ratio which leads to unique and rapidly-convergent fits. Moreover, the success of this approach is also deeply associated with the fact that it allows easy chemical interpretations of the molecular electron densities. In fact, due to the intrinsic atomic partitioning in the pseudoatom definition, the charge density of a molecule can be simply seen as the sum of atomic electron distributions that are deformed by the presence of bonding and non-bonding interactions. This important peculiarity, along with the plausible idea that chemically equivalent atoms (namely, identical atoms with a similar local connectivity) provide identical contributions to the electron density, has prompted different research groups to study the transferability of pseudoatoms from small molecules to large biomolecules with identical subunits. [START_REF] Brock | Transferability of Deformation Densities Among Related Molecules: Atomic Multipole Parameters from Perylene for Improved Estimation of Molecular Vibrations in Naphthalene and Anthracene[END_REF][START_REF] Pichon-Pesme | On Building a Data Bank of Transferable Experimental Electron Density Parameters: Application to Polypeptides[END_REF] The results of these investigations have led to the construction of pseudoatoms libraries, obtained either from experimental electron densities [START_REF] Pichon-Pesme | On Building a Data Bank of Transferable Experimental Electron Density Parameters: Application to Polypeptides[END_REF][START_REF] Jelsch | Transferability of Multipole Charge-Density Parameters: Application to Very High Resolution Oligopeptide and Protein Structures[END_REF][START_REF] Zarychta | On the Application of an Experimental Multipolar Pseudo-Atom Library for Accurate Refinement of Small-Molecule and Protein Crystal Structures[END_REF] or from theoretical electron distributions, [START_REF] Koritsanszky | Aspherical-Atom Scattering Factors from Molecular Wave Functions. 1. Transferability and Conformation Dependence of Atomic Electron Densitiesof Peptides within the Multipole Formalism[END_REF][START_REF] Volkov | Ab Initio Quality Electrostatic Atomic and Molecular Properties Including Intermolecular Energies from a Transferable Theoretical Pseudoatom Databank[END_REF][START_REF] Dominiak | A Theoretical Databank of Transferable Aspherical Atoms and Its Application to Electrostatic Interaction Energy Calculations of Macromolecules[END_REF][START_REF] Dittrich | A Simple Approach to Nonspherical Electron Densities by Using Invarioms[END_REF][START_REF] Dittrich | The Invariom Model and Its Application: Refinement of D,L-Serine at Different Temperatures and Resolution[END_REF][START_REF] Dittrich | Introduction and Validation of an Invariom Database for Amino-Acid, Peptide and Protein Molecules[END_REF] that are currently used to refine charge densities of macromolecules and to predict electrostatic properties of large systems.

Nevertheless, notwithstanding its great success, the multipole model exhibits some drawbacks: a) the fitted electron density may be characterized by unphysical negative regions; b) the refinable coefficients in the multipole expansion may have a strong correlation with the Anisotropic Displacement Parameters (ADPs); c) the number of properties that can be computed directly from the fitted charge distribution is limited because the exact functional relation between the ground state electron density and the ground state wave function is generally unknown. All these problems can be solved by another family of strategies: the so called wavefunction-or orbital-based methods for the electron density determination and analysis. Pioneered by Clinton and coworkers, [START_REF] Clinton | Direct Determination of Pure-State Density Matrices. I. Some Simple Introductory Calculations[END_REF][START_REF] Clinton | Direct Determination of Pure-State Density Matrices. II. Construction of Constrained Idempotent One-Body Densities[END_REF][START_REF] Clinton | Direct Determination of Pure-State Density Matrices. III. Purely Theoretical Densities Via an Electrostatic-Virial Theorem[END_REF][START_REF] Clinton | Direct Determination of Pure-State Density Matrices. IV. Investigation of Another Constraint and Another Application of the P Equations[END_REF][START_REF] Clinton | Direct Determination of Pure-State Density Matrices. V. Constrained Eigenvalue Problems[END_REF][START_REF] Clinton | The Cusp Condition: Constraint on the Electron Density Matrix[END_REF][START_REF] Clinton | Determination of the Electron Density Matrix from X-Ray Diffraction Data[END_REF][START_REF] Clinton | Methods for Obtaining an Electron-Density Matrix from X-Ray Diffraction Data[END_REF] these techniques consist in fitting X-ray scattering data to an assumed form for the density matrix [START_REF] Clinton | Determination of the Electron Density Matrix from X-Ray Diffraction Data[END_REF][START_REF] Clinton | Methods for Obtaining an Electron-Density Matrix from X-Ray Diffraction Data[END_REF][START_REF] Frishberg | Idempotent Density Matrices for Correlated Systems from X-Ray Diffraction Structure Factors[END_REF][START_REF] Howard | Density Matrix Refinement for Molecular Crystals[END_REF][START_REF] Snyder | A Wavefunction and Energy of the Azide Ion in Potassium Azide Obtained by a Quantum-Mechanically Constrained Fit to X-Ray Diffraction Data[END_REF][START_REF] Hibbs | A New Orbital-Based Model for the Analysis of Experimental Molecular Charge Densities: an Application to (Z)-N-Methyl-C-Phenylnitrone[END_REF][START_REF] Waller | Novel Properties from Experimental Charge Densities: an Application to the Zwitterionic Neurotransmitter Taurine[END_REF] or the wave function [START_REF] Tanaka | X-Ray Analysis of Wavefunctions by the Least-Squares Method Incorporating Orthonormality. I. General Formalism[END_REF][START_REF] Jayatilaka | Wave Function for Beryllium from X-Ray Diffraction Data[END_REF][START_REF] Jayatilaka | Wavefunctions Derived from Experiment. I. Motivation and Theory[END_REF][START_REF] Grimwood | Wavefunctions Derived from Experiment. II. A Wavefunction for Oxalic Acid Dihydrate[END_REF][START_REF] Bytheway | Wavefunctions Derived from Experiment. III. Topological Analysis of Crystal Fragments[END_REF][START_REF] Bytheway | Wavefunctions Derived from Experiment. IV. Investigation of the Crystal Environment of Ammonia[END_REF][START_REF] Grimwood | Wavefunctions Derived from Experiment. V. Investigation of Electron Densities, Electrostatic Potentials, and Electron Localization Functions for Noncentrosymmetric Crystals[END_REF][START_REF] Hudák | X-Ray Constrained Unrestricted Hartree-Fock and Douglas-Kroll-Hess Wavefunctions[END_REF] of the system under investigation and provide a two-center orbital productbased representation of the electron distribution. Therefore, unlike the multipole model, the wave function-based strategies allow to obtain electron densities with intrinsic quantum mechanical features, but, unfortunately, this results into the loss of the easy chemical interpretability of the charge distribution, which is one of the strengths of the multipole model.

In this letter we present a new technique for the determination of molecular electron densities,

technique in which we combine the quantum mechanical rigor of the orbital-based approaches with the ease of chemical interpretation typical of the multipole model. In particular, the new strategy aims at extracting from X-Ray diffraction data a single Slater determinant built up with Molecular Orbitals strictly localized (namely, without tails) on molecular fragments (e.g., atoms, bonds or functional groups) and it can be seen as a combination of the experimental wave function approach proposed by Jayatilaka [START_REF] Jayatilaka | Wave Function for Beryllium from X-Ray Diffraction Data[END_REF][START_REF] Jayatilaka | Wavefunctions Derived from Experiment. I. Motivation and Theory[END_REF][START_REF] Grimwood | Wavefunctions Derived from Experiment. II. A Wavefunction for Oxalic Acid Dihydrate[END_REF][START_REF] Bytheway | Wavefunctions Derived from Experiment. III. Topological Analysis of Crystal Fragments[END_REF][START_REF] Bytheway | Wavefunctions Derived from Experiment. IV. Investigation of the Crystal Environment of Ammonia[END_REF][START_REF] Grimwood | Wavefunctions Derived from Experiment. V. Investigation of Electron Densities, Electrostatic Potentials, and Electron Localization Functions for Noncentrosymmetric Crystals[END_REF][START_REF] Hudák | X-Ray Constrained Unrestricted Hartree-Fock and Douglas-Kroll-Hess Wavefunctions[END_REF] with the method developed by Stoll for the a priori determination of Extremely Localized Molecular Orbitals (ELMOs). [START_REF] Stoll | On the Use of Local Basis Sets for Localized Molecular Orbitals[END_REF] As in the Jayatilaka method, we assume that each crystal unit does not interact with the other ones and we write the wave function for the whole crystal like this:

Y cryst = ' k Y k (1)
where all the unit wave functions Y k are formally identical and related to each other through the crystal symmetry operations.

n c k iµ ↵ o M k i µ=1
, the a-th ELMO of the fragment can be expressed as

j k ia ↵ = M k i  µ=1 C k iµ,ia |c k iµ ↵ (2)
and the ELMO wavefunction for the k-th crystal unit (assuming it as a closed-shell system) as

Y k = Âh j k 11 j k 11 . . . j k 1n 1 j k 1n 1 . . . j k f 1 j k f 1 . . . j k f n f j k f n f i ( 3
)
where  is the antisymmetrizer, n i is the number of occupied ELMOs for the i-th fragment and j k ia is a spinorbital with spatial part j k ia and spin part b .

At this stage we consider only the case in which each non-interacting unit corresponds to a symmetry-unique portion of the crystal unit-cell. This allows to write the unit-cell charge density as a sum of N m crystal unit electron densities r j (r) that are related to a reference electron distribution r 0 (r) by the unit-cell symmetry operations n R j , r j o :

r cell (r) = N m  j=1 r j (r) = N m  j=1 r 0 (R 1 j (r r j )) (4)
Eq. ( 4) is exact provided that r 0 (r) is not obtained through an isolated crystal unit calculation. Nevertheless, in our new method this condition is not fulfilled because r 0 (r) is the electron density associated with the ELMO wave function Y 0 for the reference unit, wave function which is obtained looking for those Extremely Localized Molecular Orbitals that not only minimize the ELMO energy of the reference crystal unit, but that also reproduce as much as possible the set of collected experimental structure factors amplitudes

n F exp h o .
In other words, we look for those

J h j j j i = D Y 0 Ĥ0 Y 0 E + l ⇣ c 2 D ⌘ ( 5
)
where Ĥ0 is the non-relativistic Hamiltonian operator for the reference crystal unit, l is the Lagrange multiplier associated with the constraint, D is the desired agreement between theoretical and experimental values and c 2 is the measure of the fitting accuracy between the calculated and the observed structure factor amplitudes, namely:

c 2 = 1 N r N p  h ⇣ h F calc h F exp h ⌘ 2 s 2 h ( 6
)
with N r as the number of experimental scattering data, N p as the number of adjustable parameters (in our case only the Lagrange multiplier l ), h as the tern of Miller indexes labeling the reflection, s h as the error associated with each measure and h as an overall h-independent scale factor that is determined to minimize the c 2 value. Due to the presence of experimental errors in the collected data, it does not make sense to obtain c 2 equal to zero. Therefore, an optimal value for the desired agreement D in Eq. ( 5) is one, so that the calculated structure factors are on average within one standard deviation of the experimental values. Furthermore, for the sake of completeness, it is worth mentioning that for a proper calculation of the structure factors amplitudes

n F calc h o ,
including a proper correction for the thermal vibrational effects, we use the same equations proposed by Jayatilaka and coworkers [START_REF] Jayatilaka | Wave Function for Beryllium from X-Ray Diffraction Data[END_REF][START_REF] Jayatilaka | Wavefunctions Derived from Experiment. I. Motivation and Theory[END_REF] (see Supporting Information for further details), with the only difference that for the Fourier transforms of the basis functions pairs we have implemented an Obara-Saika scheme [START_REF] Obara | Efficient Recursive Computation of Molecular Integrals over Cartesian Gaussian Functions[END_REF][START_REF] Obara | General Recurrence Formulas for Molecular Integrals over Cartesian Gaussian Functions[END_REF] exploiting both vertical and horizontal recurrence relations. [START_REF] Head-Gordon | A Method for Two-Electron Gaussian Integral and Integral Derivative Evaluation Using Recurrence Relations[END_REF] Now, for the sake of clarity, let us omit subscripts and superscripts corresponding to the reference crystal unit and, in order to search the ELMOs that minimize the functional defined by Eq.

(5), let us consider the arbitrary variation of J ⇥ j j j ⇤ with respect to the occupied ELMO j jb ↵ (from now indicated as d ( jb ) J ). Introducing the global and the local density operators, which, due to the non-orthogonality of the ELMOs, are respectively given by r

= f  i, j=1 n i  a=1 n j  b =1 ⇥ S 1 ⇤ jb ,ia |j jb ihj ia | (7) and ri = f  j=1 n i  a=1 n j  b =1 ⇥ S 1 ⇤ jb ,ia |j jb ihj ia | (8)
with S as the overlap matrix of the occupied ELMOs, and defining the structure factor operator

Îh = N m  j=1 e i2p(R j r+r j )•(Bh) = Îh,R + i Îh,C (9)
where both Îh,R and Îh,C (real and imaginary part of Îh , respectively) are hermitian, we obtain:

d ( jb ) J = 4 ⇢ ⌦ d j jb (1 r) F(1 r + r j ) j jb ↵ + + l  h K h Re n F calc h o ⌦ d j jb (1 r) Îh,R (1 r + r j ) j jb ↵ + + l  h K h Im n F calc h o ⌦ d j jb (1 r) Îh,C (1 r + r j ) j jb ↵ (10)
with

K h = 2h N r N p h F calc h F exp h s 2 h F calc h (11)
Since the lowest value of the functional is achieved if d ( jb ) J vanishes for all j and b , we have that the ELMOs that minimize J ("experimental" ELMOs) are the ones that satisfy the following equation for each fragment:

(

(1 r) F(1 r + r j ) + l  h K h Re n F calc h o (1 r) Îh,R (1 r + r j ) + + l  h K h Im n F calc h o (1 r) Îh,C (1 r + r j ) ) j jb ↵ = 0 (12)
Now, adding the following quantities to both hand sides of Eq. ( 12)

Q1 ↵ = r † j F(1 r + r j ) j jb ↵ (13) Q2 ↵ = l  h K h  Re n F calc h o r † j Îh,R (1 r + r j ) + Im n F calc h o r † j Îh,C (1 r + r j ) j jb ↵ (14)
and defining the hermitian operators

F j = (1 r + r † j ) F (1 r + r j ) (15)
Î j h,R = K h Re n F calc h o (1 r + r † j ) Îh,R (1 r + r j ) (16)
Î j h,C = K h Im n F calc h o (1 r + r † j ) Îh,C (1 r + r j ) (17)
we have

F j,exp j jb ↵ = n j  g=1 ⌦ j jg F j,exp j jb ↵ j jg ↵ (18)
where F j,exp , which is the modified Fock operator for the j-th fragment, is given by:

F j,exp = F j + l  h Πj h,R + Πj h,C (19)
Finally, applying a unitary transformation that only mixes the occupied ELMOs of each fragment among themselves and given the invariance of F j,exp to that transformation, Eq. ( 18) becomes:

F j,exp j jb ↵ = e jb j jb ↵ (20)
Therefore, the "experimental" ELMOs are obtained solving Eq. ( 20) self-consistently for each fragment using the proper local basis-set deriving from the localization scheme defined a priori. It is important to observe that, although solved separately, the equations associated with the different subunits are coupled because each F j,exp operator depends on the global density operator, as it can be seen from Eqs. ( 15), ( 16) and (17).

Unfortunately, as in the case of the "theoretical" ELMOs, convergence problems in the resolution of modified Hartree-Fock equations (in our case Eqs. ( 20)) might arise near the functional minimum [START_REF] Stoll | On the Use of Local Basis Sets for Localized Molecular Orbitals[END_REF][START_REF] Fornili | Determination of Extremely Localized Molecular Orbitals and Their Application to Quantum Mechanics/Molecular Mechanics Methods and to the Study of Intramolecular Hydrogen Bonding[END_REF] and, therefore, to overcome this drawback, we have devised an algorithm which uses the information provided by exact first and approximate second derivatives of J with respect to the ELMOs coefficients (see Supporting Information for their expressions). In particular, following a strategy successfully used by Fornili et al. for the calculation of the "theoretical" ELMOs, [START_REF] Fornili | Determination of Extremely Localized Molecular Orbitals and Their Application to Quantum Mechanics/Molecular Mechanics Methods and to the Study of Intramolecular Hydrogen Bonding[END_REF] we have adopted a quasi-Newton procedure where an approximate analytic Hessian is calculated only at the first iteration and, afterwards, it is updated exploiting the Broyden-Fletcher-Goldfarb-Shanno formula. [START_REF] Press | Numerical Recipes: the Art of Scientific Computing[END_REF] Starting from an existing code for the ELMOs calculation, both the self-consistent resolution of Eqs. ( 20) and the quasi-Newton algorithm just described above have been implemented modifying the version 8 of the GAMESS-UK package, [START_REF] Guest | The GAMESS-UK Electronic Structure Package: Algorithms, Developments and Applications[END_REF] which has been also used to perform all the other calculations that will be mentioned below.

The capabilities of the new technique have been afterwards tested considering the crystal structure of the L-cysteyne (orthorhombic phase I) determined at 30 K. [START_REF] Moggach | L-Cysteine-I at 30[END_REF] To accomplish this task we have performed single point unconstrained calculations at the Restricted Hartree-Fock (RHF), B3LYP and ELMO levels using the Dunning cc-pVDZ basis-set and the molecular geometry obtained from the X-Ray diffraction experiment (see Figure 1). In particular, for the ELMO computation we have adopted a localization scheme almost corresponding to the Lewis structure of the molecule. In fact, we have defined both atomic fragments, which describe the core electron and the lone pairs associated with each atom, and bond subunits, which describe each "electronic couple" between two nuclei, with the only exception of a three-atom fragment for the s and p electrons of the carboxylic group (see Figure 1). We have used the same geometry, basis-set and localization pattern for all the "experimental" ELMOs computations, for which we have also considered the unit cell parameters, the ADPs and the experimental structure factors amplitudes deposited by Moggach and coworkers with their paper. [START_REF] Moggach | L-Cysteine-I at 30[END_REF] It is important to observe that experimental values characterized by F exp h < 3s h have been rejected and this resulted in the selection lence electron density and, therefore, optimizing the "core" ELMOs against the experimental data, as in the ELMO-XC/OPT calculations, does not significantly improve the agreement statistics. Finally, as already mentioned above, it is meaningless to force the fitting much beyond c 2 = 1.0 and, in the following analyses, unless otherwise stated, we will refer to the ELMO-XC/OPT and ELMO-XC/FR wave functions as those constrained ELMO wave functions obtained for l = 0.45 and l = 0.5, respectively.

In Table 2 we have reported the c 2 and energy values corresponding to all the performed calculations and, as expected, the unconstrained ELMO wave function provides the worst statistics agreement with the experimental data. This agreement improves using the RHF and B3LYP methods, even if the c 2 values are still far from an acceptable experimental error range, which, as already observed, is reached through the ELMO-XC/OPT and the ELMO-XC/FRZ wave functions. Before analyzing in detail the fitting effects on the electron density, it is worth noting that 2 the energies associated with the constrained ELMO wave functions are higher than all the other ones, especially than the one corresponding to the unconstrained ELMO wave function. This is in line with what has been already observed by Jayatilaka and coworkers, [START_REF] Grimwood | Wavefunctions Derived from Experiment. II. A Wavefunction for Oxalic Acid Dihydrate[END_REF][START_REF] Hudák | X-Ray Constrained Unrestricted Hartree-Fock and Douglas-Kroll-Hess Wavefunctions[END_REF] and it can be simply explained as follows. In a "theoretical" ELMO calculation we look for those ELMOs coefficients that correspond to a minimum point on the energy hypersurface. When we extract an "experimental" ELMO wave function from X-ray scattering data, we introduce an additional constraint (see Eq. ( 5)) without providing a new variational parameter and, therefore, the coefficients of the "experimental" ELMOs are the ones that correspond to the minimum point on the functional hypersurface, which is different from the minimum energy point.

To further study the effects of the wave function fitting, we have compared the unconstrained ELMO electron density with two ELMO-XC/OPT charge distributions: the converged one and the intermediate one obtained for l = 0.20 (indicated as ELMO-XC/OPT 0.20 ), which shows how the charge density changes during the fitting procedure. Furthermore, to better distinguish the effects of the crystal environment, we have also considered the comparison of the unconstrained ELMO electron density with the RHF and B3LYP ones. From Figures 3 and4 (see also Supporting

Information for two-dimensional plots relative to other molecular planes) it is possible to note that one of the main fitting consequences consists in a large redistribution of the electron density around the nuclei, such as the evident depletion of charge in the lone-pairs region of the oxygen atoms. This decrease in electron density is common to the RHF and B3LYP methods, but in the ELMO-XC cases it is more distinct and it is also combined with a simultaneous increase in the charge distribution close to the same oxygen atoms. Another important unique effect of the fitting to the crystal environment is represented by the reduction of charge in correspondence of all the C-H and N-H bonds. A shift of charge from the nitrogen atom to the C a atom and, in analogous way, from the carbon atom to the oxygen atoms in the carboxylic group are also noteworthy and it is possible to observe that they are just sketched in the B3LYP plots. For the sake of completeness, it is important to note that in proximity of the sulfur atom the ELMO-XC charge redisributions are characterized by an increase in electron density, which is completely missing in the RHF case and which is even replaced by a charge depletion if we consider the B3LYP method.

Finally, in order to have a more quantitative picture of the fitting effects, we have computed the values of some similarity indexes between the unconstrained ELMO charge distribution and the other electron densities just mentioned above. In particular, we have considered two indexes that allow to compare point-by-point two electron distributions in the real space: the Real-Space R value (RSR), which is defined as

RSR(r x , r y ) = 100 n p  i=1 r x (r i ) r y (r i ) n p  i=1 r x (r i ) + r y (r i ) (21)
where n p is the number of electron density grid points, and the Walker-Mezey similarity index L(r x , r y , a, a 0 ), which compares two charge densities in the space bound by the isosurfaces characterized by the values a and a 0 (see Supporting Information for further details). [START_REF] Walker | Ab Initio Quality Electron Densities for Proteins: A MEDLA Approach[END_REF] For the sake of clarity, it is important to note that a complete similarity corresponds to values of RSR and L equal to zero and 100, respectively. From Table 2 it is easy to observe that, in agreement with Figures 3 and4, the electron distributions associated with the constrained ELMO wave functions are the least similar to the unconstrained ELMO charge density. Furthermore, since the Walker-Mezey index allows to study the electron densities similarity in different regions simply changing the isosurface values a and a 0 , it is important to note that the discrepancies between the constrained and unconstrained ELMO charge distributions are more pronounced far from the nuclei (see the index L(0.001, 0.01) in Table 2) than in the region close to nuclei (see the index L(0.1, 10)). [START_REF] Genoni | Optimal Virtual Orbitals to Relax Wave Functions Built Up with Transferred Extremely Localized Molecular Orbitals[END_REF][START_REF] Genoni | Extremely Localized Molecular Orbitals: Theory and Applications[END_REF][49][START_REF] Sironi | DENPOL: a New Program to Determine Electron Densities of Polypeptides Using Extremely Localized Molecular Orbitals[END_REF] the construction of "experimental" ELMOs databases can be seen as a possible alternative to the existing pseudoatoms libraries.

Acknowledgement

The author thanks Enrique Espinosa, Emmanuel Aubert, Antonio Monari and Toru Shiozaki for helpful discussions. Maurizio Sironi is also gratefully acknowledged for providing the ELMO code.

Figure 1 :

 1 Figure 1: Molecular geometry obtained from the X-ray scattering experiment. The red frame shows the three-atom fragment defined for the s and p electrons of the carboxylic group.

Figure 2 :

 2 Figure 2: Variation of the c 2 agreement statistics in function of the Lagrange multiplier l for (A) the ELMO-XC/OPT and the (B) ELMO-XC/FRZ calculations.

Figure 3 :

 3 Figure 3: Two-dimensional plots of electron density differences in the plane of the carboxylic group: (A) r(RHF)-r(ELMO), (B) r(B3LYP)-r(ELMO), (C) r(ELMO-XC/OPT 0.20 )r(ELMO), (D) r(ELMO-XC/OPT)-r(ELMO). The contours are at linear increments of 0.01 a.u. from -0.1 a.u. to 0.1 a.u., with negative contours in blue and positive contours in red.

Figure 4 :

 4 Figure 4: Three-dimensional plots of electron density differences: (A) r(RHF)-r(ELMO), (B) r(B3LYP)-r(ELMO), (C) r(ELMO-XC/OPT 0.20 )-r(ELMO), (D) r(ELMO-XC/OPT)r(ELMO). The isosurface value is set to 0.025 a.u., with negative isosurfaces in blue and positive isosurfaces in red.

Table 1 :

 1 c 2 agreement statistics and energy values for the analyzed unconstrained and constrained wave functions.

	Method	c 2	Energy (a.u.)
	ELMO	4.12	-719.10
	RHF	3.87	-719.23
	B3LYP	3.81	-721.77
	ELMO-XC/OPT	0.99	-718.84
	ELMO-XC/FRZ	0.98	-718.82
	in Table

Table 2 :

 2 Values of the similarity indexes associated with the comparison of the ELMO electron density with the RHF, B3LYP, ELMO-XC/OPT 0.20 and ELMO-XC/OPT charge distributions.To summarize, in this letter we have proposed a new wavefunction-based technique for the determination and the analysis of molecular electron densities from high-resolution X-ray scattering data. Our preliminary tests on a crystal structure of the cysteine at 30 K have shown that to obtain Molecular Orbitals strictly localized on small molecular fragments (e.g., atoms, bonds or functional groups) reproducing a set of measured structure factors amplitudes within a predefined accuracy is straightforward. The determination of Extremely Localized Molecular Orbitals from X-Ray diffraction data allows to easily analyze the molecular charge distributions in terms of usual chemical concepts (as in the multipole model) and, at the same time, to obtain quantum mechanically rigorous electron densities (as in the orbital-based approaches) that can be used to compute many important properties. Furthermore, given the reliable transferability of the Extremely Localized Molecular Orbitals,

	Comparison Vs. ELMO	RSR	L(0.001, 10) L(0.1, 10) L(0.01, 0.1) L(0.001, 0.01)
	RHF	0.55	96.85	98.85	97.55	95.80
	B3LYP	1.04	94.25	97.74	95.67	92.31
	ELMO-XC/OPT 0.20	1.65	91.30	96.23	93.89	88.08
	ELMO-XC/OPT	2.15	89.12	95.11	92.16	85.26

Supporting Information Available

Details about the computation of the structure factors, expressions of the functional first and second derivatives with respect to the ELMOs coefficients and details about the Walker-Mezey similarity index. Figure S1: two-dimensional plots of electron density differences in the N-C a -C b plane.


ESEME Vadim S Nikolayev
email: vadim.nikolayev@espci.fr

Sergey L Gavrilyuk
email: sergey.gavrilyuk@polytech.univ-mrs.fr

Henri Gouin
email: henri.gouin@univ-cezanne.fr

Modelling of the moving deformed triple contact line: influence of the fluid inertia

Keywords: Wetting, spreading, pinning, wetting hysteresis

For partial wetting, motion of the triple liquid-gas-solid contact line is influenced by heterogeneities of the solid surface. This influence can be strong in the case of inertial (e.g. oscillation) flows where the line can be pinned or move intermittently. A model that takes into account both surface defects and fluid inertia is proposed. The viscous dissipation in the bulk of the fluid is assumed to be negligible as compared to the dissipation in the vicinity of the contact line. The equations of motion and the boundary condition at the contact line are derived from Hamilton's principle. The rapid capillary rise along a vertical inhomogeneous wall is treated as an example.

Introduction

In equilibrium, the wetting properties of a liquid in contact with a solid are well defined by the static contact angle θ eq [START_REF] De Gennes | Wetting: statics and dynamics[END_REF]. In the case of complete wetting (θ eq = 0), a thin prewetting film exists in front of the bulk of the liquid and its dynamics is well understood too. For partial wetting (θ eq = 0), a difficulty arises when the triple solid-liquid-gas contact line is moving. Huh and Scriven [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF] showed that the viscous flow in contact line vicinity cannot be described by the no-slip condition, i.e. zero fluid velocity at the solid wall. Indeed, since the contact line belongs at the same time to the liquid and to the solid, the contact line velocity is ambiguous. This leads to a nonphysical divergence of the hydrodynamic pressure and of the viscous dissipation at the contact line. In reality, the dissipation is anomalously large but finite.

A motion of the straight contact line is studied in the large majority of works. However, wetting dynamics experiments at partial wetting are almost inevitably submitted to an influence of surface heterogeneities that we call "defects" for the sake of brevity. On small scale, they lead to the contact line deformation; on macroscopic scale, to the contact angle hysteresis, i.e. to the difference between the advancing and receding contact angles. The contact line elasticity approach was proposed to describe the static deformation of the contact line by a localized defect [START_REF] De Gennes | Wetting: statics and dynamics[END_REF]. It resulted in the logarithmic contact line shape, which was shown to describe correctly experimental data in the intermediate range of distances from the defect center. Obviously, the logarithmic shape fails to describe the contact line both very close to the defect and far from it where the contact line deformation needs to be finite. To obtain the correct description in the whole range, the influence of gravity [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF] or the fluid mass conservation [START_REF] Iliev | Quasi-static relaxation of arbitrarily shaped sessile drops[END_REF] needs to be taken into account.

While the static hysteresis is understood relatively well [START_REF] De Gennes | Wetting: statics and dynamics[END_REF][START_REF] Schwartz | Contact angle hysteresis on heterogeneous surfaces[END_REF][START_REF] Gouin | The wetting problem of fluids on solid surfaces. Part 1: the dynamics of contact lines and Part 2 : the contact angle hysteresis[END_REF], the contact line dynamics in the presence of defects is under active discussion, see [START_REF] Moulinet | Dissipation in the dynamics of a moving contact line: effect of the substrate disorder[END_REF] for related references. The contact line speed is often presented as a function of the dynamic contact angle θ and of the static contact angle. However, an ambiguity appears in definition of this static value which should lie somewhere between the advancing and receding values of the static contact angle.

Among theoretical studies of the dynamics in the presence of defects (that we call dynamics of the "deformed contact line") one can distinguish the papers [START_REF] Joanny | Motion of a contact line on a heterogeneous surface[END_REF][START_REF] Golestanian | Roughening transition in a moving contact line[END_REF]. Joanny & Robbins [START_REF] Joanny | Motion of a contact line on a heterogeneous surface[END_REF] proposed a general formalism for an arbitrary contact line shape but solved it only for the straight contact line. Golestanian and Raphaël [START_REF] Golestanian | Roughening transition in a moving contact line[END_REF] used the quasi-static dynamics in conjunction with the contact line elasticity model to study the random deformation of the contact line.

To describe the effect of defects, a single phenomenological parameter ξ, the dissipation coefficient, was introduced [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF][START_REF] Nikolayev | Relaxation of non-spherical sessile drops towards equilibrium[END_REF] following [START_REF] De Ruijter | Droplet spreading: Partial wetting regime revisited[END_REF]. ξ contains all the physics of the contact line motion mechanism and has to be specific to the fluid-substrate system. To introduce it into the theory, one needs to assume that (i) the viscous dissipation in the bulk of the fluid is neglected with respect to that in the contact line vicinity, and (ii) the contact line dissipation is independent of the direction of the contact line motion (advancing or receding). The energy dissipated per unit time at the contact line is then given in the lowest order of the contact line velocity v n (measured in the direction normal to the contact line) by the expression

T = ξ v 2 n 2 dl, (1)
where the integration is performed along the contact line. The defects are modelled by the spatial variation of the surface energy of the support which, according to the Young-Laplace formula, is equivalent to the spatial variation of θ eq . It was shown in [START_REF] Iliev | Quasi-static relaxation of arbitrarily shaped sessile drops[END_REF] that in quasistatic approximation (i.e. when the fluid surface is assumed to be at equilibrium shape at any time moment) the expression (1) leads to the equation valid for any contact line (or fluid surface) geometry,

v n = σ ξ (cos θ eq -cos θ). ( 2
)
where σ is the surface tension. Therefore, in any contact line motion model that results in (2), the coefficient ξ can be interpreted as the contact line dissipation per its unit length. For instance, ξ can be directly assimilated to the friction coefficient of Blake and Haynes [START_REF] Blake | Kinetics of liquid/liquid displacement[END_REF]. The dissipation coefficient theory allows the contact line dynamics and the dynamic advancing and receding angles to be predicted once the surface heterogeneity and the initial conditions are defined [START_REF] Nikolayev | Dynamics and depinning of the triple contact line in the presence of periodic surface defects[END_REF].

The relationship (2) has been verified against multiple experiments (see [START_REF] Blake | Dynamic contact angles and wetting kinetics[END_REF] for a review). Using (2), the ξ value can be obtained from them. According to various experimental data, ξ is much larger than the shear viscosity η, which corroborates the assumption (i) above. However, there is still no certainty even about the order of value of ξ/η for a particular liquid-solid system. The experimental values range from 300 [START_REF] De Ruijter | Droplet spreading: Partial wetting regime revisited[END_REF] to 10 7 [START_REF] Andrieu | Coalescence of sessile drops[END_REF] and depend on the conditions of experiment, the state of the substrate, etc.

It is argued recently [START_REF] Moulinet | Dissipation in the dynamics of a moving contact line: effect of the substrate disorder[END_REF] that v n is non-linear in cos θ which is equivalent to the statement that ξ varies strongly with v n . At large v n , the observed deviations from the behavior (2) can be explained by the model [START_REF] Blake | Dynamic contact angles and wetting kinetics[END_REF]. However, one needs to treat the experimental ξ value at small contact speeds with caution because it is usually determined [START_REF] Moulinet | Dissipation in the dynamics of a moving contact line: effect of the substrate disorder[END_REF][START_REF] Andrieu | Coalescence of sessile drops[END_REF][START_REF] Stokes | Harmonic generation as a probe of dissipation at a moving contact line[END_REF][START_REF] Hamraoui | Can a dynamic contact angle be understood in terms of a friction coefficient?[END_REF] from the slope of the dependence of v n on cos θ (where the angle brackets mean averaging along the contact line). It is shown in [START_REF] Nikolayev | Dynamics and depinning of the triple contact line in the presence of periodic surface defects[END_REF] that due to the collective effect of the defects, the dependence (2) with constant ξ can result in a highly non-linear dependence of v n on cos θ near the pinning threshold where v n is small. The inferred from the slope ξ value turns out to be larger than its actual value. The influence of defects can thus lead to averaged v n (θ) equations different from (2) at small contact line speeds. To separate the effect of defects from the physics mechanism of the contact line motion, a description of the dynamics of the contact line deformed by the defects is necessary.

Most theoretical approaches aim to describe the slow fluid motion and thus neglect the fluid inertia. However, the inertial effects can be important in many common situations. One can mention the early stages of drop spreading [START_REF] Biance | First steps in the spreading of a liquid droplet[END_REF], possibly after its impact with a support [START_REF] Roisman | Normal impact of a liquid drop on a dry surface: model for spreading and receding[END_REF], the capillary rise [START_REF] Clanet | Onset of menisci[END_REF][START_REF] Dreyer | Capillary rise of liquid between parallel plates under microgravity[END_REF][START_REF] Stange | Capillary driven flow in circular cylindrical tubes[END_REF], and oscillating flows in various geometries: those of sessile or pendant drops [START_REF] Narhe | Contact line dynamics in drop coalescence and spreading[END_REF][START_REF] Noblin | Vibrated sessile drops: Transition between pinned and mobile contact line oscillations[END_REF][START_REF] Bahr | Oscillations of sessile drops of surfractant solutions on solid substrates with differing hydrophobicity[END_REF][START_REF] Wilkes | Forced oscillations of pendant (sessile) drops[END_REF][START_REF] Lyubimov | Behavior of a drop on an oscillating solid plate[END_REF], liquid in containers [START_REF] Hocking | The damping of capillary-gravity waves at a rigid boundary[END_REF][START_REF] Ting | Boundary conditions in the vicinity of the contact line at a vertically oscillating plate: an experimental investigation[END_REF][START_REF] Jiang | Contact line dynamics and damping for oscillating free surface flows[END_REF][START_REF] Nicolás | Effects of static contact angles on inviscid gravity-capillary waves[END_REF], or capillary bridges [START_REF] Borkar | Boundary layer analysis of the dynamics of axisymmetric capillary bridges[END_REF].

In the partial wetting regime all these situations involve a period of time where v n is small and then the defect impact on dynamics is large. Indeed, v n necessarily passes through zero during oscillations. As for relaxation processes, the influence of defects becomes important at the end of the evolution.

This problem is especially important for oscillation flows where the choice of the boundary conditions at the solid walls is a longstanding problem [START_REF] Hocking | The damping of capillary-gravity waves at a rigid boundary[END_REF]. Evidently, when the amplitude of oscillations is small enough, the contact line is pinned [START_REF] Noblin | Vibrated sessile drops: Transition between pinned and mobile contact line oscillations[END_REF], so that the fixed position boundary condition need to be applied. When the amplitude is larger, the contact line can be pinned during a part of the oscillation period [START_REF] Ting | Boundary conditions in the vicinity of the contact line at a vertically oscillating plate: an experimental investigation[END_REF][START_REF] Jiang | Contact line dynamics and damping for oscillating free surface flows[END_REF]. At large amplitude, it moves almost all the time. In the relatively well studied [START_REF] Wilkes | Forced oscillations of pendant (sessile) drops[END_REF][START_REF] Nicolás | Effects of static contact angles on inviscid gravity-capillary waves[END_REF] first regime, there is no singularity at the contact line. The transition from the first to the second regime and especially the complicated dynamics in the second regime did not yet find its explanation. In this paper we develop a model suitable to solve these problems.

The equations of motion are derived in the "shallow water" approximation in sections 2 and 3. This approximation permits to solve analytically a model problem (section 4) in order to validate the model and find a first-order "inertial" correction to the quasi-static behavior. The results are summarized in the section 5.

Equations of motion

The oscillating flows are considered often in the inviscid fluid approximation [START_REF] Lyubimov | Behavior of a drop on an oscillating solid plate[END_REF][START_REF] Nicolás | Effects of static contact angles on inviscid gravity-capillary waves[END_REF]. This is justified by the thinness of the viscous boundary model at moderate Reynolds numbers where the flow is not yet turbulent. For the contact line problems, an additional argument related to the anomalously strong contact line dissipation can be put forward. Indeed the experiments [START_REF] Noblin | Vibrated sessile drops: Transition between pinned and mobile contact line oscillations[END_REF][START_REF] Hocking | The damping of capillary-gravity waves at a rigid boundary[END_REF][START_REF] Ting | Boundary conditions in the vicinity of the contact line at a vertically oscillating plate: an experimental investigation[END_REF][START_REF] Jiang | Contact line dynamics and damping for oscillating free surface flows[END_REF] show that its contribution can be larger than that of the viscosity in the liquid bulk. Therefore, we will neglect the viscous dissipation in the bulk by assuming that the liquid is inviscid. The only source of dissipation is assumed to be at the contact line.

A Cartesian reference frame Oxyz is used to describe the fluid motion in the domain x > 0, -∞ < y < ∞, 0 < z < h(t, x, y) (Fig. 1). The free surface is described by the equation z = h(t, x, y) and is assumed to be weakly deformed

|∇h| ≪ 1, (3)
where the gradient is taken with respect to the variables (x, y). The vertical wall (x = 0) may move in z direction. The gravity is directed downward. The contact line is described by the equation

z = h w (t, y) = h(t, 0, y). (4)
For convenience, the function h(t, x, y) is assumed to be periodic (with the period 2L) in the y-direction. This is not a restrictive assumption since one might put L → ∞ in the resulting equations.

To derive the governing equations and boundary conditions in the absence of the volume dissipation, we use the Hamilton principle of stationary action.

The idea is to formulate the Lagrangian, and to present the variation of the Hamilton action as a linear functional of virtual displacements allowing to find not only the governing equations but also the natural boundary conditions [START_REF] Sedov | Mechanics of Continuous Media[END_REF][START_REF] Berdichevsky | Variational principles of continuum mechanics[END_REF].

In the present work, we will use the "shallow water" approximation which will allow the analytical results to be obtained and analyzed. The Hamilton action a and the Lagrangian L (which is the difference between the kinetic and the potential energy of the system) read

a = t 2 t 1 L dt, L = 1 2L ∞ 0 dx L -L ρh | u| 2 2 - ρgh 2 2 dy - σ 2L ∞ 0 dx L -L 1 + |∇h| 2 -1 dy + σ 2L L -L dy hw(t,y) 0 c(y, z) dz. ( 5
)
Here ρ is the density of the fluid assumed incompressible, u = (u x , u y ) is the fluid velocity field. The first integral in the Lagrangian is a difference between the kinetic and the gravitational potential energy of the system written in the shallow water approximation. The second integral is the liquid-gas interface energy, and the last integral is the energy of the fluid-solid surface [START_REF] Schwartz | Contact angle hysteresis on heterogeneous surfaces[END_REF], σc(y, z) being the difference between the interfacial tension of gas-solid and liquid-solid interfaces. According to the Young expression, c(y, z) = cos(θ eq ), where θ eq is the value of the equilibrium contact angle. The variation of the latter along the substrate (i. e. along the y -z plane) reflects the surface heterogeneity (surface defects).

The variation of the action is taken under the mass conservation constraint

∂h ∂t + ∇ • (h u) = 0 (6)
and the non-penetration boundary condition at the vertical wall

u x | x=0 = 0. (7)
Search of the extremum for the action a (see Appendix A) results in the following fluid motion equation:

∂ (ρh u) ∂t + ∇ • (ρh u ⊗ u) + ρgh∇h -σh∇(∆h) = 0. ( 8
)
Its back substitution into the Lagrangian (5) reduces the latter to the form (see also (A.3))

δL = σ 2L L -L ∂h ∂x x=0 + c(y, h w ) δh w dy. ( 9
)
The dissipation is introduced [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF] in the usual way through the dissipation functional, [START_REF] De Gennes | Wetting: statics and dynamics[END_REF]. For this, we need to find an expression for the normal velocity v n of the contact line. When the vertical wall moves vertically with respect to the reference system with the velocity v, the velocity of the contact line with respect to the wall reads

v n = ∂h w ∂t + v   1 + ∂h w ∂y 2   -1/2 ≃ ∂h w ∂t + v,
where v is assumed to be positive for the downward plate motion. Since the volume dissipation inside the fluid is neglected, the dissipation functional does not depend on the bulk fluid velocity u and the dissipation account does not influence the governing equation [START_REF] Joanny | Motion of a contact line on a heterogeneous surface[END_REF].

Following [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF], one can write now the generalized Euler-Lagrange equation which takes into account the dissipation term

d dt δL δ ḣw - δL δh w = - δT δ ḣw , (10)
where δ . . . /δ . . . means functional derivative and ḣw denotes the time derivative. The substitution of ( 9) into [START_REF] Nikolayev | Relaxation of non-spherical sessile drops towards equilibrium[END_REF] results in the equation

σ ∂h ∂x x=0 + c(y, h w ) = ξ ∂h w ∂t + v (11)
that defines the contact line position h w . Eq. ( 11) plays a role of the boundary condition for the momentum equation ( 8).

Let us note that the shallow water approximation is not necessary to obtain the main result [START_REF] De Ruijter | Droplet spreading: Partial wetting regime revisited[END_REF]. It is known [START_REF] Berdichevsky | Variational principles of continuum mechanics[END_REF] that if the full Lagrangian (involving the integration over z of the kinetic energy) is used instead of its "shallow water" counterpart (5), the classical Euler equation of motion for the inviscid fluid result is obtained instead of the equation ( 8). The back substitution into the Lagrangian would lead to the same expression [START_REF] Golestanian | Roughening transition in a moving contact line[END_REF], which is necessary to obtain Eq. ( 11).

Reduction of governing equations

The boundary conditions for the governing equations (6, 8) should be completed by the condition at infinity

h| x→∞ = d, ( 12
)
where d is the constant depth. The following initial conditions will be used,

u| t=0 = 0, h| t=0 = d. (13)
The governing equations written for the deviation from the initial state ĥ = h -d can be linearized:

∂ ĥ ∂t + d∇ • u = 0, ρ ∂ u ∂t + ρg∇ ĥ -σ∇∆ ĥ = 0. ( 14
)
The velocity u is eliminated from these equations by applying the time derivative to the first of them and subtracting the divergence of the second. The result reads ρ ∂ 2 ĥ ∂t 2 -ρgd∆ ĥ + σd∆ 2 ĥ = 0. (15a)

Evidently, the first time-dependent term of this equation is inertial. Two other terms correspond to the gravity and to the surface tension contributions respectively. Without the first term, Eq. (15a) would reduce to the quasistatic linearized equation for the interface shape [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF].

The boundary and initial conditions for (15a) are obtained from (7, 12-13) by using [START_REF] Blake | Dynamic contact angles and wetting kinetics[END_REF]:

ĥ| x→∞ = 0, (15b) ∂ ∂x ρg ĥ -σ∆ ĥ x=0 = 0, ( 15c
)
ĥ| t=0 = 0, (15d)
∂ ĥ

∂t t=0 = 0, (15e)
The nonlinear equation [START_REF] De Ruijter | Droplet spreading: Partial wetting regime revisited[END_REF] rewritten for ĥ

∂ ĥw ∂t + v = σ ξ ∂ ĥ ∂x x=0 + c(y, h w ) , (15f)
where ĥw = h w -d, closes the problem statement.

Only the solution for h w (t, y) is of interest. However, it can only be obtained by analyzing the full problem (15) for h. The set of equations ( 15) allows h w to be obtained for arbitrary distribution of the surface defects c(y, z).

Solution of a model problem

To demonstrate the applicability of the developed formalism, let us solve a problem of the capillary rise on an immobile (v = 0) wall with a single stripe defect (see Fig. 3 below)

c(y) =      c d , |y| ≤ w, c s , |y| > w, (16)
where c d , c s and w are constants, i.e. for the simplest case in which c is independent of z. The values of c d and c s are the cosines of the equilibrium contact angles inside and outside the stripe defect respectively, w being its half-width. This problem has already been solved in quasi-static approximation [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF] (in which the fluid motion is neglected). We choose the same setup to show the role of the fluid inertia by comparing the quasi-static and inertial solutions.

Since z dependence in ( 16) is absent, the problem (15) becomes linear. Following [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF], h and h w are broken into parts:

h = d + h 0 (t, x) + h 1 (t, x, y), h w = d + h w0 (t) + h w1 (t, y), (17)
where h 0 and h w0 are the averaged over the y-direction values of ĥ and ĥw respectively. The problems for h 0 and h 1 are decoupled and can be solved separately.

Solution for h w0

Consider first the problem for h 0 that can be obtained from (15) by accounting for the fact ∂h 0 /∂y ≡ 0. By applying the Laplace transform

h0 (x) = ∞ 0 h 0 (t, x) exp(-st) dt, (18)
one obtains from (15a,d,e) the ordinary differential equation

d 4 h0 dx 4 -2α d 2 h0 dx 2 + β h0 = 0, (19)
where α = ρg/(2σ), β = ρs 2 /(σd). The general solution of ( 19) has a form

h0 (x) = a 1 exp(-p 1 x) + a 2 exp(-p 2 x) + a 3 exp(p 1 x) + a 4 exp(p 2 x), (20)
where

p 1,2 = α ± α 2 -β, Re p 1,2 ≥ 0. ( 21
)
The boundary conditions for h0 have the same form (15b,c) as for h 0 ; (15f) transforms into

s hw0 = σ ξ d h0 dx x=0 + c s s , ( 22
)
where hw0 is the Laplace transform of h w0 . Because of (15b),

a 3 = a 4 = 0, (23)
and (15c) results in the following condition

d h0 dx x=0 = -hw0 (p 1 + p 2 ) = -hw0 2 α + β . (24)
The final form for hw0 can be obtained by the substitution of ( 24) into ( 22):

hw0 = c s σ ξ s s + τ -1 r s √ 2τ i + 1 -1 , ( 25
)
where τ r = l c ξ/σ is a characteristic time scale for the quasi-static relaxation [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF] and τ i = l c / √ gd is the characteristic inertial time, l c = σ/(ρg) being the capillary length. The inverse Laplace transform of ( 25) can be found rigorously as described in Appendix. It reads

h w0 (t) = c s l c 1 -erfc t 2τ i + 1 2 √ 1 + r 2 exp t τ r r + √ 1 + r 2 erfc t 2τ i r + √ 1 + r 2 - exp t τ r r - √ 1 + r 2 erfc t 2τ i r - √ 1 + r 2 , ( 26
)
where erfc(•) is the complementary error function [START_REF] Bateman | Tables of Integral transforms[END_REF], and r = τ i /τ r characterizes a relative importance of the fluid inertia. It is related to the Weber number which can be obtained from (15a) where both the terms of inertia and of surface tension are present. Taking τ r and l c as the scales of time and length, one can obtain for these terms ρh/τ 2 r and σdh/l 4 c . Their ratio is the Weber number, W e = r 2 .

In the limit r → 0, [START_REF] Wilkes | Forced oscillations of pendant (sessile) drops[END_REF] reduces to the quasi-static result [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF]

h w0 (t)| r→0 = c s l c [1 -exp(-t/τ r )]. ( 27
)
The difference between these two dependencies is transparent from Fig. 2: the fluid inertia slightly slows down the contact line relaxation. For r ≤ 0.01, the difference between the functions ( 26) and ( 27) is very small so that in practice the exponential solution ( 27) is valid.

Solution for h w1

The statement of the problem for h 1 is given by ( 15) with an additional condition of zero y-averaged value of h 1 . For a single defect, this condition reduces to h 1 (y → ±∞) = 0, which allows the Fourier transform

h1 (t, x) = ∞ -∞ h 1 (t, x, y) exp(-iky) dy (28)
to be applied. The Laplace transform [START_REF] Biance | First steps in the spreading of a liquid droplet[END_REF] can further be applied to h1 (t, x) and results in h1 (x).

The further procedure is fully analogous to that developed in the previous section. It leads to the following equation for hw1

s hw1 = σ ξ d h1 dx x=0 + c1 s , ( 29
)
where c1 is the Fourier transform of c(y

) -c s , c1 = 2δc sin kw k , (30)
and δc = c d -c s . The h1 derivative is related to hw1 through

d h1 dx x=0 = -hw1 p 1 p 2 (p 1 + p 2 ) k 2 + p 1 p 2 , (31)
where p 1,2 are defined by [START_REF] Dreyer | Capillary rise of liquid between parallel plates under microgravity[END_REF], in which

α = ρg 2σ + k 2 , β = ρ σd s 2 + ρg σ k 2 + k 4 (32)
should now be used. The substitution of ( 31) into (29) leads to an expression, for which the inverse Laplace transform is difficult to find. To overcome this difficulty, one can make use of the smallness of the parameter r, i.e. assume the smallness of the contribution of the inertial effects. In this limit, (31) reduces to

d h1 dx x=0 = -hw1 k 2 + l -2 c + B(rτ r s) 2 , (33)
where the coefficient B is

B = k 2 + l -2 c |k| - k 2 + 1.5l -2 c k 2 + l -2 c . (34)
Note that B is positive and diverges at small k like B ∼ |k| -1 . The equation for hw1 becomes

hw1 = c1 s k 2 + l -2 c + l -1 c τ r s + B(rτ r s) 2 -1 . ( 35
)
This expression is rational in s and its inverse Laplace transform can be found using conventional methods:

hw1 (t) = c1    1 k 2 + l -2 c + 1 B(rτ r ) 2 (s 1 -s 2 ) exp(s 1 t) s 1 - exp(s 2 t) s 2    . (36)
Here

s 1,2 = 1 ± 1 -4Bl c r 2 k 2 l 2 c + 1 2Bl c r 2 τ r . (37)
To obtain the time-varying contact line shape h w1 (t, y), one needs to apply the inverse Fourier transform to [START_REF] Quéré | Inertial capillarity[END_REF]. It can be done numerically. The result is shown in Fig. 3.

At r → 0, one recovers the quasi-static result [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF] hw1 (t)| r→0 = c1

k 2 + l -2 c 1 -exp - t τ r k 2 l 2 c + 1 , (38)
which corresponds to a non-zero value at k → 0. However, by considering the small k limit of [START_REF] Quéré | Inertial capillarity[END_REF], one arrives to the result 1

hw1 (t)| k→0 = c1 t 2 2Bτ 2 i → 0. ( 39
)
This small k ambiguity leads to a small but regular numerical error in h w1 (t, y) when the inverse Fourier transform of ( 36) is found at very small r. In practice, the expression (38) should be used instead for r ≤ 0.05. 3. The contact line distortion h w1 calculated for r = 0.5. The defect area (half-width of which is w = 0.2l c ) with cos θ eq = c d is shadowed. The rest of the plane has cos θ eq = c s . The parameter of the curves is the time t in the τ r units. In the insert, the result of the inertial theory for r = 0.5 (solid curve) is compared to the quasi-static (r → 0) result (dotted curve). Both curves are calculated for t = 0.1τ r .

Similarly to the case of h w0 , the inertial slowing down manifests itself in the case of h w1 . The insert in Fig. 3 shows that the topmost point of the contact line rises slower when the hydrodynamic motion is taken into account. However, the hydrodynamic effects are not limited to the simple slowing down, they change the shape of the contact line deformation. In the beginning of the capillary rise, the joint effect of the inertia and fluid mass conservation creates two wells in the h 1 shape (Fig. 3) from which the fluid flows off to form a bump in the middle. Later on, these cavities become more and more shallow and wide until they disappear in the equilibrium static contact line shape [START_REF] Nikolayev | Equation of motion of the triple contact line along an inhomogeneous surface[END_REF].

Results and discussion

The main result of the present paper is a model that allowed to describe the inertial regimes of the contact line motion whenever the impact of surface defects is important in the partial wetting regime, e.g. the pinning-depinning of the contact line during its oscillation. The equations of motion were derived for the plane vertical wall geometry and can be applied for the arbitrary defect pattern and both for the spontaneous and forced contact line motion. The hydrodynamic shallow water approximation was used to obtain the solution in the closed analytical form. However, this approach can be easily generalized to any geometry and fluid depth.

Indeed, let us have a closer look to the equation [START_REF] De Ruijter | Droplet spreading: Partial wetting regime revisited[END_REF]. Its most interesting feature is its invariance. In the hydrodynamic theory, it has exactly the same form as in the quasi-static theory [START_REF] Nikolayev | Dynamics and depinning of the triple contact line in the presence of periodic surface defects[END_REF]. Since c(y, z) = cos[θ eq (y, z)] and ∂h/∂x| x=0 ≃ -cos θ under the assumption (3), one notices that [START_REF] De Ruijter | Droplet spreading: Partial wetting regime revisited[END_REF] is nothing else than the equation [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF]. It was derived under the assumption of the small deformation of the straight contact line. By analogy with the quasistatic result [START_REF] Iliev | Quasi-static relaxation of arbitrarily shaped sessile drops[END_REF], it is however reasonable to assume that (2) is valid for any contact line deformation. Any contact line inertial problem can thus be solved by applying the Euler set of equations for the inviscid fluid with the boundary condition [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF] where θ eq (y, z) models the surface heterogeneity.

The asymptotics of the capillary rise at small times which follows from the present model can be compared to the available experimental data. It stems from the equation ( 2) that v n is constant at t = 0 because θ(t = 0) is fixed by the initial condition for h. Since h(t = 0) = 0 has been chosen, θ(t = 0) = π/2 and h(t) ∼ t at small times independently of the geometry of the problem. This asymptotics agrees with the results [START_REF] Dreyer | Capillary rise of liquid between parallel plates under microgravity[END_REF] on the capillary rise between flat walls. The asymptotics h(t) ∼ t 3/2 obtained for the complete wetting case [START_REF] Clanet | Onset of menisci[END_REF][START_REF] Quéré | Inertial capillarity[END_REF] cannot be compared to our model restricted to the partial wetting. As a matter of fact, due to existence of the prewetting film [START_REF] De Gennes | Wetting: statics and dynamics[END_REF], the contact line dissipation anomaly is much weaker for the complete wetting. The kinetics of the capillary rise is thus defined by the balance of the surface tension and inertia [START_REF] Clanet | Onset of menisci[END_REF] which results in h(t) ∼ t 3/2 . The effect of inertia on the contact line relaxation is two-fold. First, it is the slowing down of the average contact line motion. When the contact line speed is inhomogeneous, the inertia amplifies its variation along the contact line.

The impact of the inertia on the contact line relaxation can be measured by the value r = τ i /τ r of the ratio of two characteristic times. The first of them is the inertial time and the second is the quasi-static relaxation time. The quasistatic result is recovered in the limit r → 0. Note that the Weber number (the ratio of the inertial and surface tension terms) is related to r: W e = r 2 .

The simplicity of the developed approach is its main advantage. It uses a single phenomenological parameter (the contact line dissipation coefficient ξ) and allows complicated three dimensional problems to be treated. In this article we analyzed the contact line motion in presence of a stripe defect. This model can be directly verified against experiments [START_REF] Nadkarni | An investigation of microscopic aspects of contact angle hysteresis: Pinning of the contact line on a single defect[END_REF][START_REF] Marsh | Dynamics of contact line depinning from a single defect[END_REF][START_REF] Paterson | Wetting on heterogeneous surfaces: Experiments in an imperfect Hele-Shaw cell[END_REF]. However, the information presented in these papers is insufficient for a direct comparison and additional measurements are necessary.

A Derivation of the governing equation ( 8) Assumption (3) permits us to rewrite (5) in a simpler form

L = 1 2L ∞ 0 dx L -L ρh | u| 2 2 - ρgh 2 2 dy - σ 2L ∞ 0 dx L -L |∇h| 2 2 dy + σ 2L L -L dy hw(t,y) 0 c(y, z) dz (A.1)
The variation of the Lagrangian (A.1) submitted to the constraint ( 6) is taken in several steps. First, we introduce a smooth one-parameter family of virtual motions x = Φ(t, X, ε), x = (x, y) ( X stands for the Lagrangian coordinates, ε is a small parameter in the vicinity of zero and x = ϕ(t, X) = Φ(t, X, 0) is the real motion). Then, we define the virtual displacements as functions of (t, X):

δ x(t, X) = ∂ Φ(t, X, ε) ∂ε ε=0 .
Later, we consider the virtual displacements as functions of (t, x) due to the possibility to invert X through X = ϕ -1 (t, x). The variations of the free surface h(t, x) and the velocity field u(t, x) compatible with the mass conservation law [START_REF] Gouin | The wetting problem of fluids on solid surfaces. Part 1: the dynamics of contact lines and Part 2 : the contact angle hysteresis[END_REF] are given in the motion space (t, x) by (see, for example, [START_REF] Berdichevsky | Variational principles of continuum mechanics[END_REF] or [START_REF] Gavrilyuk | Hyperbolic models of homogeneous two-fluid mixtures[END_REF]) The variations δ x and δh w are independent and the term containing δ x is absent in the dissipation function because of neglect of the viscosity in the bulk of the fluid. Therefore, the integral containing δ x in (A.3) has to be equal to zero, which results in the governing equation [START_REF] Joanny | Motion of a contact line on a heterogeneous surface[END_REF].

δh
B Laplace-inversion of (25)

The expression [START_REF] Bahr | Oscillations of sessile drops of surfractant solutions on solid substrates with differing hydrophobicity[END_REF] for which the inverse Laplace transform is sought. The inverse of

f (s -m 2 ) = 1 2m 2 b 1/2 √ s -m + 1/2 √ s + m - u 2 (u 2 -u 1 ) -1 √ s -u 1 - u 1 (u 1 -u 2 ) -1 √ s -u 2 , (B.2) where u 1,2 = -b ± √ b 2 + m 2 , (B.
3) is found using the inverse of ( √ s + m) -1 listed in the tables [START_REF] Bateman | Tables of Integral transforms[END_REF]. To find the inverse of f (s), the original of f (s -m 2 ) should be multiplied by exp(-m 2 t) according to the rules of the Laplace transform.

Fig. 1 .

 1 Fig. 1. Sketch of the liquid layer near the vertical (possibly moving with the speed v) wall with surface defects.

Fig. 2 .

 2 Fig. 2. The time dependence of the y-averaged hight of the contact line rise h w0 . The exponential quasi-static law (27) (solid line) is compared to the inertial behavior (26) (dotted line) calculated for r = 0.5.

  Fig.3. The contact line distortion h w1 calculated for r = 0.5. The defect area (half-width of which is w = 0.2l c ) with cos θ eq = c d is shadowed. The rest of the plane has cos θ eq = c s . The parameter of the curves is the time t in the τ r units. In the insert, the result of the inertial theory for r = 0.5 (solid curve) is compared to the quasi-static (r → 0) result (dotted curve). Both curves are calculated for t = 0.1τ r .

  (t, x) = -∇ • (hδ x), δ u(t, x) time derivative. Finally, the expression of the variation of a in terms of the virtual displacements reads , h w )δh w dy. By integrating by parts, accounting for the periodicity in y-direction, and the boundary condition[START_REF] Moulinet | Dissipation in the dynamics of a moving contact line: effect of the substrate disorder[END_REF] in the form δ x • n| x=0 = 0, where n = (1, 0) is the unit normal vector to the plane x = 0(ρh u ⊗ u) -h∇(ρgh) + σh∇(∆h) • δ xdy + ∇ • (hδ x)∇h] dy + σ , y, h w )δh w dy, where ⊗ denotes the direct product of the vectors. The definition of δh w δh w = -∇ • (hδ x)| x=0 gives us the final expression for the variation of the Hamilton action: u) ∂t -∇ • (ρh u ⊗ u) -h∇(ρgh) +σh∇(∆h) • δ x dy + σ • n + c(y, h w )]δh w dy. (A.3)

  has a form f (s) = [s(s + 2b √ s + m 2 )] -1 , (B.1)

Notice that although s 1,2 become complex at small k, hw1 remains real.

Acknowledgements

VN is indebted to D. Beysens for stimulating discussions.