Joseph Lehec

Regularization in L 1 for the Ornstein-Uhlenbeck semigroup

Let γn be the standard Gaussian measure on R n and let (Qt) be the Ornstein-Ulhenbeck semigroup. Eldan and Lee recently established that for every non-negative function f of integral 1 and any time t the following tail inequality holds true:

where Ct is a constant depending on t but not on the dimension. The purpose of the present paper is to simplify parts of their argument and to remove the (log log r) 4 factor.

Résumé

Soit γn la mesure Gaussienne standard sur R n et soit (Qt) le semigroupe d'Ornstein-Ulhenbeck. Eldan et Lee ont montré récemment que pour toute fonction positive f d'intégrale 1 et pour temps t la queue de distribution de Qtf vérifie γn({Qtf > r}) ≤ Ct (log log r) 4 r √ log r , ∀r > 1 où Ct est une constante dépendant seulement de t et pas de la dimension. L'objet de cet article est de simplifier en partie leur démonstration et d'éliminer le facteur (log log r) 4 .

Introduction

Let γ n be the standard Gaussian measure on R n and let (Q t ) be the Ornstein-Ulhenbeck semigroup: for every test function f

Q t f (x) = R n
f e -t x + 1 -e -2t y γ n (dy).

(1) * CEREMADE (UMR CNRS 7534) Université Paris-Dauphine.

1

Nelson [START_REF] Nelson | The free Markoff field[END_REF] established that if p > 1 and t > 0 then Q t is a contraction from L p (γ n ) to L q (γ n ) for some q > p, namely for q = 1 + e 2t (p -1).

The semigroup (Q t ) is said to be hypercontractive. This turns out to be equivalent to the logarithmic Sobolev inequality (see the classical article by Gross [START_REF] Gross | Logarithmic Sobolev inequalities[END_REF]).

In this paper we establish a regularity property of Q t f assuming only that f is in L 1 (γ n ).

Let f be a non-negative function satisfying

R n f dγ n = 1,
and let t > 0. Since Q t f ≥ 0 and

R n Q t f dγ n = R n f dγ n = 1,
Markov inequality gives

γ n ({Q t f ≥ r}) ≤ 1 r ,
for all r ≥ 1. Now Markov inequality is only sharp for indicator functions and Q t f cannot be an indicator function, so it may be the case that this inequality can be improved. More precisely one might conjecture that for any fixed t > 0 (or at least for t large enough) there exists a function α satisfying lim r→+∞ α(r) = 0 and

γ n ({Q t f ≥ r}) ≤ α(r) r , (2)
for every r ≥ 1 and for every non-negative function f of integral 1. The function α should be independent of the dimension n, just as the hypercontractivity result stated above. Such a phenomenon was actually conjectured by Talagrand in [START_REF] Talagrand | A conjecture on convolution operators, and a non-Dunford-Pettis operator on L 1[END_REF] in a slightly different context. He conjectured that the same inequality holds true when γ n is replaced by the uniform measure on the discrete cube {-1, 1} n and the Orstein-Uhlenbeck semigroup is replaced by the semigroup associated to the random walk on the discrete cube. The Gaussian version of the conjecture would follow from Talagrand's discrete version by the central limit theorem. In this paper we will only focus on the Gaussian case.

In [START_REF] Ball | L 1 -smoothing for the Ornstein-Uhlenbeck semigroup[END_REF], Ball, Barthe, Bednorz, Oleszkiewicz and Wolff showed that in dimension 1 the inequality (2) holds with decay

α(r) = C √ log r ,
where the constant C depends on the time parameter t. Moreover the authors provide an example showing that the 1/ √ log r decay is sharp. They also have a result in higher dimension but they loose a factor log log r and, more importantly, their constant C then tends to +∞ (actually exponentially fast) with the dimension. The deadlock was broken recently by Eldan and Lee who showed in [START_REF] Eldan | Regularization under diffusion and anti-concentration of temperature[END_REF] that (2) holds with function

α(r) = C (log log r) 4 √ log r ,
with a constant C that is independent of the dimension. Again up to the log log factor the result is optimal.

In this article we revisit the argument of Eldan and Lee. We shall simplify some steps of their proof and short cut some others. As a result, we are able to remove the extra log log factor. We would like to make clear though that this note does not really contain any new idea and that the core of our argument is all Eldan and Lee's.

Main results

Recall that γ n is the standard Gaussian measure and that (Q t ) is the Ornstein-Uhlenbeck semigroup, defined by [START_REF] Ball | L 1 -smoothing for the Ornstein-Uhlenbeck semigroup[END_REF]. Here is our main result.

Theorem 1. Let f be a non-negative function on R n satisfying R n f dγ n = 1 and let t > 0. Then for every r > 1

γ n ({Q t f > r}) ≤ C max(1, t -1 ) r √ log r ,
where C is a universal constant.

As in Eldan and Lee's paper, the Ornstein-Ulhenbeck semigroup only plays a rôle through the following lemma. Lemma 2. Let f : R n → R + . For every t > 0, we have

∇ 2 log(Q t f ) ≥ - 1 2t id, pointwise.
Proof. This is really straightforward: observe that (1) can be rewritten as

Q t f (x) = (f * g 1-ρ )(ρ x),
where ρ = e -t and g 1-ρ is the density of the Gaussian measure with mean 0 and covariance (1 -ρ)id. Then differentiate twice and use the Cauchy-Schwarz inequality. Details are left to the reader.

What we actually prove is the following, where Q t does not appear anymore.

Theorem 3. Let f be a positive function on R n satisfying R n f dγ n = 1. Assume that f is smooth and satisfies

∇ 2 log f ≥ -β id (3)
pointwise, for some β ≥ 0. Then for every r > 1

γ n ({f > r}) ≤ C max(β, 1) r √ log r ,
where C is a universal constant.

Obviously, Theorem 3 and Lemma 2 altogether yield Theorem 1.

Let us comment on the optimality of Theorem 1 and Theorem 3. In dimension 1, consider the function

f α (x) = e αx-α 2 /2 .
Observe that f α ≥ 0 and that R f α dγ 1 = 1. Note that for every t ≥ 1 we have

γ 1 ([t, +∞)) ≥ c e -t 2 /2 t ,
where c is a universal constant. So if α > 0 and r ≥ e then

γ 1 ({f α ≥ r}) ≥ c exp -1 2 log r α + α 2 2 log r α + α 2 .
Choosing α = √ 2 log r we get

γ 1 ({f α ≥ r}) ≥ c r √ log r .
Since (log f α ) = 0 this shows that the dependence in r in Theorem 3 is sharp.

Actually this example also shows that the dependence in r in Theorem 1 is sharp. Indeed, it is easily seen that

Q t f α = f αe -t ,
for every α ∈ R and t > 0. This implies that f α always belongs to the image Q t . Of course, this example also works in higher dimension: just replace f α by

f u (x) = e u,x -|u| 2 /2
where u belongs to R n .

Theorem 4. Let X be a random vector having density f with respect to the Gaussian measure, and assume that f satisfies (3). Then for every r > 1

P (f (X) ∈ (r, e r]) ≤ C max(β, 1) √ log r .
Theorem 4 easily yields Theorem 3.

Proof of Theorem 3. Let G be standard Gaussian vector on R n and let X be a random vector having density f with respect to γ n . Then using Theorem 4

P[f (G) > r] = +∞ k=0 P f (G) ∈ (e k r, e k+1 r] ≤ +∞ k=0 (e k r) -1 E f (G) 1 {f (G)∈(e k r,e k+1 r]} = +∞ k=0 (e k r) -1 P f (X) ∈ (e k r, e k+1 r] ≤ +∞ k=0 (e k r) -1 C max(β, 1)
log(e k r)

≤ C e e -1 1 r max(β, 1) √ log r ,
which is the result.

The rest of the note is devoted to the proof of Theorem 4.

Preliminaries: the stochastic construction

Let µ be a probability measure on R n having density f with respect to the Gaussian measure. We shall assume that f is bounded away from 0, that f is C 2 and that ∇f and ∇ 2 f are bounded. A simple density argument shows that we do not lose generality by adding these technical assumptions.

Eldan and Lee's argument is based on a stochastic construction which we describe now. Let (B t ) be a standard n-dimensional Brownian motion and let (P t ) be the associated semigroup:

P t h(x) = E[h(x + B t )],
for all test functions h. Note that (P t ) is the heat semigroup, not the Ornstein-Ulhenbeck semigroup. Consider the stochastic differential equation

X 0 = 0 dX t = dB t + ∇ log(P 1-t f )(X t ) dt, t ∈ [0, 1]. ( 4
)
The technical assumptions made on f insure that the map

x → ∇ log P 1-t f (x)
is Lipschitz, with a Lipschitz norm that does not depend on t ∈ [0, 1]. So the equation ( 4) has a strong solution (X t ). In our previous work [START_REF] Lehec | Representation formula for the entropy and functional inequalities[END_REF] we study the process (X t ) in details and we give some applications to functional inequalities. Let us recap here some of these properties and refer to [4, section 2.5] for proofs.

Recall that if µ 1 , µ 2 are two probability measures, the relative entropy of µ 1 with respect to µ 2 is defined by

