Yann Chalopin
email: yann.chalopin@ecp.fr

Marc Hayoun

Sebastian Volz

Hichem Dammak

Surface enhanced Infrared Absorption in dielectric thin films with ultra-strong confinement effects

a By formulating a microscopic description of the non-local dielectric constant, we have investigated the mechanisms of infrared absorption in dielectrics thin films by molecular dynamics simulations. We found that light absorption in dielectric slabs does not occur predominantly at the polaritons resonances but through anomalous surface modes extremely confined in space. This demonstrate that any macroscopic description of electrodynamics in dielectrics breaks down at the nanoscale.

The study of the optical properties of nano-materials involves two main classes of physical systems: metals and dielectrics. At the nanoscale, both exhibit particularly interesting properties due to the resonances activated by the confining boundary conditions -e.g surface phononpolariton for dielectrics and surface plasmon for metals. Understanding the microscopic mechanisms of light coupling with nanomaterials forms the cornerstone of the fields of nanophotonics and nano-optics from which a large panel of applications have recently emerged such as the enhancements of Raman scattering 1 , the fluorescence imaging 2 or the control of coherent thermal radiation. When studying how size-confinement impacts optical properties at the short-scales, it has been reported that non-local effects take place in metals 4,5 such as the spill-out of electrons 6,7 : For systems smaller than 10nm, a strong increase of the intrinsic electron-phonon interaction can be observed. On the other hand, very few work reported peculiar size effects in dielectric materials 9 although there has been much effort devoted to exploring optical modes of vibration in these systems 12,21 .

In this paper, the microscopic mechanisms of IR-light absorption in a model (MgO) of confined dielectric materials of few atomic layers thickness are investigated. In these systems, it is notably demonstrated that the optical response in the THz regime exhibits anomalous size effects indicating that any macroscopic description of infra-red (IR) light/matter interaction breaks down at the atomic scale. A microscopic non-local model that allows the recovery of the frequency versus wave-vector polarizability spectra is presented. It is based on the tracking of the equilibrium polarization fluctuations. In this study, the IR absorption spectra for MgO slabs thickness ranging from 0.4 to 4nm which are shown to be size dependent. Ultra strong confined modes are unveiled and their physical origin are discussed.

To this aim, a microscopic formalism that captures the non-local dielectric constant is developed. It is adequate to equilibrium molecular dynamics (MD) simulations. Our approach is formulated on the basis of the linear response theory. The susceptibility χ of a system subjected to an external electric field E can be expressed as

p β (ω, k) = V ε 0 χ βγ (ω, k) E γ (ω, k),
with p the dipolar moment induced by E. The indexes β and γ stand for the space coordinates (x, y, z) and V is the volume of the dielectric system. For any anisotropic materials, the fluctuation-dissipation theorem 13 allows to generally express the tensor χ as 14

