Fabien Brieuc

Hichem Dammak
email: hichem.dammak@centralesupelec.fr

Marc Hayoun

Quantum Thermal Bath for Path Integral Molecular Dynamics Simulation †

The quantum thermal bath (QTB) method has been recently developed to account for the quantum nature of the nuclei by using standard molecular dynamics (MD) simulation. QTB-MD is an efficient but approximate method for dealing with strongly anharmonic systems, while path integral molecular dynamics (PIMD) gives exact results in a huge amount of computation time. The QTB and PIMD methods have been combined in order to improve the PIMD convergence or correct the failures of the QTB-MD technique. A new power spectral density of the random force within the QTB has therefore been developed. A modified centroid-virial estimator of the kinetic energy, especially adapted to QTB-PIMD, has also been proposed. The method is applied to selected systems: a one-dimensional double well system, a ferroelectric phase transition, and the position distribution of an hydrogen atom in a fuel cell material.

The advantage of the QTB-PIMD method is its ability to give exact results with a more reasonable computation time for strongly anharmonic systems.

1 Introduction 1 The quantum nature of nuclei can play a major role at low temperatures and/or in systems containing light atoms. In this case, nuclear quantum effects cannot be neglected and have to be taken into account in atomistic simulations. A standard way of including them is to use path integral methods [1][START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF][START_REF] Marx | [END_REF] such as path integral molecular dynamics (PIMD) which provide exact quantum results even for strongly anharmonic systems. Unfortunately, time correlation functions are not accessible and require the use of techniques such as ring polymer molecular dynamics 4 or centroid molecular dynamics. 5 Moreover, the computation time required for PIMD reduces its range of applicability, in particular when using a first-principle description of the interatomic forces.

Alternative techniques based on a modified Langevin equation have recently been proposed to include nuclear quantum effects in molecular dynamics (MD) simulations. 6,7 Among them the quantum thermal bath (QTB) method 6 is an approximate approach, very simple to implement in an existing MD code, and that yields accurate results for various types of systems. [8][9][10][11][12][13] In addition, it requires the same amount of computation time as standard MD and time correlation functions are directly accessible. The QTB method becomes approximate when dealing with strongly anharmonic systems and therefore suffers from the zero-point energy (ZPE) leakage. 14,15 The ZPE leakage is a known problem [16][17][18][19] in classical trajectories where energy flows from the high-frequency modes, with large zero-point energies, to lowfrequency modes with smaller zero-point energies. An empirical solution reducing the ZPE leakage within the QTB method has been attempted but with limited efficiency. 14,15 In this work, the QTB is combined with the PIMD method in order (i) to avoid the ZPE leakage within QTB for strongly anharmonic systems and (ii) to reduce the computation time of PIMD. The convergence of the PIMD is significantly accelerated since less replicas are needed for convergence. The combination of QTB with PIMD is first described and then applied to selected systems: a one-dimensional double well system, a ferroelectric phase transition, and the position distribution of an hydrogen atom in a fuel cell material.

2 Combining Path integrals with the quantum thermal bath [START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF] The goal is to combine the PIMD with the QTB method in a way similar to that developed by Ceriotti et al. 20 In our case, the idea is to replace the classical thermostat applied to each bead of the polymer by the QTB. To achieve this, it is necessary to find the appropriate power spectral density of the random force which will then depend on the number of replicas, P .

Equations of dynamics

In both the PIMD and the combined QTB-PIMD formalisms, the equation of motion of each

atom i (i = 1, • • • , N) of a replica s (s = 1, • • • , P ) is given by ṗi,s = 1 P f i,s -m i ω 2 P (2r i,s -r i,s+1 -r i,s-1 ) -γ p i,s + R i,s (1)
where r i,s , p i,s , and f i,s are the atomic position, the momentum, and the force exerted by all the other atoms of the replica s. The spring constant between beads is equal to

m i ω 2 P = m i P k 2 B T 2 h2 (2)
The last two terms of eq 1 correspond to the friction and stochastic forces of the thermostat, respectively.

The power spectral density of the random force, I R i , is obtained from the fluctuation-dissipation theorem: 21,22

I R i (ω, T ) = 2 m i γ κ(ω, T ) ( 3 )
In standard PIMD, when the Langevin thermostat is used, the stochastic force is a white noise and κ(ω, T ) = k B T . For the combined QTB-PIMD, the power spectral density is ω-dependent and corresponds to the energy of the oscillator ω, which matches the energy θ(ω, T ) of the quantum harmonic oscillator when P = 1. κ(ω, T ) will be determined in section 2.2. The correlation function of the random force satisfies the Wiener-Khinchin theorem:

R i,s,α (t)R i,s,α (t + τ ) = +∞ -∞ I R i (ω, T ) exp[-iω τ ] dω 2π (4)
The random force R i,s,α (t) is computed using the numerical technique described in the Appendix of Ref. 23

Derivation of the power spectral density

x 2 = 1 m ω 2 θ(ω, T ) = h 2 m ω coth β hω 2 (5)
where β = 1/(k B T ). In the PIMD scheme, one can transform the coordinates of the P

replicas, x 1 , • • • , x s , • • • , x P , into normal modes, q 0 , • • • , q k , • • • , q P -1 with pulsations ω 2 k = ω 2 P + 4 ω 2 P sin 2 k π P (6)
and the mean square fluctuation is then obtained according to

x 2 = 1 P P -1 k=0 q 2 k (7)
The mean potential energy of the normal modes is equal to 1 2 κ(ω k , T ) and then

q 2 k = 1 m ω 2 k κ(ω k , T ) ( 8 )
Now, let us determine κ when performing QTB-PIMD. The function κ(ω, T ) must allow one to recover the expected position fluctuation given by eq 5

1 P P -1 k=0 1 m ω 2 k κ(ω k , T ) = h 2 m ω coth β hω 2 (9)
Defining the dimensionless quantities

u = β hω 2 (10)
h(u) = u coth(u) (11) f (0)
P (u) = β P κ 2u β h ( 12
)
eq 9 becomes

P -1 k=0 u 2 u 2 k f (0) P (u k ) = h(u) (13)
where u k is the reduced pulsation according to equations 6 and 10

u 2 k = u 2 P + P sin 2 k π P (14)
In the definition of κ through eq 9 all the normal modes are treated in the same way.

There is an alternative definition 24 in which the normal mode at k = 0 (centroid of the ring polymer) is classically considered, i. e. κ(ω 0 , T ) = k B T . In this case, eq 9 becomes

1 P k B T m ω 2 0 + 1 P P -1 k=1 1 m ω 2 k κ(ω k , T ) = h 2 m ω coth β hω 2 (15)
Using the same dimensionless quantities eq 15 leads to a new equation to be solved

P -1 k=1 u 2 u 2 k f (1)
P (u k ) = h(u) -1 (16)
Equations 13 and 16 can be solved by using the self-consistent iterative technique of Ceriotti et al. 20 The numerical calculation of f (0) P [START_REF]The f P function is different from the g P function of Ceriotti et al. 20 because the authors used another formulation of the PIMD (at temperature P × T ). Indeed[END_REF] and f

P are reported in Appendix A.

2.3

The estimation of macroscopic properties 4 In the PIMD method, the potential energy, U, is calculated using the expression

U = 1 P s V (r 1,s , • • • , r N,s ) (17)
whereas two expressions for the estimator of the kinetic energy are usually used; the primitive estimator given by

K prim = i,s p 2 i,s 2 m i - i,s 1 2 m i ω 2 P (r i,s -r i,s+1 ) 2 (18)
and the centroid-virial estimator given by

K Cvir = 3N 2 k B T - 1 2P i,s (r i,s -r c i ) • f i,s (19)
where r c i = s r i,s /P is the centroid of the ring polymer i. The last estimator is known to exhibit weaker fluctuations, which are insensitive to P compared to the primitive estimator for which fluctuations grow with P . [START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF] Combining the QTB and PIMD methods includes a part of the quantum fluctuations in the momenta, i.e. i,s

p 2 i,s 2 m i > 3NP k B T 2 ( 20
)
thus expression 19 underestimates the kinetic energy for the QTB-PIMD method and it is more suitable to replace the classical energy 3 2 Nk B T by P times the kinetic energy of the N normal modes centroids as follows

K mCvir = P i (p c i ) 2 2m i - 1 2P i,s (r i,s -r c i ) • f i,s (21)
where p c i is the momentum of the centroid i. The factor P takes into account the factor √ P between the centroid coordinate and the normal mode coordinate q i,0 = s r i,s / √ P of the ring i.

The two expressions, 19 and 21, are equivalent for any interatomic potential in the case of the original PIMD, but only for harmonic potentials in the case of QTB-PIMD using the f

P function. In the case of an anharmonic potential, each centroid i, i.e. the normal mode q i,0 , is coupled to the other normal modes of the corresponding polymer. Consequently, the centroid temperature is different from the thermostat one (QTB), since a part of the quantum effects is included in the momenta of the internal modes q i,k>0 of the polymer. On the other hand, only expression 21 may be used for the QTB-PIMD using the f

P function, because in this case the dynamics of the centroid also includes quantum effects.

Next, we want to compare the two formulations of the QTB-PIMD method using either the f (0) P or f

(1) P functions (eq 13,16

) and to choose the adequate kinetic energy estimator among expressions 18, 19, and 21. For this, let us consider a system for which the pressure is zero and thus the kinetic and potential contributions cancel each other out. This allows one to express the kinetic energy using the virial estimator, 26,27 K vir , from the general expression of the pressure

p V = 2 3 K + 1 3P i,s r i,s • f i,s (22)
K vir = - 1 2P i,s r i,s • f i,s (23)
This estimator is used as a reference to validate the relevance of the other estimators. The test is performed on two cases of the Morse potential characterized by the dimensionless parameter λ defined as

1 λ 2 = h2 α 2 /(2mD) (24)
where D and α -1 are the depth and width of the well, and m is the reduced mass. The anharmonicity of the potential can be expressed as the relative shift of the ground state energy with respect to the harmonic approximation value, (4λ 2 ) -1 . The first case is weakly anharmonic (λ -2 = 0.0015 as the case of the HCl molecule) while for the second case the anharmonicity is chosen to be more important (1/λ 2 = 0.024).

Figure 1 shows that the QTB-PIMD method, either using f