H(µ 1 | µ 2 ) = log dµ 1 dµ 2 dµ 1 , if µ 1 is
v t = ∇ log P 1-t f (X t ).
Using Itô's formula it is easily seen that

d log P 1-t f (X t ) = v t , dB t + 1 2 |v t | 2 dt.
Therefore, for every t ∈ [0, 1]

P 1-t f (X t ) = exp t 0 v s , dB s + 1 2 t 0 |v s | 2 ds . (5)
Combining this with the Girsanov change of measure theorem one can show that the random vector X 1 has law µ (again we refer to [START_REF] Lehec | Representation formula for the entropy and functional inequalities[END_REF] for details). Moreover we have the equality

H(µ | γ n ) = 1 2 E 1 0 |v s | 2 ds . (6)
Also, if (u t ) is any drift and if ν is the law of

B 1 + 1 0 u t dt, then H(ν | γ n ) ≤ 1 2 E 1 0 |u s | 2 ds . (7)
So the drift (v t ) is in some sense optimal. Lastly, and this will play a crucial rôle in the sequel, the process (v t ) is a martingale. Eldan and Lee introduce a perturbed version of the process (X t ), which we now describe. From now on we fix r > 1 and we let

T = inf{t ∈ [0, 1], P 1-t f (X t ) > r} ∧ 1
be the first time the process (P 1-t f (X t )) hits the value r, with the convention that T = 1 if it does not ever reach r (note that our definition of T differs from the one of Eldan and Lee). Now given δ > 0 we let (X δ t ) be the process defined by

X δ t = X t + δ T ∧t 0 v s ds.
Since T is a stopping time, this perturbed process is still of the form Brownian motion plus drift:

X δ t = B t + t 0 (1 + δ1 {s≤T } )v s ds.
So letting µ δ be the law of X δ 1 and using [START_REF] Talagrand | A conjecture on convolution operators, and a non-Dunford-Pettis operator on L 1[END_REF] we get

H(µ δ | γ n ) ≤ 1 2 E 1 0 (1 + δ 1 {s≤T } ) 2 |v s | 2 ds = 1 2 E 1 0 |v s | 2 ds + δ + δ 2 2 E T 0 |v s | 2 ds . ( 8
)
4 Proof of the main result

The proof can be decomposed into two steps. Recall that r is fixed from the beginning and that X δ 1 actually depends on r through the stopping time T . The first step is to prove that if δ is small then µ and µ δ are not too different.

Proposition 5. Assuming (3), we have

d T V (µ, µ δ ) ≤ δ (β + 1) log r,
for every δ > 0, where d T V denotes the total variation distance.

The second step is to argue that f (X δ 1 ) tends to be bigger than f (X 1 ). An intuition for this property is that the difference between X δ 1 and X 1 is somehow in the direction of ∇f (X 1 ). Proposition 6. Assuming (3), we have

P f (X δ 1 ) ≤ r 1+2δ e -4 ≤ P(f (X 1 ) ≤ r) + (β + 4)δ 2 log(r),
for all δ > 0.

Remark. Note that both propositions use the convexity hypothesis (3).

It is now very easy to prove Theorem 4. Since X 1 has law µ, all we need to prove is

P(f (X 1 ) ∈ (r, er]) ≤ C max(β, 1) √ log r .
We choose δ = 5 2 log r .

almost surely. We shall use this inequality several times in the sequel. Recall that X 1 has law µ and that X δ 1 has law µ δ . Taking expectation in the previous inequality and using (10) we get

H(µ δ | µ) ≤ H(µ δ | γ n ) -H(µ | γ n ) -δ E T 0 v 1 , v s ds + βδ 2 2 E T 0 |v s | 2 ds .
Together with ( 6) and ( 8) we obtain

H(µ δ | µ) ≤ -δ E T 0 v 1 -v s , v s ds + (1 + β)δ 2 2 E T 0 |v s | 2 ds .
Now since (v t ) is a martingale and T a stopping time we have

E v 1 , v s 1 {s≤T } = E |v s | 2 1 {s≤T }
for all time s ≤ 1. This shows that the first term in the previous inequality is 0.

To bound the second term, observe that the definition of T and the equality (5) imply that

T 0 v s , dB s + 1 2 T 0 |v s | 2 ds ≤ log r,
almost surely. Since (v t ) is a bounded drift, the process ( t 0 v s , dB s ) is a martingale. Now T is a bounded stopping time, so by the optional stopping theorem

E T 0 v s , dB s = 0.
Therefore, taking expectation in the previous inequality yields

E T 0 |v s | 2 ds ≤ 2 log r,
which concludes the proof.

Proof of Proposition 6

The goal is to prove that P f (X δ 1 ) ≤ r 1+2δ e -4 ≤ P(f (X 1 ) ≤ r) + δ 2 (β + 4) log r.

Obviously P f (X δ 1 ) ≤ r 1+2δ e -4 ≤ P(f (X 1 ) ≤ r) + P f (X δ 1 ) ≤ r 1+2δ e -4 ; f (X 1 ) > r
Now recall the inequality (11) coming for the weak convexity of log f and rewrite it as log f (X δ 1 ) ≥ K 1 + 2δK T + Y where (K t ) is the process defined by

K t = log(P 1-t )(f )(X t ) = t 0 v s , dB s + 1 2 t 0 |v s | 2 ds,
and Y is the random variable

Y = -2δ T 0 v s , dB s + δ T 0 v 1 -v s , v s ds - βδ 2 2 T 0 |v s | 2 ds.
Recall that the stopping time T is the first time the process (K t ) exceeds the value log r if it ever does, and T = 1 otherwise. In particular, if

K 1 = log f (X 1 ) > log r then K T = log r. So if f (X 1 ) > r then f (X δ 1 ) > r 1+2δ e Y .
Therefore P f (X δ 1 ) ≤ r 1+2δ e -4 ; f (X 1 ) > r ≤ P(Y ≤ -4). So we are done if we can prove that P(Y ≤ -4) ≤ (β + 4)δ 2 log r.

(12)

There are three terms in the definition of Y . The problematic one is

δ T 0 v 1 -v s , v s ds.
We know from the previous section that it has expectation 0. A natural way to get a deviation bound would be to estimate its second moment but it is not clear to us how to do this. Instead we make a complicated detour.

Lemma 8. Let Z be an integrable random variable satisfying E[e Z ] ≤ 1. Then

P(Z ≤ -2) ≤ -E[Z].
Remark. Note that E[Z] ≤ 0 by Jensen's inequality.

Proof. Simply write

E e Z ≥ E e Z 1 {Z>-2} ≥ E (Z + 1) 1 {Z>-2} = E[Z] -E Z 1 {Z≤-2} + 1 -P(Z ≤ -2) ≥ E[Z] + P(Z ≤ -2) + 1. So if E[e Z ] ≤ 1 then P(Z ≤ -2) ≤ -E[Z].
Lemma 9. Let Z be the variable

Z = -δ T 0 v s , dB s + δ T 0 v 1 -v s , v s ds - (β + 1)δ 2 2 T 0 |v s | 2 ds.
Then P(Z ≤ -2) ≤ δ 2 (β + 1) log r.

Proof. As we have seen before the first two terms in the definition of Z have expectation 0 and

E[Z] = - (β + 1)δ 2 2 E T 0 |v s | 2 ds ≥ -δ 2 (β + 1) log r.
By Lemma 8 it is enough to show that E[e Z ] ≤ 1. To do so, we use the Girsanov change of measure formula. The process (X δ t ) is of the form Brownian motion plus drift:

X δ t = X t + δ T ∧t 0 v s ds = B t + t 0 (1 + δ1 {s≤T } )v s ds,
Note also that the drift term is bounded. Therefore, Girsanov's formula applies, see for instance [5, chapter 6] (beware that the authors oddly use the letter M to denote expectation). The process (D δ t ) defined by

D δ t = exp - t 0 (1 + δ 1 {s≤T } ) v s , dB s - 1 2 t 0 (1 + δ 1 {s≤T } )v s 2 ds
is a non-negative martingale of expectation 1 and under the measure Q δ defined by dQ δ = D δ 1 dP the process (X δ t ) is a standard Brownian motion. In particular

E[f (X δ 1 )D δ 1 ] = E[f (B 1 )] = 1.
Now use inequality (11) once again and combine it with equality (5) written at time t = 1. This gives exactly

f (X δ 1 ) D δ 1 ≥ e Z .
Therefore E[e Z ] ≤ 1, which concludes the proof.

We now prove inequality (12). The idea being that the annoying term in Y is handled by the previous lemma. Observe that Putting everything together we get (12), which concludes the proof.

Y = Z -δ

T 0 v 2 T 0 |v s | 2 0 v s , dB s ≥ 1 + P δ 2 2 T 0 |v s | 2 ds ≥ 1 .v s , dB s 2  = δ 2 E T 0 |v s | 2 δ 2 2 T 0 |v s | 2

 0202012021202202 s , dB s -δ 2 ds. So P(Y ≤ -4) ≤ P(Z ≤ -2) + P δ T Recall that T 0 v s , dB s has mean 0 and observe that ds ≤ 2δ 2 log r.So by Tchebychev inequalityP δ T 0 v s , dB s ≥ 1 ≤ 2δ 2 log r.Similarly by Markov inequality P ds ≥ 1 ≤ δ 2 log r.

  absolutely continuous with respect to µ 2 (and H(µ 1 | µ 2 ) = +∞ otherwise). Also in the sequel we call drift any process (u t ) taking values in R n which is adapted to the natural filtration of (B t ) (this means that u t depends only on (B s ) s≤t ) and satisfies