χ βγ (ω, k) = 1 V ε 0 k B T [ p β (0, k)p γ (0, k) + iω ∞ 0 e iωt p β (0, k)p γ (t (1
) where p is expressed in the reciprocal space as

p β (t, k) = j q j • r jβ(t) e ik•r 0 j , (2)
r 0 j is the equilibrium position of atom j caring the charge q j . Hence, the optical spectra of any dielectric material are extracted from the fluctuations of the system at thermal equilibrium and the non-local description of their dielectric properties corresponds to the Fourier transform of the charges fluctuations. For a bulk material, the dielectric constant ε is complex, ε = ε ′ + iε ′′ , and one has:

ε ′ = Re{χ} + 1 (3) ε ′′ = Im{χ} (4)
For a dielectric slab:

V χ = α = α ′ + iα ′′ (5)
α is the polarizability. MD simulations were performed by using an empirical potential 15 which has been demonstrated to be valid for bulk and surface properties 16,17 . We have first studied the MgO bulk with a crystal of 8000 atoms which corresponds to a 10 × 10 × 10 cubic cell ( 10 3 points in the first Brillouin zone) with periodic boundary conditions.

We have validated our approach by recovering the dielectric constant at k = 0 to allow a fair comparison with the available experimental data 18 (see Fig. 1). We have calculated the non-local imaginary part of the dielectric constant, ǫ ′′ (ω, k), with k ranging along the high symmetry lines of the first Brillouin zone (Fig. 1). From these two plots, we can see that the bulk absorption at the long wavelength limit corresponds to the TO mode at the frequency of f T O = 12 THz. In contrast, we can see that for k = 0, the polar LO (25THz) is absent (Fig. 1.a).

We now consider the case of a crystal slab of finite thickness L = 3a in the z direction [001] and periodic in the two lateral ones x,y ([100] and [010]). The in-plane modes are defined by the wave-vector components (k x , k y ) and the cross plane one by k z as α βγ (k x , k y , k z , ω). The absorption at long wavelength limit (k x = k y = k z = 0) is given on Figure 2 a). It can be seen that IR absorption occurs through mainly two bands at 12 and 25 THz.

To elucidate the physical origin of these IR resonances which are not observed in the bulk, the spatial mapping of the absorption processes has been performed at each resonant frequency (Figure 2 middle). To this aim, each atom i of charge q i has been subjected to a polarized electric force q i E and the corresponding local absorption, a i , has been traced back with the local dipolar moment

q i x i as a ,⊥ i (ω) = ω k B T ∞ 0 q 2 i x ,⊥ i (0)x ,⊥ i (t) e iωt dt, (6) where x
,⊥ i refers to the displacement of atom i in the cross plane (⊥) or in the in-plane directions ( ) of the slab.

The lowest frequency band at 12 THz corresponds to in-plane polarization ([100], [010]) whereas absorption at the highest frequency band 25 THz corresponds to cross plane polarization ([001]). It can be seen that the polar mode at the frequency f 1 involves surface atoms whereas f 2 and f 3 corresponds to volume modes. This clearly indicates that in this slab and more generally suggests that in dielectric thin films, IR light is absorbed at the bulk TO and LO 24 frequencies in addition to a surface mode. Interestingly, the absorption occurs predominantly through this surface mode (f 1 ) which is completely localized at the topmost layers. The modes f 2 and f 3 are the transverse and longitudinal polaritons, respectively. The f 1 and f 3 modes originate from the finite boundary conditions of the slab. The appearance of the f 3 mode is a clear evidence of the Berreman effect 25 . We think that a signature of this mode can be possibly observed experimentally [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF] even while measuring the bulk dielectric constant (see Figure 1 at 28 THz ) as the mea-surements always involve a dielectric/air interface at a surface of the crystal sample.

The non-local polarizability calculated in the first Brillouin zone (Fig. 2 bottom) shows that the LO and TO polaritons have dispersion relations similar to those of the bulk LO and TO branches. On another hand, it can be seen that the surface mode (f 1 ) is non-dispersive along the k z direction. In contrast to the bulk non local properties, extra polar modes are activated at lowest frequencies from the K point to the Γ point.

To investigate the size effects, the optical spectra have been studied at different film thicknesses ranging from L = 1a to 10a. The corresponding absorption spectra are presented on top of Figure 3. For the thinner thickness (L = 1a ; two planes), it is observed two peaks at frequencies of 11 and 23 THZ, respectively. These two peaks correspond to polaritons modes involving the displacements of all the atoms parallelly and perpendicularly to the slab, respectively. Let us first consider the size dependence of the 11 THz peak : When increasing the layer thickness from L=2a, we observe a splitting of this mode into two contributions identified as f 1 and f 2 . Such a splitting occurs because the topmost atomic layers undergo interatomic force constants weaker than that of the volume layers where the atoms are isotropically bonded. Below L=5a, the absorption by the surface mode, f 1 , predominates over the absorption of the TO polariton mode, f 2 . Recalling that while increasing L, the number of surface atoms get decreased with respect to the volume one, it can be noticed that even for a ratio of 20 % at L = 5a, the absorption due to the topmost layers remains as important as that of the volume ones.

We observe two additional peaks between 20 and 25 THz with very much smaller amplitudes. Their associated frequencies and amplitudes are size dependent. They progressively disappear at the thicker thicknesses. By mapping the local polarisability (Equ.6) for L = 3a, we found that the absorption peak around 21 THz involves the atoms of the planes underneath the outermost planes of the slab (not displayed). The peak around 23-25 THz originates from the atoms of the third deepest surface planes. Thus, these two absorptions are very localized. The frequencies indicates that the physical origin of the additional peaks does not correspond to the modes with a wavelength λ that fits the thickness L as L = nλ/2 11,21 . Indeed, the number of modes does not increase as the layer thickness increases. We attribute this effect to the dielectric screening by the topmost layers [START_REF] Chalopin | [END_REF] .

By looking at the size dependence of the non-local polarisability (fig. 3-bottom) in the [001] direction (i.e along the slab thickness), we see that the TO absorption modes in the range 11-15 THz corresponds to a competition between the polariton branch and the surface branch strongly localized at the topmost atomic layers. The surface branch remains non-dispersive whereas the dispersion of the polariton branch tends to that observed for the bulk (Fig. 1 bottom) as the thickness increases. In a similar fashion, the LO branch at higher frequen-cies (20-25THz) corresponds to the polariton mode and its dispersion converges to the bulk one at thicker thicknesses. In addition, its intensity (i.e the polarisability) is weaker than that of the TO modes for the various thicknesses. These discrepancies can be attributed to the effect of dielectric screening by the surface layers as it has been observed for the case of MgO nanoparticles of few nm size 9 .

To conclude, we have demonstrated that IR optical absorption in dielectric slabs of few nanometers is driven by strongly localized surface processes which are not predicted by a macroscopic description of electrodynamics. These peculiar resonances originating from the vibrational dynamics at the surfaces of finite thickness slabs enhance the radiation absorption in an important way and are predominant among other usual absorption processes (polaritons). To unveil the physical origin of these effects, we have presented a microscopic non-local model based on molecular dynamics simulations and which provides a convenient framework for investigating optical properties of matter at the nanoscale.

FIG. 1 .FIG. 2 .

 12 FIG. 1. Top: frequency dependent bulk dielectric constant ((ǫxx + ǫyy + ǫzz)/3 ) compared to experimental data 18 ; Bottom: the corresponding non-local spectra (T r[ǫ(ω, k]) in log scale) along the main directions of the MgO first Brillouin zone.

FIG. 3 .

 3 FIG. 3. Size dependence of the imaginary part of the polarizability in MgO slabs ranging from L = 1a to L = 10a. Top: Absorption of the slabs (k = 0); Bottom: α ′′ (ω, kz) in log scale. a a Applied Physics Letter, 104(2014)011905. DOI: 10.1063/1.4860989


Grégory Geneste

Hichem Dammak

Marc Hayoun

Mickael Thiercelin

Low-temperature anharmonicity of barium titanate: a path-integral molecular dynamics study

We investigate the influence of quantum effects on the dielectric and piezoelectric properties of barium titanate in its (low-temperature) rhombohedral phase, and show the strongly anharmonic character of this system even at low temperature. For this purpose, we perform path-integral molecular-dynamics simulations under fixed pressure and fixed temperature, using an efficient Langevin thermostat-barostat, and an effective hamiltonian derived from first-principles calculations. The quantum fluctuations are shown to significantly enhance the static dielectric susceptibility (≈ by a factor 2) and the piezoelectric constants, reflecting the strong anharmonicity of this ferroelectric system even at very low temperature. The slow temperatureevolution of the dielectric properties observed below ≈ 100 K is attributed (i) to zero-point energy contributions and (ii) to harmonic behavior if quantum effects are

turned off.

I. INTRODUCTION

It is a common result of solid state physics that, below its Debye temperature, θ D , a solid exhibits a behavior that significantly deviates from the predictions of classical mechanics, because of quantum fluctuations associated to atomic motions: in such conditions, one or several vibration modes of the system do not behave classically, i.e. the energy quantum hω separating their eigenstates is significantly higher than the thermal energy k B T . In harmonic, or mildly anharmonic systems, quantum effects can be theoretically treated in the framework of the harmonic or quasiharmonic approximation, using the standard quantization of atomic vibrations through the introduction of the normal coordinates, and many solids can be considered as harmonic below a certain temperature. Such harmonic systems usually exhibit weak dependence of their dielectric properties with temperature, and their dielectric permittivity is not sensitive to quantum effects 1 . However, the influence of quantum effects is more complex in systems in which the microscopic degrees of freedom move inside a strongly anharmonic energy landscape. For example, in a number of crystals exhibiting a form of ferroic order related to atomic displacements (ferroelectric, ferroelastic, antiferrodistortive, ferrotorroidic, etc), most of the degrees of freedom typically evolve throughout a multiple-well energy surface. If the system is considered as classical (quantum effects neglected), there always exists a temperature below which such system exhibits harmonic behavior (it is finally trapped in one single well of the energy surface), but in a quantum system, if the quantum fluctuations on the displacements in the ground state extend beyond the harmonic region associated to each well, the standard treatment of quantum effects by means of the harmonic approximation becomes irrelevant.

In such cases, where quantum fluctuations of the ground state cause the system to probe the anharmonic part of the potential, deviations from the harmonic approximation are expected down to zero K.

Ferroelectric systems are typically characterized by complex multiple-well energy landscapes affecting the polar degrees of freedom, and the interplay between quantum fluctuations and polar modes at low temperature is strongly system-dependent. In the case of a deep multiple-well surface, the effects might be limited to weak quantum delocalization. But Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 in the case of a shallow multiple-well surface, dramatic consequences can be observed, up to strong tunneling effect that can be able to fight against the order parameter, and even make it disappear. One of the most spectacular interaction between quantum fluctuations and polar degrees of freedom is indeed the suppression of the ferroelectric polarization due to quantum zero point motions in the so-called "quantum paraelectric" materials such as SrTiO 3 2,3 and KTaO 3 , 4 two systems in which the existence of an underlying ferroelectric instability (thus associated to a shallow multiple-well energy surface) is experimentally infered from the very large and saturated values of the dielectric permittivity just above zero K, and suggested from density-functional theory (DFT) calculations in the case of SrTiO 3 5 . Similar effects are suggested in another family of compounds, the "high-temperature" quantum paraelectrics, such as CaTiO 3 6,7 , La 1/2 Na 1/2 TiO 3 8,9 or in a more general way RE 1/2 Na 1/2 TiO 3 , 10 in which RE features a rare-earth element. However, even in "standard" ferroelectric crystals such as BaTiO 3 (BTO), the quantum effects have been shown by path-integral Monte-Carlo technique to significantly decrease the phase transition temperatures by ≈ 30-50 K 3 and to strongly modify the shape of the pressure-temperature phase diagram 11 up to room temperature. Quantum effects can thus strongly influence the structural and dielectric properties of ferroelectric (FE) systems, not only at low temperature but also at room temperature and beyond.

Anharmonicities in BTO are impacting the physics even at low temperature, as shown by its lowest phase transition (rhombohedral-orthorhombic -expt: 183 K 12 , ≈ 190 K with the hamiltonian used in the present work [START_REF] Zhong | [END_REF] ). In a simple picture, each phase transition in BTO corresponds, upon heating, to the temperature at which the local modes (dipoles) get out of the potential well(s) in which they were confined, and come to visit new potential energy minima, giving rise to a new value and direction of the macroscopic polarization.

In other words, phase transitions are the points at which anharmonicities in the potential energy surface strongly manifest and impact the physics. This simple microscopic picture of BTO has been theorized more than 40 years ago through the well-known "eight-site" order-disorder model of Comes et al. 15 . In this model, the local dipoles evolve among eight off-center satellite sites located along the <111>-type directions. In the paraelectric phase, all the sites are visited with equal probability (resulting in a zero polarization), while in the FE phases, only a subset of these sites (the same at all cells due to strong intersite correlations) is visited, giving rise to a non-zero polarization. Modern calculations using an Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113

effective hamiltonian fitted on first-principles calculations [START_REF] Zhong | [END_REF] show that, in the paraelectric phase, the local density of probability is in fact quasi-isotropic (with slight maxima along the <111> directions) since the dipoles also spend a significant part of their time between these off-center sites 16,17 . At any rate, this simple model confirms the strongly anharmonic character of this system, at least above its first phase transition and thus, the complete impossibility to describe its phase transitions in the harmonic approximation. Let us mention that the signature of anharmonicities is seen also below the first phase transition: in the rhombohedral phase, order-disorder mechanims associated with local reversal of the dipoles have been observed, both experimentally 18 and by calculations, 16,17 suggesting anharmonic behavior also in the ground state phase of BTO.

On the other hand, the Debye temperature of BTO is commonly placed about 150 K above room temperature (θ D ≈ 480 K 14 ), showing that the dipole dynamics should be impacted by quantum effects up to such high temperature. Therefore, since the structural phase transitions occur well below θ D , subtle interplay between quantum fluctuations and anharmonicities are expected in BTO. Let us also mention the fact that the so-called "temperature rescaling" method 19 does not apply to BTO since the effective temperature calculated from the phonon density of states obtained at T = 0 K in the rhombohedral phase is about 250 K, which falls above the two first phase transitions of BTO, namely in the tetragonal phase.

In the present work, we investigate the anharmonicities of barium titanate at low temperature (in its rhombohedral phase), that manifest through the quantum fluctuations of its ground state. We describe the influence of quantum effects on the order parameter, on the dielectric permittivity and on the piezoelectric constants of BTO, by using path-integral molecular-dynamics simulations.

II. COMPUTATIONAL DETAILS AND THEORETICAL BACKGROUND

A. Path-integral formalism

The quantum effects related to atomic motions are accounted for by using the path integral (PI) formalism. In this formulation of quantum statistical mechanics, the canonical partition function Z is written as a discretized imaginary time path integral. For a quantum system containing N (discernable) particles of mass m, Z can be expressed according to:

Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113

Z = lim P →∞ ( 2πmP k B T h 2 ) 3NP/2 × U (1)
...

U (P )
e -βV ef f ({ U (1) ... U (P ) }) d U (1) ...d U (P ) .

(1)

This integral involves P (Trotter number) replicas of the system labeled by the integer s, each replica consisting of a set of N positions

U (s) = ( u (s) 1 ... u (s)
N ). These replicas characterize the discretization of the PI in imaginary time (imaginary time slices). β is the statistical temperature, β = 1 k B T . The effective potential V ef f ({ U (1) ... U (P ) }) couples the positions of the P slices through:

V ef f ({ U (1) ... U (P ) }) = P s=1 [ N i=1 1 2 k( u (s) i -u (s+1) i ) 2 + 1 P Φ({ u (s) i })]. ( 2
)
The harmonic term of this effective potential involves a spring constant k

= mP k 2 B T 2 h2
. Φ is the physical potential energy computed inside each imaginary time slice s. Each particle (i, s) is thus interacting by harmonic forces with the particles (i, s + 1) and (i, s -1), forming a ring that closes onto itself ( u

(P +1) i = u (1) i ).
In the limit of infinite Trotter number, Eq. (1) tends to a functional integral (imaginarytime path integral): that respectively depend on the momenta d U dτ and positions U (τ ). The discretized expansion of Eq. 1 is at the root of a formal analogy 20 ("classical isomorphism" 21 ) between any quantum system and an equivalent classical system made of P images of the initial set of N particles, because the multidimensional integral of Eq. 1 can be viewed as the canonical partition function of this equivalent classical system. Each quantum particle is associated to a ring polymer of P classical particles, these classical particles interacting with each other through the "true" physical forces (divided by P ) inside each slice, and through harmonic forces (acting between each particle of slice (s) and the corresponding particles of slice (s -1) and (s + 1)). Each set of harmonic interactions is assumed to close onto itself, forming a closed ring (P +1 → 1). In the limit where the Trotter number P → ∞, this equivalent classical system has exactly the same partition function as Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 that of the quantum system under study. The extension to the isothermal-isobaric ensemble (NPT) is straightforward 22 .

Z = D Ue -1 h R βh 0 [T ( d U dτ )+Φ( U (τ ))]dτ , (3)
As a consequence, classical simulation techniques such as Monte Carlo (MC) or molecular dynamics (MD) can be applied to the classical equivalent to estimate numerically the thermodynamic properties of the quantum system. The corresponding methods are respectively called path-integral Monte Carlo and path-integral molecular dynamics (PIMD). The estimated properties should be converged with the Trotter number, and if this condition is fullfilled, it is possible to compute thermodynamic quantities that exactly include all the quantum dispersion effects. However, the path-integral expansion of Eq. 1 assumes distinguishible particles. The physical properties computed by the present PI formalism thus do not include exchange between particles, which are assumed to obey Boltzmann statistics 23 .

Practical application of PIMD raises technical problems, related to the fact that for high Trotter number, the forces acting inside the classical equivalent are mainly harmonic. In such harmonic system, obtaining ergodic trajectories by using the standard algorithms of MD such as Nosé-Hoover thermostat is difficult 25 . In order to recover ergodicity, two different approaches can be employed:

(i) a deterministic approach using efficient schemes based on thermostat chains 23,24,25 ;

(ii) a stochastic approach based on the use of the Langevin thermostat 26,27,28,29 .

In the present case, we have found very efficient to use the Langevin dynamics, which is extremely powerful to produce an ergodic exploration of phase space by introducing at each time step a random force that mimics the "noise" observed in the motion of a brownian particle. We have successfully tested the scheme on simple systems (1D and 3D harmonic oscillator, double well potential, quartic potential) and found an excellent agreement between PIMD and the exact result (obtained by analytic formulae or by a numerical solution of the Schrödinger equation).

Another practical difficulty of PIMD is the existence of modes evolving on very different time scales 24 . To circumvent this difficulty, a specific coordinate transformation ("normal mode" or "staging") that diagonalizes the harmonic parts of the PIMD forces could be employed. An appropriate choice of fictitious masses gives the same time scale for the dynamics of each degree of freedom. Instead of using such a coordinate transformation, we performed long enough trajectories (see Sec. II D).

Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113

B. Hamiltonian

We use the effective hamiltonian of Zhong et al. [START_REF] Zhong | [END_REF] , which is derived from first-principles density-functional calculations and has been shown to provide an excellent description of the thermodynamics of BTO, especially its complex sequence of phase transitions: rhombohedral(R) -orthorhombic(O) -tetragonal(T) -cubic(C) -and the (first) order of its phase transitions. In particular, although the Curie temperature T c is predicted at ≈ 300 K, i.e. about 100 K too low with this hamiltonian, it provides a good value for the R-O phase transition temperature (≈ 190 K without inclusion of quantum effects [START_REF] Zhong | [END_REF] , and ≈ 160 K after inclusion of quantum effects, to be compared with the experimental value of 183 K [START_REF] Mitsui | Landolt-Bornstein Numerical Data and Functional Relationships in Science and Technology[END_REF] ). The degrees of freedom of this hamiltonian are the local modes { u i }, the mechanical displacement modes { v i } and the (homogeneous) strain tensor {η l }. u i is, roughly speaking, the local polar displacement inside cell i, related to the local dipolar moment d i through an effective charge

Z * : d i = Z * u i .
For simplicity, the mechanical displacement modes -that allow appearance of an inhomogeneous component of the total strain and have been shown to be of very weak influence in this material [START_REF] Zhong | [END_REF] -are not accounted for in the present study. The hamiltonian Φ({ u i }, {η l }) is thus a function of the local modes { u i }, and of the (homogeneous) strain tensor components {η l }. It consists of a (local) "onsite" part, a long-range dipole-dipole interaction term, a term describing short-range interactions between neighboring local modes (up to 3 rd neighbor), a (local) term that couples the local mode to the strain and an elastic energy [START_REF] Zhong | [END_REF] . In the presence of an external electric field E, a termi Z * u i . E is added.

In this work, we use in the discussion the mean local mode < u> as order parameter (this is a displacement). An unambiguous relationship with the spontaneous polarization P can be made by P = Z * < u> Ω , where Ω is the unit cell volume and Z * the effective charge associated with the local mode (9.956 e). The components of < u> are expressed in a 0 units, where a 0 is the theoretical lattice constant of BTO [START_REF] Zhong | [END_REF] Since Langevin dynamics is very efficient to recover ergodicity within the PI formalism in the canonical ensemble, we wish to use this method, not only under fixed temperature but also under fixed pressure. The extension of the Langevin method to the isothermal-isobaric ensemble has precisely been achieved by Quigley and Probert 30,31 , giving rise to an algorithm in which random and friction forces are applied, not only on the atomic coordinates, but also on the supercell vectors. We have thus implemented this "Langevin barostat" within the PI formalism. In what follows, second-rank tensors are written in bold. The equations of motion on local mode i of slice (s) (with mass m) using the Langevin barostat are:

d p (s) i dt = f (s) i -γ p (s) i + L (s) i - p G W g p (s) i - 1 N f . T r(p G ) W g p (s) i , (4) with f (s) i = -1 P ∇ u (s) i Φ( u (s) 1 , ..., u (s) N ) -k(T, P )(2 u (s) i -u (s+1) i -u (s-1) i
) the PIMD force that includes the quantum kinetic energy contribution, which takes the form of an harmonic force with spring constant k(T, P ) = mP k 2 B T 2 /h 2 . The term -γ p

(s) i dt = p (s) i m + p G W g u (s) i , (5)
while the matrix of the supercell vectors h and its conjugate momentum p G evolve according to

dh dt = p G h W g , (6)
and

dp G dt = V (t)(X -P ext Id) + 1 N f i,s p (s) 2 i m Id -γ G p G + L G , ( 7
)
in which V (t) is the supercell volume (that evolves with time), W g is the "mass" associated to the barostat, N f is the number of degrees of freedom, P ext is the external pressure, Id is the identity tensor and X is the internal pressure tensor 30,32 . This random force on the barostat is symmetrized at each time step to avoid global rotation of the supercell during the simulation. Let us precise that the same supercell h is common to all the imaginary time slices (there is no replication of the supercell).

The Langevin algorithm is very sensitive to the quality of the random number generator.

In this study, we have adopted the routines of R. Chandler and P. Northrup 33, 34 . Finally, the algorithm is very stable and thermalization is fully achieved within ≈ 20000 MD time steps.

D. Molecular dynamics

The MD simulations are performed using a 12×12×12 supercell with periodic boundary conditions. The time step is δt=1.0×10 Equilibrium trajectories of typically 3×10 5 (low T) up to 5×10 5 (high T) steps are generated at each temperature after an equilibration time that can be long at low temperature, due to a very slow dynamics. Thus, at low temperature (T ≤ 60 K), an electric field along [111] is applied during the equilibration procedure to help the system reach faster its equilibrium state. The friction coefficient of the Langevin thermostat is γ = 0.5 THz.

Dielectric and piezoelectric tensors have been computed using a finite differences method, by directly applying a finite electric field along [001]. Equilibrated trajectories of PIMD steps are generated under static electric field and fixed pressure (-4.8 GPa). Different values of the electric field, in the range [-5, +5] × 10 6 V/m, are chosen for each temperature. The field is weak enough to induce a linear dielectric response, in both the polarization and the strain.

In order to ensure the sampling accuracy, we compare our PIMD averages with the one obtained by PIMD under the staging mode transformation. Since the difficulty to obtain such a good sampling could occur with high Trotter number, let us consider the most unfavourable case (P = 64 at T = 30K). A long NP T trajectory has thus been followed by a ≈ 100000 step-NV T trajectory performed by fixing the strain tensor components to the mean values deduced from the NP T trajectory. This second trajectory was computed by using both the primitive coordinates, and the staging coordinates. This allows to demonstrate that long trajectories, even using primitive coordinates, provide very well converged results, as shown in Tab. I. It is important to mention that the primitive averages converge to the staging ones after ≈ 20000-30000 times steps, whereas our NP T trajectories contain more than 300000 time steps.

III. SPONTANEOUS POLARIZATION

As a first step, we investigate the phase sequence of BTO by both classical and quantum simulations in order (i) to evaluate the importance of the quantum contributions as a function of temperature, and (ii) to examine whether the low-temperature spontaneous polarization is impacted by quantum fluctuations, providing a first idea of low-temperature anharmonicities.

Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 NV T trajectory with staging coordinates. The strain fixed for the NV T cases is the one averaged along the NP T trajectory. In the classical case (P = 1), the temperature evolution of polarization is displayed on In the quantum case, the transition temperatures (see Fig. 1) are lowered compared to the classical case by ≈ 10-15%, and agree very well to the ones obtain by path-integral Monte

E x = E y = E z =
Carlo by Zhong and Vanderbilt 3 . This behavior is expected since the quantum fluctuations destabilize the ferroelectric order. We also observe an important decrease in the spontaneous polarization in the three ferroelectric phases (Fig. 1). At T = 30 K, we find in the quantum case u x = u y = u z = 0.0216 a 0 (0.221 C/m 2 ), whereas in the classical case, we have 0.0257 a 0 (0.263 C/m 2 ). This difference reflects the fact that the mean polar displacement in the ground state is not at the minimum of the potential energy. This is due to the quantum fluctuations extending to the asymmetric (and thus anharmonic) region of the potential energy surface of the system. Indeed, an harmonic system would have its ground state symmetric with respect to the minimum energy point. Despite the multidimensional character of this energy surface, it is possible to have an idea of its anharmonicities by plotting the energy as a function of polarization P , with all the local modes fixed to the same value in all the unit cells of the system: We now examine the influence of quantum fluctuations. In both the classical and quantum-mechanical cases, our calculations show that the dielectric constants slowly evolve with temperature for temperatures below ≈ 100 K (Fig. 3). This common feature is explained differently according to the case. In the classical approach, the system, as explained previously, eventually becomes harmonic at sufficiently low temperature, leading to the slow temperature-evolution in the dielectric response. In the quantum case, the system rapidly reaches its ground state upon cooling, generating a freezing (saturation) of all the physical quantities. The slow evolution of the dielectric response at low temperature is therefore attributed to the quantum zero-point effects. For higher temperatures (> 100 K) the dielectric constants increase and become very large as approaching from below the R-O phase transition. We systematically find that the inclusion of quantum fluctuations enhances the dielectric response, approximately by a factor 2 for both 11 and 33 . Such differences between the quantum and classical descriptions are the signatures of a strong anharmonicity in the potential energy landscape of the rhombohedral phase of barium titanate, since in an harmonic system, the static dielectric tensor would not depend on the inclusion of quantum effects and would be independent on the temperature 1 .

u 1 = u 2 = ... = u N = u.
Finally, we have also performed measurements of the dielectric constant, * 33 , on a [001] oriented single crystal in the 10-300 K temperature range by using an impedance analyzer.

These results are given in Fig. 4 Hence, in the R phase of BTO the direction of enhanced piezoelectricity is different from the direction of the polarization. This feature was previously studied in giant-piezoelectric single crystals such as PMN-PT 41,42 or PZN-PT 43 . Such piezoelectric properties in [001] oriented crystal has been attributed to the very large value of the d 24 shear coefficient. This property apparently specific to morphotropic compounds is thus observed in a simple perovskite like BTO. Indeed, in the case of BTO the d 24 /d 33 ratio is around 21 and close to the ratio of 22 found for PMN-PT in its R phase.

VI. DISCUSSION

It is useful to compare qualitatively the case of BTO to other ferroelectric systems. Thus we now discuss the importance of quantum fluctuations and their interplay with anharmonicities in such systems. We limit the discussion to standard FE systems, that do not exhibit any other order parameter than polarization, and to their lowest-temperature phase.

Let us denote by E zp the zero-point energy associated to atomic motions, and by ΔV the typical (free) energy barrier to overcome to reverse the polarization (as would be given for instance by a phenomenological treatment within Landau theory). Different prototypical cases can be described.

(i) FE crystals with very deep double-well free energy landscape exhibit large barrier height (ΔV ) compared to E zp and thus a Curie temperature T c relatively high. The typical example of such systems is PbTiO 3 . Upon cooling, for T << T c , the crystal eventually behaves as an harmonic system. Therefore, its dielectric permittivity is expected to be rather flat at low temperature, and not sensitive to quantum effects. state is expected to extend over the different energy minima. In such cases, the ferroelectric transition might be suppressed by quantum zero-point fluctuations and the system remains paraelectric for all temperatures, despite the existence of a FE instability. This is typically the case of the quantum paraelectric crystals, such as KTaO 3 . 4 In such systems, the dielectric permittivity increases upon cooling but, since the phase transition does not take place, it eventually saturates to a large value down to zero K.

(iii) An intermediate situation corresponds to a zero-point energy E zp < ΔV , but large enough anyway so that the ground state extends up to the anharmonic region of the potential (while staying confined in a single energy minimum). In such a case, all the physical quantities exhibit therefore anharmonic behavior down to T = 0 K. The dielectric permittivity, in particular, saturates at low temperature, as a property of the ground state. In the present work we have demonstrated that BTO belongs to this last case, and is thus an anharmonic system down to zero Kelvin. It is important to point out that in the case of the classical treatment, harmonic behavior is observed at sufficiently low T , leading here again to a plateau in the low-temperature evolution of the dielectric permittivity. However, the difference of origin between the quantum and classical plateaus is reflected by their different values (the quantum plateau is higher than the classical one).

VII. CONCLUSION

In this work, the PIMD simulations have been performed to study low-temperature (R phase) dielectric and piezoelectric properties of BTO: spontaneous polarization, dielectric susceptibility and piezoelectric constants. We have shown that these three properties are different whether they are treated classically or quantum-mechanically. More precisely, significant enhancement of the dielectric tensor components and of the piezoelectric constants are observed as a consequence of the inclusion of quantum fluctuations. Such enhancement is attributed to the anharmonic contributions to the ground state, in which the system saturates at low temperature.

By contrast, without including the quantum effects, the system eventually becomes harmonic at low temperature. In BTO, the polarization quantum-fluctuations in the ground state therefore extend over a region in which the potential energy surface is strongly an-Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 harmonic. This is corroborated by the significant difference between the low-temperature spontaneous polarizations computed classically and quantum-mechanically.

BTO is a strongly anharmonic system down to zero Kelvin, the anharmonicity being retained at low temperature by the quantum zero-point effects. The anharmonicity should thus be accounted for to achieve realistic predictions in this system, even in its low-temperature rhombohedral phase.

APPENDIX A: PIEZOELECTRIC-COEFFICIENT CALCULATIONS

We choose the electric field E and the stresses σ as independent variables. The electromechanical properties are therefore described by the dielectric constants at constant stress σ ij , and the piezoelectric constants d iα . Since there is no ambiguity, the superscripts σ will be omitted for simplicity in the following.

  the integral being over all paths [τ ∈ [0; βh] → U (τ ) ∈ R 3N ] with the cyclic condition [ U (0) = U (βh)], T and Φ in the Euclidean action being the kinetic and potential energies,

  computed in the framework of the localdensity approximation to DFT: a 0 = 7.46 Bohrs. For each temperature, the macroscopic order parameter < u> and the homogeneous strain are obtained by averaging over unit cells, (real) time steps, and imaginary time slices, allowing to determine the symmetry of the phase (R, O, T or C). Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 C. Langevin barostat BTO is simulated under fixed (hydrostatic) pressure and fixed temperature conditions.

  friction force of the Langevin thermostat and L (s) i is the so-called Langevin force, which is randomly drawn at each time step in a gaussian of variance 2γmk B T δt , δt being the time step. The momentum p

  -15 s. The mass associated to the local mode (not important in the classical case for the computation of ensemble averages, but crucial in the quantum case) is 39.0 atomic mass units, determined from the force constant matrix eigenvector defining the local mode. This value is the same as that used in Ref.35. The external pressure is fixed at -4.8 GPa, as in Ref.13, a negative value that compensates the underestimation of the lattice constant within the local-density approximation to DFT. The Langevin-PIMD equations of motion are integrated within the Verlet algorithm: at each step, the new positions at time t + δt are computed from the positions at t and tδt, and also from the velocities at t, as required by the equations of motion of the Langevin thermostat-barostat. These velocities are calculated self-consistently (a short internal loop is performed at each step of the Verlet algorithm) starting from an estimation of the velocity taken from Ref.36.Convergence with the number of imaginary time slices has been carefully studied. Our tests show that various properties (polarization, strain and phase transition temperatures) of the system from T = 120 to 300 K do not exhibit significant change from P = 8 to P = 16. Thus the computation of the dielectric and piezoelectric tensors is achieved at low temperature (from T = 30K to T = 137K) by maintaining P × T = 120 × 16 = 1920, leading to the use of Trotter numbers as large as P = 64 at the lowest temperature studied Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 (T = 30 K). Note that classical mechanics is recovered by setting the Trotter number P to 1.

Fig. 1 .

 1 Fig.1. We find the expected sequence of phase transitions: R -O -T -C, as experimentally observed[START_REF] Lines | Principles and Applications of Ferroelectrics and Related Materials[END_REF] and in excellent agreement with the classical Monte carlo calculations of Zhong et. al.[START_REF] Zhong | [END_REF] using the same hamiltonian. In particular, the transition temperatures are very close to those obtained by these classical Monte Carlo simulations[START_REF] Zhong | [END_REF] : the Curie temperature is found at ≈ 295 K, while the R-O and O-T transition points are obtained at ≈ 190 K and ≈ 230 K, respectively. These preliminary calculations illustrate the accuracy of the Langevin barostat to sample the (NP T ) ensemble. Note that these transition points are subject to uncertainties of ≈ 10-15 K related to hysteretic phenomena common to the simulation of first-order phase transitions. A more accurate determination of the phase transition temperatures would require to achieve a finite-size scaling analysis[START_REF] Challa | [END_REF] which is beyond the

  This simplified static energy landscape is displayed in Fig.2along three different directions and help understanding the effect of quantum fluctuations: since the energy surface is below (resp. above) the harmonic approximation along the [111] direction when the polarization is decreased (resp. increased) with respect to its value in the energy minimum, the quantum fluctuations decrease the value of the spontaneous polarization, as found in the calculations. This significant difference between the classical and quantum systems is a signature of the importance of anharmonicities in the rhombohedral phase of BTO. In other words, quantum fluctuations lower the spontaneous polarization at low temperature by 15-20%. This manifestation of quantum mechanics takes place through the strong anharmonicity of the potential energy surface of BTO.IV. DIELECTRIC PERMITTIVITYWe now focus on the rhombohedral phase of BTO below 140 K and compute the dielectric and piezoelectric tensors in the presence of quantum fluctuations. This calculation is achieved, as explained in the computational part, by applying a finite external electric field E along the [001] direction. Fig.3shows the temperature evolution of the transverse Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 and longitudinal components of the static dielectric tensor, computed classically (P =1) and quantum-mechanically (P × T = 1920). These components are systematically given, and discussed, in the rhombohedral reference system (see Appendix).First, we note that, independently from the inclusion of quantum fluctuations, the transverse component, 11 , is much larger than the longitudinal component, 33 . This behavior is expected and in good agreement with the density functional perturbation theory (DFPT) calculations of Wu, Vanderbilt and Hamann 39 , who obtained 11 = 265 and 33 = 50 at T = 0 K. The ratio between these two components can be explained through the curvature of the potential energy surface around the rhombohedral minimum. Fig.2shows how the energy landscape varies as a function of the polarization around this minimum. The[ 110] and[ 11 2] directions give insight into 11 , the [111] direction into 33 . The much sharper increase observed along the longitudinal [111] direction (Fig.2) shows that 33 is smaller (smaller polarization fluctuations) than 11 along the transverse ones.

  Figures 5 and 6 display the temperature evolution of the longitudinal, d 33 , tranverse, d 31 , and shear, d 22 and d 24 , piezoelectric constants in the R phase of BTO, computed in the classical and quantum-mechanical descriptions. The same types of effects are observed as for the dielectric constants. We note that, independently from the inclusion of quantum fluctuations, d 24 > d 22 > d 33 > d 31 and the shear piezoelectric constants are one order of magnitude larger than the longitudinal and transverse constants. The values obtained at T = 30 K are in good agreement with the DFPT calculations of Wu, Vanderbilt and Hamann 39 (performed at T = 0 K). Note that the comparison with T = 0 K results, reported in Tab. II, is relevant since the properties slowly evolve with temperature below 100 K. The components of the piezoelectric tensor have not been experimentally determined in the R phase, to the best of our knowledge. Nevertheless, the longitudinal constant along [001], d * 33 , was derived from the slope of the strain versus electric field curves measured at T = 173 K 40 (≈ 10 K below the experimental R-O transition temperature). d * 33 (see Appendix) calculated at T = 137 K (≈ 20 K below the R-O transition temperature) is rather close to this experimental value (see Tab. II).It is important to mention that the longitudinal coefficient d * 33 , along [001], is higher than d 33 , along[111]. This strong anisotropy is due to the contribution of the super large shear piezoelectric constants, d 22 and d 24 . Hence, it is interesting to calculate the orientation dependence of d * 33 in order to identify the direction along which its value is the more enhanced. Figure7shows the orientation dependence of d * 33 in the ( 110) plane along which the maximum value is obtained. It corresponds to a direction close to [001] and this behavior is clearly temperature-independent.

d 24 -2d 22 -d 22 d 22 0 d 24 0 0 d 31 d 31 d

 2231 The dielectric and piezoelectric matrix, ij and d iα , of the rhombohedral single domain state using the [110], [112] and[111] axis and according to the R3m symmetry is as below 33

FIG. 1 :FIG. 2 :FIG. 3 :FIG. 4 :FIG. 5 :FIG. 6 :FIG. 7 :

 1234567 FIG. 1: Temperature evolution of the mean local mode components (proportional to the polarization) as obtained by standard MD and PIMD.

  . In the right member of Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113Eq. 7, one recognizes a friction force on the supercell -γ G p G (γ G is a friction coefficient for the barostat) and a random force L G , a 3 × 3 matrix whose components are randomly drawn at each time step in a gaussian with variance 2γ G Wgk B T δt

TABLE I :

 I Averaged local modes and stress at T = 30 K and for P = 64 obtained over (i) a NP T trajectory with primitive coordinates, (ii) a NV T trajectory with primitive coordinates, (iii) a

TABLE II :

 II Piezoelectic constants (pC/N) of barium titanate obtained by calculations. T (K) d 31 d 33 d 22 d 24 d *

			33 [001]
	present work 30 7.5 11 75 236 137
	DFPT 39	0 6.8 15 70 243 137
	exp. 40	173	275
	present work 137	315 ± 20

  Physical Review B 87, (2013) 014113. DOI: 10.1103/PhysRevB.87.014113 (ii) At the opposite, if E zp > ΔV or if ΔV and E zp have similar values, the ground

  * ij and d * iα obtained by axis rotation, i. e. in the referential [100], [010] and [001], are as below

		⎡				⎤
		⎢ ⎢ ⎢ ⎣	* 11 * 31	* 31 * 11	* 31 * 31	⎥ ⎥ ⎥ ⎦
			* 31	* 31	* 11
	⎡
	⎢ ⎢ ⎢ ⎣	d * 33 d * 31 d * 31 d * 36 d * 34 d * 34 d * 31 d * 33 d * 31 d 34 d * 36 d * 34
		d * 31 d * 31 d * 33 d * 34 d * 34 d * 36

We consider, at a given temperature, the rhombohedral phase with a spontaneous polarization, P 0 where P 0 1 = P 0 2 = P 0 3 and spontaneous homogenous strains, η 0 α where η 0 1 = η 0 2 = η 0 3 and η 0 4 = η 0 5 = η 0 6 . When an electric field E is applied along the [001] direction, the polarization and homogenous strain variations are given by

From equilibrium trajectories for different values of the electric field, it is easy to derive the six dielectric and piezoelectric constants, ij and d iα , by using equations (A1-A8). * Corresponding author: gregory.geneste@cea.fr 1 In the harmonic approximation, the static dielectric constant writes

where τ runs over the vibration modes at the Γ point, S (τ ) is the mode-oscillator strength tensor of mode τ , Ω is the unit cell volume and ∞ αβ is the electronic dielectric tensor. See for instance: X. Gonze, C. Lee, Phys. Rev. B 55, 10355 (1997).


Pascale Henner

Christophe Schwartz

Eric Lichtfouse
email: eric.lichtfouse@dijon.inra.fr

Éric Lichtfouse

Pipette Pasteur extraction: a fast, cheap, convenient, exhaustive and environmentally friendly method for the extraction of solid samples

come     Pipette Pasteur extraction : a fast, convenient, exhaustive and environmentally friendly method for the extraction of solid samples

INTRODUCTION

Although extraction methods are rarely described in details in the literature, it is one the key steps of relevant biogeochemical investigations. Soxhlet extraction usually give high yields [START_REF] Lindhardt | Comparison of soxhlet and shake extraction of polycyclic aromatic hydrocarbons from coal tar polluted soils sampled in the field[END_REF][START_REF] Collister | An isotopic biogeochemical study of the Green River oil shale[END_REF], but since this system is operated with a boiling solvent during long time periods, organic compounds may be either removed by evaporation or degraded by oxidation, photooxidation and mineral catalysis. Shake extraction at low temperature [START_REF] Radke | Removal of soluble organic matter from rock samples with a flow-through extraction cell[END_REF][START_REF] Lichtfouse | A molecular and isotopic study of the organic matter from the Paris basin, France[END_REF] should lower the degradation of organic molecules, but losses by evaporation are still possible, especially for low molecular weight substances. Ideally, those methods should be conducted under inert atmosphere in closed systems, but from E.L. experience in several laboratories worldwide this is rarely done. Furthermore, mixing organic compounds with minerals using either high energy appartus, e.g. ultrasonic [START_REF] Lichtfouse | Molecular, 13 C, and 14 C evidence for the allochthonous and ancient origin of C 16 -C 18 n-alkanes in modern soils[END_REF], microwave [START_REF] Budzinski | Extraction assistée par chauffage micro-ondes focalisées (MOF) à pression ambiante des composés organiques dans les matrices naturelles : application à l'analyse des composés aromatiques[END_REF] , or even low temperature, e.g. shake extraction, should favor organic/organic and organo/mineral reactions. Indeed, clay-catalysed reactions can occur under very mild conditions. For example, the treatment at room temperature for 1 h. of choslet-5-ene in cyclohexane in presence of montmorillonite clay gives quantitatively rearranged cholest-13(17)-enes (Sieskind and Albrecht, 1985). Under similar conditions, cholesterol reacts with benzene to give phenylcholestenes with 50% yield in 30 min (Siekind and Albrecht, 1993). Therefore, ideally, the following conditions should be lowered for extraction : O 2 contact, high temperature, high energy appartus (sonication, micro-wave), extraction time, transfert and concentration steps, costly materials and solvent volume (health hazard).

Here, we report a simple method which fulfil the following requirements : it is fast to shorten contact with minerals and atmosphere; yields are higher than shake and ultrasonic extraction; the absence transfert steps and the extraction at room temperature minimizes degradation and losses by evaporation; the use of few small materials and solvent volumes makes it cheap, environnementally friendly, and reduces possible contamination of the sample, e.g. by plasticizers and atmospheric fallout. It requires less than 1 m 2 of laboratory space under a ventilated hood and may be also operated in field studies.

EXPERIMENTAL Material preparation

Silica gel 60 (Merck, particle size 0.063-0.200 mm, 70-230 mesh ASTM) and cotton wool are washed three times with methylene chloride either by soxhlet extraction or by shake extraction in a beaker. The cotton wool is allowed to dry over an aluminium foil under a ventilated hood overnight. The silica gel is allowed to dry 2 days under a ventilated hood then reactivated in a beaker in a ventilated oven 4 h at 60°C, 4 h at 90°C, 4 h at 110°C, then kept in a dessicator with blue silica gel dessicant. Methylene chloride (CH 2 Cl 2 ) is either purchased at high purity (GC, HPLC grades) or distillated. Pasteur pipette columns are made by cutting carefully the narrow end of 150 mm Pasteur pipettes (length 95 mm, top external diameter 7 mm, see Figure 1). A small piece of cotton wool is introduced into the pipette and push firmly on the bottom. The pipette is washed with CH 2 Cl 2 (3 x 500 µl) using a 1 ml syringe, then kept in a clean box until extraction. Screw top vials (Varian) are washed twice with 1 ml methylene chloride.

Pasteur pipette extraction

The extraction appartus is drawn on Figure 1. Typically, 5 kg of soil from an ancient gaswork are well-mixed. 50 g are dried under a ventilated hood at room temperature, sieved to 2 mm then finely grounded. 300 mg are weighted in a 10 ml beaker then mixed with 300 mg of silicagel to easy solvent elution in later extraction. A small amount of silica gel is loaded onto the modified pipette Pasteur column above the cotton wool (silica height ~ 3 mm) in order to prevent fine mineral particles from falling into the 2 ml vial during extraction. The soil/silica gel mixture is loaded using a glass funnel (external diameters : 8 mm, 57 mm) plugged in a teflon tube (length 2 cm, diameter 8 mm). The column is gently tapped with a fingernail (3x) then loaded with 1000 µl CH 2 Cl 2 using a 1 ml glass syringue (Hamilton # 1001). Elution is allowed until the solvent surface reaches the sol/silicagel top. A low pressure is applied using a pear-shaped rubber (taken from a soap gas flowmeter) filled with cotton wool to prevent rubber particles from falling into the Pasteur pipette column, and fitted with a 5 cm teflon tube (length 5 cm, diameter 8 mm). Care should be taken to apply a low pressure even when removing the pear after elution, otherwise vacuum may pull the soil/silica gel upward. A dark-brown colored extract typical of heavy fossil fuel products is observed. The elution is repeated twice with 500 µl CH 2 Cl 2 : the second load is allowed to elute until the last drop of CH 2 Cl 2 (colorless). The elution time is about 5 min. The elution volume (1439 µl) is calculated using extract and solvent weight (1.9 g) and CH 2 Cl 2 density (1.32). At this point, the sample may be directly analysed by gas chromatography coupled to mass spectrometry in order to prevent losses of low-molecular weight substances by evaporation. Internal standards may also be introduced for further analytical steps. The mixture in the 2 ml vial is allowed to concentrate to dryness under a ventilated hood, then rapidly weighted. Alternatively, a low pressure of nitrogen may be applied provided that care is taken to stop the N 2 flow when dryness is reached, otherwise low molecular weight substances will be rapidly lost by evaporation.

Shake extraction

The same solvent volume/soil weight ratio was used in all extractions : pipette Pasteur, shake, ultrasonic. 1000 mg of the same soil were extracted three times 10 min. at room temperature with 2220 µl CH 2 Cl 2 in a 10 ml beaker using a magnetic stirrer. The supernatants were filtered using a glass funnel and a filter paper, collected in a 50 ml round-bottom flask, then concentrated under reduced pressure at 20°C to near dryness. The extract was then transfered with 4 x 400 µl CH 2 Cl 2 into a 2 ml vial using two 1 ml syringes (one for the sample, one for the solvent), then allowed to concentrated as described above.

Ultrasonic extraction 1000 mg of the same soil were extracted 10 min. three times with 2220 µl CH 2 Cl 2 in a 10 ml beaker using an ultrasonic bath. The solution is then treated as in the shake extraction (see above).

RESULTS AND DISCUSSION

Nine samples (3 replicates) of the same soil were extracted giving 3.0, 3.0 and 3.0 mg/g averaging at 3.0 mg/g for the Pasteur pipette method, 1.8, 1.7 and 1.8 mg/g averaging at 1.8 mg/g for the shake method, and 1.9, 1.9 and 1.9 mg/g averaging at 1.9 mg/g for the ultrasonic method. Obviously, the Pasteur pipette method shows the highest yields. Several explanations may explain this difference. First, shake and ultrasonic extraction are prone to losses by evaporation because of the vigorous mixing and of the higher contact surface with atmosphere during extraction and filtration. Second, the additional step of transfert should also increase experimental handling losses. Third, reactions such as degradation and polymerisation are more likely to occur due to atmosphere contact (O 2 ) during stirring with minerals. Whatever the case, the pipette Pasteur is clearly more exhaustive, faster, more convenient, cheaper, and may be applied to any kind of solid samples such as soil, sediment, plant, food and biological samples. Moreover, this method using a small amount of sample (300 mg) is particulary well suited for the preparation of samples for current analytical tools such as gas chromatography coupled to mass spectrometry because detection levels are now very low, e.g. in the range of 1 pg to 1 ng.

Figure 1 .

 1 Figure 1. Apparatus for Pasteur pipette extraction of solid samples. Lengths in mm


Florent Mousque
email: jlore111.111ousque@e11siace1jt-abstracf

Marianne Boix

Stéphane Égny"

Ludovic Montastruc'

Louis Genty'

Serge Domenech''

Optimal on-grid hybrid po\-ver system for eco industrial parks planning and influence of geographical position

Keywords: eco-industrial parks, renewable energies, multi-objective optimization, environmental assessment

Nowadays, many Eco-lndustrial Parks (Efi)) are emerging ail around the world. In these eco-induslrial parks, industries gaLher lo share several sLreams like energies, water, resourccs, and wastcs in ordcr to rcalisc cconomic bencfits whilc lowcring their environmental impact. However, most of the time, these EIP are historically constituted and therefore sub-optimal. To really observe tbe sougbt gains, exchange net w orks need to be optimally designed. This article aims to develop a methodology for a grassroots design of hybrid power system (HPS) based on mu Ili-objective optimization to supply an EIP. Overhang of this sludy is multi-objective (MO) approach which aims to minimize energy cost for project lifespan and to minimize environmental impact of energy production LO supply ETP demand. This lasl criterion is evaluated by a life cycle assessment (LCA) on each sources. Another advantage is that geographical innuence on HPS design is taken into accom1t to evaluate climate effect on RE production and LCA impact depending on selected position for the EIP. A case study compares single objective and MO optimization and geograpbical intlue□ce on designed HPS for the same ETP, in two countrics, Malaysia and in German y. Rcsults shows that a big reduclion on environmental impact (above 70%) lhanks to a little invest (cost increased about 10%) can be reached. Also, compared wilh different geographic positions, dcsigned solution is really contrasted according to cost, environmental impact and sclcctcd RE sources.

Introduction

ln recent years, al climate conference of the UNO (COP23 in 2017) important measures have bccn takcn to fight against climatc change and global wanning. and to rcducc human footprint on his cnvironmcnl. ln this uim, to rcducc extraction of ruw matcrials and to limit amount of generated wastes, industrial ecology proposes to industries to gather within eco-industrial parks (ErP) 10 exchange energies and rcsources. Thereby, the main opportunity for compc111ics is to rcducc thcir cnvironmcnlal impact whilc incrcasing thcir compctitivcncss by makmg cconomic gains. ln tbcir rcv1cw [START_REF] Boix | 01,1imization me-1hods a1>1>liéd l<) lhc design of ero--induS-lrial parks: A litc1'::tlure rcview[END_REF] have analysed methods applied to design and to optimize networks and streams in EIP. J>rcvious s-tudies have becn muinly focuscd on \VIiier and matcrials oetworks, whilc there is already a little number of 1 >ublic.ation dealing with energy networks for EJ'f>. Howcvcr 10 supply cm. :rgy dcmand whilc roducing its cnvirOmTiental impact an El P could produce ils own local encrgy by using a renewablc. encrgy hybrid power systt '. lll (1-IPS), An HPS ,s an e.nergy neiwork where several 1"C11e\,;ablc energy sources are The ê1im of Lhis paper is to <lcvt.:lop an on-grid HPS modc1 to optimizc cnergy planning to suppl y an El f>. ? , ;fQ optimization mcthod is driven by both antagonist critcria: economic and environmental. To me.}Sure geog.raphical impac. t of the location, Lhe same HPS is optimiz<.• d in diffcrcnt countrics (J , .falaysia. G<. • nnany) using the samc EIP demand, reg.ardless ofHPS technologies. Moreover, this model for the grassroots design of an EU) allows 10 Oplimizc eoc.i: gy .Slôrage size, aod .1: eoe\v(1bk� eoergy powtr rating: (for wind Turbines, photovoltaic panels and biomass. combustion) in orde.r to provide known energy de,nand of EIP's processes.