(0) P or f (1)
P , allows a faster convergence with the number of beads than the PIMD one, as expected. For both methods the convergence is slower in the case of the more anharmonic potential. For the weakly anharmonic model, the potential energy (Figure 1 For the two potential models, it is worth noting that K prim and K Cvir systematically give the highest and the lowest values for the kinetic energy, respectively (Figures 1(b) and 1(d)). In contrast, K mCvir provides values very close to those obtained with K vir which is considered as a reference in this example. Figure 2 shows the influence of the effective friction coefficient γ on the kinetic energy estimators. It is clearly shown that K mCvir (eq 21) provides the better estimation and is, in particular, insensitive to γ when using the f

P function, whereas the primitive estimator needs the use of a low value of γ leading to the increase of the computation time.

This example shows that both definitions of the f P functions (eq 13 and eq 16) allow similar convergences of the potential and kinetic energies with the number of beads. The best estimator for the kinetic energy is the modified centroid-virial one given by eq 21. In the following sections it is also shown that position distributions of atoms obtained by using P functions respectively. The relative deviation from the virial estimator values, (K -K vir )/K vir , is plotted. γ is normalized by the pulsation ω min which is equal to that of the centroid (normal mode k = 0 of the ring polymer) in the harmonic approximation: ω min = ω 0 = ω/ √ P (eq 6) . Case of the Morse potential where 1/λ 2 = 0.00145 using P = 4 and T = 0.02 D/k B . 5 Let us consider a particle of mass m in a double well potential

Applications

Position distribution in a one-dimensional double well potential

V (x) = V 0 (x/a) 2 -1 2 (25)
Using reduced units

y = x/a, = E/V 0 ( 26
)
for the position and the energy respectively, the equation for the stationary wave functions φ writes

-C d 2 φ dy 2 + (y 2 -1) 2 φ = φ ( 27
)
where C is a parameter depending on the barrier height, V 0 , and the distance between the two wells, 2a

C = h2 /(2ma 2 V 0 ) (28)
The numerical resolution of eq 27 shows that there exists a critical value for the parameter, C 0 = 0.731778, when the ground state energy equals V 0 . The eigenvalues rely on C and we find in particular that the energy of the ground state, 0 , is lower than the height of the energy barrier, i. e. 0 < 1, when C < C 0 . The motivation behind using such a one-dimensional quartic double-well potential is to provide a simple model to check the ability of the QTB-PIMD method to recover the tunnel effect in contrast to the original QTB approach. 28,29 We investigate the position distribution of a particle in this double-well potential at a reduced

d P = +∞ -∞ (ρ -ρ 0 ) 2 dy +∞ -∞ ρ 2 0 dy (29)
similar to the reliability factor used in the Rietveld method [START_REF] Young | The Rietveld Method[END_REF] to refine a theoretical line profile until it matches an experimental profile.

In all cases, the evolution of the divergence factor and the distributions as a function of P obtained with the QTB-PIMD method are similar when using either the f (0)

P or the f

(1) P function.

Figure 3 shows the evolution of the divergence factor as a function of the number of beads for the three investigated cases. In the first case, the value of the C parameter , C = 1, is greater than the critical value C 0 (eq 28). As shown by the divergence factor in Figure 3 The QTB position distribution is also plotted (dotted line) for comparison purposes. The potential energy curves are superimposed (gray full line) and the horizontal gray dotted lines show the total energy of the system in each case; the values of the barrier height V 0 are deduced from C through eq 28 for the hydrogen atom and the distance between minima is fixed at 2a = 0.8 Å. one, within an error of about d P = 2% at P = 16 and P = 21, respectively. Hence, the full convergence is especially difficult to reach in this case. Indeed, the exact distribution exhibits only one maximum located at y = 0 (Figure 3(d)), while a poor convergence leads to a distribution with two maxima. One of these situations is illustrated by the QTB method which dramatically fails with an error of d P = 40%. In the two other cases, C is lower than C 0 and all the position distributions show two maxima (Figures 3(e) and 3(f)). For the QTB-PIMD simulation, the number of beads required to converge whithin an error of about d P = 2% decreases when C decreases; P equals 7, and 4 for C = 0.3, and C = 0.1, respectively (Figures 3(b) and 3(c)). It is worth noting that, for the lowest value of C, the distribution obtained with the QTB method is in good agreement with the exact one. This case corresponds to high values of either the barrier height, V 0 , or the distance, 2a, between the two wells, or the particle mass, m. In other words, the QTB method succeeds when the tunnel effect is not predominant.

In summary, the combination of the QTB and the PIMD methods allows a better convergence with the number of beads. In fact, for a high accuracy in the position distribution the gain obtained by QTB-PIMD with respect to PIMD is not so important and is lower than a factor of two. In contrast, the advantage of the QTB-PIMD is more substantial -a gain of a factor of three -with respect to the convergence of the total energy. 6 We now investigate the BaTiO 3 (BTO) ferroelectric crystal which is characterized by a complex multiple-well energy landscape affecting the polar degrees of freedom. In this compound the quantum effects significantly decrease the phase transition temperatures by ≈ 30-50 K 31 and strongly modify the shape of the pressure-temperature phase diagram 32 even at room temperature.

Ferroelectric-paraelectric phase transition

The ferroelectric properties of BTO are modeled by using the effective hamiltonian of Zhong et al., 33 which was derived from first-principles density-functional calculations. This approach yields an excellent description of its complex sequence of phase transitions: rhom-bohedral(R) -orthorhombic (O) -tetragonal (T) -cubic (C). The T-C phase transition temperature, T C , is shifted from 300 K to 260 K when using the standard MD and PIMD methods, respectively. 34 In this work, we investigate the convergence of the value of the phase transition temperature (T C ) as a function of the number of beads by performing either QTB-PIMD or PIMD simulation in the isothermal-isobaric ensemble.

Figure 4 displays the polarization as a function of the temperature obtained by QTB-PIMD compared to the converged PIMD result. Favorably, the QTB-PIMD method with a number of beads equal to 2 gives a transition temperature T C = 257 K very close to that obtained by PIMD with P = 16 (259 K). In addition, the inset in Figure 4 shows that the convergence of the polarization of the ferroelectric phase is faster when using the f P referenced as (0) and (1), respectively. Vertical dashed lines show the quick convergence of the Curie temperature obtained by QTB-PIMD: 202 K (P = 1), 257 K (P = 2), and 259 K (P = 3). The inset provides the convergence of the polarization at T = 240 K with the number of beads for the three methods.

Proton disorder in a fuel-cell oxide BaZrO 3

7 We now investigate the proton disorder in the cubic phase of BaZrO 3 (BZO). 35,36 Indeed, it has been shown that perovskite type oxides (ABO 3 ) can exhibit high protonic conductivity 37 and that quantum effects play an important role. Among these oxides, doped BZO exhibits one of the highest conductivity combined with a good chemical stability 37,[START_REF]Proton Conductors: Solids, membranes and gels -materials and devices[END_REF] The interactions between the atoms are modeled by an ab initio based force field [START_REF] Raiteri | [END_REF] that reproduces the ab initio computed activation energies for the transfer and reorientation mechanisms. [START_REF] Raiteri | [END_REF]40 The reactivity of the proton is included through the empirical valence bond model. No dopant is included in our simulations, so the distributions presented here illustrate the situation in a region far from the dopant atom. To ensure the electrical neutrality in the computations, a compensating uniform background charge is added. The MD simulations were performed within the canonical ensemble on a 3 × 3 × 3 simulation box containing 136 atoms.

The 3D position distribution of one proton in BZO is computed at T = 300 K using QTB-PIMD and PIMD. As an illustration, the classical MD proton position distribution in the O-Zr-O plane is displayed in the top part of Figure 5. From this distribution one can note that the hydrogen atom remains bonded to the same oxygen atom, most of the time.

One type of reorientation is essentially seen where the hydrogen atom performs a rotation around the nearest oxygen atom and remains in the same O-Zr-O plane. Figure 5 shows the evolution of the position distribution with the number of beads as obtained by QTB-PIMD and PIMD. Quantum effects are important at this temperature since the classical distribution corresponding to PIMD with P = 1 is clearly different from the quantum distribution obtained by PIMD or QTB-PIMD with high P values. In the classical case the distribution exhibits two peaks, whereas in the quantum case the proton freely rotates around the nearest oxygen atom. Anyway, the transfer mechanism is rarely observed at T = 300 K. For the QTB-PIMD simulation, Figure 5 shows that the two peaks disappear from P = 3, and that the position distribution already converges to the exact one for P = 4 . In contrast, the PIMD method requires a value of at least P = 8 for convergence. A more accurate comparison of the two methods is given in Figure 6 in which the corresponding divergence factors d P (a generalization of eq 29 to 3D distributions) are plotted as a function of P .

The PIMD calculations are converged within an error of d P = 8% for P = 6, while the QTB-PIMD method requires a smaller number of beads, P = 4. We conclude as expected that QTB-PIMD is more efficient than the standard PIMD, since less replicas are needed for convergence.

4 Conclusion 8 We have combined the path integral MD (PIMD) method with the quantum thermal bath (QTB) approach. This combination can be seen in two ways. In the first one, the QTB is used as a thermostat applied to standard PIMD in order to improve the PIMD convergence.

Compared to standard PIMD this combination needs less replicas to converge -the gain is a factor of two or three -since a part of the quantum effects is included directly in the dynamics of the beads through the QTB. In the second way, a small number of replicas of the path integral allows the correction of the failures of the QTB-MD. Indeed, the QTB-MD technique gives rather satisfactory macroscopic properties for weakly anharmonic systems, especially away from phase transitions and in the case of non predominant tunnel effects. Thus, QTB- MD stays an alternative method for who is not interested in high accuracy computations.

Small failures of the QTB can be fixed by only two replicas in QTB-PIMD. In contrast, when QTB-MD dramatically fails, QTB-PIMD can address the problem at the price of high computation cost. The advantage of the QTB-PIMD method is its ability to give exact results with a more reasonable computation time for strongly anharmonic systems including the tunnel effect between two wells. Unfortunately, as in PIMD, the major drawback of the new method is that the time dependent correlation functions are not directly accessible.

The combination with PIMD requires the modification of the power spectral density of the random force within the QTB. This spectral density is proportional to a reduced function (f P ) which can be defined in two ways. In the first one, random forces are applied to each bead of the ring polymer (f (0) P ). An alternative way is to apply random forces to the normal modes of the ring polymer (f The combination procedure is similar to the one presented by Ceriotti et al. 20 The iterative algorithm described in Appendix A to determine the f P function is the same as the one used to establish the g P function of Ref. 20. The difference lies in the choice of the thermostat used to include the quantum effects. In the colored-noise thermostat (GLE) of Ceriotti et al., 7 the quantum effects are introduced through a dispersive friction coefficient while in the QTB case the quantum effects are included through the power spectral density of the random force. The two methods are basically equivalent but the QTB-MD method is easier to implement in a PIMD code. Knowing the f P function, the random forces can be directly generated. 23 In contrast, the GLE method requires to carefully optimize different parameters in the equations of motion (see eq. 8 of Ref. 20) in order to recover the quantum fluctuations.

The modified QTB is easy to include in any PIMD code and the implementation does not increase its complexity. The combination of the QTB and the PIMD methods within first principle descriptions is even more interesting. Such an implementation in the ABINIT code 41,42 is in progress. Moreover, the possibility of combining the modified QTB with methods such as ring polymer molecular dynamics (RPMD) 4 or centroid molecular dynamics (CMD) 5 will allow one to go beyond PIMD. In this case, time dependent correlation functions could be computed. Since the combination with the QTB consists in replacing the classical thermostat by the quantum thermal bath, the combination with RPMD (thermostated RPMD has been recently proposed 43 ) or CMD would have the same complexity as the one of the QTB-PIMD. Let us point out that the combination requires the use of physical bead masses instead of fictitious ones.

A Self-consistent resolution of equations 13 and 16 9 The f (0) P and f (1) P functions are determined through a self-consistent resolution of the equations.

In the case of the f (0) P function, eq 13 can be reformulated by isolating the k = 0 term for which

u 0 = u √ P (30)
f P ( u √ P ) = 1 P h(u) - P -1 k=1 f P (u k ) (u 2 k /u 2 ) (31)
The superscripts (0) is omitted for simplicity. Before solving, eq 31 is rewritten using the function

F P (u) = f p (u/ √ P ) (32)
F P (u) = 1 P h(u) -

P -1 k=1 F P (u k √ P ) (u 2 k /u 2 ) ( 33
)
We choose an initial solution with a good asymptotic behavior for eq 31:

P (u) ( 35
)
where α is a weighting parameter whose value giving the best convergence is close to 1/P .

The random force of the QTB is generated from the power spectral density in the pulsation space in a range [ω min ; ω max ]. ω max and ω min are related to the MD time step δt and the time-length of a MD trajectory, respectively. It results in the F P (u) function needing to be determined in the correlated range [u min ; u max ]. Unfortunately, eq 35 shows that F P must be calculated for the value u k √ P which reaches a value, u 2 max + P 2 , greater than u max .

To overcome this problem, the values of F P for u greater than u max are linearly extrapolated from the last 20% of the u range. The F P function is obtained with enough accuracy over f

P functions

I R i = 2m i γ P β f P (βhω/2) ( 43
)
corresponds to that of the random forces R i,k,α to be applied to normal modes (i, k) of the ring polymer (i). In the case of the f (0) P function, one can show that the power spectral densities of the random forces applied to beads or to normal modes are equal. Hence, the random forces R i,s,α can be directly generated according to eq 4. In contrast, for the f

  (a)) remains almost constant as a function of P showing that the QTB approach (P = 1) already provides a value very close to the exact one. The kinetic energy behaves similarly (Figure 1(b)) but depends on the estimator.

FFigure 1 :Figure 2 :

 12 Figure 1: Convergence of the potential, U, and kinetic, K, energies as a function of the bead number for PIMD and QTB-PIMD simulations using either f (0) P or f (1) P , herein referenced as (0) and (1), respectively. Case of two Morse potentials, λ -2 = 0.00145 for (a) and (b), and λ -2 = 0.024 for (c) and (d). The energies are normalized by the well depth, D, and calculations were carried out at T = 0.02 D/k B . K values were obtained using either primitive (prim), virial (vir), centroid-virial (Cvir) or modified centroid-virial (mCvir) estimators given by equations 18, 23, 19 and 21, respectively. The exact values were derived by numerically solving the Schrödinger equation and occupying the resulting eigenstates with the proper Boltzmann factor.

  temperature T * = k B T / 0 = 0.4 for three values of the C parameter: 1, 0.3, and 0.1. The position probability density, ρ(y), obtained by QTB-PIMD simulation is compared to the exact one, ρ 0 (y), obtained by numerically solving the Schrödinger equation and occupying the resulting eigenstates with the proper Boltzmann factor. The convergence of the density is evaluated by calculating the divergence factor d P

3 DFigure 3 :

 33 Figure 3: Divergence factor, d P , as a function of the number of beads for the double-well potential given by eq 27. Three cases of this potential are investigated: (a) C = 1.0, (b) C = 0.3, and (c) C = 0.1. Calculations have been carried out at a reduced temperature T * of 0.4 for PIMD, QTB, and QTB-PIMD using the f (1) P function. The arrows indicate the smallest numbers of beads obtained for d P < 2%. The calculated position distributions corresponding to the QTB-PIMD case are shown (full line) in (d) P = 16, (e) P = 7, and (f) P = 4, together with the exact solution (dashed line).The QTB position distribution is also plotted (dotted line) for comparison purposes. The potential energy curves are superimposed (gray full line) and the horizontal gray dotted lines show the total energy of the system in each case; the values of the barrier height V 0 are deduced from C through eq 28 for the hydrogen atom and the distance between minima is fixed at 2a = 0.8 Å.

  the QTB-PIMD method. In contrast, the use of the QTB method (QTB-PIMD with P = 1) strongly underestimates the value of T C by ∼ 55 K. It is worth noting that the error of the QTB method is fixed solely by two replicas when combining QTB with PIMD. This means that the effects of the zero-point energy leakage associated to the use of the QTB have been suppressed by the combination.

Figure 4 :

 4 Figure 4: Temperature evolution of the reduced polarization associated to the tetragonal (T) to cubic (C) ferroelectric transition as obtained by PIMD (P = 16) and QTB-PIMD simulations using either f (0) P or f (1)

  which makes it a potential candidate as electrolyte materials for proton conducting fuel cells. The long range migration of a proton in a BZO crystal is a combination of transfer and reorientation mechanisms. During the transfer step the hydrogen atom jumps between the two neighboring oxygen atoms while, during the reorientation stage the hydrogen atom performs a rotation around the nearest oxygen atom. There are two possible rotations: one around the Zr-O-Zr axis and one around an axis orthogonal to the Zr-O-Zr axis with the hydrogen atom remaining in the same O-Zr-O plane.

Figure 5 :

 5 Figure 5: Position distribution of the proton in the O-Zr-O plane at T = 300 K as obtained by PIMD and QTB-PIMD simulations. The top image displays an example of the distribution in the unit cell (PIMD with P = 1) where one can distinguish the eight equivalent positions for the proton. The oxygen and zirconium atoms are shown by red and yellow full circles, respectively. The other images are enlargements of the distribution including contour lines to bring out the increase of the proton disorder by rotation around the neighbor oxygen atom as the number of beads increases.

Figure 6 :

 6 Figure 6: Divergence factor, d P , (generalisation of eq 29 in 3D) as a function of the number of beads for the proton distribution in BZO at T = 300 K obtained by PIMD and QTB-PIMD simulations using either f (0) P or f (1)P , herein referenced as (0) and (1), respectively. The horizontal dotted line allows one to determine the number of beads corresponding to d P = 0.08.

( 1 )

 1 P ) except for the centroid mode dealt with a Langevin thermostat. It is shown that the f (1) P function gives a better convergence of the macroscopic quantities with the number of beads than f (0) P . Considering that a part of the quantum fluctuations is included in the momenta through the QTB contribution, a modified centroidvirial estimator of the kinetic energy is proposed. This estimator is accurate and insensitive to the effective friction coefficient when using the f (1) P function.