1 0 |u t | 2 ds < +∞, almost surely. Let (v t ) be the drift

Acknowledgments

We would like to thank Ronen Eldan and James Lee for introducing us to the problem and communicating their work to us. We are also grateful to Michel Ledoux for some fruitful discussions at the final stage of this work.

For this value of δ, Proposition 6 gives P(f (X δ 1 ) ≤ e r) ≤ P(f (X 1 ) ≤ r) + 25 4

β + 4 log r ,

whereas Proposition 5 yields

.

Combining the two inequalities we obtain

which is the result.

Remark. We actually prove the slightly stronger statement:

.

Proof of the total variation estimate

We actually bound the relative entropy of µ δ with respect to µ. Recall that log f is assumed to be weakly convex: there exists β ≥ 0 such that

pointwise.

Proposition 7. Assuming (9), we have

for all δ > 0.

This yields Proposition 5 by Pinsker's inequality.

Proof. Observe that


H Perron

C Domain

J Roques
email: roques@ipno.in2p3.fr

R Drot

E Simoni

H Catalette

Optimisation of Accurate Rutile TiO 2 (110), (100), (101) and (001) Surface Models from Periodic DFT Calculations Authors

Keywords: Rutile, TiO 2, DFT, surface

) and (001) faces were computed. The calculated surface energies are consistent with the natural rutile powder composition, even if a dependence on the number of layers of the slab used to model the surface was identified. Internal constraints, consisting in freezing some internal layers of the slab to atomic bulk positions, were thus added to mimic the bulk hardness in order to stabilize the computed surface energies for thinner systems. In parallel, the influence of pseudopotentials was studied and it appears that four valence electrons for titanium atoms are sufficient. The aim of this study was to optimise accurate rutile TiO 2 surfaces models that will be used in further calculations to investigate water and uranyl ion sorption mechanisms.

Introduction

For years, the titanium dioxide, TiO 2 , has been widely used as a white pigment and opacifier.

Its recent application in catalysis and its photocatalysis properties make it a substrate of great interest from an experimental point of view as well as on a theoretical one. In the nature, this compound can be found under three crystallographic phases which are in order of decreasing abundance: rutile, anatase and brookite. A lot of microscopic studies carried out under ultrahigh vacuum (UHV) conditions and with preliminary treatments (Ar + -ion bombardment, irradiation, high temperature) have been performed especially considering some selected faces of the rutile phase, allowing a detailed knowledge at the atomic level in terms of relaxation and reconstruction [1][2][3][4][5][6][7]. Then, many metals and metal oxides overlayers growth [8][9][10][11][12][13], organic and inorganic molecules adsorption [14][15][16][17][18][19][20][21], have been also studied and the surface chemistry on this phase has been abundantly investigated. Among the low index faces naturally present in the rutile phase powder, the (110) face was found as the most stable one and thus has been much more studied than the others. In the case of the anatase phase, several applications were found, for example in nanostructured solar cells. All these experimental investigations have besides been summarized by Diebold [22].

From a theoretical point of view, the rutile and anatase bulk phases and low index faces have been also investigated (see Diebold review article [22] and reference therein). Calculated intrinsic bulk parameters such as lattices parameters or bulk moduli were well reproduced [23]. The computed surface energies are also in agreement with the natural repartition of the different faces [24]. Adsorption studies also have been performed, especially on the rutile (110) face, with small organic molecules (H 2 O, CO, NH 3 , …) [25][26][27][28][29], atomic ions (Na + , K + , Ca 2+ , …) [30][31][32] or metallic atoms (Au, Cu, Ag, …) [33][34][35], but some works have given opposite results. As an example, a disagreement arised for the adsorption of the water molecule on the rutile (110) face which is still a matter of controversy. Some authors found a spontaneous dissociation [36][37] while others agreed with a molecular adsorption mechanism [26,38]. This behaviour seems coming from the different models used (thickness and methods) showing the real importance of an accurate description of these surfaces.

Experimental data suggests that the two forms may coexist simultaneously, and the relative proportions depend on several factors such as the initial oxygen vacancies density or pressure [39][40][41].

The regular increasing of the modern supercomputers power and the improvements of the methods allow to better describe the real systems. However, before studying solid/liquid interfacial reactions, a detailed knowledge of the surface, in terms of chemical and physical properties, should be beforehand established. Thus, it is of first importance to previously define well-characterised titanium dioxide TiO 2 surface models that could be then used as accurate models for studying complex adsorption processes. Among the previous ab initio calculations on TiO 2 rutile low index faces, several approaches have been used (Hartree-Fock or DFT in different formalisms) with multiple surface models. Nevertheless, to our knowledge, no systematic comparison of the different surfaces as a function of the slab thickness with the same formalism has been performed.

In this paper, the (110), (100), ( 101) and (001) TiO 2 rutile surface properties have been characterised within the same DFT formalism based on the plane waves methodology. The ultrasoft pseudopotential (USPP) and projector augmented wave (PAW) formalisms are compared, as well as the effect of the number of valence electrons, on bulk properties and surface energies. A comparison of unconstrained and constrained surface models and the effect of the number of layers needed to reach convergence are also presented. Then, the surfaces relative stabilities and their atomic relaxations are correlated to their relative unsaturation and their structures. Finally, accurate (110), (100), ( 101) and (001) constrained surface models have been optimised to study sorption processes on large systems with a limited number of layers in order to speed up calculations.

Computational Details

All DFT periodic calculations were performed using the Vienna ab initio Simulation Package, VASP 4.6 [42][43][44][45], in local density approximation (LDA) and generalized gradient approximation (GGA as defined by Perdew and Wang [46]) formalisms, for the exchangecorrelation energy evaluation. All atoms were described with pseudopotentials taken from the VASP library and developed on plane waves basis sets. Different pseudopotentials were used.

The first ones are based on Vanderbilt-type [47][48] pseudopotentials also known as ultrasoft pseudopotentials (USPP). Titanium atoms were described by ten valence electrons (3p 6 4s 2 3d 2 ) and oxygen ones by six electrons (2s 2 2p 4 ). This first set of pseudopotentials will be referred as US10 in the following. The second set is generated with the Projector Augmented Wave (PAW) method [49][50]. Two pseudopotentials were used for describing titanium atoms: the first one describes explicitly ten electrons (3p 6 4s 2 3d 2 ) and will be referred as PAW10; the second one, takes into account only four valence electrons (3d 2 4s 2 ), and will be noted PAW4.

For oxygen atoms, a PAW pseudopotential with six valence electrons (2s 2 2p 4 ) was used.

These two titanium PAW pseudopotentials were used in order to determine if the 3p electrons have to be included as semi-core in the pseudopotential to correctly describe the crystal properties. The Brillouin zone was integrated using the Monkhorst-Pack sets of k-points [51] (centered at the Γ point) depending on the supercell dimensions and of the number of atoms.

Results for bulk relaxations were checked for convergence with respect to the number of kpoints as well as energy cutoff. The density of states (DOS) calculations on the bulk were performed at the equilibrium volume using the tetrahedron method with Blöchl corrections for accuracy [52]. Except for calculating bulk parameters, all atomic relaxations were performed at constant volume (at the bulk equilibrium lattice parameters) by using the conjugate gradient optimisation scheme. All faces were built from direct bulk cleavage and thus exhibit undercoordinated atoms relative to their bulk structure. For each face, a sufficient vacuum thickness has been introduced between the slabs allowing to neglect their interaction.

Bulk Properties

This study started by optimising the rutile bulk parameters in order to build the faces and to determine the accuracy of the modelling by comparing calculated parameters to experimental ones. These first calculations were performed using different sets of k-points and energy cutoff to evaluate their effects on the bulk parameters and to optimise them. Relaxations were performed in two steps: the first one, with constant volume to optimise atom positions in the lattice, and the second one, with unconstrained volume, to determine equilibrium bulk parameters and reference energies. In addition, the anatase phase lattice parameters were also determined in order to compare cohesive energies with the rutile phase.

The bulk rutile unit cell is tetragonal (a = b = 4.587 Å, c = 2.954 Å, internal parameter x = 0.305 and c/a = 0.644 [23]), space group No 136. Titanium atoms are in 2a positions and O atoms in 4f positions using Wyckoff's notations. The primitive cell contains two TiO 2 units.

For the LDA calculations, the optimised parameters (see Table 1) were obtained with an optimised 5×5×5 k-point mesh and a 400 eV cutoff, whereas for the GGA calculations the energy cutoff was optimised to 350 eV with the same k-point sampling. For a given functional (LDA or GGA), the results obtained with US10 and PAW10 pseudopotentials are similar.

Moreover, only small differences can be noticed between PAW10 and PAW4 pseudopotentials showing that the four valence electrons are sufficient to correctly describe the titanium atoms and 3p semi-core electrons do not change significantly the results. Furthermore, this smaller number of electrons explicitly treated will allow to consider larger surfaces. By comparing to the experimental values, the LDA leads to a slight underestimation of the bulk parameters (around 0.04 and 0.03 Å on a and c parameters, respectively) while GGA gives an overestimation (0.04 and 0.02 Å on a and c parameters, respectively). Despite these little differences, there is a good agreement between calculated and experimental bulk parameters (1% error). Previous works from other groups gave similar results (see Table 1).