Methodology and model formulation

The (irst S(cp of tbis mclh<.H.lulogy is Lo deline a mo<lcl rcpn. : �en1ing HPS. To this cm. l, u su 1 >e1'Slruct11re ha� beeo developOO. in whic:h are repr-esented links betwee.n sourùt?S and power demand (cach company owns a nmnbcr of pI• oce::-s unîts rcpresonting sinks). Scheme of this superstructure is illustrated in figure 1. The EIP consis1s in 4 i11dus1ries composed of AC or OC proccss 10 supply. While it îs a theoretical El P, energy d,mumds arc thosc of real proccss industries, in order to rcprcscm a complcx dcmand, with many variations ovcr tirnc. lndccd in 1hc alm't!' IÎ1lium casting industry works as a batch production (most of cncrgy dcmand is bctwccn 7AM and 6PM) white the other industries works continuously. The cement plant is the biggcst consumer (abovc 75% or the consumption) duc to proccsses like raw-mill driver, ccmcnt mill driver and convcyor system. Mor<.. -ovcr, rcnewable cncrgy production ca,1 be s1ored in OC cncrgics st0ragc which is lcad 4 acid bancries technology duc to thcir bc-.:;.t pricc ratio to makc profitable s1oragc use.

The main mathematic:al constraints are energy balances at power source which is supplied dircctly to the 1>0wcr demand or charged into storage. This cons1r. aint is availablc for each powc1• source and for cach type of clec1rical currem. Co11cerniog the dcmands. the 1>owcr demand must be satislied by cqui\ 1 alent power of sources or baueries 1nuhiplied by 1>01ential energy loss if il has to be convened or discharged. Any loss corresponds to chnrging cfficicncy (90%) and convertcr cfficicncy (95%).

Encrgy bnhmcc for the storagc system is also considcrcd. i.e. the total amount of <:ncrgy stored at the end of a time slice is the sum of the initial amount of energy stored and the net ditî crcncc of charging and discharging. Addi1ional cons1rain1s arc added in order to limit converter and batteries capacitics to ma. ximal values observed throughout the duration time of the simulation. Funhermore, the geographical innuence is considered through input modcl paramctcrs. Thcsc paramc1crs arc outsourccd clcc1rici1y (OE) pricc, wind and solar 1>ower, and environmenlal impact of each renewable energy sources. The 1>ararnctcrs rclatcd to googrnphical 1>0sîtîon arc cxtcrnal clectrîcîty pricc, îmr,ac, of sources. wind speed and solar irradimion.

The variables are the J)()WCr rate of each rcnewable energ_y source, 1>0wer capacity of convcncr (invc. rtcr and rectifier). 1hc storage cap.acity and 1he amount of OR purchased.

A multi 4 objective optimization model is solved with an economic and an environme.ntal criterin. The cconomic cri1cri1: 1 to minimizc is dcfined as the touai Net Prcscnt Value of HPS owr ils lifelime, including capital investment, maintenance cost, OJ>er.:llional ex1>enditurc for dcsigncd dccîsion variable which arc rencwablc cncrgy sources. lcad ucid batteries, and convcrtors. Finnlly cost of purchased OE fr om extcmal nct,"-' Ork is 1ake1l into accounl. The e11vironmen1al objec1ive is based on a LCA approach (c1' adle-10gravc), for cach source (rcnew1:1blc encrgy and outsourccd clectricity) producing I kWh, the impact is quantified through SIMAPRO® software. TI1is value depcnds < > n the country whcrc cnergy is produce<I. The func1ion to minimize is calculatcd by ;;1dding total impact produced by each source.

The resulting 01>timi;, .ation model is an LP model able to find a global optimal solution and il has bcen solved wi1h IOM ILOG C J> LEX Optimizcr.

Case study

ln this case study, Malaysian data arc takcn from a prcvious work so that the modcl can be validatcd (Theo et al., 20 1 6). The second location. i.e. Gennany, is chose. n to be vcry dîflCrcnt from the lirst one bc<:ausc thcy have adoptcd a policy to promotc rcncwablc energy sources, favournble to rcncwablc cnergy HPS implantation. Gcrmany is a windy country with ail in11>0t1.: uu 1>01cmial foi• wi11d trn•bi11es and low solar J )Otcn1ial while This section highlights the diff crcnces betwccn mono and muhi-objcclivc op1imi1. .ation and also invcstigatcs the i11flucncc of gcog.raphic.al position. Figure 3 shows the Parc10 fron1 solu1ions ob1aincd for Malaysia (MY) and Gcrrnany (GE). ln addition, dc1ailcd results for dcsigned HPS arc prescntcd in Table l whcrc the gcographic innuc-nce be1wccn Malaysia and Germany can be clcarly highligh1cd.