FF

  exact solution in the case P = 1. Using this initial solution and following Ceriotti et al.20 the equation is iteratively solved as

Acknowledgement 10

 10 the centroid mode (k = 0) is not dealt in the same way as the other normal modes (k > 0), one has to first generate the random forces R i,k,α and then to transform them into R i,s,α by using the orthogonal transformation matrixC sk R i,s,α = k C sk R i,k,αThis work was performed using HPC resources from the "mesocentre" computing center of CentraleSupélec, which is supported by CentraleSupélec and CNRS. 10 J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

Considering a one-dimensional (1D) harmonic potential energy: V (x) = 1 2 mω 2 x 2 , the quantum mean square fluctuation of the position x at temperature T is given by

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI: 10.1021/acs.jctc.5b01146

J. Chem. Theory Comput. 12, (2016) 1351-1359. DOI:

10.1021/acs.jctc.5b01146

about 30 iterations in eq 35. Results show that the function f P (u) = F P (u √ P ), directly related to the spectral energy κ(ω, T ), tends to 1/P and u/(P √ P ) when u tends to 0 and ∞, respectively. Figure 7 presents the curves of the f P (u) solutions obtained for different values of P .

The superscripts (1) is omitted for simplicity. Before solving, eq 37 is rewritten using the function

We choose an initial solution with a good asymptotic behavior for eq 37:

Using this initial solution, the equation is iteratively solved as

where α is a weighting parameter whose value giving the best convergence is close to 1/P .

The F P (u) function is determined following the above mentioned procedure. Results show that the function

directly related to the spectral energy κ(ω, T ), tends to u/((P -1) √ P ) when u tends to ∞.

Graphical TOC Entry


Keywords: efficacité énergétique, routage économe en énergie, réseaux mobiles ad hoc, réseaux de capteurs, OLSR, protocole de routage, durée de vie du réseau, relais multipoint

Etude comparative de différentes stratégies de routage efficaces en énergie et basées sur OLSR Résumé : L'efficacicité énergétique est un problème majeur des réseaux sans fil ad hoc et des réseaux de capteurs. Plusieurs axes ont été explorés pour maximiser la durée de vie d'un réseau, parmi eux, le routage économe en énergie. Dans ce rapport, nous montrons comment étendre le protocole de routage OLSR, standardisé à l'IETF, pour améliorer son efficacité énergétique. Pour prendre en compte l'énergie résiduelle des noeuds, l'algorithme natif de sélection des relais multipoint dans OLSR doit être modifié. Nous évaluons trois variantes et choisissons la meilleure. L' extension d'OLSR que nous proposons sélectionne le chemin minimisant l'énergie consommée par la transmission de bout-en-bout d'un paquet du flux et évite les noeuds avec une faible énergie résiduelle. Comme par ailleurs, il a été montré que router un flux sur deux routes préserve l'énergie, nous comparons cette extension avec un routage par la source maintenant deux routes par flux (les routes sont à liens disjoints ou à noeuds disjoints). Une évaluation de performances nous permet de conclure que notre proposition maximise la durée de vie du réseau d'une part et la quantité de données remises d'autre part.