Therefore, these different parameters have been used in all the coming constant volume relaxations.

The anatase phase is also tetragonal (a = b = 3.79 Å, c = 9.51 Å, internal parameter z = 0.208, c/a = 2.509 [23]), space group No 141. Titanium atoms are in 4a positions and oxygen ones in 8e. The elementary cell contains four TiO 2 units. Using the same sets of k-points and energy cutoff as for rutile, the anatase cell was optimised (see Table 2). Here, calculated values were still in good agreement with experimental ones but it can be noticed that the c parameter is slightly overestimated in GGA that implies a largest error on the c/a ratio. This deviation has already been observed in the previous works reported in Table 2. Nevertheless, this error does not exceed 3% relatively to the experimental values.

The computed bulk moduli for anatase and rutile phases are also reported in Table 1 and2. The bulk modulus measures the response in pressure, or resistance to a uniform compression, due to a change in the volume relative to the equilibrium one:

2 2 V E V V P V B ∂ ∂ = ∂ ∂ - = , ( 1
)
where E is the total energy of the supercell as a function of its volume V , P is the pressure and B the bulk modulus evaluated at the minimum of E .

In agreement with experimental results, the rutile's bulk modulus was always found higher than the anatase one whatever the functional (LDA or GGA) or the pseudopotential (US10, PAW10 or PAW4) used. For a given phase and functional, the US10 and PAW10 pseudopotentials gave similar results whereas the PAW4 one gave slightly higher values. The LDA calculations always lead to 30-40 GPa higher bulk moduli than GGA ones which were in better agreement with experimental values (rutile phase: 210 GPa [53][54], anatase phase:

180 GPa [55][56]).

The cohesive energy of the crystal can be calculated according to Eq. 2:

2 TiO O O Ti Ti bulk coh N E N E N E E - - = , (2)
where Ti E is the reference energy for one titanium atom in the hexagonal compact crystal, O E is the reference energy for oxygen atom in an isolated dioxygen molecule in gas phase, Ti N , O N and 2

TiO

N are respectively the numbers of titanium atoms, oxygen atoms and TiO 2 units contained in the supercell and, finally, bulk E is the total energy of the supercell. Using this definition, a negative value is favourable. By comparing cohesive energies, the relative stability of these two phases can be calculated (see Table 3). The PAW4 pseudopotential gave about 10 meV energies lower than the two others pseudopotentials (PAW10 and US10). The relative energy stability is quite different when it is evaluated from LDA or GGA calculations.

This disagreement has already been observed and these differences come from the exchangecorrelation functionals used [57][58]. However, ab initio calculations always predict the rutile phase as the most stable one.

The energy gaps (difference between HOMO and LUMO) of the rutile and the anatase phases have been experimentally evaluated at 3.0 and 3.2 eV respectively [59]. The values calculated here, 1.6-1.9 eV and 1.9-2.3 eV for rutile and anatase respectively, are always lower (around 1-1.3 eV) relative to the experimental ones and depend strongly on the pseudopotential used.

However, this result was predictable since it is known that DFT always underestimate gap energies because it is based on the description of the fundamental state and thus is not the most appropriate formalism to described excited states. Comparing the values obtained for the bulk with these three pseudopotentials, the set of pseudopotentials noted PAW4 (Ti: 4s 2 3d 2 , O: 2s 2 2p 4 ) in this study gives accurate results.

Surfaces: Effect of the Slab Thickness

According to Jones et al. [60][61], it is known that rutile powder present three major crystallographic faces: (110), ( 101) and (100) with respectively a ratio of 60% / 20% / 20%.

Surface energies were evaluated using Eq. 3:

S N E E E TiO TiO slab surf 2 2 2 - = , (3)
where slab E is the total energy of the supercell, 2 TiO E the reference energy for a TiO 2 unit in bulk phase, 2 TiO N the number of TiO 2 unit in the supercell and S the surface area of one side of the slab depending on the considered face (see Table 4). For the four faces, a layer (L) has been defined as a plane containing titanium atoms. All atomic positions were able to relax during calculations.

Rutile (110) Face

The TiO 2 (110) face (Fig. 1-a) exhibits atoms with different environments. First, a pentacoordinated titanium atom, noted Ti (5), which is undercoordinated relative to its bulk structure. There are also two kinds of oxygen atoms, the first one is localised in the surface plane and is threefold coordinated (noted Os for "surface" oxygen); the second is prominent from the surface around 1 Å and is only doubly coordinated (noted Ob for "bridging" oxygen). All (110) surface energies were obtained using a 3×1×5 k-point grid (the system

being c vaccuum nL a × + × ) ( 2
). In Fig. 2, surface energies for systems from 1 to 13 layers are represented from LDA and GGA calculations with the three sets of pseudopotentials previously defined (US10, PAW10 and PAW4). It can be noted that the surface energy oscillates with the number of layers as already noted by Bates et al. [62] and Ramamoorthy et al. [24]. In a recent paper, Bredow et al. [63] have showed that the energy gap and the interlayer distances also oscillate with nL. However, the amplitude of these fluctuations decreases as the number of layers increases. Convergence was reached (within 0.02 J/m² for the three pseudopotentials) from the 10-12 layer systems. Referring to the converged surface energies, even nL gave smaller values and odd nL led to bigger ones. The differences between the LDA calculated surface energies and the GGA ones were important whatever the pseudopotential used (about 0.4 J/m²) but, for a given formalism, they gave similar results.

Atom displacements due to the surface relaxation were of various amplitude according to nL but remain in agreement with those experimentally observed [1][2] and calculated [24,26,62,64]. On the surface, the Ti-O bond lengths were modified relative to the bulk ones (two at 1.97 and one at 2.00 Å in GGA(PAW4)): d(Ti(6)-Ob) = 1.85 Å, d(Ti(6)-Os) = 2.05 Å. The Ti(5) atom falls down into the surface leading to a shortening of the remaining apical bond length at 1.83 Å. The calculated surface energies reported in the literature are always larger in GGA (0.73 [62], 0.81 [65] and 0.84 [66], in J/m²), while a better agreement was observed in LDA (0.89 [24], 0.84 [67] and 0.66 [68] in J/m²).

In agreement with the results on the bulk, the PAW4 pseudopotential gives similar results than the two other ones. Then, for the three following faces, (100), ( 101) and ( 001), the systematic study with the number of layers was performed only with the PAW4 pseudopotential.

Rutile (100) Face

The natural rutile powder is composed of about 20% of (100) face (Fig. 1-b). Although its "sawtooth" orientation, the (100) face exhibits exactly the same three kinds of atoms as the (110) face (see § IV.A). As for the (110) face, the surface energy was calculated with different number of layers (Fig. 3) with a 1×5×5 k-point mesh (with supercell dimensions:

c a vacuum nL × × + ) (
). The (100) surface energy also oscillates with nL and the convergence was reached within 0.01 J/m² only from 8-10 layer systems (0.69 J/m² and 1.18 J/m² in GGA and LDA with PAW4 respectively) even if the amplitude of the oscillations is smaller than for the (110) face. These results are in agreement with values reported in the literature: 1.12 [24],

1.38 [69] and 1.30 [70] J/m² in LDA and 0.83 [69] J/m² in GGA. The (100) surface energy was found to be higher than the (110) one which is in agreement with the fact that the (110) face is the most stable one and with previous theoretical works [24,69]. The bond lengths on the surface are also modified similarly to those observed on the (110) face. As previously noted for the (110) face, the calculated LDA surface energies were always larger than the GGA ones around 0.5 J/m² with the PAW4 pseudopotential. These calculations were also performed with the PAW10 pseudopotential and gave similar results.

Rutile (101) Face

The (101) face (Fig. 1-c) is the third face occurring naturally in rutile powder for about 20%.

This face looks like the (100) one, it exhibits two kinds of atoms, first, a pentacoordinated titanium, noted Ti(5); second, a twofold oxygen atom with two different Ti-O bond lengths, noted O(2). The surface energy as a function of the number of layers was also calculated and reveals a similar behaviour as for the other faces (Fig. 3), the k-point mesh used was 4×5×1

(with system dimensions:

) ( ² ² vacuum nL a c a + × × +
). The converged surface energy (1.03 J/m² in GGA and 1.47 in LDA with PAW4) was reached, within 0.01 J/m², from 6-8 layer systems. The (101) surface is also less table than the (110) one. These values are in agreement with Ramamoorthy et al.

[21] reporting 1.39 J/m² from LDA calculations and with the faces repartition in natural rutile powder. Surface energies evaluated from LDA calculations are still larger than the GGA ones around 0.4 J/m², as already noted for the (110) and the (100) faces.

The optimised Ti-O bond lengths on this surface are very different than those of the bulk (two at 1.97 and one at 2.00 Å in GGA(PAW4)): 1.84 and 1.89 Å for the twofold oxygen atoms while, for the threefold ones, 2.07 and 2.12 Å bond lengths were found. The PAW10 pseudopotential still gave similar results.

Rutile (001) Face