According to thcsc rcsulls, il can be no1iccd 1ha1 1hcsc critcria arc antago11ist (whcn cos1 increases, environmental impact (ENV) decreases) whatjustifies the choice to use a MO ,.,. 12" . . ,. ,,. •I "' "°" ''"' . , .. . . Tnblc 1 1.X:taikd resul1s for t' .ach designcd so1ulion

Conclusion

ln suppon of a fr. uncwork for enc. rgy sup1>ly of an cco-industrial park, wc have presented a generic me1hodology 10 0 1 Hiinî1.e 1he reocwnble energy soutces with a MO oplimizmion approach. Usint this mcthodology, wc can dclcrminc, economic cost and environmental impact of such a solution blll also the corre.s 1 xmding size for storage batteries. h also allows 10 dctcrminc the rnost profitable RE sources for a gîvcn geographic position. An interesling perspeclive, aflc.r this design of 1he EIP de 1 >ending Oil i1s gc.!ogtaphic im1)lan1.11îo1l is 10 achicvc a scnsîbility �rnalysis. To go funher flnurc improvcments for this mode) can be to allow it to select hctwccn scveral technologie.al altcma1ives for RI; sources or s-toragc. duc to its gencricily it can be donc <. •asily by 1ransforming il into an MILP. Finally becausc an EIP is conslantly cvo1ving. by using a multi•1>eriocl approach, füture research may make possible to modify nctwork ovcr time 10 adapt Îl to its moving contcxt.

  availabk also callcd RE•HPS. Lc-c et al. (2014) dcvclopcd an LP optimization modcl Lo minimi1.e. energy loss or a RE-111'$. More recen1ly, Theo et al. (2016) have completed füis Sludy by introJucic�g an l'vilLP model lo mini.mite opcrnting und invcstrncn1 cosls and to sek�t storage tcchnologles.