Motivation

The diversity of the applications supported by wireless sensor networks explain the success of this type of network. These applications concern as various domains as environmental monitoring, wildlife protection, emergency rescue, home monitoring, target tracking, exploration mission in hostile environments... Sensor nodes are characterized by a small size, a low cost, an advanced communication technology, but also a limited amount of energy. This energy can be very expensive, difficult or even impossible to renew. That is why, energy efficient strategies are required in such networks in order to maximize network lifetime.

Network lifetime can be defined as:

Definition D1: the time to the first node failure due to battery outage;

Definition D2: the time to the unavailability of an application functionnnality. For instance, a vital parameter of a patient is no longer controlled, a target is no longer tracked. Definition D2 differs from D1 insofar as redundancy is provided: if for instance any area to be monitored is covered by k sensors, the failure of k -1 sensors is perfectly tolerated.

Definition D3: the time to the first network partitioning. As soon as the network is no longer connected, vital information can no longer be transferred to its destination.

In the absence of knowledge of the application supported by the network, definitions D1 and D3 are the most useful ones to compare different energy efficent strategies. In the following, we use definition D3 to evaluate the network lifetime.

Solutions to maximize network lifetime can be classified into four categories. Some adjust the transmission power of wireless nodes like [START_REF] Ingelrest | Optimal Transmission Radius for Energy Efficient Broadcasting Protocols in Ad Hoc Networks[END_REF][START_REF] Cardei | Topology Control in Ad hoc Wireless Networks with Hitch-hiking[END_REF]; others reduce the volume of information transferred by means of aggregation like [START_REF] Kalpakis | Maximum Lifetime Data Gathering and Aggregation in Wireless Sensor Networks[END_REF][START_REF] Tang | Extending Network Lifetime for Precision-Constrained Data Aggregation in Wireless Sensor Networks[END_REF][START_REF] Kulik | Adaptive Protocols for Information Dissemination in Wireless Sensor Networks[END_REF]; other too make nodes sleep in order to spare energy like [START_REF] Cardei | Energy-efficient target coverage in wireless sensor networks[END_REF][START_REF] Carle | Energy-Efficient Area Monitoring for Sensor Networks[END_REF][START_REF] Cardei | Improving wireless sensor network lifetime through power aware organization[END_REF]; finally the last ones use energy efficient routing, [START_REF] Heinzelman | Energy-efficient communication protocol for wireless microsensor networks[END_REF][START_REF] Xia | Near-optimal node clustering in wireless sensor networks for environment monitoring[END_REF][START_REF] Shah | Energy Aware Routing for Low Energy Ad Hoc Sensor Networks[END_REF][START_REF] Ganesan | Highly-resilient,energyefficient multipath routing in wireless sensor networks[END_REF][START_REF] Srinivas | Minimum Energy Disjoint Path Routing in Wireless Ad-Hoc Networks[END_REF][START_REF] Senouci | Energy efficient routing in wireless ad hoc networks[END_REF][START_REF] Kwon | Energy-Efficient Interference-Based Routing for Multi-hop Wireless Networks[END_REF][START_REF] Hassanein | Reliable Energy Aware Routing In Wireless Sensor networks[END_REF][START_REF] Shresta | Reception Awarness for Energy Conservation in Ad Hoc Networks[END_REF]. In this paper, we focus on solutions belonging to this last category, for multiple reasons: a multihop transmission is energy consuming and reducing the energy spent in the transmission of a packet from its source to its destination would increase network lifetime; avoiding nodes with a low residual energy would also contribute to prolong network lifetime; optimizing network flooding would also reduce the number of transmissions needed to reach all nodes in the network and therefore spare node energy; avoiding nodes that already have a high traffic load would reduce medium access contention, collisions if the medium access type is CSMA-CA and then spare energy lost in useless transmissions.

INRIA

The OLSR routing protocol has been standardized at IETF [START_REF] Adjih | Optimized Link State Routing Protocol[END_REF]. However, this protocol does not take into account the energy. We propose to extend this protocol in order to make it energy efficient. The results that are shown in this paper are valid for any wireless network meeting the following constraints: Assumption 1: Power conditions: P idle < P receive and P idle < P transmit , where P idle (resp. P transmit , P receive ) is the power used by a wireless node in the idle (resp. transmit, receive) state.

An IEEE 802.11 network meets those conditions. The performance evaluation reported in Sections 4.2 and 4.3 are based on values taken from Lucent Wavelan silver cards.

For simplicity reasons, we also assume that: Assumption 2: Radio interferences are limited to two hops. This assumption is generally admitted in ZigBee and IEEE 802.11 networks. This paper is organized as follows. In Section 2, we present a brief state of the art dealing with energy efficiency. We then focus more particularly on the OLSR routing protocol and more generally on energy efficient routing. In Section 3, we define our energy consumption model. In Section 4, we show how to extend OLSR to make it energy efficient. This extension, called RE selects the path minimizing the energy consumed in the end-to-end transmission of a flow packet and avoids nodes with low residual energy. To take into account residual node energy, the native selection of multipoint relays of OLSR is changed. Three selection algorithms based on the minimum residual energy are evaluated and the best one, M 1E is chosen. As it has been shown that two-path routing is energy efficient, we compare the RE extension with two two-path source routing strategies (DL with different links and DN with different nodes). We show that RE outperforms both DN and DL in terms of network lifetime and delivery rate. Finally, we conclude in Section 5 and give some perspectives for our further work.

State of the art

As said in Section 1, energy efficient routing is a way to maximize network lifetime. Before proposing an extension of OLSR routing to make it energy efficient, we will first recall the main principles of OLSR and then briefly present a state of the art related to energy efficient routing.

The OLSR routing protocol

OLSR (Optimized Link State Routing), [START_REF] Adjih | Optimized Link State Routing Protocol[END_REF], is a proactive routing protocol where nodes periodically exchange topology information in order to establish a route to any destination in the network. It is an optimization of a pure link state routing protocol, based on the concept of multipoint relays (MPRs). First, using multipoint relays reduces the size of the control messages: rather than declaring all its links in the network, a node declares only the set of links with its neighbors that have selected it as "multipoint relay". The use of MPRs also minimizes flooding of control traffic. Indeed only multipoint relays forward control messages. This technique significantly reduces the number of retransmissions of broadcast messages.

OLSR consists in two main functionalities: Neighborhood discovery. Each node acquires the knowledge of its one-hop and two-hop neighborhood by periodic Hello messages. It independently selects its own set of multipoint relays (MPRs), among its one-hop neighbors in such a way that its MPRs cover (in terms of radio range) all its two-hop neighbors.

Topology dissemination. Each node also maintains topological information about the network obtained by TC (Topology Control) messages broadcast by MPR nodes.

Each node computes its routing table by the Dijkstra algorithm. This table provides the shortest route (i.e. the route with the smallest hop number) to any destination in the network. In [START_REF] Plesse | OLSR performance measurement in a military mobile ad-hoc network[END_REF], we reported performance evaluation results showing that a MANET with OLSR routing achieves very satisfying performances.

Complexity results

Because of radio interferences, while a node N is transmitting a frame to a node D, no other node in the transmission range of N can simultaneously receive another frame. Similarly, while node D is receiving this frame, no other node in its transmission range can simultaneously transmit another frame. To make it simple, no node in the interference area of a transmitter can meanwhile transmit or receive another frame. As already said, we assume that interferences are limited to two-hops from the transmitter. The two following complexity results related to radio interferences have been established:

Result 1: Because of radio interferences, the selection of a unicast path, between a source and a destination, meeting the requested bandwidth is NP-hard : see [START_REF] Allard | Bandwidth Reservation in Multihop Wireless Networks: complexity, heuristics and mechanisms[END_REF].

Result 2: Because of radio interferences, the selection of a unicast path, between a source and a destination,ensuring that each node has sufficient residual energy is NP-hard : see [START_REF] Mans | On the Complexity of Reducing the Energy Drain in Multihop Ad Hoc Networks[END_REF] where it is also shown that this result is still valid, even if interferences are limited to a single hop.

Energy efficient routing protocols

Some routing protocols organize wireless nodes into clusters, such as Leach [START_REF] Heinzelman | Energy-efficient communication protocol for wireless microsensor networks[END_REF]. The author of [START_REF] Xia | Near-optimal node clustering in wireless sensor networks for environment monitoring[END_REF] has established the conditions under which such protocols are energy efficient and determined the optimal radius of a cluster. In the following, we assume that such conditions are not met.

INRIA

Existing energy efficient routing protocols can be first distinguished by the number of paths maintained to a destination: a single path or multiple paths. Multipath routing protocols, such as [START_REF] Shah | Energy Aware Routing for Low Energy Ad Hoc Sensor Networks[END_REF] or [START_REF] Ganesan | Highly-resilient,energyefficient multipath routing in wireless sensor networks[END_REF], have the advantage of sharing load of any flow on several paths, leading to a lesser consumption on the nodes of the selected paths. It has been shown in [START_REF] Srinivas | Minimum Energy Disjoint Path Routing in Wireless Ad-Hoc Networks[END_REF] that two paths with different links are generally sufficient. In Section 4, we will evaluate the performance obtained by multipath routing, when two paths are used with two variants: different links or different nodes.

We can distinguish three families of energy efficient routing protocols: the protocols selecting the path consuming the minimum energy. The advantage is that each transmission of a packet from its source to its destination minimizes the energy consumed. We can cite for example [START_REF] Senouci | Energy efficient routing in wireless ad hoc networks[END_REF] and a more sophisticated protocol [START_REF] Kwon | Energy-Efficient Interference-Based Routing for Multi-hop Wireless Networks[END_REF] where the selected path minimizes the additional energy dissipated by the routing of the new flow, taking into account the SINR and the energy lost in interferences. However, such protocols use always the same nodes (those minimizing the energy consumed) without any consideration on their residual energy. Consequently, these nodes will exhaust their battery more quickly than the others and the network lifetime is not maximized.