The dry (001) face (Fig. 1-d) does not belong to the three major faces present in powder however some experimental [22,70] and theoretical [24,[71][72][73][74] works performed on this face are available in the literature. Then, this fourth face was investigated in order to make comparison with the three previous ones. This face also exhibits two kinds of atoms: a tetracoordinated titanium atom, noted Ti(4), with a high unsaturated valence and thus very acidic; and a twofold oxygen atom noted O(2). As above, the evolution of the surface energy with an increasing number of layers was studied (Fig. 3) with a 5×5×1 k-point mesh (the

supercell being ) ( vacuum nL a a + × ×
). In this case, the space between two layers is very small (only 2 / c ≈ 1.5 Å) and, relative to the three other faces (around 3 Å), for an equivalent thickness, the number of layers is about twice than the three other ones. This explains why for the same number of layers, as for the three other faces, the surface energy is not converged and thus requires a higher number of layers. The fluctuations observed were as important as for the (110) one. Convergence within 0.1 J/m² was reached from 13-15 layer systems (1.25 J/m² and 1.76 J/m² in GGA and LDA with PAW4 respectively). The surface energy is larger than for the three previous faces in agreement with the fact that the (001) orientation has a minor contribution in the natural powder. This is also in agreement with the previous Ramamoorthy et al. [24] theoretical work giving in J/m²: 1.65 for the (001), 0.89 for the (110) and 1.12 for the (100). On this face, short Ti-O bond are also reported. The Ti(4) atom is surrounded by two oxygen atoms at 1.91 Å (O(2) on Fig. 1-d) and two others at 1.81 Å. As for the three previous faces, the surface energies calculated from LDA calculations are larger than for the GGA ones around 0.5 J/m². Calculations with PAW10 are also in agreement with results obtained with PAW4.

As a last test, the influence of the three pseudopotentials on the computed GGA surface energies (from the biggest and fully relaxed systems) was studied, results are summarized in Table 5. For the four faces, the calculated surface energies, with the three pseudopotentials, were always in the same order: ] directions. Therefore, the (100) face is thus, in a first approximation, equivalent to the (110) one with a 2 factor that is also found in the ratio of these two surface energies (see Table 5). Following these observations and those done previously for the bulk, it can be concluded that the PAW4 pseudopotential can correctly describe these systems. It has thus been used for all coming calculations.

Internal Constraints

Effects on the Surface Energies

Since the surface energies oscillate with nL, the models need to be described by slabs composed of a large number of layers to reach convergence. However, to study large surfaces, the number of layers needs to be the smallest as possible. In this way, internal constraints were added in the slabs. The surface energies and the atomic relaxations have been used to find the best compromise between size and accuracy. These constraints consist in freezing some of the most internal layers to their bulk atomic positions. They introduce the hardness of the bulk that isolates the relaxation of each side of the slab from one to the others and thus decrease size effects. There are several possibilities to add these constraints. As an example, for a five layers system, if only the most internal layer is frozen, the two external layers on each side of the slab will be able to relax, while if the three most internal layers are frozen, only the most external layer on each side will be able to relax. These systems are noted nL_m with n the number of layers and m the number of layers able to relax on each side. The surface energies of the totally constrained systems (noted Unrelaxed) were also calculated for comparison. Similar studies on partially constrained systems have already been performed by Thompson et al. [75] and Hameeuw et al. [76], on the (110) face only and focus more particularly on geometrical parameters. All the following calculations were performed in GGA with the PAW4 pseudopotential.

Surface energies of the constrained (110) face are reported in Fig. 4-a. For the most constrained systems (Unrelaxed), the surface energies are larger around 0.8-0.9 J/m² relative to the totally unconstrained ones (Fully relaxed), but the parity of nL has smaller effects. By unconstraining slightly these systems (nL_1), the surface energies decrease significantly around 0.6-0.7 J/m² and the oscillation with nL parity becomes smaller. By again unconstraining these systems (nL_2), the surface energies still decrease (around 0.1 J/m²) and stay quasi-stable with nL. For the last systems (nL_3), there is no significant difference relative to the nL_2 ones. For the three nL_m systems, surface energy convergence is here reached within 0.02 J/m² for 5-8 layer systems, even if these converged energies are larger than the one of the totally unconstrained systems.

For the (100) face (see Fig. 4-b), the energies of the unrelaxed systems are larger by more than 0.8 J/m² as for the (110) face. The progressive relaxation of the external layers has significant effects on the surface energies and, the 3L_1 case excepted, the surface energies are linear with nL for the 3 constrained systems. Since the surface energy of the 5L_2 system (0.71 J/m²) is close the converged one (0.69 J/m²), this constrained system is a good compromise to model the (100) surface.

The internal constraints were then added to the (101) face slabs (see Fig. 4-c). The effects were the same as for the two previous faces. The energy of the unrelaxed surface is larger by 0.6 J/m² relative to the converged one. For a 5L_2 system, the calculated surface energy is 1.05 J/m² which is very close to the converged one for the fully relaxed system with 1.03 J/m². As observed for the (100) face, only the surface energies of the fully relaxed systems oscillate while those of the constrained ones are linear with the number of layers.

Finally, these internal constraints were added to the (001) face (see Fig. 4-d). The effects on the surface energies were exactly the same as for the (110) face: i) high surface energies for the unrelaxed systems (about 2.3 J/m²); ii) accurate stability with nL but larger surface energy for the second most constrained systems (nL_2); iii) net decrease of the surface energy for the nL_4 and nL_6 systems; iv) the surface energy is then very close to that of the fully relaxed systems. By adding the constraints, the convergence was reached with 8-10 layer systems.

Four surface models have been optimised: five layer slabs with their most internal layer frozen to bulk positions for the (110), (100) and (101) face and a nine layer slab with the same constraint for the (001) face. These constrained models have been used in the following calculations.

Surface Unsaturation and Relative Energies

The relative stabilities of these four surfaces can be explained by their respective unsaturation.

In Table 6 are reported the unsaturated valence densities, from GGA(PAW4) calculations, for titanium and oxygen atoms per nm², defined as the number of "broken bonds", relative to their bulk coordination (6 for titanium atoms and 3 for oxygen ones). The (110) face has the lowest number of unsaturated valence per nm² and is also the most stable face, followed by the (100) and then the (101) ones. The (001) face is the less stable one with a higher surface energy relative to the (110) one. The relaxed surface energies were found to increase with the surface densities of unsaturated sites (considering Ti and O simultaneously). The relaxed surface energy is thus directly linked to its unsaturation. The more the surface is unsaturated, the less it is stable. This correlation was also observed for the converged surface energies whatever the functional or the pseudopotential used (see Table 5).

In Table 6, the relaxed surface bond lengths are also reported and compared to the calculated bulk ones. These deviations were correlated to the unsaturation of the surface because the more it is unsaturated, the more the bond lengths should be modified. The smallest deviations were observed for the (110) and the (100) faces. For the threefold oxygen atoms, these deviations were more important for the (101) and then for the (001) faces. Regarding the twofold oxygen atoms, only bond shortening were observed with similar amplitude for the four faces. Considering the threefold oxygen atoms, a bond lengthening was observed in the case of the (110), (100) and (101) faces while a shortening occurred for the (001) face. This behaviour can be explained by the fact that the titanium atoms are fivefold coordinated for the (110), (100) and (101) faces while they are only fourfold coordinated on the (001) face.

Conclusion

First-principles calculations on the titanium dioxide crystals and on the TiO 2 rutile surfaces were performed in order to build accurate surface models which could be used to investigate surface defects, reconstruction or adsorption processes.

First, the lattice bulk parameters, the bulk moduli, the gap energies and the relative energies of rutile and anatase phases were calculated. The computed lattice bulk parameters and bulk moduli were found to be in good agreement with experimental data and previous theoretical works. The rutile phase was always found as the most stable one. It was found that the calculated energy gaps underestimate the experimental values as already mentioned for other systems from DFT calculations. Nevertheless, it appears that four valence electrons for titanium atoms and six for oxygen ones with the PAW pseudopotential formalism are well adapted to correctly describe bulk properties.

Secondly, rutile low index faces were investigated. In agreement with the natural rutile powder composition, the (110) face was found as the most stable one, followed by the (100)

and the (101) ones. The (001) face was also studied and its calculated surface energy is much higher than the three previous faces. A nearly linear correlation between the four surface energies and their respective unsaturation and relaxation amplitude was highlighted. In addition, it was observed that the surface energies of these four faces depend on the number of layers used in the slab to model the system. Internal constraints were thus added in the slabs for the four faces to stabilise their surface energies. These constrained systems allowed to build stable surface models for smaller slab thicknesses than the unconstrained ones. For the (110), (100) and the (101) faces, five layer models with their most internal layer frozen to bulk positions were found to be a good compromise to model these faces. Finally, for the (001) face, a nine layer slab with its most internal layer frozen to bulk positions was determined as an accurate model. Comparing the effect of the pseudopotentials (US10, PAW10 and PAW4), it was found that the one noted PAW4 in this study can be used to properly describe these systems.