Figure 1 :

 1 Figure 1: Supcrstruclurc or tbe HPS model ln this HPS, the runcwabk ',j llcrgy Altcrm1tive Currcnl (AC) power sources arc Wind Turbines, Biomass and Direct Current (DC) one is Photovoltaics panels. lt is also possible 10 buy AC OUIS<.lurcéd clectricîty (OCT) i'i •o,n grid. To COl'I\ I Cl'l AC or DC electric cum:nt, rc(.:titkr and invcrtt�r arc availablc. Duc to Che intcrmittcncy of rcncwable energy production. il is n dynam1c model, in which tnnt: : is discreflzed by hour, .n any hour or tbè horizon lime.

  PhotovoltaicWind tutbines

  photovohaic pands) (cf. figur1: 2). Simulating in �1 more sunny country and including govcrnmenl he11ls should promote this tcchnology. Biomass can be choscn. depcnding on eountrics, ind(, " éd, in Gcnnany îts impac1 is. really low, Al las1, duc 10 ils low impac1 on cnvironmcnt, wind turbines appe--Jr to be the most us<.. -d solution in this design.W,nd

  Malaysia is sunnicr. Furthcm1ore, in Gcrmany, impacl of RE sources is globally really lower than in Malaysia (Figure2). To gel reliable dma, horiwn lime of this simulmion is sc1 at 4 days, one per scason, whîch allows to considcr annual climatc changing. lndccd, whilc Malaysia is ncar from e< 1 uator with low scasonal impact, in Gennany huge variations can be notcd bclwcen scasons.The diff crcnt sources takcn into account in this mode) and thcir impact arc shown in ligure 2. for Malaysia and Gennany. Values corne from SIMAPRO® sonware. using Ecoinvcnt® databasc and Impact 2002+ mc1hod. lt allows to propose a.n aggrcg_at<.-d indicator, gathering 4 damage categories: Human health, Ecosystem qualily, Climate change and Rcsourccs. The unity of thîs unîquc indicator îs the mcan environmcntal impact of a Europe.an during a year.

	Solar information, arc obtaincd with t > VGIS® tool from the Éuro1>ean commissio1t On
	1hc othcr hand wind SJ> eèd for Malaysia cornes from Theo et al. (2016) and for Gennany
	from Fraunhofer ISË Ëuro1> Can rcscarch institutc (2014). Outsourccd clcctricity pricc
	uscd in this rnodel is the average for industries in Malaysia: 0.078 USD/kWh (Theo et
	al., 2016) and in Gcrmany: 0.103 USD/kWh according 10 Clcan Encrgy \\lire (2015).

  •S3"•I"°' .. ,,. . ,.,. �• ,,. . 1 3 • . .

	Min lN\f M1,l,,ys.ia (MY) G.ermanv IGE) c.,,,p,-GE/MY (1')	124,1 149,4	9.9 '• ' 42,8 0 16S78 3S,l 0 14717	0	' '
	TOP$1S


Hichem Dammak

Ekaterina Antoshchenkova

Marc Hayoun

Fabio Finocchi
email: fabio.finocchi@upmc.fr

Isotope effects in lithium hydride and lithium deuteride crystals by molecular dynamics simulations ‡

Keywords: numbers: 05.10.-a, 88.30.rd, 65.40.De, 67.63-r

Molecular dynamics (MD) simulations have been carried out to study isotope effects in lithium hydride and lithium deuteride crystals. Quantum effects on nuclear motion have been included through a quantum thermal bath (QTB). The interatomic forces were described either within the density-functional theory (DFT) in the generalized gradient approximation (GGA) or by the phenomenological approach using the shell model. For both models, the isotopic shift in the lattice parameter can be successfully predicted by QTB-MD simulations. The slope of the experimental isotopic shift in pressure is satisfactorily reproduced by QTB-MD within DFT-GGA, in contrast to both density-functional perturbation theory and QTB-MD with the shell model. We have analyzed the reasons of these discrepancies through the vibrational densities of states and the isotopic shifts in bulk modulus. The results illustrate the importance of anharmonic contributions to vibrations and to the isotopic pressure shift between LiH and LiD.

Introduction

A considerable amount of experimental data and theoretical studies evidence that isotope substitution influences thermal, elastic and dynamical properties of crystals [1]. Isotope effects are especially pronounced whenever the reduced mass changes considerably. This happens in passing from LiH to LiD, two materials that have recently attracted much attention, as they are used in nuclear fusion [START_REF] Atzeni | The physics of inertial fusion[END_REF] and are involved in hydrogen-based energy production [START_REF] Grochala | [END_REF].

From the theoretical point of view, 7 LiH and 7 LiD are relevant test cases for any method that includes the description of quantum effects [4]. In general, the crystal with lighter isotopes has larger lattice parameters. This purely quantum effect is due to the dependence of phonon frequencies on atomic masses and on anharmonic interactions. Isotope effects are important below Debye temperature, θ D , and progressively weaken at higher temperatures.

Isotopic shift in the lattice parameter and in the equation of state of many crystals can be treated through the quasi harmonic approximation (QHA) [5]. However, this approach is not well suited for disordered systems, such as liquids, amorphous solids or crystals with defects. Moreover, many important dynamical phenomena, such as diffusion or density fluctuations, are outside the range of applicability of QHA.

Another technique that gives full access to statistical properties of materials, independently of their phase, is molecular dynamics (MD). However, in most cases, nuclei are treated as classical particles, and statistical averages are derived from classical partition functions. As quantum effects are missed, standard MD is suited at temperatures beyond θ D . Nevertheless, some authors used standard MD to describe isotope effects in lithium hydrides, at the expense of introducing two distinct empirical interatomic potentials, one fitted for 7 LiH and the other for 7 LiD [6]. This procedure is questionable from the theoretical point of view, because isotope effects are primarily of statistical nature and are not a consequence of changes of interatomic interactions. Moreover, these potentials are not necessarily transferable to complex lithium hydride phases and might fail to describe isotope effects in a wide temperature range.

Quantum effects can be included in MD simulation through the well-known path-integral formalism of Feynmann [7].

This approach provides averages of physical quantities by taking into account the exact quantum statistical distributions. Unfortunately, it is very time consuming, which limits applications to complex or large systems. Recent alternatives consist in coupling the system to a quantum thermal bath (QTB) [8] or to a colored-noise thermostat [9]. Both approaches were designed to describe quantum nuclear dynamics and save at least two orders of magnitude of computer time with respect to path-integral MD technique [10,11].

In this work, we combine QTB-MD with first-principles description of the potential energy surfaces, within the density functional theory (DFT). On one side, QTB-MD accounts for the quantum nature of the light H, D and 7 Li elements; on the other side, DFT provides a full description of the Li-H bond, in a large range of applied pressures.

We also present QTB-MD simulations based on an empirical interatomic potential [6]. The structure of the paper is as follows. The implementation of the QTB in the ABINIT code and the computational details are presented in section II. Section III reports the results of MD simulations for the isotope effects on the lattice parameter, on phonons and on the equation of state.

Computational Details

We employed the generalized gradient approximation (GGA) for describing the exchange-correlation energy, in the PBE form [12]. Previous numerical studies on LiH and LiD within the DFT employed pseudopotentials in order to describe the interaction between core (1s) electrons of Li and valence electrons from both Li and H. However, it appears that some structural and dynamical properties are sensitive to the relaxation of the Li core electrons [13]; therefore, we generated a norm-conserving pseudopotential for Li, within the Troullier-Martins scheme [14], by including the (1s) electrons in the valence. The reference atomic configuration was (1s) 2 (2s) 0 (2p) 0.25 . r s = 1.45 a 0 and r p = r d = 1.3 a 0 cutoff radii were adopted for the s, p and d components, respectively. The H pseudopotential has been generated in the Troullier-Martins scheme, with standard parameters [15]. Our calculations are thus all-electron, although Kohn-Sham orbitals were pseudized in the space region around the nuclei. Within this choice, static calculations, including no quantum effect, yielded converged lattice parameters for LiH already at 70 Ryd energy cutoff for the plane-waves basis set.

As a starting point, we computed the equilibrium lattice parameter, a 0 , the bulk modulus, B 0 , and its pressure derivative, B 0 , within the static approximation, that is, neglecting any quantum motion. These quantities are given in table 1 and agree well with previous calculations [16,17].

Table 1. Results of static calculations without taking into account the zero-point motion for LiH: lattice parameter a 0 , bulk modulus B 0 and pressure derivative B 0 of bulk modulus. All calculations adopted the GGA for the exchange-correlation energy functional. We also specify the type of equation of state (EOS) used for fitting the numerical data. present

Ref. [33] Ref. [17] Ref. [ We implemented the Quantum Thermal Bath (QTB) in the parallel version of the ABINIT code [18,19]. We started from ABINIT 6.0.4 and included the QTB by modifying subroutine moldyn. The reversible algorithm by Martyna and coworkers [22] was used to integrate the equations of motion. This algorithm allows to simulate either the canonical (NV T ) or the isobaric-isothermal (NP T ) ensembles. In both cases, the classical thermostat of the original algorithm was replaced by the QTB [8]. The ABINIT code was very slightly modified, as the original atomic forces that are derived from the electronic structure are unchanged. Two additional contributions to atomic forces have been included: a dissipative force with an effective frictional coefficient and a Gaussian random force having the power spectral density given by the quantum fluctuationdissipation theorem, thus different from a white noise as in classical Langevin dynamics. Before each run, random forces were generated according to the method described in Ref. [8] and by using a high-frequency cutoff (ω max = 100 THz, above the highest frequency of the physical LiH system), in order to avoid artifacts due to divergence of the energy spectral density at high frequencies [23]. Random forces were stored in a file that is read at each time step along the simulation.

In the equations of motion, the effective frictional coefficient was set at γ = 10 -4 au -1 , that is γ = 0.414 × 10 13 s -1 . This value is smaller than the full width at half maximum of the simulated infared absorption peak within the GGA. The mass of the barostat was fixed at 5 × 10 6 m a 2 0 , where m is the electron mass and a 0 the Bohr radius. DFT simulations were performed for samples consisting of 64 atoms, in a cubic cell. Starting from the face-centered cubic structure, the system was equilibrated for several ps, for different temperatures ranging from 10K to 300K. Afterwards, statistics on the lattice constant in the (NP T ) ensemble was taken for about 10 ps. We also calculated the pressure for 7 LiH and 7 LiD as a function of the volume in the (NV T ) ensemble, along QTB-MD runs about 6-7 ps as long, after equilibration.

For comparison purposes, we also performed QTB-MD simulations based on an empirical interionic potential. We adopted, for both 7 LiH and 7 LiD, the same incompressible shell model developed by Haque and Islam [6] for 7 LiH. This model allows to simulate electronic polarization in a phenomenological way. Periodic boundary conditions have been employed and the time step was equal to 5 × 10 -16 s. For each simulation, the system of 1000 ions was equilibrated during 50 ps before the computation of average values, which were obtained along QTB-MD trajectories 50 ps as long.

We computed the density of states (DOS) of phonons by two distinct methods: (i) from MD runs through the Fourier transform of the velocity autocorrelation function [24] and (ii) from lattice dynamical properties within density functional perturbation theory (DFPT) by interpolating the interatomic force constants on a 40 × 40 × 40 mesh in the Brillouin Zone. Concerning (i), the equilibrium configuration were generated in the (NV E) ensemble at a kinetic energy corresponding to that found within the QTB-MD simulation at 300 K. This temperature rescaling was used in order to recover the mean quantum vibrational energy [25]. For instance, at T = 300 K within the QTB, LiH shows a classical mean kinetic energy of 515 K. We chose to compute DOS along (NV E) trajectories, as, within QTB or Langevin dynamics, a spurious high-frequency tail would appear, because of the frictional coefficient γ.

Results and Discussion

Lattice parameter versus temperature

We choose the lattice parameter of 7 LiD at T = 0 K (hereafter noted a LiD 0 ), as a reference value to define the isotopic shift. The experimental value (4.049 Å) [26] is reproduced within 1.2% by QTB-MD in the GGA (4.097 Å) and within 2.5% by QTB-MD employing the shell model [6] (4.152 Å). Our DFT value fairly compares with the QHA calculation in the GGA (4.077 Å) by Yu and coworkers [17]. We ascribe the tiny discrepancy with the latter value to the different pseudo-potentials used in the two calculations.

Here we focus on the isotopic shift in lattice parameter between 7 LiH and 7 LiD, as a function of temperature T . The shifts are collected in figure 1. The experimental shifts are close to 0.016 and 0.014 Å at 0 K and 300 K, respectively. The isotopic shifts obtained from QTB-MD simulations within the GGA amount to 0.019 and 0.016 Å at 0 K and 300 K, respectively. QTB-MD simulations with the shell model yields 0.010 and 0.008 Å at 0 K and 300 K, slightly underestimating the experimental isotopic shift in the whole 0-300 K range.

As shown in figure 1, in both cases, the experimental trend in the isotopic shift is correctly reproduced by QTB-MD, although its precision depends on the representation of interatomic forces. We point out that the divergence between the lattice parameter that was computed neglecting quantum effects and that obtained by including them through the quantum thermal bath is much bigger than the isotopic shift itself. For instance, at 0 K, the inclusion of zero-point energy (ZPE) increases the 7 LiD lattice parameter by 0.069 Å when using the GGA and by 0.043 Å within the shell model.

Phonons

Isotope effects between vibrational spectra of 7 LiH and 7 LiD are mainly due to the difference in atomic masses. Indirect effects on interatomic forces due to the small difference in lattice constant can be safely neglected in the following discussion. figure 2 shows that vibrational DOS are similar in the acoustic part of the spectrum, at low frequencies, while optical contributions at higher frequencies differ substantially. This behavior fully agrees with Raman spectroscopy [27] and previous DFT calculations [13]. The shift in the optical region of the spectrum also involves larger zero-point effects for 7 LiH than 7 LiD, which implies larger equilibrium volumes for the lighter isotope.

For instance, the experimental transverse optical frequency ν T O , at ambient conditions, shifts from 17.7 THz ( 7 LiH) down to 13.4 THz ( 7 LiD) [28,29]. In DFPT within the GGA, we chose 7 LiH and 7 LiD equilibrium volumes as determined by QTB-MD runs at zero pressure and T = 300 K, in order to compute the transversal optical frequencies. Both ν T O ( 7 LiH) = 15.5 THz and ν T O ( 7 LiD) = 11.9 THz were underestimated by about 10%, mainly as a consequence of the use of the GGA. The shell model severely underestimated these frequencies, ν T O ( 7 LiH) = 13.2 THz and ν T O ( 7 LiD) = 10.0 THz, which were computed through QTB-MD runs at constant volume. The ratio of transverse optical frequencies ν T O (LiH)/ν T O (LiD) can be taken as a measure of the isotopic shift. Both experimental and simulated ratios give the expected inverse ratio between the square root of reduced masses, that is, μ( 7 LiD)/μ( 7 LiH) 1.3.

Moreover, we found that anharmonic contributions to vibrations may be relevant. In order to estimate the importance of these contributions for thermal and mechanical properties of 7 LiH and 7 LiD, we compared their phonon DOS, obtained through MD simulations in the (NV E) ensemble, as specified in Section 2, which correspond to distinct mean kinetic energies. At T = 1 K, in particular, the system is in the harmonic regime.

Boronat et al. claimed that anharmonic effects are negligible at 20 K [30]; however, as shown in figure 3, we found that they are relevant in the vibrational spectra. For instance, for the lighter and more abundant 7 LiH isotope, the gap in the phonon DOS around ν 15 THz in the harmonic regime was partially filled. Moreover, in the range of optical phonons, all peaks are smoothed; for ν ≥ 32 THz, phonons appear at frequency overtones, which were absent in the harmonic regime. Both characteristics are signatures of phonon-phonon interactions that are revealed by quantum effects, which cause LiH to probe the anharmonic part of the interatomic interactions. In the 7 LiD case, differences are less visible because of the smaller frequency range and the absence of the gap between acoustic and optical frequencies (figure 3).

Isotopic shift in pressure

The behavior of 7 LiH and 7 LiD at high pressures has been the object of considerable attention in past years. Experimental data were obtained by neutron [31] or x-ray diffraction [32]. The latter work showed that 7 LiH remains in a rocksalt structure up to 36 GPa and 7 LiD up to 94 GPa. Isotope effects on the equation of state (EOS) were discussed.

The EOS of 7 LiH and 7 LiD were fitted to ab initio QTB-MD data at T =300 K through the Vinet equation [21]. The resulting parameters, a(P = 0), B and B are listed and compared with experimental values in table 2. Therefore, within the statistical error, we did not find any isotopic shift in bulk modulus. This is consistent with the fact that acoustic phonons are pretty much the same in the two isotopes, as discussed In all simulations, the lattice parameter was set, for each isotope, to its ambienttemperature value including quantum effects, as shown in figure 1.

in the previous section. A similar behavior was reported [1] for other elastic constants, within the experimental uncertainty.

Table 2. Lattice parameter a, bulk modulus B, and its pressure derivative B , for both 7 LiD and 7 LiH at T = 300 K and P = 0. All calculations adopted the GGA for the exchange-correlation energy functional.

MD-QTB Exp.

QHA this work

Ref. [32] Ref. [33] a(LiD) ( Å) We computed the isotopic shift in pressure, ΔP , between 7 LiH and 7 LiD versus the applied pressure P on 7 LiD from QTB-MD simulations at 300K. Our results are close to the experimental data [32], in the whole 0-20 GPa range of applied pressures, both in value, in spite of being overestimated, and increasing trend (figure 4). Therefore, QTB-MD within the GGA reproduced the slope of isotopic shift in pressure accurately; in particular, the higher P , the smaller the volume, the more relevant quantum effects on nuclear motion become. Open circles are experimental results [32] and dotted line is DFPT calculation of Ref. [13]. We mention that ΔP = 0 in standard MD, whatever the atomic forces.

The pressure on lighter 7 LiH thus increases faster than on heavier 7 LiD. QTB-MD simulations using the 7 LiH shell model of Haque and Islam [6] yielded constant ΔP (figure 4) close to the experimental value obtained for zero pressure on 7 LiD. Therefore, this interatomic potential is not suitable to study lithium hydride at high pressures.

QHA within the DFPT [13] provided a fair estimate of the isotopic shift ΔP at low pressure on 7 LiD (figure 4). However, the slope of the curve ΔP significantly differed from the experimental one, especially at high pressures. Although the authors of ref. [13] used the LDA and different pseudopotentials from ours, we attribute this discrepancy to a failure of the quasi-harmonic approximation to represent the isotopic shift in pressure faithfully within the linear response regime. As pressure increases and volume shrinks, interatomic force constants get bigger and optical frequencies may substantially increase. Anharmonic contributions at high frequencies, which are not the same in 7 LiH and 7 LiD, are not accounted for by the QHA ( §3.2 and figure 3). We guess that this is the main reason why ΔP from the QHA diverges from experimental data at high pressure.

In order to discuss the behavior of isotopic shift in pressure ΔP at low P , we rely on the Murnaghan EOS [20]. From a Taylor expansion of the Murnaghan EOS in terms of the equilibrium volumes V 0 , bulk moduli B 0 and their pressure derivatives B 0 for 7 LiH and 7 LiD, respectively, one obtains:

∂ΔP ∂P B 0 ΔV 0 V 0 + ΔB 0 B 0 (1)
where ΔV 0 and ΔB 0 are the isotopic shifts in equilibrium volume and bulk modulus, respectively; other quantities are referred to 7 LiD. At P 0, experimental ∂ΔP ∂P 2.5; in our QTB-MD simulation within GGA, the first term (4.1) in equation ( 1) is much bigger than the second one (-1.3), yielding ∂ΔP ∂P 2.8. Indeed, the isotopic shift in bulk modulus is almost null within DFT-GGA. QTB-MD simulations using the shell model provided a non negligible shift ΔB = -1 GPa [START_REF]At T = 300K, fitted B 0 are 44.5 GPa and 45.5 GPa for 7 LiH and 7 LiD, respectively, when using the shell model[END_REF]. As a consequence, the two terms on the right side in equation ( 1) have opposite signs and almost compensate each other (1.9,-2.2), yielding a roughly constant ΔP , with null slope at the origin. Possibly, an improved shell model, fitted on high-pressure data, might reproduce the experimental isotopic shift in pressure.

Conclusion

We studied isotope effects in 7 LiH and 7 LiD by using a quantum thermal bath (QTB) that accounts for quantum effects on nuclear motion in molecular dynamics (MD) simulations. Two non-equivalent descriptions of interatomic forces have been considered: the ab initio DFT-GGA and a phenomenological approach using the shell model. In both cases, the experimental isotopic shift in the lattice parameter is reproduced, showing the adequacy of the QTB for describing such effects.

Acoustic vibrational DOS are similar for 7 LiH and 7 LiD, while optical phonons are shifted to higher frequencies for 7 LiH. The comparison of the phonon DOS obtained either in the harmonic regime or at kinetic energies that are representative of system dynamics including quantum effects illustrates the importance of anharmonic contributions.

QTB-MD within DFT-GGA reproduces the experimental isotopic shift in pressure with the correct slope, in contrast to the quasi-harmonic approximation, which fails at high pressures. The discrepancy between the two methods again demonstrates the importance of anharmonic contributions. Besides, QTB-MD by using the 7 LiH shell model of Haque et al. [6] does not provide the experimental slope of the isotopic shift in pressure. We attribute this discrepancy to the (likely erroneous) isotopic shift in bulk modulus as found within the shell model. Indeed, in both experiments and GGA-based QTB-MD simulations, bulk moduli of 7 LiH and 7 LiD are equal within the statistical uncertainty whereas, for the shell model, the isotopic shift in bulk modulus is significant. These results might provide a route to derive suitable shell models for reliable simulations of light elements at high pressures. A new potential could be obtained by fitting the shell-model parameters to the Vinet equation of state derived from static DFT-GGA calculations without any quantum effects. Its validity could be checked by taking into account isotopic shifts in bulk modulus and pressure.

Figure 1 .for both 7

 17 Figure 1. Isotopic shift in lattice parameter, Δa = aa LiD 0 for both 7 LiH and 7 LiD, versus T . (a) Experimental results [26]; QTB-MD simulations by using (b) DFT within the GGA or (c) the shell model. Fits (dotted lines) were obtained through the expression Δa = A + B/[exp(T D /T ) -1] where T D is analoguous to a Debye temperature. In order to point out quantum effects, panel (d) shows also the results of standard MD with the shell model up to T = 600 K, which leads to same results for 7 LiH and 7 LiD.

Figure 2 .

 2 Figure 2. Vibrational DOS for 7 LiH and 7 LiD obtained through DFPT within the GGA for the optical and acoustic bands.

Figure 3 .

 3 Figure 3.Vibrational DOS for7 LiH and 7 LiD. They have been computed via the Fourier transform of the velocity autocorrelation function within (NV E) MD, at 1 K (dashed line), or at 515 K ( 7 LiH) and 435 K ( 7 LiD) mean kinetic energies (full lines). In all simulations, the lattice parameter was set, for each isotope, to its ambienttemperature value including quantum effects, as shown in figure1.

Figure 4 .

 4 Figure 4. Isotopic shift, ΔP , defined as the difference in pressure between7 LiH and 7 LiD at a given volume, versus pressure P in 7 LiD. Filled circles and squares with a linear fit are QTB-MD calculations within DFT-GGA and shell model, respectively. Open circles are experimental results[32] and dotted line is DFPT calculation of Ref.[13]. We mention that ΔP = 0 in standard MD, whatever the atomic forces.

‡ This work was performed using HPC resources from the "mesocentre" computing center of Ecole Centrale Paris (ECP), which is supported by ECP and CNRS.