the protocols selecting the path visiting the nodes with the highest residual energy, such as [START_REF] Hassanein | Reliable Energy Aware Routing In Wireless Sensor networks[END_REF]. Each flow is ensured to have enough energy on the selected path: depleted nodes are avoided. However, the path selected does not minimize the energy needed to transmit a flow packet from its source to its destination. Hence, the network lifetime may not be maximized.

the hybrid protocols selecting the path with the minimum cost, where the cost takes into account the resiudual energy of each visited node (and possibly its neighbors) and the energy consumption of a packet on this path. These protocols avoid the problems encountered by the protocols of the two previous categories by weighing the factors used in the cost computation. We can cite for instance [START_REF] Shresta | Reception Awarness for Energy Conservation in Ad Hoc Networks[END_REF]. In Section 4, we will present different OLSR extensions belonging to this category.

Energy consumption model

A wireless node's radio can be in one of the following four states: T ransmit, Receive, Idle or Sleep and each of which consumes different levels of energy.

Transmit: node is transmitting a frame with transmission power P transmit ;

Receive: node is receiving a frame with reception power P receive . This frame can be decoded by this node or not, it can be intended to this node or not; Idle (listening): Even when no messages are being transmitted over the medium, the nodes stay idle and keep listening the medium with P idle ; Sleep: when the radio is turned off and the node is not capable of detecting signals: no communication is possible. The node uses P sleep that is largely smaller than any other power.

As we are interested in the additional energy spent during the transmission of a flow packet from its source to its destination with reference to the idle state, the values used for transmission power and reception power are calculated as follows:

P trans = P transmit -P idle , P rcv = P receive -P idle .

In Table 1, we report the reference values of P transmit , P receive , and P idle taken from a Lucent silver wavelan PC card. These values are used in the performance evaluation reported in Section 4. Now, we can determine the energy dissipated in transmitting (E trans ) or receiving (E rcv ) one packet (this energy is the additional energy consumed when transmitting or receiving packet with reference to the idle state). Let Duration denote the transmission duration of a packet. We have:

State

E trans = P trans * Duration, E rcv = P rcv * Duration.
When a transmitter transmits one packet to next hop, because of the shared nature of wireless medium, all its neighbors receive this packet even it is intended to only one of them. Moreover, each node situated between transmitter range and interference range receives this packet but it cannot decode it. These two problems generate loss of energy. So to compute the energy dissipated by one transmission, we must take into account these losses as follows [START_REF] Allard | Evaluation of the energy consumption in MANET Adhoc-Now 2006[END_REF]:

cost transmission (i) = E trans + n * E rcv ,
where n represents the number of non-sleeping nodes belonging to the interference zone of the transmitter i.

We now compute the energy cost of a flow dissipated on its path P .

cost(f low) = i∈sender(f low) cost transmission (i),
where i is a sender of f low on its path P . This cost indicates the quantity of energy consumed by one packet of the flow to reach its destination. Our idea to maximize network lifetime, is to minimize the energy consumed by the flow in selecting the best path with minimum energy dissipated. For that, we will use this cost to compute routes for flows. Indeed, instead of using the number of hop between source and destination to select the best route (i.e.; every link has a cost of one), as done in OLSR, we will use cost(f low) as the criterion to choose the best path, where every link i → j has a cost equal to cost transmission (i).

INRIA 4 Routing strategies

Goal and architecture

In this section, we show how to extend the OLSR routing protocol to make it energy efficient. Hence, our goal is to maximize the network lifetime by: Minimizing the energy consumed by a packet transmission from its source to its destination;

Avoiding the nodes with a low residual energy;

Avoiding the nodes with a small available bandwidth;

Reducing the overhead.

Hence, energy must be taken into account in route computation. Indeed, energy must be a criterion of route selection. We can notice that with OLSR, the intermediate nodes in a path are MPRs. As a consequence, the selection of MPRs must also take energy into account. For this purpose, energy information must be included in the topology information exchanged in OLSR. Hence, we get the following architecture, illustrated in Figure 1. Notice that the energy efficient routing we propose can be combined with a solution allowing nodes to sleep: the very short naps of nodes at the MAC level are kept transparent to the routing protocol.

Energy Routing

MAC with Energy

Signaling

Energy efficient selection of MPRs

Presentation of EMPR selection policies

We now focus on an energy efficient selection of MPRs: the multipoint relays are selected according to the residual energy and are denoted EMPRs. We study three variants depending on the nearby nodes considered for the computation of the minimum residual energy, where E R (M ) denotes the residual energy of node M :

The E policy considers only E R (M ), the residual energy of the EMPR candidate, M ;

The M 1E policy considers the weighted residual energy of the EMPR candidate M and its 1-hop neighbors:

min( E R (M ) Ptrans +Prcv , min D∈1hop(M ) ( E R (D) 2 * Prcv )).
The weights of ER(M ) and ER(D) take into account the role played by the nodes M and D in a transmission from N , the node performing the EMPR selection, to D, one of its two-hop neighbors, via the node M . It represents the maximum transmission duration that can be sustained;

The M 2E policy considers the weighted residual energy of the EMPR candidate M and its 1-hop and 2-hop neighbors:

min( E R (M ) Ptrans+Prcv , min D∈1hop(M ) ( E R (D) 2 * Prcv ), min D∈2hop(M ) ( E R (D) Prcv )).
The EMPR selection algorithm, performed by any node N , is the following:

Node N puts in the set U ncovered all its two-hop neighbors.

Node N sorts its 1-hop neighbors by decreasing order of the selection criterion (e.g.; E R (M ) for E). Let N 1 be this ordered set.

Node N selects M the first node in N 1 and removes the nodes covered by M from the set U ncovered.

If the successor of M in N 1 covers nodes that are in U ncovered, this node is selected as EMPR, and all nodes covered by it are removed from the set U ncovered;

and so on until U ncovered becomes empty.

Control messages in the OLSR extension

As previously said, no additional message is required in the OLSR energy extension we propose. In order to select the EMPRs, the Hello messages include the following additional information:

the residual energy of the sending node, in case of E, the residual energy of the sending node, and the minimum residual energy of the one-hop nodes, in case of M 1E, the residual energy of the sending node, the minimum residual energy of the one-hop nodes and the minimum residual energy of the two-hop nodes, in case of M 2E.

INRIA

In order to compute the energy cost of a flow, we need to know the number of nodes up to two-hop of the node considered (see Section 3). Hence, the T C messages include the number of nodes belonging to the interference area of the T C originator.

Performance evaluation

When a QoS (Quality of Service) criterion, such as energy, local available bandwidth or medium access delay, is used to select QoS MPRs (i.e.; EMPRs when the energy criterion is used), we have the following result:

Result 3: the average number of neighbors selected as QoS MPRs is in O(n 1/3 log(n)), where n is the average number of neighbor per node, whereas the average number of MPRs selected according to the native procedure in OLSR is in O(n 1/3 ): see the proof in [START_REF] Nguyen | Anlysis of MPR selection in the OLSR protocol[END_REF]. The average number of QoS MPRs selected per node is increased by a log(n) factor.

Simulation parameters: In the following simulations, nodes are uniformly distributed in the network area. The transmission range is equal to 250m. Interferences are limited to 500m. The initial energy of nodes is uniformly distributed in the interval [20J, 60J]. The powers used are those given in Table 1. Results reported in this section are the average of five simulations. First, we evaluate the average number of EMPRs per node as a function of network density, with the three selection variants. The number of nodes is set to 100. Similarly, the network density (i.e.; the average number of one-hop neighbors per node) being fixed to 10, we study the impact of the node number on the average number of EMPRs per node. In both cases, the native MPR selection is used as a reference. Simulation results are illustrated on Figures 2 and3. Simulation results confirm Result 3: the number of EMPR per node is higher than the number of MPR per node. We can also notice that the difference increases with network density and does not depend on node number. Furthermore, with the power values taken, the M 1E and M 2E selections give the same results. The reason is that the criterion used to sort the nodes leads to the same value with both M 1E and M 2E. Indeed the third term in formula 2 (see Section 3) is negligible compared to the first two terms. So the use of M 1E is more interesting than M 2E because it is less complicated to compute and needs less information from the network. OLSR forwarding rule: a node forwards once a broadcast message with a non-null time-to-live only if it has received this message for the first time from a node that has selected it as MPR.

Figure 4 depicts the number of retransmissions per TC, as a function of network density. Similarly, Figure 5 depicts the number of TCs received per node in a TC period, as a function of network density. These numbers are a good indication of the overhead induced by the extension of the OLSR routing protocol. As expected, the number of retransmissions of a flooded message is higher when EMPRs are used, instead of MPRs. Similarly, the number of TCs received per node in a TC period is higher with EMPRs. As we want to reduce the overhead, we recommend to use two types of multipoint relays:

INRIA

EMPRs to build energy efficient routes: all intermediate nodes on a path to a destination are EMPRs. Moreover, as in our case, both M 1E and M 2E provide the same results and the criterion used in M 1E is simpler to compute than in M 2E and requires less additional information in the Hello messages, we recommend the use of M 1E. This will be confirmed in the next section.

Energy efficient routing

In this section, we present three energy efficient routing strategies based on the OLSR protocol and evaluate their performance in terms of network lifetime, delivery rate and residual energy distribution. For the three strategies, the cost of a route is equal to the energy consumed by the transmission of a packet from its source to its destination as computed in Section 3. Each node computes its route of minimum energy towards any other node in the network, using the Dijkstra algorithm, with the energy cost.

Presentation of different routing policies

The first policy we will study is one hop-by-hop energy efficient routing, denoted RE, where each node forwards the received packet toward the next hop on the minimum cost route. In order to avoid frequent route changes, the selection of EMPRs is changed only when the topology changes or the ratio

E(new EM P R)-E(old EM P R) E(old EM P R) > 10%.
Another policy that can be used to spare energy is to share the load among the nodes. Multipath routing has been introduced for that purpose. In this paper, we will also study this policy, assuming multipath source routing. The source is in charge of computing the two paths used by its flow. For each flow packet, the source selects one of the paths inversely proportionally to its cost. The selected path will be encapsulated in the packet. Hence, each visited node forwards the packet to the next hop on the path selected by the source. We distinguish two variants: two paths with disjoint links, denoted DL. The source of the flow computes a first path minimizing the cost, with Dijkstra algorithm. It then removes all the links used to compute the second path minimizing the cost with this new topology. Notice however that a same node (see node 2 in Figure 6) can belong to both paths and then deplete early.