Consequently, these four constrained surface models with PAW pseudopotentials (Ti: 3s a The surface normal is along the axe with one k-point. a There is only one unsaturated Ti per (001) unit area but it has 2 unsaturated valences instead of only one for the three other faces.

Figure 1 .

 1 Figure 1. The dry rutile (a) (110), (b) (100), (c) (101) and (d) (001) faces.

Figure 2 .

 2 Figure 2. Rutile (110) surface energy as a function of the slab thickness for 1 to 13 layer

Figure 3 .

 3 Figure 3. Rutile (100), (101) and (001) surface energy as a function of the slab thickness with

Figure 4 .

 4 Figure 4. Effect of internal constraints on surface energies (in J/m²) from GGA(PAW4)

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Figure 4 .

 4 Figure 4.

Table 1 .

 1 2 4p 2 and O:2s 2 2p 4 ) could be used to investigate, with a good accuracy, adsorbates-(TiO 2 surface) interactions or surface reconstruction phenomenon. The (110) face constrained model have been used to study the co-adsorption of molecular and dissociated water molecules[START_REF] Perron | Surface Science[END_REF] and to investigate the uranyl ion sorption [78] that are compared to experimental data[16,79]. Optimised rutile bulk parameters (in Å), bulk moduli (in GPa) and energy gaps (in eV). Data from previous studies: GGA and LDA using plane waves, and HF using linear combination of atomic orbitals, are also reported.

		a	c	x	c/a	B	Gap
	Exp.	4.587	2.954	0.305	0.644	210	3.0 b
	GGA(US10)	4.639	2.983	0.305	0.643	218	1.90
	GGA(PAW10) 4.642	2.973	0.305	0.640	217	1.80
	GGA(PAW4) 4.649	2.972	0.304	0.640	225	1.65
	LDA(US10)	4.554	2.937	0.304	0.645	254	1.94
	LDA(PAW10) 4.554	2.927	0.304	0.643	256	1.81
	LDA(PAW4) 4.561	2.929	0.303	0.642	263	1.66
	GGA a	4.651	2.964	0.307	0.637		
	LDA a	4.574	2.927	0.304	0.640		
	HF a	4.575	2.999	0.306	0.656	239	

a

[23]

. b

[59]

.

Table 2 .

 2 Optimised anatase bulk parameters (in Å), bulk moduli (in GPa) and energy gaps (in eV). Data from previous studies: GGA and LDA using plane waves, and HF using linear combination of atomic orbitals, are also reported.

		a	c	z	c/a	B	Gap
	Exp.	3.782	9.502	0.208	2.512	180	3.2 b
	GGA(US10)	3.794	9.822	0.205	2.589	191	2.30
	GGA(PAW10) 3.804	9.724	0.206	2.557	189	2.08
	GGA(PAW4)	3.810	9.726	0.206	2.552	199	1.93
	LDA(US10)	3.751	9.498	0.208	2.532	219	2.22
	LDA(PAW10) 3.747	9.488	0.208	2.532	219	2.04
	LDA(PAW4)	3.756	9.497	0.208	2.528	226	1.94
	GGA a	3.792	9.714	0.206	2.562		
	LDA a	3.758	9.495	0.208	2.527		
	HF a	3.781	9.735	0.203	2.575	202	

a

[23]

. b

[59]

.

Table 3 .

 3 Relative stabilities between rutile and anatase phases (in meV).

	Refs

a

[23]

. b

[66]

.

Table 4 .

 4 Rutile face's unit surface area and the corresponding k-point grid.

		Unit surface area (Ų) k-point grid a	Space between layers
	(110)	a × 2	c		3×1×5	(a	2	)	/	2
	(100)	a ×	c			1×5×5	a	/	2
	(101)	a	×	a	² c +	²	4×5×1	c
	(001)	a ×	a			5×5×1	c	/	2

Table 5 .

 5 Converged surface energies (in J/m²) calculated from the biggest systems in GGA for the three pseudopotentials.

	Face	E	10 surf US	E	PAW surf	4	E	PAW surf	10
	(110)	0.54	0.50		0.48
	(100)	0.76	0.69		0.67
	(101)	1.08	1.03		1.01
	(001)	1.32	1.25		1.21

Table 6 .

 6 Surface energies for the four constrained surface models (in J/m²) and their unsaturation (per nm²) from GGA(PAW4) calculations. Effects on the relaxed Ti-O surface bond lengths (in Å) relative to the bulk ones (1.97 and 2.00 Å, deviation in parenthesis in Å).

	E	4 constraine PAW	d	Unsaturation Twofold oxygen atoms	Threefold oxygen atoms
	(110) 0.60		10.23	1.85 (-0.12)	2.06 (+0.06)
	(100) 0.71		14.47	1.86 (-0.11)	2.06 (+0.06)
	(101) 1.05		15.59	1.83 (-0.14)	2.05 (+0.08)
					1.90 (-0.10)	2.12 (+0.12)
	(001) 1.38		18.50 a	1.88 (-0.09)	1.82 (-0.15)
					1.90 (-0.10)

Acknowledgments

This research has been done using the CEA CCRT supercomputers, in the framework of an EDF-CEA contract, as well as the cluster of CNRS-IN2P3 in Lyon.


Rémi Delassus

Romain Giot

Raphael Cherrier

Gabriele Barbieri

Guy Mélançon

Broken Bikes Detection Using CitiBike Bikeshare System Open Data

It seems necessary to detect a broken bike rooted at a station in near realtime as the number of bikes within bikeshare systems has reached more than a million in 2015. Indeed, a bike that cannot be moved is not cost effective in terms of number of trips. This brings frustration to users who were expecting to find a bike at that station without knowing that it is actually defective. We thus propose a methodology from feature extraction to anomaly detection on a distributed cloud infrastructure in order to detect bicycles requiring a repair. Through a first step of K-means clustering, and a second step consisting of spotting samples that do not clearly belong to any cluster, we separate anomalies from normal behaviors. The proposal is validated on a publicly available dataset provided by Motivate, the operator of the New-York bikeshare system. The number of distinct bikes that have been classified by this algorithm as broken at least once during a month is close to the number of repairs given in monthly reports of Motivate.

I. INTRODUCTION

Bikeshare systems generally refer to a set of stations, maintained by an operator, within a city in which bikes can be picked up and dropped off. As of 2015, more than one million bikes have been shared within more than 700 bikeshare systems worldwide [START_REF] Fishman | Bikeshare: A review of recent literature[END_REF]. Each one of those bikes has been heavily used, traveling up to an average of ten times a day in Barcelona as reported by Bikesharing Napoli [START_REF] Napoli | Bike share boom: 7 cities doing it right[END_REF]. The high number of bikes as well as the high usage level make breakdowns very likely to happen. Hence, many bikes are broken, waiting to be repaired, immobilized in a station. This is so common that bikeshare users have picked up the habit of turning the saddle backwards when encountering a broken bike, thus signaling to others that it is not worth renting [START_REF] Susi | How to be a good citi bike citizen[END_REF]. Thanks to this personal investment, other users avoid waisting time on renting and returning unusable bikes.

In the city of New York, although CitiBike's maintenance staff regularly performs roundtrips to inspect each bike's state, an alert system has been setup. Each dock has a button that can be pressed in order to alert that the docked bike is broken [START_REF] Susi | How to be a good citi bike citizen[END_REF]. The need for such a system proves that rounds are not efficient enough and operators are aware of it.

Bikes are repaired on site whenever possible (for example when the tire is deflated). Otherwise, the maintenance staff picks up the faulty bike and brings it to a specialized workshop that will fix it.

In a previous work, we have shown that as a bikeshare system sees its number of stations doubled, the number of bike trips is multiplied by three [START_REF] Bonnotte | Real-time data analytics and optimization of shared-vehicles networks[END_REF]. Therefore the number of trips in an expanding bikeshare system grows fast, implying that the number of failures and the number of users wanting to rent a bike will also quickly grow. This leads to a greater number of frustrated users that are in a situation where the only available bikes are broken.

To overcome this problem, we propose a model to classify, for each day, the status of a bike, deciding if it needs repairing or not. Our classification would let maintenance staff know whether a bike needs to be repaired so they could optimize their rounds. Moreover, it could be used within mobile applications where users can plan their trips, displaying the number of bikes not only available, but also usable.

The originality of our study relies on the fact that it is bike centered. This means for example that if a broken bike stays immobile a long time in a station, it is more relevant to study the fact that the bike is not moving rather than the fact that the station is not emptying. In order to follow each bike we work with the New York trips data [START_REF] Motivate | System data[END_REF] shared by Motivate, the company operating the CitiBike bikeshare system. Most of bikeshare systems operators have a website displaying the number of bikes and the number of free docks available at every station, which is used by the users to plan their trips. Few operators share their trips data (station and time of departure and arrival of identified bikes) and none of them share their maintenance data. Motivate exposes nearly three years of trips [START_REF] Motivate | System data[END_REF], along with monthly reports stating, among others informations, the number of repaired bikes in their workshops [6]. We propose a machine learning process to detect broken bikes. First we extract features from the trips dataset, and then perform an anomaly detection on the extracted features. At that point we use the data from the monthly reports as labels. Those labels are an aggregation of the repairs over a month. Since there is not a one-to-one correspondence between the features and the labels, but a many-to-one correspondence between the features and aggregated labels, we face a problem of Aggregated Output Learning [START_REF] Musicant | Supervised learning by training on aggregate outputs[END_REF]. This paper is organized as follows. Section II presents the existing work in anomaly detection applied to bikeshare systems. Section III describes the proposed method, including the feature extraction and the anomaly detection modules. Section IV depicts the protocol, such as the parameters selection and the infrastructure used. Finally sections V and VI respectively details the results of our method and discusses its limits.