two path with disjoint intermediate nodes, denoted DN . As previously, the source of the flow computes a first path minimizing the cost, with Dijkstra algorithm. It then removes all the intermediate nodes used and then computes the second path minimizing the cost with this new topology.

INRIA

Evaluation of RE with different EMPR selections

We first compare the network lifetime obtained with RE, using E, M 1E and M 2E EMPR selection policies. Results are shown on Figure 7. First, we observe that M 2E and M 1E EMPR selection policies provide the same result. We notice that RE obtains the best performance when used with M 1E or M 2E. That is because, unlike E, M 1E takes into acount the energy dissipated in packet transmission and reception on the 1-hop neighbors of the transmitter. As M 1E is simpler than M 2E, the EMPR selection is done according to M 1E, for all the energy efficient routing strategies evaluated. In this second series of experiments, we compare the three energy efficient routing strategies RE, DN and DL with regard to the network lifetime, when the number of nodes ranges from 50 to 200, the network density is set to 10. The OLSR routing protocol that does never take energy into account and always selects the shortest path (i.e.; the path with the minimum hop number), is used as a reference. Simulation results are illustrated in Figure 8. As expected, OLSR provides the smallest network lifetime. This shows that the selection of the shortest path is not sufficient to save energy. Concerning the two 2-path source routing strategies, DN provides better results than DL. This is not surprising insofar as energy is dissipated per nodes and not per wireless link. Hence, DL that allows common nodes in the two paths can exhaust INRIA the energy of these common nodes more quickly. The main conclusion of these experiments is that RE significantly outperforms DN and DL whatever the number of nodes. Moreover, the gain is increasing with the network size. RE prolonges the network lifetime of 50% compared with OLSR for a network of 200 nodes. Notice that in the same conditions, DN prolonges the network lifetime of only 10%. Indeed, the two paths chosen by the source of the flow are used for all flow packets independently of the residual energy of these nodes. So the intermediates nodes exhaust their energy more quickly.

Evaluation of the amount of data delivered with RE, DN and DL versus OLSR Notice that the selection of an energy efficient routing protocol maximizing network lifetime, would provide no advantage to the application, if this increase in network lifetime was not followed by an increase in the amount of user data delivered. That is why, in this third series of experiments, we evaluate the delivery rate with the three routing strategies RE, DN and DL, when the number of nodes ranges from 50 to 200, the network density being set to 10. This delivery rate is compared with this provided by OLSR. Results are given in Figure 9. As previously, RE provides the best delivery rate, followed by DN , DL and finally OLSR. In conclusion, RE allows a real benefit to the application by delivering a higher amount of data with the same initial energy.

Distribution of residual energy

In this fourth series of experiments, we measure the residual energy at each node at a given time, (time = 400s corresponding to the network lifetime obtained with OLSR), for a network of 200 nodes with a network density of 10. We recall that the initial energy of a node is randomly selected in the interval [START_REF] Allard | Bandwidth Reservation in Multihop Wireless Networks: complexity, heuristics and mechanisms[END_REF]60] Joules and is represented in violet on Figure 10. We compare the values obtained by OLSR and RE. We observe that at time 400s, 86 nodes have a residual energy less than 20 Joules. With RE, half of the nodes have an energy in [START_REF] Allard | Bandwidth Reservation in Multihop Wireless Networks: complexity, heuristics and mechanisms[END_REF]40] Joules. This can be explained by: the EMPR selection that takes into account the residual energy of nodes.

the selection of the route with the minimum energy cost. As RE outperforms the other routing policies with regard to both network lifetime and delivery rate, we decide to select it, with M 1E, the best variant of EMPR selection. That is the solution we recommend for an energy efficient routing based on OLSR.

Conclusion

Wireless ad hoc and sensor networks are faced with the problem of energy efficiency in order to maximize network lifetime. The aim of this paper was to extend the OLSR routing protocol to make it energy efficient. We have defined an energy model and express its validity conditions. We have studied different variants of multipoint relay selection based on residual energy. The variant M 1E takes into account the energy dissipated in transmission and reception up to one-hop from the transmitter. M 1E presents the best tradeoff between the energy consumed and the overhead induced. Moreover, we recommend to keep the native MPR selection to optimize network flooding and to use the M 1E selection to build energy efficient routes.

We have also compared different routing strategies: RE, a hop-by-hop adaptive strategy selecting the route with minimum energy cost while avoiding nodes with a low residual energy, INRIA DL, a two-path routing strategy with different links, DN , a two-path routing strategy with different nodes.

Simulation results show that RE outperforms the other two strategies with regard to both network lifetime and delivery rate. Moreover, RE tends to maintain a higher residual energy level at each node. For these reasons, we recommend to use this strategy. In a further work, we will see how to use this strategy in a network where nodes are allowed to switch to the sleeping state for short periods.

Figure 1 :

 1 Figure 1: An energy efficient architecture In order to support the energy efficient extension, no additional message is required. Only additional energy information is included in both Hello and T C messages, as shown in Section 4.2.2.

Figure 2 :

 2 Figure 2: Number of multipoint relays per node, function of network density.

Figure 3 :

 3 Figure 3: Number of multipoint relays per node, function of node number.

Figure 4 :

 4 Figure 4: Number of retransmissions per TC, function of network density.

Figure 5 :

 5 Figure 5: Number of TCs received per node, function of network density.

Figure 6 :

 6 Figure 6: Multipath routing

Figure 7 :

 7 Figure 7: Impact of the EMPR selection policy on the network lifetime Evaluation of network lifetime with RE, DN and DL versus OLSR

Figure 8 :

 8 Figure 8: Comparison of network lifetime with RE, DN and DL versus OLSR

Figure 9 :

 9 Figure 9: Comparison of delivery rate with RE, DN and DL versus OLSR

Figure 10 :

 10 Figure 10: Distribution of residual energy at time 200s

  The OLSR routing protocol . . . . . . . . . . . . . . . . . . . . . 2.2 Complexity results . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Energy efficient routing protocols . . . . . . . . . . . . . . . . . . Goal and architecture . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Energy efficient selection of MPRs . . . . . . . . . . . . . . . . . 4.2.1 Presentation of EMPR selection policies . . . . . . . . . . 4.2.2 Control messages in the OLSR extension . . . . . . . . . . 4.2.3 Performance evaluation . . . . . . . . . . . . . . . . . . . 4.3 Energy efficient routing . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Presentation of different routing policies . . . . . . . . . . 4.3.2 Performance evaluation . . . . . . . . . . . . . . . . . . .

	Contents
	1 Motivation
	2 State of the art
	2.1 3 Energy consumption model
	4 Routing strategies
	4.1

Table 1 :

 1 Power value in each radio state

		Power value	Power Increment
	Transmit P transmit = 1.3W P trans = 0.56W
	Receive	P receive = 0.9W	P rcv = 0.16W
	Idle	P idle = 0.74W
	Sleep	P sleep = 0.047W

RR n 6374


Florentin Millour
email: fmillour@oca.eu

Denis Mourard

Julien Woillez

Philippe Berio

Andrea Chiavassa

Orlagh Creevey

Eric Lagadec

Marc-Antoine A Martinod

Anthony Meilland

Nicolas Nardetto

Karine Perraut

Philippe Stee

Perspectives of a visible instrument on the VLTI

Keywords: high angular resolution, visible interferometry, fundamental stellar parameters, stellar imaging

Introduction

The last decade of Optical Interferometry has seen an important expansion of the astrophysical results coming mainly from the VLTI [START_REF] Glindemann | VLTI technical advances: present and future[END_REF] and the CHARA arrays [START_REF] Ten Brummelaar | First Results from the CHARA Array. II. A Description of the Instrument[END_REF]. At the same time, the expansion of the user's community has progressed with the preparation of the second generation instruments GRAVITY [START_REF] Gravity Collaboration | [END_REF] and MATISSE [14] on one side, and the opening of CHARA to the wider community on the other side. Both arrays will be equipped soon with adaptive optics systems on all their telescopes.

While science programs have continued to be well diversified, it is also clear that significant efforts have been devoted to the determination of fundamental stellar parameters, with an important complementarity with asteroseismology, spectroscopy, photometry, and now astrometry [START_REF] Perryman | GAIA: Composition, formation and evolution of the Galaxy[END_REF]. In a recent work [START_REF] Stee | Science cases for a visible interferometer[END_REF], Stee et al. produced a review of the potential science cases for optical interferometry at visible wavelengths. With the launch of new space missions dedicated to stars and planets (TESS, CHEOPS, PLATO), optical interferometry has an important role to play in bringing direct measurements of stellar parameters. They will permit one to correctly disentangle the planet parameters from the potential stellar activity. It is clear however that not only the best possible angular resolution will be necessary to correctly sample the selected stars, but also good imaging quality will be required to correctly characterize the stellar surface features. The long baselines of the CHARA Array and the dense configurations of the VLTI are thus very complementary.

In [START_REF] Mourard | SPICA, Stellar Parameters and Images with a Cophased Array: a 6T visible combiner for the CHARA array[END_REF] [START_REF] Mourard | SPICA, Stellar Parameters and Images with a Cophased Array: a 6T visible combiner for the CHARA array[END_REF], we presented our perspective for a new 6telescope beam combiner for the CHARA Array. In the present paper, we describe our ideas to develop a visible instrument for the VLTI: Sec. 2 presents the most promising science cases, Sec. 3 details the main instrumental requirements that will be necessary to achieve our science goals, and Sec. [START_REF] Chiavassa | Radiative hydrodynamics simulations of red supergiant stars. II. Simulations of convection on Betelgeuse match interferometric observations[END_REF] presents a study of the potential samples of targets in the different classes.

Unique science that can be done at VLTI in the visible

In 2012, O. Chesneau started a prospective work in the Action Spécifique Haute Résolution Angulaire (ASHRA) from the Centre National de la Recherche Scientifique (CNRS). The initial idea was to make a past, present and future appraisal of the type of science that can be addressed with visible interferometry. This work culminated with the publication of the visible interferometry white book [START_REF] Stee | Science cases for a visible interferometer[END_REF].

High spatial resolution is the key to understanding various astrophysical phenomena. But even with the future ELT, single dish instruments are limited to a spatial resolution of about 4 mas in the visible. As a reference point, most stellar photospheres in our Galaxy remain smaller than 1 mas.

We started exploiting this whitebook with the SPICA project on the CHARA array. SPICA will provide many improvements compared to now: 6 telescopes, full use of the new adaptive optics, use optical fibers to perform spatial filtering, fringe tracking, and the use of new low noise EMCCD cameras. It will focus on two main science cases: an angular diameter survey of all observable stars up to the magnitude m R = 10, and an imaging survey program to produce surfaces images in the visible of giant and supergiant stars. This instrument will benefit from two unique features of CHARA : the longest operating baselines in the world (330 m), and 6 telescopes at full throttle. It is also considered as a testbench for the technologies of a future visible VLTI visible instrument.

Indeed, the VLTI has very complementary characteristics compared to CHARA, which makes it very appealing for hosting a visible instrument: the telescopes are larger (1.8 m and 8 m), there are shorter baselines (starting at 8 m), allowing a better sampling of small spatial frequencies, and the array is denser (due to the movable auxiliary telescopes). The Paranal site has also an excellent median seeing of 0.8" in the visible domain, and clear skies all over the year. In addition, the arrival of the second-generation instruments of the VLTI triggered a thorough assessment and improvement of the VLTI infrastructure performance [START_REF] Woillez | The 2nd generation VLTI path to performance[END_REF][START_REF] Woillez | VLTI status update: Three years into the second generation[END_REF], making relevant the perspective of a visible addition to the VLTI instruments suite.

In summary, while CHARA reaches a higher angular resolution, the VLTI has an advantage for imaging larger targets (AT and short baselines), and fainter targets (ATs are bigger than CHARA telescopes, and UTs even more). In the current paper, we focus on the most promising science cases for development on the VLTI, taking into account the overall characteristics of the interferometric array at Cerro Paranal.

Mass loss from evolved stars and galactic chemical enrichment

The chemical enrichment of our Galaxy is strongly influenced by Asymptotic Giant Branch (AGB) and Red Supergiant (RSG) stars, which supply heavy elements into the ISM through a mass-loss mechanism not yet fully understood. This mass loss is certainly the result of a two-step process. First, the dust is formed behind shocks in an atmosphere extended by pulsation and/or convection. Then, radiation pressure on dust and/or molecules triggers the mass loss, in a complex mechanism that can involve magnetic fields.

Their surface is covered by a few convective cells [START_REF] Chiavassa | Radiative hydrodynamics simulations of red supergiant stars. II. Simulations of convection on Betelgeuse match interferometric observations[END_REF][START_REF] Montargès | The convective surface of the red supergiant Antares. VLTI/PIONIER interferometry in the near infrared[END_REF], evolving on timescales of several weeks to years [START_REF] Chiavassa | Radiative hydrodynamic simulations of red supergiant stars. III. Spectro-photocentric variability, photometric variability, and consequences on Gaia measurements[END_REF], and causing large temperature inhomogeneities. Their velocity can levitate the gas and contribute to the mass loss [START_REF] Chiavassa | Radiative hydrodynamics simulations of red supergiant stars. IV. Gray versus non-gray opacities[END_REF]. Imaging of the surface of these stars is a key aspect to understand their dynamics and the link to the mass-loss process [START_REF] Paladini | The VLTI/MIDI view on the inner mass loss of evolved stars from the Herschel MESS sample[END_REF], with the possibility to resolve spatially and spectroscopically lines features in the photosphere.

The visible domain is still not fully explored at the high angular and spectral resolutions available with interferometers. Getting an interferometer equipped with a spectral resolution of about 30,000 on a wide bandwidth would enable studying large scale movements at the surface of the star, temperature stratification due to convection, and the kinematics of the cells. In addition, the comprehensive interpretation of their surface dynamics, induced by pulsation and convection, requires multi-epochs observations. With interferometry, one can determine the geometries of these stars, and the VLTI will clearly excel at it with its movable telescopes and imaging capabilities on ≈ 10 mas targets. More and more observations indicate that giant stars and their descendant are not spherical. Determining the geometry of a substantial number of AGB stars is crucial to obtain accurate radiative transfer models and thus estimations of the mass-loss rate. This will will have to be a vital part of a community-large long-term multi-wavelength observing project aiming at determining the quantitative impact of AGB stars on the enrichment of the ISM.

Finally, a visible instrument will enable finding companions to AGB stars, and the VLTI is again very well suited for that with its 4 movable telescopes. Key binary parameters (mass and separation) can be extracted to quantitatively determine how binaries affect the mass loss from AGB stars and the Galaxy's chemical enrichment.

A & B stars binarity below 30 au

De Rosa et al. [START_REF] De Rosa | The VAST Survey -IV. A wide brown dwarf companion to the A3V star ζ Delphini[END_REF] have conducted a large volume-limited (D < 75 pc) multiplicity survey of A type stars, sensitive to companions beyond 30 au, by using a combination of adaptive optics imaging and a multi-epoch common proper motion search. They reach the conclusion that 69.9 ± 7.0% of A-type stars have companions. Studies done on OB stars using AO and interferometry [START_REF] Sana | Southern Massive Stars at High Angular Resolution: Observational Campaign and Companion Detection[END_REF] conclude to similar results.

Knowing the importance of hot stars for extragalactic stellar astrophysics (for example through early type eclipsing binaries in distant galaxies [START_REF] Riess | A 2.4% Determination of the Local Value of the Hubble Constant[END_REF][START_REF] Vilardell | The distance to the Andromeda galaxy from eclipsing binaries[END_REF]), and the fact that A-type stars have not yet been probed extensively with interferometry, we propose to go below this limit of 30 au with the VLTI (down to 1 au) to search for companions to A-type stars with the large baselines available at the VLTI array. This will permit us to correctly calibrate surfacebrightness-color relationships of early-type stars, on which relies the distance determination of early-type eclipsing binaries [START_REF] Challouf | Improving the surface brightness-color relation for earlytype stars using optical interferometry[END_REF], and to remove potential biases due to the multiplicity of objects.

The main instrumental requirements for this program will be to reach magnitude up to m R = 10, long baselines (from 8 to 130 m) and to be able to do fast snapshot observations with frequent revisits for the determination of the orbits. Therefore, the VLTI is uniquely suited to reach faint magnitudes and long baselines. Combining this survey with the Gaia data [START_REF] Perryman | GAIA: Composition, formation and evolution of the Galaxy[END_REF], we will be also able to determine an important number of stellar masses. Combining a spectrograph with an interferometer was one of the innovation of visible interferometry in the 70's [START_REF] Labeyrie | Interference fringes obtained on VEGA with two optical telescopes[END_REF]. This combination has proven to be extremely efficient at characterizing circumstellar environments where the gas is excited by a central ionizing source [START_REF] Mourard | The rotating envelope of the hot star Gamma Cassiopeiae resolved by optical interferometry[END_REF], or masking underlying continuum emission [START_REF] Ohnaka | Vigorous atmospheric motion in the red supergiant star Antares[END_REF]. Indeed, spectral lines analysis with an interferometer give access to the geometry, kinematics, chemical and physical conditions in the circumstellar medium at scales of a fraction of millisecond of arc.

In young stellar objects (YSOs), like Herbig or T Tauri stars, this corresponds to potential planet-forming regions. While optical direct imaging and sub-millimetric images revealed complex structures around YSOs (10 auscale spirals, gaps, holes), reaching the au-scale regions, where key star-disk-protoplanet(s) interactions processes occur, requires sub-mas angular resolution, only provided by visible interferometry. Improving the sensitivity of visible interferometry and the ability to spatially resolve emission and absorption lines would be of strong interest, as one would gain in angular resolution and probe lines complementary to the infrared hydrogen Brγ line. We can cite the hydrogen recombination lines (Hα, Hβ, ...) that are considered as accretion tracers, He I 1083 nm that provides a unique diagnostic of kinematic motion in the regions close to the star, or forbidden lines (for instance [OI], [NII], [FeII], and [SII]) that trace low-density gas that is located on relatively large spatial scales.

In Cepheid stars, this corresponds to where lies the still-misunderstood circumstellar environments (CSEs). Indeed, a few of them observed by longbaseline interferometry show faint (few %) and close (3 R * ) CSEs in the infrared. Characterizing these CSEs, and in particular their radial brightness profile in several spectral bands will enable one to unbias the period-luminosity relation. Recently, δ Cep was observed with the VEGA/CHARA instrument and an unexpected visible CSE contributing to 7% to the total flux was discovered [START_REF] Nardetto | VEGA/CHARA interferometric observations of Cepheids. I. A resolved structure around the prototype classical Cepheid δ Cep in the visible spectral range[END_REF]. Many more visible CSEs can be discovered with a more sensitive visible instrument. Combining visible and infrared observations of Cepheids is probably one of the keys to characterize their environment.

In evolved emission-line stars, like Be stars or B[e] stars, this corresponds to where the disk is formed. This excretion disk is still a puzzle, as its exact formation process is not yet well understood. Non-radial pulsations, as well as the possible presence of a companion star, could be the culprits to trigger the ejection of matter around these stars. Visible interferometry can not-only fully characterize the gaz disk, but also detect unseen companion stars, and resolve the stellar disk itself. Spectrally dispersed interferometry can be also applied to fast-rotating stars. The detection and characterisation of fast and/or non-solid rotation of these objects could provide decisive information on the inner structure of these stars, in which the internal redistribution of matter and evolution is strongly influenced by rotation-induced processes (e.g. gravity darkening, non-spherical mass loss, and internal circulation of mass, chemical elements, and angular momentum)

In O, B or A-type hot stars, this is where their radiative winds form. Emission lines seen in the visible can be resolved with the longest VLTI baselines, giving a direct access to previously-unconstrained physical parameters of the stellar wind. It is worth noting that, today, the VLTI is the only interferometer able to reach the southern sky and it exists many emission-line star studies in that region of the sky.

Stellar masses and fundamental parameters of Main Sequence stars and Giants

A direct determination of the radius of the star is an important measurement complementing effective temperature and luminosity. These three properties are fundamental, and with only two independent, the third provides a method to test systematic errors or to disentangle less-well-known variables that play a role in their determination e.g. line-of-sight absorption, distance, bolometric correction, surface gravity. Having accurate determinations of these properties along with better constraints on their correlations results directly in a much more accurate and precise determination of the star's mass and evolution stage.

In particular, complementary information from, for example, asteroseismology, allows one to provide a direct estimate of the absolute stellar age, because the mass can be derived with much higher accuracy.