II. EXISTING WORK

Although Bikeshare systems have been heavily studied [START_REF] Fishman | Bikeshare: A review of recent literature[END_REF], there is a lack of work related to anomaly detection. In an example of such a proposition [START_REF] Kaspi | Detection of unusable bicycles in bike-sharing systems[END_REF] Kaspi et al. use the same trips we are using to simulate the behavior of a station. This study is thus station centered. A stochastic model is then described. This model is parametrized by the probability of failure of bikes which must be known a priori. It is supposed that users do not rent bikes that are broken, which is only true if there are other possibilities. This model has not been confronted to the reality so we cannot know if the assumptions made to make a model were detrimental.

Bertens et al. defined a co-occurrence anomaly [START_REF] Bertens | Beauty and brains: Detecting anomalous pattern co-occurrences[END_REF] as two normal events occurring together where they are usually exclusive. They tried to detect such anomalies in a bikeshare system. They used rental data as well as weather data. Unfortunately, the only co-occurrence anomaly they found was a very rare co-occurrence of low real temperature and high perceived temperature.

As far as we know, no machine learning algorithm has been applied to bikeshare systems with the intent of detecting broken bikes, which is the focus of our work. One explanation for this lack of studies is that no label is available to train machine learning models. There is no operator that shares detailed maintenance data, such as the dates and nature of reparations of each bike.

However some operators, such as Motivate, share the number of repairs that occur each month. Musicant et al. [START_REF] Musicant | Supervised learning by training on aggregate outputs[END_REF] described machine learning models that can use such an aggregated output to be trained instead of usual labels.

Also, Fanaee et al. [START_REF] Fanaee-T | Event labeling combining ensemble detectors and background knowledge[END_REF] explain how without any prior label, they were able to label events with bikeshare data. But the events labeled did not occur within the bikesare system. The bikeshare system is rather used as an event detector and background knowledge retrieved from Google is used to provide labels.

III. PROPOSED METHOD

The method can be applied in two independent steps. The first step consists to extract features from the raw data, while the next one uses this enhanced and reshaped data to detect anomalies. A diagram recapitulating the process is presented in Figure 1.

A. Features extraction

The features extraction step is strongly linked to the available source of information. So far of our knowledge, only Motivate proposes data that could be used for fault detection. This data is an aggregation of the desired labels, which leads us an Aggregate Output Classification problem.

1) Raw data: We have at our disposal two kinds of data.

(i) The first one is massive and accurate; it describes the trips done with the CitiBike bikeshare system which can be downloaded on CitiBike's website [START_REF] Motivate | System data[END_REF]. Each trip is characterized by heterogeneous data described in Table I. We refer to that table as the trips dataset.

(ii) The second one is rather small and imprecise. It consists of a set of monthly reports [6] describing the activities of both the operator and the users during the month. It provides various informations such as the average number of daily trips per bike, the number of users, the respect of the service level agreement and, the one we use to optimize our model, the number of bicycle repairs. We refer to that set as the ground truth dataset. Note that we do not know what happens when a bike has been repaired twice in a month (does it counts for one or two bicycle repairs ?). We do not know how bad the breakage has to be to lead to the bike being picked up by the maintenance team.

We use the trips dataset to build a set of features and the ground truth dataset to train and evaluate our model.

2) Extracted Features: This study is not trip centered so it does not characterize a trip to know if it is a broken bike's trip or not. Indeed the lack of trip can be a sign of abnormality. Instead, it builds a sample per day and per bike, then determine if at a given day, a given bike was broken or not. For each of those samples, the tailored features are statistics on the life of the bike during a given interval of T days. They characterize the last of those days. By sliding a window of T days for each bike, several millions of samples are computed. If during the duration of a window (of period T ) a given bike has not been used at all, then this sample is not created.

The features extracted are as follow: number of trips made, number of distinct departure stations, number of distinct arrival stations, number of loop trips (same arrival and departure station), average distance traveled over the trips, maximal distance traveled over the trips, minimal distance traveled over the trips, average duration over the trips, maximal duration over the trips, minimal duration over the trips.

That means that each sample exists in a 10 dimensions space, each of those features being one of those dimensions.

Those features are chosen based on three criteria. Firstly, they are an easy way to aggregate several days of trip for a bike. Secondly, once the classification has been made, it is possible to define heuristics to check the coherence of the classification. For example, a broken bike would not make a lot of long trips across the city compared to other bikes. Lastly, the shape of the distribution of each feature is either Gaussian, or at least with a long tail of outliers events, see Figure 2. This shape is desirable for anomaly detection since it allows us to easily separate rarely occurring events from normal ones.

B. Anomaly detection

After the extraction of the features for each dichotomy of day and bike, it is necessary to classify them into one of the two categories: normal if the bike is running smoothly, abnormal if the bike has to be picked up by the maintenance.

Musicant et al. [START_REF] Musicant | Supervised learning by training on aggregate outputs[END_REF] proposed algorithms to perform machine learning on aggregated outputs such as these. However there is a major difference with their data that makes their solution not directly applicable. In their case the aggregated output is the sum of the output they would normally use. For example, given a set of five balls, classify each ball as a blue ball or

Parameters set (K, delta) seletion

Average Error For each Parameters

Training sets

Fig. 1. Process overview. We can see here the two main distinct processes; first the features extraction creates a features set from raw data. Then the cross validation process that trains and evaluates the model. During the grid search an error is computed for each parameters set and each fold. An averaged error over the folds is then computed for each set of parameters and the one with the smallest error is selected. It is then used to perform a classification on the validation set which are features that have never been used to train the model. The error obtained on the validation set is an indactor of how well our model can generalize to new data. a red ball, given that there are 3 blue balls and 2 red balls. The aggregated output is the sum of balls of a given color. However, in this work, the aggregated ouptut is the sum of the repairs that occurred during each month. Thus their method would be directly relevant if we were trying to classify each day of a bike's life as being one where a repair occurred. In this case, each day classified as a repair day would be counted and the resulting sum should approach the aggregated output.

But this is not the classification that is performed because the aim of this work is to classify each day of a bike's life as a day where it was broken or not. Since there is no information if a bike can be repaired several times in a month, or how that would impact the given number of repairs, it is assumed here that a bike is only repaired once a month. The aggregation function used is thus a count of the distinct bikes that have been classified as broken each month. In consequence, the algorithms proposed by Musicant bike, (ii) then the second steps consists in identifying samples that do not clearly belong to a cluster. They are the one that are too far from the closest cluster center.

1) Clustering: First, we suppose that a bike has a different behavior depending of its situation. For example, a bike docked near a green area such as a park or a water plan is not used in the same way than a bike parked near a residential area. The former is more likely to be used for a stroll, finishing in the same station it started. The latter might be used early in the morning to go to a working area and in the evening, to go back to a living area. This leads us to perform a clustering on our samples to separate different normal behaviors. To make this clustering, we use the K-Means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. K, the number of clusters and thus the number of relevant normal behaviors, being a parameter of our model. Each sample belongs to a cluster k, with clusters indexed by k ∈

[1; K] of size M k . Each sample of a cluster k is designated by S k i , i ∈ [1; M k ].
Each cluster possesses a centroid, C k , which is the average position of all the points of the cluster k. We can place a centroid in our multidimensional space by defining each of its features {C k j }, j ∈ [1; n] as the average of the values of this feature for the samples belonging to this cluster;

C k j = M k i=1 S k ij M k
. 2) Distance from cluster center: Once the samples have been assigned to a cluster, we still have to determine which ones are not normal. Since our clusters are supposed to represent different normal behaviors, we assume that the abnormal ones are those that are not clearly in a cluster. We consider that a sample is not clearly in a cluster when its distance to the centroid is significantly larger than the one of the other members of the cluster to the same centroid.

This distance is defined as follow:

dist(S k i , C k ) = n j=1 S k ij -C j k 2 (1)
Then we can define the function that is used to label sample as normal or abnormal as

anom(S i ) = 1 if dist(S k i , C k ) > δ * σ k 0 else (2)
where anom is a function that returns 1 if the sample S i is abnormal (the bike is supposed to be broken) and 0 if it is normal (the bike is not broken), σ k is the standard deviation of the distances to the centroid, δ is another parameter of our model, it is used to modulate the distance from the centroid at which a sample is considered abnormal.

IV. EXPERIMENTAL PROTOCOL

We describe here the protocol followed to train and evaluate our model.

1) Parameters optimization: We are left with two parameters to optimize; K the number of clusters (i.e., the number of relevant behaviors) and δ the factor of the standard deviation of the distance from the centroid at which a sample is classified as an anomaly. T , the number of days used to build the features is not optimized here but fixed to seven.

The parameter selection is controlled by a process of cross validation. The selection has been made over a year and a half, start from June 2013 down to June 2014. Then once the parameters have been fixed, we performed a classification over a period of 15 months, start from January 2015 down to March 2016, which is our validation set. During the selection phase, the features dataset has been separated five times in a training set and a testing set (5-folds cross validation). Each time all months from June 2013 to December 2014 are used, 80% of them for the training and 20% of them to compute an error (the result of the cost function described in Equation 3) on the testing set. As for now we are using the Root Mean Square Error (RMSE). The error of a set of parameters being the average testing error of the five different folds. Then the parameters set with smallest average error on the testing set is selected and the classification can be run on the validation set.

We find the best set of parameters thanks to a grid search since the space search is little enough to perform an exhaustive search. K is searched between 1 and 5; 1 means that there is no need for a clustering step, 5 is high enough so several behaviors may appear and low enough so the computation time is reasonable. If the error would have presented a tendency to decrease as K increases we would have had to push this limit farther but it has not been the case. δ is searched between 1 and 3, with steps of 0.25. That's because many of our features have a Gaussian-like shape, so we know that with σ < 1 we have too many abnormal samples and with σ > 3 we have too few of them. In order to judge if a set of parameters is better than an other one, we define a cost function (equation 3) that the models have to minimize. It measures the distance between the number of distinct bikes that have a sample detected as abnormal during a month and the number of repairs reported for this month.

1 N N i=1 (Anomalies K,δ (i) -Repairs(i)) 2 (3)
with Anomalies K,δ (m) the number of distinct bikes that have been classified as broken at least once during the month i and N the number of months used. This particular cost function is called Root Mean Square Error (RMSE)

The domains of the extracted features are different which could lead the computed distance between two vectors of this multi-dimensional space to depend almost exclusively of the difference between one of their dimensions. It is necessary to reduce and center them, so every feature has the same average and a similar value range.

2) Execution platform: This procedure has been implemented with Spark [12] 2.6.0. We ran the feature extraction and the parameters selection on both a desktop machine with an 8 cores i7 CPU, 32GB of RAM, and 100GB SSD hard drive and a cluster of 15 machines, each one with 50 GB of RAM and 48 cores, using the data stored in Hadoop Distributed File System. The configuration allocates 50 workers with 5 cores per worker (there are more than one worker per machine) and the possibility to use 5GB of RAM per worker.

V. RESULTS

This section presents for the results obtained with the optimal set of parameters.

This optimal set of parameters is K = 1 and δ = 2.25 as illustrated in Figure 4. These parameters provide us a model that performs better with the chosen features and cost function on unseen data.

Equation 3 presents to minimize. It implies that we are trying to minimize the distance between two values; the given number of repairs for each month and the number of distinct bikes that have been detected as broken at least once during said month. Figure 3 presents the evolution of those values. The training period where we performed cross validation ends in December 2014. The validation period starts in January 2015. During the latter we perform a classification on data that we have never seen before. Using the root mean square error (RMSE), we compute an error depending of the distance between the two curves.

On the training set the RMSE is 852 and on the validation set 1304. Since these two RMSE values are of the same order, there is no overfitting.

We computed a normalized cross-correlation (Equation 4) between the two curves during the training period and the validation period. The closer it is to 1, the more our classification and the ground truth follow the same variations The normalized cross-correlation is defined as follows:

1 n x (f (x) -f )(t(x) -t) σ f σ t (4)
where t(x) is the value of index x in the vector containing the values of the green curve, f (x) is the value of index x in the vector containing the values of the red curve, n is the number of months in t(x) and f (x), f is the average of f and σ f is the standard deviation of f . The normalized cross-correlation during the training period is 0.47 and during the validation period 0.35. Having low values for both the RMSE and the correlation metrics indicates that the selected δ is helpful to fit the average number of repairs, but does not enable us to predict the variability.

The implementation on a cloud architecture is justified by the size of the data. On a regular machine the features extraction requires 25 minutes of processing while the optimisation of the parameters is done in 2 hours without cross validation, and can not be performed with cross validation executed with a parallel Spark implementation: after several hours of computation, the disk is full and the process stops. On the cluster the features extraction is performed in 2 minutes and 30 seconds, the optimization without cross validation in 17 minutes and with the cross validation 2 hours.

VI. DISCUSSION

Since K = 1 is the best parameter, we failed to separate multiple normal behaviors. Three explanations are possible. Firstly, it is maybe not optimal to try to separate these behaviors and no solution with k > 1 can be better. Secondly, the most probable explanation is that the chosen features do not allow the behaviors to be separated. Indeed, we have not been able to visualize a separation between the clusters obtained with the Kmeans algorithm: scatter plots on each pair of dimensions do not separate them. We have to find other features that will lead to a meaningful clustering.

We could use t-sne [START_REF] Maaten | Visualizing data using t-sne[END_REF] or other advanced techniques for high dimensional data visualization in order to distinguish them. If it is not possible we will train a linear model on our features with our classification as target in order to analyze the parameters of such a model. We will then gain some insight into what makes a feature more susceptible to separate clusters. The third possibility is that the cost function does not favor such a separation. Several clues indicate that this function is currently suboptimal. For example, if it is not unusual to see a bike being repaired several times a month, and if we are sure that every broken bike is repaired during the month, then the distinct number of broken bikes our model outputs should always be lower than the registered number of repairs. We will thus heavily penalize the opposite case in our cost function.

On the other hand, the correlation is expected to be high and will be rewarded. Finally, a bike being detected as broken only one day or all days of the month makes no difference. A regularization term will help maintaining the overall number of anomalies lower. Another limitation of our approach is that we cannot yet figure out if the detected anomalies match broken bikes. The curve which represents the performance of our classification (green curve) can fit very closely to the curve representing the real number of repairs (red curve), without representing broken bikes.

One impactful thing to do can be using a classic Aggregate Output Classification algorithm to try to detect for each bike / day pair if this is a day at which the bike has been taken by the maintenance team for a repair. The count of these occurrences is supposed to be exactly the given aggregated labels. Thus we know that those Aggregated Output Classification algorithms, designed for this very case, will have better performances. This also unambiguously handles the case when a bike is picked up more than once. So we will be trying to fit a number of repairs to the given number of repairs, which will be less prone to error of modeling or error of interpretation. Then if that new model is good enough, we can use its output as unaggregated new labels. It will be necessary to have a new way of measuring the correctness of this model so we will be able to use it as a new "ground truth" giving for each bike the pick-up time. Using this new information, we will have a better way to determine how our current model is performing. All the anomalies we detect should occur just before the bike has to be picked up.

Finally, the factor of distance determining if a sample is too far to its centroid to be classified as normal, δ, can be unique for each cluster. The advantage would be that we might have whole clusters of anomalies. Also, the parameter K will be used to control the complexity of the model. As for now we cannot allow this because it would make our model very complex. With so many parameters to optimize both the computation time and the probability of overfitting are very high.

VII. CONCLUSION

Bike sharing systems suffer of broken bike and this paper present a methodology which help to find these broken bikes.

Using the CitiBike trips data as well as the given number of monthly repairs, we propose a method to classify each day of each bike's life as a day where this bike is not used, broken or normal. We do so by building features through the computation of statistics on the bike usage during the previous days, applying a clustering algorithm on the obtained samples and then detecting outliers.

We find that with the features designed and the cost function used, the clustering is not needed. By classifying samples too far from the average of all samples as broken bikes, we obtain a good fit to the curve representing the number of repairs each month. However, a small normalized cross correlation indicates that the number of distinct bikes classified at least once as broken during the month does not follow the same variations than the number of repairs. That is a clue to the fact that anomalies detected do not correspond to broken bikes.

Future work will aim at improving the model by using a better cost function, tailoring a new set of "raw data" via the use of standard Aggregate Output Classification methods, and finally factoring features that will let us take full advantage of the clustering algorithm.

Fig. 2 .

 2 Fig. 2. Probability distribution of all tailored features

Fig. 4 .

 4 Fig. 4. Best parameters. This represents the average error on the 5 testings sets with each set of parameters. The lighter the color, the smaller the error. The minimum error is reached at K = 1 and δ = 2.25.

TABLE I

 I

		RAW DATA DESCRIPTION
	Feature Name	Description	Frequency
	Departure location Departure station coordinates Always
	Arrival location	Arrival station coordinates	Always
	Departure station	Departure station name	Always
	Arrival station	Arrival station name	Always
	Departure date	Departure timestamp	Always
	Arrival date	Arrival timestamp	Always
	Bike	Bike unique identifier	always
	Birth	User birth year	When the user is a subscriber
	Gender	User gender	When the user is a subscriber
	Repairs	Reported repairs number	Once a month

  Classification versus ground truth. The red dotted curve shows the number of repairs that have been reported each month. The green curve represents the number of distinct bikes that have been classified as abnormal at least once this month. The training period goes from the beginning to December 2014. The validation set starts on January 2015. They are separated by a black vertical line.

					6000 7000							number of distinct bikes classifies as anormal reported number of repairs limit between training and validation sets
					5000
				count	3000 4000
					2000
					1000
					0	Sep 2013		Jan 2014		May 2014	Sep 2014	Jan 2015 months	May 2015	Sep 2015	Jan 2016	May 2016
	Fig. 3. 1.0	1.25	1.5	1.75	2.0	2.25	2.5	2.75	3.0	3.25	3.5

ACKNOWLEDGMENT

The authors would like to thank Motivate for their open data, since they are one of the few bikeshare operators who are sharing their trips. This work was supported by the CIFRE n o 1303/2014 of the French Association Nationale de la Recherche et de la Technologie.