Combining interferometry and asteroseismology has been possible for some stars using data for example from CoRoT, Kepler and K2 [START_REF] Kervella | The diameter and evolutionary state of Procyon A. Multi-technique modeling using asteroseismic and interferometric constraints[END_REF][1] [START_REF] White | Beyond the Kepler/K2 bright limit: variability in the seven brightest members of the Pleiades[END_REF]. In a near future, TESS (NASA) and Plato (ESA) will provide seismic data on brighter sources and for a much larger number.

Focusing on the brighter stars allows effectively to determine an interferometric diameter. In addition, some of the stars whose diameters will be constrained will have an orbiting exoplanet detected by TESS. The constraints provided by interferometry will give the diameter and the age of the exoplanet.

On this particular science case, the VLTI will complete the SPICA survey (Northern hemisphere) with Southern stars, and especially those close to the galactic center.

Instrumental requirements

The high-level instrumental requirements that could be set with respect to the science goals presented above are the limiting magnitude to reach representative samples of stellar classes, the overall efficiency for rapid imaging and survey programs, the possibility to reach dense (u,v) coverage for imaging, and finally high to very high spectral resolution. Recent examples of stellar programs [START_REF] Vedova | VLTI/PIONIER images the Achernar disk swell[END_REF][START_REF] Kluska | SPARCO : a semi-parametric approach for image reconstruction of chromatic objects. Application to young stellar objects[END_REF] are driving our considerations. The instrument is proposed as a dispersed fringes system with low (R=300) and medium (R=3000) spectral resolutions. A high (R=30000) spectral resolution Echelle spectrograph adapted to interferometric data is also considered. The instrument will take benefit of the installation of adaptive optics on all the telescopes and the light of each beam will be coupled into single mode fibers. This design is based on the experience gained on different instruments and on a prototype called FRIEND that we have installed on the CHARA Array (Martinod et al., 2018, in press).

Two main milestones are focusing our efforts. The first one is the optimization of the injection into single mode fibers at visible wavelengths with partial correction by adaptive optics. Coupling efficiency of at least 25% should be reached and the stability of the coupling should be guarantied. Indeed combining 3, 4 or 6 beams for interferometric measurements suppose that the stability of the coupling on each beam is very high if one wants to reach a good interferometric efficiency. Photometric drops to 0 should be avoided, and optimized command-laws of the deformable mirrors should be implemented in this respect. The second milestone is the fringe stabilization to permit long exposure and thus to reach fainter magnitudes or lower fringe visibilities. The recent successes of the Gravity instrument [START_REF] Gravity Collaboration | [END_REF] are very encouraging. As a consequence, we study a fringe sensor similar to the Gravity one, adapted to shorter wavelengths.

A high spectral resolution mode has been already used in the VEGA/CHARA instrument [START_REF] Mourard | VEGA: Visible spEctroGraph and polArimeter for the CHARA array: principle and performance[END_REF][START_REF] Mourard | Spatio-spectral encoding of fringes in optical long-baseline interferometry. Example of the 3T and 4T recombining mode of VEGA/CHARA[END_REF] but with a limitation on the spectral bandwidth due to the sampling characteristics of this instrument. With a single-mode instrument in mind, we are considering an extension of the spectral band to a large fraction of the visible domain, by using the principle of Echelle spectroscopy coupled with the necessary anamorphic system required for the fringe sampling. The goal is to reach a spectral resolution R=30000, while keeping a band ∆λ = 300 nm centered at 700 nm. The resolution element is δλ = 0.02 nm and the sampling in the spectral direction should then be 0.01 nm. A 6-telescope instrument generates 15 frequencies (in a linear non-redundant configuration) and each baseline will generate 4 fringes sampled over a minimum of 4 × 3 pixels. With all these numbers at hand, we can deduce that 17 orders of an Echelle spectrograph could fit in a 4K x 4K detector, as presented in Fig. 1. It is clear however that such a detector could not be read at the usual frequency in optical interferometry (of the order of the coherence time of the atmosphere) and thus long exposures of the order of 1s will be mandatory and are possible. This means that this mode should only be considered with the assistance of an active fringe tracker, which will put even more stringent requirements on the adaptive optics system.

The current VLTI is optimised for the infrared and therefore not ready to support a visible instrument yet. The adaptive optics beam splitters are dichroic mirrors optimised to only send the infrared light beyond 1 µm to the VLTI instruments, setting aside visible light for the wavefront sensors of the adaptive optics systems. Permanently replacing the existing dichroic to send more visible light to the VLTI would affect the limiting magnitude of the adaptive optics and impact the red science carried out by the existing infrared instruments. The use of a reconfigurable dichroic in place of the current static one seems like the most promising solution. The installation of such a system in the existing star separators may turn out challenging, due to the crowded environment of the Coudé focus on both UTs and ATs.

The limited order of correction of the existing adaptive optics systems is an even greater challenge. Neither MACAO on the UTs nor NAOMI on the ATs have sub-apertures small enough to deliver a reasonable Strehl in the visible. Supporting our visible instrument with an infrared fringe tracker could make a low order adaptive optics system sufficient for fringe tracking, but the transmission to the visible instrument would remain dramatically low. On the ATs, the NAOMI adaptive optics system has the particularity of using an off-the-shelf mirror with an actuator count much larger than its wavefront sensor sub-aperture count. The actuator spacing of ∼16 cm is well adapted to the visible. However, in order to make full use of this deformable mirror, one would need to design a matching wavefront sensor. This could be done using the now decommissioned single-field Relay Optics Structure (ROS), the interchangeable interface between the telescopes and the tunnels to the delay lines. The dual-field ROS is the one currently in operation at VLTI; it contains the star separators and the NAOMI wavefront sensor. The nominal design of a visible compatible ROS would be based on a 11x11 wavefront sensor and require an upgrade to the real-time controller to support the correction order and bandwidth increase. It is worth mentioning that such a high order system on a small telescope in addition to supporting visible operation, would deliver an extremely high Strehl in the near-infrared, and especially at L band, in support to the high contrast Hi-5 endeavour [START_REF] Defrère | The path towards high-contrast imaging with the VLTI: the Hi-5 project[END_REF]. For the 8 m UTs, the upgrade would be more extensive and also include the deformable mirror itself, since the MACAO deformable mirror, and its ∼1 m inter-actuator spacing, is not adapted to the visible.

Study of the potential targets

In order to study the number of potential targets accessible in the sky of Paranal (-90 • < δ < +20 • ), we made use of the SIMBAD database with different criteria in terms of angular diameters and magnitude in the R band.

A visible instrument addition will enable the VLTI to measure the angular diameter of several hundreds or thousands of stars. Indeed, an interferometer of 150 m baseline at 550 nm wavelength can fully resolve a star (i.e. attain the first zero of the Bessel function of a uniform disk) of 0.9 m diameter. Stars smaller than this diameter can get their diameters measured, depending on the accuracy of the interferometer. We consider, based on the VEGA experience, that stars down to 0.5 m can get their diameters measured. Crossing this information (and declination ≤ 20 • ) with the JSDC catalog of predicted angular diameters (which contains itself 465,877 stars diameters), provides the number of ≈ 21000 stars that can be resolved in the visible with the VLTI, given it can reach a limiting magnitude of m R < 8. However, this huge number has to be reduced depending on the sensitivity of the considered instrument. A m R < 5-limited instrument can resolve ≈ 2000 stars (i.e. almost all naked-eye visible stars), while a m R < 3-limited instrument can resolve ≈ 300 stars.

The results are presented in Fig. 2 in the top line. The figure titles give numbers for each stellar type subset. This study shows that several hundreds of main sequence FGKM stars can get their diameters measured with a m R < 5limited instrument on the VLTI already, with an over-representativity of Kand M-types stars. A m R < 8-limited instrument will outperform this by an order of magnitude, and the limit will rather come from available observing time than instrumental capacities.

On the other hand, imaging capabilities need the full angular resolution of a star. Therefore, the available sample is purposely limited to stars larger than 2 m angular diameter (i.e. have at least 2 resolution elements on the star disk). As can be seen in Fig. 2 in the bottom line, this restricts the stellar types to mainly K and M stars, i.e. mostly AGB and RSG stars. To dig a bit further, we selected known RSG and AGB stars, and splitted them as a function of possible spectral resolution that can be used to observe them (the higher the spectral resolution, the worse the sensitivity for a given exposure time). The result of this study is presented in Fig. 3. It shows that ≈ 250 RSG and AGB stars can be imaged with medium spectral resolution (≈ 1000), and that a minimum of ≈ 56 RSG and AGB stars can be observed in high spectral resolution (≈ 30000).

Conclusions

In this paper we have shortly presented our preliminary study of the potential for a visible interferometric instrument for the VLTI. Apart from additional considerations on the VLTI infrastructure for visible wavelengths, it is clear that the science goals presented above will require the assistance of a fringe sensor aiming at stabilizing the fringes in the visible. More detailed studies of the science programs will also be necessary in the future in order to correctly assess the need for the different spectral resolutions presented here. With our preparatory work on the CHARA Array, we continue to study critical aspects of this kind of instrument: optimization of the injection in the fibres, compensation of the birefringence, and preliminary study of a H-Band fringe sensor inspired from the Gravity one. We also plan to progress in the coming years on the Echelle spectrograph (SPINE) design in order to offer a wide spectral band at very high spectral resolution. This will clearly open new frontiers in stellar astrophysics.

2. 3

 3 Gaseous circumstellar environments : YSOs, Be and B[e] stars, OBA winds, Cepheids, ...

Fig. 1

 1 Fig. 1 Numerical simulation of a 4k detector fed with an Echelle spectrograph accommodating 17 orders of dispersion, each order being anamorphosed in order to correctly sample the high frequencies generated by the interferometric recombination of 6 pupils.

Fig. 2

 2 Fig. 2 Top row: HR diagram for stars with 0.5 < θ(mas) < 10 and m R < 8/5 from left to right. Bottom row: same but for 2 < θ < 10. The range of declination is [-90 • ; +20 • ].

Fig. 3

 3 Fig. 3 Red giants and AGB stars that could be imaged with the visible instrument

Acknowledgements This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n • 730890). It was ade possible thanks to the OPTICON project, the Optical Infrared Coordination Network for Astronomy.

We thank A. Domiciano de Souza for reading the paper and suggesting changes to the content.