T Hmidi

A Mantile
email: andrea.mantile@univ-rennes1.fr

F Nier

Time-dependent delta-interactions for 1D Schrödinger Hamiltonians

Keywords: MSC (2000 37B55, 35B65, 35B30, 35Q45, keywords: Point interactions, solvable models in Quantum Mechanics, non autonomous Cauchy problems

The non autonomous Cauchy problem i∂tu = -∂ 2 xx u+α(t)δ0u with ut=0 = u0 is considered in L 2 (R) . The regularity assumptions for α are accurately analyzed and show that the general results for non autonomous linear evolution equations in Banach spaces are far from being optimal. In the mean time, this article shows an unexpected application of paraproduct techniques, initiated by J.M. Bony for nonlinear partial differential equations, to a classical linear problem.

Introduction

This work is concerned with the dynamics generated by the particular class of non-autonomous quantum Hamiltonians: H α(t) = -d 2 dx 2 + α(t)δ, defining the time dependent delta shaped perurbations of the 1D Laplacian. Quantum hamiltonians with point interactions were first introduced by physicists as a computational tool to study the scattering of quantum particles with small range forces. Since then, the subject has been widely developed both in the theoretical framework as well as in the applications (we refer to [START_REF] Albeverio | Solvable Models in Quantum Mechanics[END_REF] for an extensive presentation). For real values of the coupling parameter α, the rigorous definition of: H α = -d 2 dx 2 + αδ arises from the Krein's theory of selfadjoint extentions. In particular, H α identifies with the selfadjoint extension of the symmetric operator: H 0 = -d 2 dx 2 , D(H 0 ) = C ∞ 0 (R\ {0}) defined through the boundary conditions

ψ ′ (0 + ) -ψ ′ (0 -) = αψ(0) ψ(0 + ) -ψ(0 -) = 0 (1.1)
ψ(0 ± ) denoting the right and left limit values of ψ(x) as x → 0 [START_REF] Albeverio | Solvable Models in Quantum Mechanics[END_REF]. Explicitely, one has

D(H α ) = ψ ∈ H 2 (R\ {0}) ∩ H 1 (R) ψ ′ (0 + ) -ψ ′ (0 -) = αψ(0) (1.
2)

H α ψ = - d 2 dx 2 ψ in R\ {0} . (1.3)
When α(t) is assigned as a real valued function of time, the domain D(H α(t) ) changes in time with the boundary condition (1.1), while the form domain is given by H 1 (R). The quantum evolution associated to the family of operators H α(t) is defined by the solutions to the equation

i d dt u = H α(t) u u |t=0 = u 0 . (1.4)
The mild solutions are the solutions to the associated integral equation u(t) = e it∆ u 0 -i t 0 e i(t-s)∆ q(s)δ 0 ds (1.5) with q(s) = α(s)u(s, 0) . The questions are about:

• the regularity assumptions on t → α(t) for which (1.4) defines a unitary strongly continuous dynamical system on L 2 (R)

• the meaning of the differential equation (1.4), according to the regularity of t → α(t) .

General conditions for the solution of this class of problems have been long time investigated. In the framework of evolution equations in Banach spaces, Kato was the first who obtained a result for the Cauchy problem

d dt u = A(t)u u t=0 = u 0 , (1.6)
when t → A(t) is an unbounded operator valued function [START_REF] Kato | Integration of the equation of evolution in a Banach spac[END_REF]. This result, which applies to the quantum dynamical case for A(t) = -iH(t), requires the strong differentiability of the function t → A(t) and the time independence of the domain D(A(t)). Afterwards, a huge literature was devoted to this problem in the main attempt of relaxing the above conditions (e.g. in [START_REF] Simon | Quantum Machanics for Hamiltonians defined as quadratic forms[END_REF] and [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF]; a rather large bibliography and an extensive presentation of the subject are also given in [START_REF] Fattorini | The Cauchy problem[END_REF]).

In particular, the time dependent domain case was explicitely treated in [START_REF] Kisyński | Sur les opérateurs de Green des problèmes de Cauchy abstraits[END_REF] by Kisyński using coercivity assumptions and C 2 loc -regularity of t → A(t). The regularity conditions in time were substantially relaxed into a later work of Kato [START_REF] Kato | Linear evolution equations of 'hyperbolic' type[END_REF] who proves weak and strong existence results, for the solutions to (1.6), when {A(t)} forms a stable (a notion defined in [START_REF] Kato | Linear evolution equations of 'hyperbolic' type[END_REF]-Definition 2.1 expressing uniform bounds for the norms of resolvents) family of generators of contraction semigroups leaving invariant a dense set Y of the Banach space X, and the map t → A(t) is norm-continuous in L(Y, X) .

Due to the particular structure, point interaction Hamiltonians allows rather explicitely energy and resolvent estimates, so that most of the techniques employed in the analysis of non-autonomous Hamiltonians can be used to deal with the equation (1.4), provided that α(t) is regular enough. At this concern, the Yafaev's works [START_REF] Yafaev | On "eigenfunctions" of a time-dependent Schrödinger equation[END_REF], [START_REF] Yafaev | Scattering theory for time-dependent zero-range potentials[END_REF] and [START_REF] Sayapova | The evolution operator for time dependent potentials of zero radius[END_REF] (with M. Sayapova) on the scattering problems for time dependent delta interactions in the 3D setting are to be recalled: There the condition α ∈ C 2 loc (t 0 , +∞) is used to ensure the existence of a strongly differentiable time propagator for the quantum evolution. Such a condition, however, could be considerably relaxed. In our case, for instance, a first appraoch consist in adapting the strategy of [START_REF] Kisyński | Sur les opérateurs de Green des problèmes de Cauchy abstraits[END_REF] by constructing a family of unitary maps V t,t0 such that: V t,t0 H α(t) V * t,t0 has a constant domain; then, it is possible to solve the evolution problem for the deformed operator by using results from [START_REF] Kato | Linear evolution equations of 'hyperbolic' type[END_REF]. To fix the idea, let y t be the time dependent vector field defined by ẏt = g(y t , t)

y t0 = x (1.7) with g(•, t) ∈ C 0 (0, T ; C ∞ 0 (R)), supp g(•, t) ⊂ (-1, 1
), and g(0, t) = 0 for each t. Under these assumptions, (1.7) allows an unique solution depending continuously from time and Cauchy data {t 0 , x}. Using the notation: y t = F (t, t 0 , x), one has

∂ x F (t, t 0 , x) = e R t t 0 ∂1g(ys,s)ds > 0 ∀x ∈ R , (1.8)
∂ 1 g(•, s) denoting the derivative w.r.t. the first variable. This condition allows to consider the map of x → F (t, t 0 , x) as a time-dependent local dilation and one can construct the family of time dependent unitary transformation associated to it

(V t0,t u) (x) = (∂ x F (t, t 0 , x)) 1 2 u(F (t, t 0 , x)) V t,t0 = V -1 t0,t
Under the action of V t,t0 , the equation (1.4) reads as

i d dt v = V t,t0 H α(t) V t0,t v -i 1 2 [∂ y g](y, t) + g(y, t)∂ y v v t=0 = u 0 (1.9) with V t,t0 H α(t) V t0,t = -∂ y b 2 ∂ y + a 2 -(∂ y ab) + b 0 α(t)δ b(y, t) = e R t t 0 ∂1g(F (s,t,y),s) ds ; b 0 = b(0, t) a(y, t) = 1 2 t t0 ∂ 2 1 g (F (s, t, y), s) e R t t 0 ∂1g(F (s ′ ,t,y),s ′ ) ds ′ ds and V t,t0 u = v. Set: A(t) = iV t,t0 H α(t) V t0,t + 1 2 ∂ y g(y, t) + g(y, t)∂ y , the domain D(A(t)) is the subspace of H 2 (R\ {0}) ∩ H 1 (R) identified by the boundary condition b(0, t) u ′ (0 + ) -u ′ (0 -) = α(t)u(0)
For α ∈ W 1,1 (0, T ) and sign(α(t)) = const., one can determine (infinitely many) g(y, t) such that: 1 k b(0, t) = α(t) for a fixed constant k. With this choice, the operator's domain is constant, D(A(t)) = Y . Moreover, under the same conditions on α, one can shows that A(t) defines a stable family of skew-adjoint operators t-continuous in L(H 1 , H -1 )-operator norm. Thus, one can use the Theorem 5.2 and Remark 5.3 in [START_REF] Kato | Linear evolution equations of 'hyperbolic' type[END_REF] to get strongly solutions for the related evolution problem. This is summarized in the following Proposition.

Proposition 1.1 Let α ∈ W 1,1 (0, T ), sign(α(t)) = const. and u 0 ∈ D(H α(0)
). There exists an unique solution u t of the problem (1.4), with:

u t ∈ D(H α(t) ) for each t and u ′ t ∈ C 0 (0, T ; L 2 (R)).
In spite of those improved results with already known general tools, our aim is to prove that they are far from being optimal. An additional structure can lead to the same conclusions with weaker regularity assumptions. This is the question that we propose to explore with one dimensional δ-interactions which allow direct and explicit computations. The main result of this paper is the following.

Theorem 1.2 1) Assume α ∈ H 1 4 loc (R) , then for any u 0 ∈ H 1 (R) the integral equation (1.5) admits a unique solution u ∈ C(R; H 1 (R)) with i∂ t u -α(t)u(t, 0)δ 0 ∈ C(R; H -1 (R)) and (1.4) is weakly well-posed.
2) With the same assumption, (1.5) defines a unitary strongly continuous dynamical system U (t, s)

on L 2 (R) . 3) If additionally α ∈ H 3/4 loc (R), then for any u 0 ∈ D(H α(0) ) the solution u of (1.5) belongs to the space C 1 (R; L 2 (R)), with u(t) ∈ D(H α(t) ) for every t ∈ R.
Remark. The problem of defining the quantum evolution for 1D time dependent delta interactions has also been considered in a nonlinear setting [START_REF] Adami | A class of nonlinear Schrödinger Equations with concentrated nonlinearitie[END_REF] where α is assigned as a function of the particle's state, in our notation: α(t) = γ |u t (0)| 2σ , with γ ∈ R, σ ∈ R + . In this framework, the authors prove that solutions to the nonlinear evolution problem exist locally in time for u 0 ∈ H ρ with ρ >

i ∂ x = F -1 • (ξ×) • F
, where F is the Fourier transform in position (or in time) normalized according to

F ϕ(ξ) = R e -iξx ϕ(x) dx .
The Sobolev spaces are denoted by H s (R), s ∈ R, their local version by H s loc (R). The notation u ∈ H s+0 (R) means that there exists ε > 0 such that u ∈ H s+ε (R) (inductive limit) and its local version u ∈ H s+0 loc (R) allows ε R > 0 to depend on R > 0 while considering the interval [-R, R]. More generally the Besov spaces are defined through dyadic decomposition: For (p, r) ∈ [1, +∞] 2 and s ∈ R, the space B s p,r is the set of tempered distributions u such that

u B s p,r := 2 qs ∆ q u L p ℓ r < +∞ ,
where ∆ q = ϕ(2 -q D), q ∈ N, is a cut-off in the Fourier variable supported in C -1 2 q ≤ |ξ| ≤ C2 q . Details are given in Appendix A. Finally, the notation ' ', appearing in many of the following proofs, denotes the inequality: '≤ C', being C a suitable positive constant.

Proof of Theorem 1.2

This theorem is a consequence of simple remarks, explicit calculations and standard applications of paraproduct estimates in 1D Sobolev spaces. Let us start with some elementary rewriting of the Cauchy problem (1.4).

Preliminary remarks

• First of all equation (1.4) or its integral version (1.5) are local problems in time so that t 0 = 0, t ∈ [-T, T ] for some T > 0 and even supp α ⊂ [-T /2, T /2] can be assumed after replacing α with α T (s) = α(s)χ( s T ) for some fixed χ ∈ C ∞ 0 ((-1/2, 1/2)) and χ ≡ 1 near s = 0. The dependence of H s -norms of α T with respect to T will be discussed when necessary.

• The equation (1.4) or its integral version (1.5) makes sense in S ′ (R x ) as soon as u(t, 0) is well defined for almost all t ∈ [-T, T ] and q(t) = α T (t)u(t, 0) is locally integrable. Then it can be written after applying the Fourier transform as a local problem in ξ ∈ R

i∂ t u(t, ξ) = |ξ| 2 u(t, ξ) + q(t) u(0, ξ) = u 0 (ξ) (2.1) with q(t) = α T (t)u(t, 0) = α T (t) R u(t, ξ) dξ . (2.2)
This is equivalent to the integral form

u(t, ξ) = e -it|ξ| 2 u 0 (ξ) -i t 0 e -i(t-s)|ξ| 2 q(s) ds (2.3) with q(t) = q 0 (t) -iα T (t) t 0 R e -i(t-s)|ξ| 2 q(s) dsdξ (2.4)
by setting q 0 (t) = α T (t)[e it∆ u 0 ](0) . The assumption

u 0 ∈ H s (R), s > 1/2, (resp. u 0 ∈ L 1 (R)) ensures that [e it∆ u 0 ](0) ∈ C 0 ([-T, T ]) (resp. t 1/2 [e it∆ u 0 ](0) ∈ C 0 ([-T, T ]))
. Such an assumption as well as looking for u(t) ∈ H 1 (R) ensures that the quantities q 0 (t) and q(t) make sense for almost all t ∈ [-T, T ].

• With the support assumption supp α T ⊂ [-T /2, T /2], the convolution equation (2.3) can be written

q(t) = q 0 (t) -iα T (t) t 0 R 1 [-T,T ] (t -s)e -i(t-s)|ξ| 2 1 [-T,T ] (s)q(s) dsdξ in D ′ (R)
:= q 0 (t) -iα T (t)Lq(t) := q 0 (t) + L α q(t) .

(2.5)

• Once q is known after solving (2.5), equation (2.3) with t ∈ R reads simply

u(t, ξ) = e -it|ξ| 2 u 0 (ξ) -ie -it|ξ| 2 F q1 [0,t] (-|ξ| 2 ) . (2.6)
• When u 0 and q are regular enough the time-derivative of the quantity (2.3) gives

i∂ t (∂ t u)(t, ξ) = |ξ| 2 ∂ t u(t, ξ) + q ′ (t).
By Duhamel formula, this implies

∂ t u(t, ξ) = e -it|ξ| 2 ∂ t u(0, ξ) -i t 0 e -i(t-s)|ξ| 2 q ′ (s)ds , while (2.1) says for t = 0 ∂ t u(0, ξ) = -i|ξ| 2 u 0 (ξ) -iq(0) .
Therefore we obtain for t ∈ R

i∂ t u(t, ξ) = e -it|ξ| 2 |ξ| 2 u 0 (ξ) + q(0) + e -it|ξ| 2 F (∂ s q)1 [0,t] (-|ξ| 2 ) . (2.7)

Reduced scalar equation for q

Let us now study the equation (2.5) written:

q = q 0 + L α q with L α q := -iα T (t)Lq = -iα T (t) t 0 R 1 [-T,T ] (t -s)e -i(t-s)|ξ| 2 1 [-T,T ] (s)q(s) dsdξ .
Solving this fixed point equation relies on the next result.

Proposition 2.1 The estimate

Lq H s T 1 2 1 [-1,1] (D t )q L 2 + T 1 2 -θ q 1 [0,T ] H s-θ + q 1 [-T,0] H s-θ
holds for every s ∈ R and θ ∈ [0, 1 2 ] .

Proof: Owing to R e ±iλ|ξ| 2 dξ 2π = e ±i π 4 √ 4πλ for λ > 0, Lq writes as

1 √ π Lq(t) = e -i π 4 R 1 [0,T ] (s)q(s) 1 [0,T ] (t -s) (t -s) 1 2 ds + e i π 4 R 1 [-T,0] (s)q(s) 1 [-T,0] (t -s) (s -t) 1 2
ds Passing to the Fourier transform, we get

1 √ π Lq(τ ) = e -i π 4 T 0 e -itτ √ t dt F (1 [0,T ] q)(τ ) + e i π 4 T 0 e itτ √ t dt F (1 [-T,0] q)(τ )
One easily checks

T 0 e ±itτ √ t dt ≤ 2 √ T and T 0 e ±itτ √ t dt |τ | -1 2 .
This yields for every θ ∈ [0, 1 2

], τ ∈ R T 0 e ±itτ √ t dt T 1 2 1 [-1,1] (τ ) + T 1 2 -θ |τ | -θ 1 {R\[-1,1]} (τ ).
Thus we get for every s ∈ R, θ ∈ [0, 1 2 ],

Lq H s T 1 2 1 [-1,1] (D t )q L 2 + T 1 2 -θ 1 [0,T ] q H s-θ + 1 [-T,0] q H s-θ . Proposition 2.2 1. Let u 0 ∈ H s (R x ) with s > 1/2 and let α ∈ H 1 4 loc (R t ) . Then the equation (2.5) has a unique solution q ∈ H 1 4 loc (R t ). Moreover, for a fixed u 0 ∈ H s (R x ) with s > 1/2, the map α → q is locally Lipschitzian from H 1 4 (R t ) to H 1 4 (R t ) . 2. Let u 0 ∈ H 1 (R x ) and let α ∈ H 3 4 loc (R t ) . Then the equation (2.5) has a unique solution q ∈ H 3 4 loc (R t ). Moreover, for a fixed u 0 ∈ H 1 (R x ) the map α → q is locally Lipschitzian from H 3 4 (R t ) to H 3 4 (R t ).
3. Let u 0 ∈ H 1+2ε (R x ) and let α ∈ H Proof: 1) Let us first prove that q 0 ∈ H

1 4 when u 0 ∈ H 1 2 +ε . Write first (e it∆ u 0 )(0) = 1 -1 e -it|ξ| 2 u 0 (ξ) dξ 2π + +∞ 1 e -itτ u 0 ( √ τ ) + u 0 (- √ τ ) dτ 2π √ τ = I(t) + II(t) .
The first term I(t) defines a C ∞ function with

I B s ∞,∞ u 0 L 2 ,
for every s ∈ R . On the other hand, the Fourier transform of II equals

II(-τ ) = 1 [1,∞[ (τ ) u 0 ( √ τ ) + u 0 (- √ τ ) 1 √ τ .
The Sobolev regularity of the second term, II, is given by:

II 2 H ν +∞ 1 τ 2ν u 0 ( √ τ ) √ τ 2 dτ u 0 2 H 2ν-1/2 , (2.8)
for any ν ∈ R . Now, write q 0 (t) = α T (t)I(t) + α T (t)II(t).

Lemma A.2-b) applied to the first term, implies

α T I H 1 4 α T H 1 4 I B s+ε ∞,∞ α T H 1 4 u 0 L 2 .
For the second term, use Lemma A.2-a), Sobolev embeddings and (2.8)

α T II H 1 4 α T H 1 4 II B 1 2 2,∞ ∩L ∞ α T H 1 4 II H 1 2 + ε 2 α T H 1 4 u 0 H 1 2 +ε .
By combining these estimates, we get

q 0 H 1 4 α T H 1 4 u 0 H 1 2 +ε .
It remains to estimate α T . Let χ ∈ D(R) with χ ≡ 1 in [-1, 1] and set α(t) = χ(t)α(t) . By using again Lemma A.2-a), we get for 0 ≤ T ≤ 1

α T H 1 4 α H 1 4 χ(T -1 •) B 1 2 2,∞ ∩L ∞ . A change of variable in the Fourier transform F χ(T -1 .) (τ ) = T χ(T τ ) leads to χ(T -1 .) H µ ≤ T 1 2 -µ χ H µ and χ(T -1 .) B µ 2,∞ T 1 2 -µ χ B µ 2,∞ , (2.9)
for µ ≥ 0 and T ≤ 1. Hence we get

α T H 1 4 α H 1 4
.

and q 0 H 1 4 α H 1 4 u 0 H 1 2 +ε (2.10)
In order to estimate the the operator L, use Lemma A.2-a)

L α q H 1 4 α T H 1 4 Lq B 1 2 2,∞ ∩L ∞ α H 1 4
Lq H 1 2 +ε , while Proposition 2.1 says

Lq H 1 2 +ε T 1 2 q L 2 + T 1 4 -ε 1 [0,T ] q H 1 4 + 1 [-T,0] q H 1 4
.

Hence we get for 0 ≤ T ≤ 1 and by Lemma A. 2-a)

Lq H 1 2 +ε T 1 4 -ε q H 1 4
This yields

L α q H 1 4 α H 1 4 T 1 4 -ε q H 1 4 . (2.11)
This proves that L is a contracting map in H 1 4 for sufficiently small time T . The time T depends only on α H 1 4 and then we can construct globally a unique solution q ∈ H 1 4 loc (R) for the linear problem (2.5). This proves that L is a contracting map in H 3 4 for sufficiently small time T . The time T depends only on α H For the locally Lipschitz dependence with respect to α, the proof is left to the reader: it can be done easily done like for the case α ∈ H 1 4 .

3) Like in the proof of the second point 2) we get

q 0 H 3 4 +ε α T H 3 4 +ε u 0 L 2 + II H 3 4 +ε T -1 4 -ε u 0 H 1+2ε .
Reproducing the same computation as (2.12) leads to

L α q H 3 4 +ε α T L ∞ Lq H 3 4 +ε + α T H 3 4 +ε Lq L ∞ α L ∞ T 1 12 q1 [0,T ] H 1 3 +ε + q1 [-T,0] H 1 3 +ε + T -1 4 -ε α H 3 4 +ε Lq H 1 2 +ε α L ∞ T 1 12 q H 1 3 +ε + T -1 4 -ε α H 3 4 +ε T 1 3 +ε q H 1 3 +2ε T 1 12 α H 3 4 +ε q H 3 4
+ε .

With the fixed point argument we can conclude the proof.

Regularity of u

We start with the following result.

Lemma 2.3 For s ∈ R, let H s T be the closed subset of H s (R t )

H s T = {u ∈ H s (R t ), supp u ⊂ [-T, T ]} ,
endowed with the norm H s . For any T > 0 and any s ∈ R, there is a constant C T,s such that

∀f ∈ H 2s-1 4 T , F -1 F f (-|ξ| 2 ) H s ≤ C T,s f H 2s-1 4
.

Proof: It suffices to compute

R (1 + |ξ| 2 ) s f (-|ξ| 2 ) 2 dξ = +∞ 0 (1 + τ ) s 2τ 1/2 f (-τ ) 2 dτ ≤ max τ ∈[0,1] f (-τ ) 2 + ∞ 1 (1 + τ ) s-1/2 f (-τ ) 2 dτ ≤ max τ ∈[0,1] f (-τ ) 2 + f 2 H s 2 -1 4 ,
where f (τ ) = e iτ x χ(x) , f with χ ∈ D(R) with value 1 in [-T, T ]. By duality we have for ν ∈ R

sup 0≤τ ≤1 | f (τ )| ≤ f H ν sup 0≤τ ≤1 e iτ • χ(•) H -ν ≤ C 1 T,ν f H ν .
The main result of this section is the following.

Proposition 2.4 1. Let u 0 ∈ H 1 (R x ), α ∈ H 1 4 loc (R t ), then the equation (1.5) has a unique solution u ∈ C(R; H 1 (R)) . with T u v = q S q-1 u∆ q v and R(u, v) = |q ′ -q|≤1 ∆ q u∆ q ′ v .
Let us now recall the definition of Besov spaces through dyadic decomposition. For (p, r) ∈ [1, +∞] 2 and s ∈ R, the space B s p,r is the set of tempered distribution u such that u B s p,r := 2 qs ∆ q u L p ℓ r < +∞.

This definition does not depend on the choice of the dyadic decomposition. One can further remark that the Sobolev space H s coincides with B s 2,2 . Below is the Bernstein lemma that will be used for the proof of product laws and which is a straightforward application of convolution estimates and Fourier localization.

Lemma A.1 There exists a constant C such that for q, k ∈ N, 1 ≤ a ≤ b and for f ∈ L a (R), sup |α|=k ∂ α S q f L b ≤ C k 2 q(k+ 1 a -1 b ) S q f L a , C -k 2 qk ∆ q f L a ≤ sup |α|=k ∂ α ∆ q f L a ≤ C k 2 qk ∆ q f L a .
The following product laws have been used intensively in the proof of our main result.

Lemma A.2 In dimension d = 1 the product (u, v) → uv is bilinear continuous a) from H s × (B
To estimate the remainder term we use the fact

∆ q R(u, v) = j≥q-4 |j-j ′ |≤1 ∆ q (∆ j u∆ j ′ v).
According to Bernstein lemma one gets

2 qs ∆ q (R(u, v)) L 2 2 q(s+ 1 2 ) j≥q-4 |j-j ′ |≤1 ∆ j u L 2 ∆ j ′ v L 2 j≥q-4 |j-j ′ |≤1 2 (q-j)(s+ 1 2 ) 2 js ∆ j u L 2 2 j ′ 1 2 ∆ j ′ v L 2 v B 1 2 2,∞ j≥q-4 2 (q-j)(s+ 1 2 ) 2 js ∆ j u L 2 .
It suffices now to apply the convolution inequalities. b) First remark that the case s = 0 is obvious: L 2 × L ∞ → L 2 and B ε ∞,∞ ֒→ L ∞ . Hereafter we consider s > 0 . To estimate the first paraproduct, use the embedding B s+ε ∞,∞ ֒→ B s ∞,2 , for ε > 0.

T u v 2 H s q 2 2qs S q-1 u 2 L 2 ∆ q v 2 L ∞ u L 2 v B s ∞,2 u L 2 v B s+ǫ ∞,∞
For the second term we use the result obtained in the part a):

T v u 2 H s q 2 2qs S q-1 v 2 L ∞ ∆ q u 2 L 2 v 2 L ∞ u 2 H s v 2 B s+ǫ ∞,∞ u 2 H s
To estimate the remainder term we write

2 qs ∆ q (R(u, v)) L 2 2 qs j≥q-4 |j-j ′ |≤1 ∆ j u L 2 ∆ j ′ v L ∞ v L ∞ j≥q-4 |j-j ′ |≤1
2 (q-j)s 2 js ∆ j u L 2 .

Since s > 0 then we obtain by using the convolution inequalities

R(u, v) H s ≤ v L ∞ u H s .
c), d) These results are standard, see for example [START_REF] Chemin | Perfect incompressible fluids. Translated from the 1995 French original by Isabelle Gallagher and Dragos Iftimie[END_REF]. 2,∞ (R)).

B Sobolev and

Proof:

1) The Fourier transform of 1 [0,t] (s) equals F (1 [0,t] )(τ ) = e -iτ t -1 -iτ . One gets for ν ∈ [0, 1/2) Choosing judiciously λ then we obtain for |t 2 -

t 1 | ≤ 1 1 [0,t1] -1 [0,t2] H ν ≤ C ν |t 2 -t 1 | 1/2-ν .
2) We set f t (s) := 1 [0,t] (s), then have

f t 2 B 1 2 2,∞ ≤ f 2 L 2 + max q∈N 2 q 2 q ≤|τ |≤2 q+1 | f t (τ )| 2 dτ |t| + max q∈N 2 q 2 q ≤|τ |≤2 q+1 |τ | -2 dτ |t| + 1.

1

 1

2 + 1 4 loc

 24 and: U (t, s)u 0 (0) ∈ H ρ as a function of time. This corresponds to the condition: α ∈ H ν loc with ν >1 2 . Although a linear problem is considered here, the result of Theorem 1.2 improves significantly the regularity condition: α ∈ H 1/4 loc (R) is weaker than a continuity assumption. In what follows, D x denotes 1

4

 4 t ) , for some ε > 0. Then the equation (2.5) has a unique solution q ∈ H 3 +ε loc (R t ).

3 4

 3 and then we can construct globally a unique solution q ∈ H 3 4 loc (R) for the linear problem.

R ( 1 + 4 R( 1 + τ 2 )I |t 2 -t 1 | 2 λ 0 τ 2ν dτ + +∞ λ τ 2ν-2 dτ |t 2 -

 14122202 τ 2 ) ν F (1 [0,t1] ) -F(1 [0,t2] ) 2 (τ ) dτ = ν |sin(τ |t 2 -t 1 |/2t 1 | 2 λ 2ν+1 + λ 2ν-1 .

  Besov regularity of cut-offs . The map t → 1 [0,t] (s) belongs to L ∞ (R + ; B

	1
	2

Lemma B.1 1. For any ν < 1/2, t → 1 [0,t] (s) belongs to C(R + ; H ν (R)) . More precisely we have for |t -t ′ | ≤ 1 1 [0,t] -1 [0,t ′ ] H ν |t -t ′ | 1 2 -ν .

2

It remains to prove the continuity dependence of q with respect to α. Let α, ᾱ ∈ H 1 4 loc and q, q the corresponding solutions then we have q(t)q(t) = α T (t)Lq(t) -ᾱT (t)Lq(t), with ᾱT (t) = ᾱ(t)χ(t/T ).

Since L is linear on q then q(t)q(t) = α T (t)L(qq)(t) + (α T -ᾱT )(t)Lq(t) = L α (qq)(t) + L α-ᾱ q(t).

To estimate the terms of the r.h.s we use (2.11)

With the choice of T done above we get

This achieves the proof of the continuity.

2) Write again q 0 (t) = α T (t)

Since H 3 4 is an algebra the inequality

holds for T ∈ [0, 1], owing to (2.9). It follows

The second term is estimated with (2.8):

Finally we get for T ∈ [0, 1]

, Sobolev embeddings and Proposition 2.1 (with θ = 5 12 and θ = 1 6 ), gives

+ε .

(2.12)

Thus we get for 0 ≤ T ≤ 1,

loc (R t ), then the equation (1.5) has a unique solution u belonging to the space C 1 (R; L 2 (R)), with u(t) ∈ D(H α(t) ) for all t ∈ R .

Proof: 1) The solution of (1.5) is obtained via the equation (2.6)

Let us check that we have the required regularity for u. Applying Lemma 2.3 to (2.6) implies

(2.13) Lemma A.2-a) and Lemma B.1 yield

. This proves that u ∈ L ∞ loc (R; H 1 ). It remains to prove the continuity in time of u. First notice that we need for this purpose to prove only the continuity in time of v(t) := u(t) -e it∆ u 0 . This will be done in two steps. In the first one we deal with the case α ∈ H 3 4 loc . In the second one, we go back to the case α ∈ H

loc . Remark that according to Proposition 2.2-2) we can construct a unique solution q ∈ H 3 4 loc for the problem (2.5). An easy computation gives for t, t

.

Using the fact sin 2 x ≤ |x| ε , ∀ε ∈ [0, 1], and Lemma 2.3 gives

It suffices now to use Lemma A.2-a)

For the second term we use again Lemma A.2-c) combined with the proof of Lemma B.1

This concludes the proof of the time continuity of u when α ∈ H

loc . We smooth out the function α leading to a sequence of smooth functions α n that converges strongly to α in H 1 4 loc . To each α n we associate the unique solutions q n and u n . From the first step u n belongs to C(R; H 1 ). Similarly to (2.13) we get for n, m ∈ N

By Proposition 2.2-a), {q n } is a Cauchy sequence in H 1 4 and thus {u n } converges uniformly to u in L ∞ T H 1 . This gives that u ∈ C([-T, T ], H 1 ), for every T > 0. 2) Recall from (2.7) that

Since u 0 ∈ D(H α(0) ) then the first term of the r.h.s belongs to C(R; L 2 ). On the other hand we have D(H α(0 ) ⊂ H 3 2 -d for any d > 0. It follows from Proposition 2.4-1) that we can construct a unique solution q ∈ H 3 4 . Now, let w(t) := i∂ t u -e it∆ H α(0) u 0 . Then Lemma 2.3 yields

Thus we get for every t ∈ [-T, T ]

(2.14)

It follows that w ∈ L ∞ loc (R; L 2 ). To prove the continuity in time of w we use the same argument as for the first point of this proposition. We start with a smooth function α, that is α ∈ H

For the second term of the r.h.s we use Lemma A.2-c-with s = -1 4 + ε, s ′ = 1 2 -ε and Lemma B.1

This achieves the proof of the continuity of w in time for α ∈ H loc we do like the first point of the proposition: we smooth out α and we use the continuity dependence of q with respect to α stated in Proposition 2.2-2) combined with the estimate (2.14). By writing i∂ t u = H α(t) u(t) we get that for every t ∈ R, u(t) ∈ D(H α(t) ).

A Paraproducts and product laws

The aim of this section is to prove some product laws used in the proof of the main results. For this purpose we first recall some basic ingredients of the paradifferential calculus. Start with the dyadic partition of the unity: there exists two radial positive functions χ ∈ D(R) and ϕ ∈ D(R\{0}) such that χ(ξ) + q≥0 ϕ(2 -q ξ) = 1, ∀ξ ∈ R.

For every tempered distribution v ∈ S ′ , set ∆ -1 v = χ(D)v ; ∀q ∈ N, ∆ q v = ϕ(2 -q D)v and S q = q-1 j=-1 ∆ j .

For more details see for instance [4][3]. Then Bony's decomposition of the product uv is given by


Lucas Fourel

Jean-Philippe Noyel

Etienne Bossy

Xavier Kleber

Philippe Sainsot

Fabrice Ville
email: fabrice.ville@insa-lyon.fr

Towards a grain-scale modeling of crack initiation in rolling contact fatigue -Part 1: Shear stress considerations

Keywords:

Introduction

In many engineering applications, fatigue is the primary failure mechanism that governs the life of a component. During the 19th century, it was thought that the sudden failure of metal railway axles was due to the crystallization of the metal, but this lead quickly gave way to another theory: fatigue. In 1870, Wöhler summarized his work on railroad axles. He concludes that cyclic stress range is more important than peak stress and introduces the concept of endurance limit [START_REF] Wöhler | Über die festigkeitsversuche mit eisen und stahl[END_REF]. Since then, numerous experiments and studies have investigated this mechanism of failure. Basquin (1910) [START_REF] Basquin | The exponential law of endurance tests[END_REF], Palmgren and Miner (1945) [START_REF] Pålmgren | Die lebensdauer von kugellagern[END_REF][START_REF] Wilkins | Cumulative damage in fatigue[END_REF], Coffin and Manson (1954) [START_REF] Coffin | A study of the effects of cyclic thermal stresses on a ductile metal[END_REF][START_REF] Manson | Behavior of materials under conditions of thermal stress[END_REF] proposed rules to predict the number of cycles prior to fatigue failure based on stress-life or strain-life relationships.

Rolling Contact Fatigue (RCF) refers to parts that are subjected to repeated contact pressure such as rolling bearings or gears. RCF is characterized by the initiation of a subsurface micro-crack and results in irreversible surface damages such as micropitting, wear or spalling [START_REF] Olver | The mechanism of rolling contact fatigue: An update[END_REF]. In order to predict the service life of these components and to understand the mechanisms causing their failure, empirical and numerical models were developed. Many of them are described in Sadeghi's review [START_REF] Sadeghi | A Review of Rolling Contact Fatigue[END_REF]. However, there is currently no complete predictive life model, and understanding RCF failure mechanism remains a significant challenge. RCF is characterized by multiple damage stages: crack initiation, crack propagation and material failure [9]. Microstructure is known as having an important role especially in the first two stages [START_REF] Sangid | A physically based fatigue model for prediction of crack initiation from persistent slip bands in polycrystals[END_REF].

In 1981, Tanaka and Mura proposed a micromechanical model for crack initiation [START_REF] Tanaka | A Dislocation Model for Fatigue Crack Initiation[END_REF]. Recent studies have proposed implementations of the Tanaka-Mura law in mesoscopic numerical model using Voronoi tessellations. Bruckner-Foit et al. [START_REF] Brückner-Foit | Numerical simulation of microcrack initiation of martensitic steel under fatigue loading[END_REF][START_REF] Bruckner-Foit | On the determination of material parameters in crack initiation laws[END_REF] studied the fatigue crack threshold on martensitic steel and Toyoda et al. [START_REF] Toyoda | Numerical simulation of fatigue crack initiation in thinwalled high strength steel as modeled by voronoi-polygons[END_REF] on thin-walled high strength steel. Hilgendorff et al. [START_REF] Hilgendorff | Simulation of irreversible damage accumulation in the very high cycle fatigue (vhcf) regime using the boundary element method[END_REF] used boundary element method to compute stress field and Briffod et al. [START_REF] Briffod | Fatigue crack initiation simulation in pure iron polycrystalline aggregate[END_REF] investigated crack initiations in pure iron aggregate. They all use numerical methods for stress computation combined with the Tanaka-Mura model for fatigue crack initiation. However, these works only cover uniaxial loading conditions. Zhou et al. [START_REF] Zhou | Micropitting in Rolling and Sliding Contact Under Mixed Lubrication[END_REF][START_REF] Zhou | Surface topography and fatigue life of rolling contact bearing[END_REF] and Cheng et al. [START_REF] Cheng | Micromechanics Modeling of Crack Initiation Under Contact Fatigue[END_REF][START_REF] Cheng | Semi-Analytical Modeling of Crack Initiation Dominant Contact Fatigue Life for Roller Bearings[END_REF] have applied Tanaka-Mura's law to RCF, although they used analytical methods for computing the stress field which restricted the studies to Hertzian contact conditions. Other researchers used Voronoi tessellations in RCF combined with continuum damage mechanics [START_REF] Raje | A discrete damage mechanics model for high cycle fatigue in polycrystalline materials subject to rolling contact[END_REF][START_REF] Paulson | A coupled finite element ehl and continuum damage mechanics model for rolling contact fatigue[END_REF][START_REF] Lorenz | A continuum damage mechanics finite element model for investigating effects of surface roughness on rolling contact fatigue[END_REF], the Dang-Van criterion [START_REF] Bossy | Competition between surface and subsurface rolling contact fatigue failures of nitrided parts: A dang van approach[END_REF][START_REF] Cerullo | Application of dang van criterion to rolling contact fatigue in wind turbine roller bearings under elastohydrodynamic lubrication conditions[END_REF] or Fatemi-Socie criterion [START_REF] Wang | Effects of microstructure on rolling contact fatigue of a wind turbine gear based on crystal plasticity modeling[END_REF] to predict crack initiation. Some studies also considered competition between surface and subsurface initiated failure due to complex contact conditions in RCF [START_REF] Morales-Espejel | A model for rolling bearing life with surface and subsurface survival-tribological effects[END_REF][START_REF] Morales-Espejel | A model for gear life with surface and subsurface survival: Tribological effects[END_REF].

In this study, the Tanaka-Mura model is used in combination with finite element method (FEM) which allows for moving contact pressures with complex contact conditions such as a dented surface. But rough surfaces or ElastoHy-droDynamics (EHL) lubrication could also be considered. A microstructure is randomly generated by a Voronoi tessellation.

The first part of this article will describe the Tanaka-Mura model and material properties. The second part will detail a macroscopic approach which is utilized to study the influence of the different shear stresses, including orthogonal and maximum shear stresses which are often discussed. The final part will detail a mesoscopic approach which takes advantage of the polycrystalline geometry and grain size and distribution.

Nomenclature

Fatigue predictive modeling

Sadeghi et al. [START_REF] Sadeghi | A Review of Rolling Contact Fatigue[END_REF] reviewed the different probabilistic engineering life models for RCF like the Lundberg-Palmgren [START_REF] Lundberg | Dynamic capacity of rolling bearings[END_REF] and Ioannides-Harris [START_REF] Ioannides | A New Fatigue Life Model for Rolling Bearings[END_REF] theories. These early works are based on empirical results and propose formulas to calculate the fatigue life of classical components such as rolling bearings. In Sadeghi's review [START_REF] Sadeghi | A Review of Rolling Contact Fatigue[END_REF], deterministic research life models for RCF are also presented. These models are based on the physical principles of failure. Different types of models have been originally developed to compute the fatigue lives: models which focus on crack initiation and models which focus on crack propagation. The model used in this study is the Tanaka-Mura law. It is based on the micromechanics of fatigue crack initiation and can be easily implemented under its stress-life relationship form. Tanaka and Mura created their micromechanical model in 1981. It predicts the fatigue crack threshold based on the accumulation of dislocations in Persistent Slip Bands (PSBs) [START_REF] Tanaka | A Dislocation Model for Fatigue Crack Initiation[END_REF]. It uses the energy balance of the dislocation structure along the PSB which can lead to crack initiation after a certain number of loading cycles [START_REF] Sangid | The physics of fatigue crack initiation[END_REF].

The Tanaka-Mura model does not consider normal stress but only shear stress. In addition, in the energy balance some simplifications are made concerning the hysteresis loop. Improvements of the model have been suggested by Mura and Nakasone [START_REF] Mura | A theory of fatigue crack initiation in solids[END_REF][START_REF] Mura | A theory of fatigue crack initiation[END_REF] which proposed using Gibb's free energy. Other authors offered variants of this model with the use of additional parameters such as the Burgers vector [START_REF] Wu | A fatigue crack nucleation model for anisotropic materials[END_REF] or the PSB width [START_REF] Chan | A microstructure-based fatigue-crack-initiation model[END_REF]. The original model has nevertheless received a great popularity within the scientific community because it captures the essence of crack initiation via slip and accumulation of dislocations. The Tanaka-Mura model is consistent with three fundamental fatigue concepts [START_REF] Tryon | A reliability-based model to predict scatter in fatigue crack nucleation life[END_REF]:

-The Manson-Coffin [START_REF] Manson | Behavior of materials under conditions of thermal stress[END_REF][START_REF] Coffin | A study of the effects of cyclic thermal stresses on a ductile metal[END_REF] empirical equation where the fatigue life is inversely proportional to the square of the plastic strain amplitude -The Hall-Petch [START_REF] Hall | The deformation and ageing of mild steel: III discussion of results[END_REF][START_REF] Petch | The cleavage strength of polycrystals[END_REF] relation between the grain size and the high-cycle fatigue -The Palmgren-Miner [START_REF] Pålmgren | Die lebensdauer von kugellagern[END_REF][START_REF] Wilkins | Cumulative damage in fatigue[END_REF]

law of damage accumulation for loads with variable amplitudes

Several studies were able to correlate the Tanaka-Mura fatigue crack initiation model with experimental results [START_REF] Bruckner-Foit | On the determination of material parameters in crack initiation laws[END_REF][START_REF] Tryon | A reliability-based model to predict scatter in fatigue crack nucleation life[END_REF][START_REF] Jezernik | Numerical modelling of fatigue crack initiation and growth of martensitic steels[END_REF][START_REF] Kramberger | Extension of the tanaka-mura model for fatigue crack initiation in thermally cut martensitic steels[END_REF]. Cheng and Cheng compared the model to experiments in the case of contact fatigue [START_REF] Cheng | Micromechanics Modeling of Crack Initiation Under Contact Fatigue[END_REF][START_REF] Cheng | Semi-Analytical Modeling of Crack Initiation Dominant Contact Fatigue Life for Roller Bearings[END_REF].

The number of cycles required to initiate a crack is described in Equation 1.

= 8 (1 -) (Δ -2 ) 2 (1)
= (Δ -) 2 (2)
Using and as material constants:

= 8 (1 -) (3) = 2 (4)
Different values of the material properties can be found in the literature, therefore, Table 1 compiles ranges given for bearing steels [START_REF] Bruckner-Foit | On the determination of material parameters in crack initiation laws[END_REF][START_REF] Toyoda | Numerical simulation of fatigue crack initiation in thinwalled high strength steel as modeled by voronoi-polygons[END_REF][START_REF] Hilgendorff | Simulation of irreversible damage accumulation in the very high cycle fatigue (vhcf) regime using the boundary element method[END_REF][START_REF] Tryon | A reliability-based model to predict scatter in fatigue crack nucleation life[END_REF][START_REF] Jezernik | Numerical modelling of fatigue crack initiation and growth of martensitic steels[END_REF][START_REF] Kramberger | Extension of the tanaka-mura model for fatigue crack initiation in thermally cut martensitic steels[END_REF]. Therefore, parameters and which are functions of the material properties can vary according to Table 2.

Table 2

Ranges of material constants.

Parameter

Minimum Maximum The number of cycles required to initiate a fatigue crack depends on two variables: the length of the persistent slip band (which can be approximate to the grain size) and the shear stress range Δ . Equation 2 can be transformed in order to give Δ as a function of and then to plot a S-N curve that gives the shear stress range required to initiate a crack for a certain number of cycles.

Δ = + √ (5)
Figure 1 shows the influence of material constants and on the S-N fatigue crack initiation curve for an average grain size of 25

(which corresponds to a standard bearing steel).

The dashed lines correspond to the lower and upper bounds using respectively the most critical and values while the colored curves only show the influence of one parameter by setting the second to an intermediate value. has more influence at low-cycle fatigue. This parameter is related to rigidity and fracture energy of the material. A larger corresponds to a material that is strongly resistant to fatigue crack initiation.

acts as a threshold parameter. If the shear stress range Δ does not overcome , the fatigue crack threshold is infinite. Indeed, this constant is directly correlated to the frictional stress which is the resistance that a dislocation needs to overcome in order to move.

In order to implement this fatigue model into a numerical RCF model, shear stress range needs to be calculated. In this study two approaches are compared: macroscopic and mesoscopic.

Macroscopic modeling

Model

The model is designed to represent a contact between two bodies in 2D. Only one of the bodies is modeled. The action of the second one is emulated by a moving contact pressure. This method is commonly used in the literature [START_REF] Bossy | Competition between surface and subsurface rolling contact fatigue failures of nitrided parts: A dang van approach[END_REF][START_REF] Noyel | Development of a granular cohesive model for rolling contact fatigue analysis: Crystal anisotropy modeling[END_REF][START_REF] Šraml | Numerical procedure for predicting the rolling contact fatigue crack initiation[END_REF][START_REF] Raje | A discrete damage mechanics model for high cycle fatigue in polycrystalline materials subject to rolling contact[END_REF].

Hertz contact theory [START_REF] Hertz | Über die Berührung fester elastischer Körper und über die Härte[END_REF] provides tools to calculate the pressure field in contacts. The contact pressure distribution ( ), the half-width of the contact area and the maximum contact pressure 0 are described by Johnson [START_REF] Johnson | Contact mechanics and the wear of metals[END_REF]. The objective is to obtain the stress distribution in the material. In order to be computationally efficient, there is no need for a high level of accuracy in a very large area, only the area of stress concentration is of interest. Generally, in the case of a Hertzian contact, the area of maximum stress is between the surface and a depth of 2 . In the direction of the contact pressure displacement, 2.5 was found to be sufficient to cover the width to the stress concentration. Ideally, this area of analysis of 2 × 2.5 at the surface would be embedded into a semi-infinite body so that there is no effect of the boundary conditions on the stress field. With fixed perpendicular displacements at the left, right and bottom edges, a body of dimension 50 ×25 is equivalent to a semi-infinite body (the stress difference being negligible) [START_REF] Noyel | Analyse de l'initiation de fissures en fatigue de contact : Approche mésoscopique[END_REF]. Figure 2 shows model dimensions, boundary conditions and moving contact pressure positions.

The relative displacement of the bodies is modeled by the incremental displacement of the contact pressure from its initial position to its final position at each calculation step. The initial position is set to 0 = -4 for the center of the contact pressure field and its final position is set to = 4 . This loading condition is representative for perfectly smooth linear contact [START_REF] Johnson | Contact mechanics and the wear of metals[END_REF]. In RCF, surface roughness, manufacturing defects, dents caused by particle contamination or EHL lubrication strongly affect the standard Hertz contact pressure. These factors can introduce stress concentrations in the material and thus, modify the location of the critical areas [START_REF] Johnson | Contact mechanics and the wear of metals[END_REF]. Nowadays, it is well accepted that considering surface fatigue, for rolling bearing applications, surface roughness of the components being generally low, dents appear to be more severe [START_REF] Coulon | Effect of a dent on the pressure distribution in dry point contacts[END_REF][START_REF] Coulon | An abacus for predicting the rolling contact fatigue life reduction due to debris dents[END_REF][START_REF] Coulon | Pressure profiles measured within lubricated contacts in presence of dented surfaces. comparison with numerical models[END_REF]. Whereas, for gear applications, as surface roughness is important compared to film thickness, it can be the predominant driver of fatigue failure. To emphasis the impact of surface defects compared to smooth surfaces, the authors choose to consider dents in this study. However, similar analyses considering surface roughness are possible.

Coulon et al. characterized dent geometry with three parameters: dent diameter , dent depth ℎ and shoulder height is determined by [START_REF] Coulon | Effect of a dent on the pressure distribution in dry point contacts[END_REF]. The shape function of a dent on the ( , ) plane is:

( , ) = -ℎ exp - 2 + 2 4 2 cos √ 2 + 2 (6)
Since the model is two-dimensional, the dent shape is considered as the section of the shape function at = 0. Figure 3 illustrates the resulting 2D dent shape.

As a comparison between a perfectly smooth surface and a more critical case, the simulation is performed with both loading conditions.

-Hertzian contact: = ℎ -Dented contact: = ℎ + Δ Δ has been determined by Coulon et al. [START_REF] Coulon | Effect of a dent on the pressure distribution in dry point contacts[END_REF][START_REF] Coulon | An abacus for predicting the rolling contact fatigue life reduction due to debris dents[END_REF] depending on the dent characteristics. Half-width of the contact , maximum contact pressure 0 and dent characteristics used in this study are described in Table 3. Dent is symmetric and is considered to be centered at = 0. The contact pressure of the two loading cases are shown in Figure 4, the darker area represents the coordinates of the area of analysis.

Stress components , , and are obtained by FEM computation with plane deformation assumptions. Stress tensor is calculated for each incremental position of the moving contact pressure with quasi-static conditions. The simulation is performed for dry contact, in order to analyze the most severe case of damage [START_REF] Coulon | Effect of a dent on the pressure distribution in dry point contacts[END_REF].

FEM analysis is performed with isotropic and homogeneous material properties. The values that are used in this study are shown in Table 4.

Mesh size has been adapted to be fine in regions of high stress gradient (near dent), coarser in region of lower stress gradient and a lot coarser outside the area of analysis. Therefore, element size is set to:

-5µm between a depth of 0 and 0.5a inside the area of analysis, -10µm between a depth of 0.5a and 2.5a inside the area of analysis, -Progressive element size between edges of the area of analysis and left, right and bottom edges of wide body where element size is prescribed to 5mm.

Mesh sensitivity has been assessed to suppress mesh dependency. The mesh is the same for both smooth and dented surfaces.

Shear stress

The Tanaka-Mura model was originally intended for uniaxial tensile fatigue. In the case of uniaxial loading, the shear stress is maximized for the directions at 45°from the loading direction.

In rolling contact fatigue, the loading is multi-axial, therefore, several approaches can be considered. At least three different critical shear stresses have been used in different RCF models [START_REF] Cheng | Semi-Analytical Modeling of Crack Initiation Dominant Contact Fatigue Life for Roller Bearings[END_REF]:

-Orthogonal shear stress [29] -Maximum shear stress

[49] -Octahedral shear stress [START_REF] Zhou | Surface topography and fatigue life of rolling contact bearing[END_REF] In this study, a rotation of 45°of the stress tensor is performed since this is the direction that maximizes for a Hertzian contact.

is directly computed by the FEM simulation, 45 can be deduced by rotation, and are obtained using the principal coordinate system.

In continuum mechanics, the stress tensor can be diagonalized when expressed in the principal coordinate system.

Maximum shear stress (also known as Tresca shear stress) is defined by:

= 1 2 ( - ) (7)
With and , the maximum and the minimum principal stresses ( ≥ ≥ ). is oriented according to the bisector of the 1st and the 3rd principal directions.

Octahedral shear stress is defined by:

= 1 3 √ ( - ) 2 + ( - ) 2 + ( -) 2 (8)
is oriented according to the plane whose normal vector forms equal angles with the three principal directions.

In the case of Hertzian contact pressure, shears maximum values and critical depths are different depending on the shear stress (see Table 5). For dented surfaces, no analytical law provides the critical depths and the maximum values of the different shear stresses. However, these results can be obtained numerically (see Table 6). Figure 5 shows the different shear stresses for the centered contact pressure. For Hertzian and dented contact cases, stress fields are symmetric ( 45 , , ) or antisymmetric ( ), therefore, the illustrations are split between both cases.

is the most critical stress for centered static contact pressure, making it a commonly used criterion in fatigue. In the following section, shear stress ranges are calculated for a moving contact pressure to simulate the passing of a contact body.

Shear stress range

Shear stress range Δ can be defined as the difference between the maximum and the minimum of shear stress during a full stress cycle.

Figure 6 gives the shear stress ranges for the different shears. Since all stress fields are symmetric, the illustrations are split between Hertzian and dented contact cases. and depend on the principal coordinate system whose orientation can change at each simulation step (see Figure 7). This is not ideal in the case of a moving contact pressure since neither of Δ nor Δ represents the shear stress range on any fixed resolved direction. Moreover, and are continually positive values which results in smaller ranges than fixed coordinate system shear stresses. For example, > but Δ < Δ . For Hertzian contacts, Table 7 gives the depths that maximize the shear stress ranges.

Maximum values and depths of maximum values of shear stresses for the dented contact are given thanks to numerical computations (see Table 8).

It appears in Table 7 and Table 8 that Δ is the most critical shear stress of the four in these configurations. This result is also visible when looking at the variations of the different shear stresses as functions of the contact pressure position (see Figure 8). Maximum and octahedral shears are two critical stresses that are therefore more relevant in case of repeated static pressure. For contact fatigue where the pressure is moving, fixed coordinate system critical stresses should be privileged. This is the case of Δ and Δ 45 .

Fatigue crack threshold

Figure 9 shows the number of cycles resulting from the Tanaka-Mura (Equation 2) using four different stresses. These results are obtained using the parameters given Table 9. In the case of Hertzian contact pressure, the use of Δ gives the critical results (smaller , ). The nondestructive testing by Yoshioka showed that cracks first initiate at the depth of the maximum of Δ [START_REF] Yoshioka | Detection of rolling contact sub-surface fatigue cracks using acoustic emission technique[END_REF]. These results were also found in tests by Chen et al. [START_REF] Chen | Study on initiation and propagation an-gles of subsurface cracks in gcr15 bearing steel under rolling contact[END_REF][START_REF] Chen | Measurement of the critical size of inclusions initiating contact fatigue cracks and its application in bearing steel[END_REF]. These observations confirmed that Δ is more critical than Δ 45 in the case of a smooth surface.

The slip band length is considered as a constant, rep- resentative for the mean grain size. This assumption can be questioned as in a real material, some of the grains are larger than the mean grain size, and some others are smaller. This distribution has a direct effect on fatigue crack threshold. Surface defects such as a dent cause contact pressure to increase and consequently alterations in the stress field. As Figure 9 shows, the presence of a dent induces a variation of , but also a variation of the spatial distribution. The comparison between multiple shear stresses have shown that is more relevant than 45 . However, the orientation that maximize Δ could be different than 0°and 45°. Moreover, this orientation could vary depending on the position in the material.

The use of a polycrystalline geometry would open opportunities of improvement concerning the calculation of both the shear stress range Δ and the PSB length . Different shear stress directions could be calculated for each grain and an algorithm could even optimize these directions in order to consider the most critical cases. The distribution of slip band length would be much more representative of a real material as long as the polycrystalline geometry is realistic.

Mesoscopic modeling

The mesoscopic approach is similar to the macroscopic approach except for the area of analysis (see Figure 2) which is replaced by a granular geometry. The geometry of this area is generated with Neper [START_REF] Quey | Large-scale 3d random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF][START_REF] Quey | Optimal polyhedral description of 3d polycrystals: Method and application to statistical and synchrotron xray diffraction data[END_REF][START_REF] Quey | Nearly uniform sampling of crystal orientations[END_REF]. It uses Voronoi's tessellations to generate polycrystalline geometry. Grain sizes follow a log-normal distribution with a mean value of 25 and a standard deviation of 11 (see Figure 10). Each grain is meshed into quadrangle quadratic elements.

Shear stress range

In this study, mesoscopic shear stress is the shear stress calculated in the fixed direction that maximizes the shear stress range Δ for each grain. Only in-plane shear is considered. The mesoscopic shear stress range Δ is obtained by an algorithm which tests every resolved direction from 0°to 90°with an angular increment of 0.1°. For each grain, the algorithm keeps the direction that maximized the shear stress range. The four components , , and of each element of a grain are averaged to give mean values for each grain. The following results are all grain-related values in opposition of section 3 where the results were elementrelated values.

The main advantage of the mesoscopic method is that the direction of Δ corresponds to the most favorable orientation for dislocation slip is ensured. In the case of the Hertzian contact (Figure 11), there is no significant difference between (a) Δ and (c) Δ because coordinate system is already favorable for maximization of Δ . It confirms that orthogonal shear is relevant in smooth surface cases.

However, there is a 5% difference between the maximum of Δ and the maximum of Δ in a dented surface (see Figure 11 (b,d)). Mesoscopic method is more critical than the uniform resolved direction and stress field distribution is different. High stress concentration zone covers a larger area in case of mesoscopic method, in contrast with Δ where concentration region is smaller.

Orthogonal shear stress is commonly used for RCF, this comparison with mesoscopic shear stress shows that it is valid in case of Hertzian contact, however, results with a dented surface illustrate that it could be less legitimate in other cases.

Fatigue crack threshold

Fatigue crack threshold is computed thanks to Tanaka-Mura equation (see Equation 2). Same material constants and than the macroscopic approach are used (Table 9). Shear stress range Δ calculation is detailed in the previous section and is computed using the granular geometry. For each grain, is the length of the longest segment calculated in the direction of Δ . In Hertzian contact, the value of Δ is constant for a given depth, therefore, the variation of at constant depth is exclusively due to grain size variations.

Figure 12 shows fatigue crack threshold results obtained in polycrystalline geometry using three different methods. is 24082 cycles while it is only 3563 for the dented contact with the mesoscopic approach. As expected, surface conditions have a predominant influence on fatigue crack threshold. A damaged surface will highly decrease .

Computation with moving contact pressure is significantly heavier in terms of CPU time than a simple static contact pressure due to the requirement of numerous calculation steps. Therefore, the benefit of the moving contact pressure compared to a static pressure can be questioned. Indeed, it is possible to determine Δ simply by considering the stressed state (static centered contact pressure for example) and a relaxed state (stresses are nil). Repeated stress cycles correspond to alternations between the two states. It turns out that the static contact pressure is not suitable to approximate a rolling contact since the results obtained with the moving contact pressure are significantly different. Critical initiation depths can be identified Figure 12. The benefit of modeling a polycrystalline geometry is that it introduces a probabilistic aspect of crack initiation. For the Hertzian contact, the combination of Δ and distributions results in a depth range where cracks are more likely to first initiate. Using a threshold of , = 10 × , , this depth ranges from = -1.1 to = -0.2 . This range is consistent with the experimental observations [START_REF] Chen | Study on initiation and propagation an-gles of subsurface cracks in gcr15 bearing steel under rolling contact[END_REF][START_REF] Chen | Measurement of the critical size of inclusions initiating contact fatigue cracks and its application in bearing steel[END_REF]. Concerning the dented contact, cracks first appear under the dent. The critical depth range varies from = -1.0 to = 0 (using , = 10 × , ). This implies that damaged surfaces tend to originate at smaller depths and eventually at the surface itself.

Concerning the choice of the method used to calculate Δ , it should be noted that in Hertzian contact , is 0.3% higher using Δ than using Δ . For the dented contact, , is 24.1% higher using Δ than using Δ . This confirms that orthogonal shear is suitable in smooth surface contact. In the case of damaged surfaces that induce less uni-form contact pressure, the mesoscopic shear stress is more relevant and will provide more critical results.

Conclusions

The numerical model for rolling contact fatigue which is presented in this study leads to several conclusions. [START_REF] Wöhler | Über die festigkeitsversuche mit eisen und stahl[END_REF] Modeling the motion of the contact pressure gives significantly different results compared to a static pressure which can imply that static pressure is not relevant to model crack initiation in rolling contact fatigue. [START_REF] Basquin | The exponential law of endurance tests[END_REF] The orthogonal shear stress is a relevant indicator of the resolved shear stress for most favorable orientations in smooth surface cases. (3) A mesoscopic approach should be preferred in rough or damaged surface cases. The resolved direction of each grain is the one that maximizes Δ . (4) In Hertzian contact cases, cracks first initiate at the depth of maximum orthogonal shear stress of = -0.5 . ( 5) In damaged surface cases, cracks first initiate at smaller depth values and eventually surface itself.

Even though material properties are actually a drawback to obtain accurate results, the modeling can replicate the physics of fatigue crack initiation and show the most critical locations in different cases. It may be useful for engineers to know the depth at which a surface treatment of the material is beneficial. However, the model described in this article does not replicate all the physics involved in fatigue crack initiation. Firstly, the material model is assumed as homogeneous and isotropic. However, polycrystalline materials present anisotropic characteristics. Anisotropic simulation would allow to model the slip bands that form inside the polycrystalline material. Indeed, intra-granular and inter-granular slip bands are the physical phenomena that led Tanaka and Mura to develop their equation. Although slip band length is approximated by the grain size in this study, the modeling of the PSB field would provide more physically based results. All these phenomena are as many interesting aspects to pursue in order to try to approach a better modeling of crack initiation.

Δ

  Dent contribution to contact pressure [ ] Δ Shear stress range [ ] Poisson's ratio Dent diameter [ ] Normal stress in direction [ ] Normal stress in direction [ ] Normal stress in direction [ ] 45 Orthogonal shear stress rotated by 45°from direction [ ] Maximum shear stress [ ] Mesoscopic shear stress [ ] Octahedral shear stress [ ] Orthogonal shear stress in direction [ ]

-

  Shear modulus, [ ] -Fracture energy, [ ∕ 2 ] -Poisson's ratio, [-] -Slip band length, [ ] -Δ Shear stress range, [ ] -Frictional stress, [ ] Tanaka-Mura equation can be simplified as:

Figure 1 :

 1 Figure 1: S-N curves depending on (a) variation of material constant and (b) variation of material constant .

Figure 2 :

 2 Figure 2: Model dimensions and contact pressure position.

Figure 3 :

 3 Figure 3: Analytical dent shape and parameters.

Figure 4 :

 4 Figure 4: Contact pressure for load centered on (blue) = -1.3 , (red) = 0, (green) = 1.3 , for (a) Hertzian contact (b) dented contact.

Figure 5 :

 5 Figure 5: (a) Orthogonal shear stress , (b) orthogonal 45°shear stress 45 , (c) maximum shear stress , (d) octahedral shear stress ; centered contact pressure.

Figure 6 :

 6 Figure 6: (a) Orthogonal shear stress range Δ , (b) orthogonal 45°shear stress range Δ 45 , (c) maximum shear stress range Δ , (d) octahedral shear stress range Δ .

Figure 7 :

 7 Figure 7: Principal coordinate system orientation relative to fixed coordinate system at = 0 and = -0.78 .

Figure 8 :

 8 Figure 8: Variation of the different shear stresses for Hertzian contact at = 0 and = -ℎ .

Figure 9 :

 9 Figure 9: Number of cycles for crack initiation using (a) Δ , (b) Δ 45 , (c) Δ , (d) Δ .

Figure 10 :

 10 Figure 10: Polycrystalline geometry used for the area of analysis (generated with Neper [53, 54, 55]).

Figure 11 :

 11 Figure 11: Comparison between (a,b) orthogonal shear stress range Δ and (c,d) mesoscopic resolved shear stress range Δ for (a,c) hertzian contact pressure and (b,d) dented surface.

  (a) and (b) show the results for a static centered contact pressure and the use of , (c) and (d) show the results for a moving contact pressure and the use of orthogonal shear stress range Δ , and (e) and (f) shows the results for a moving contact pressure and the use of mesoscopic shear stress range Δ . For the Hertzian contact, the minimum of the number of cycles ,

Figure 12 :

 12 Figure 12: Number of cycles required to initiate fatigue crack for (a,b) static load using , (c,d) moving contact pressure using Δ , (e,f) moving contact pressure using Δ , (a,c,e) Hertzian contact pressure and (b,d,f) dented contact pressure.

Table 1

 1 Ranges of material properties.

	Property	Minimum	Maximum
	[	]	80	90
	[-]		0.27	0.31
	[ ∕ 2 ]	2.0	440
	[	]	100	500

Table 3

 3 Contact and dent characteristics.

	Characteristic		Value
	Half-width [ ]		615
	Max contact pressure 0 [	]	2100
	Dent diameter [ ]		250
	Dent depth ℎ [ ]		37
	Dent shoulder parameter	[-]	12

Table 4

 4 Material properties of the simulation.

	Property		Value
	Young's modulus [	]	200
	Poisson's ratio [-]		0.3

Table 5

 5 Max and depth at max of different shear stresses for Hertzian contact[START_REF] Johnson | Contact mechanics and the wear of metals[END_REF].

	Shear stress	Depth at max	Max
		0.50	0.25 0
	45	0.78	0.30 0
		0.78	0.30 0
		0.67	0.26 0

Table 6

 6 Max and depth at max of different shear stresses for dented contact.

	Shear stress	Depth at max	Max
		0.069	0.74 0
	45	0.078	0.71 0
		0.062	0.78 0
		0.061	0.70 0

Table 7

 7 Max and depth at max of different shear stress ranges for Hertzian contact.

	Shear stress range	Depth at max	Max
	Δ	0.50	0.50 0
	Δ 45	0.63	0.37 0
	Δ	0.78	0.30 0
	Δ	0.67	0.26 0

Table 8

 8 Max and depth at max of different shear stress ranges for dented contact.

	Shear stress range	Depth at max	Max
	Δ	0.021	1.2 0
	Δ 45	0.068	0.82 0
	Δ	0.061	0.80 0
	Δ	0.059	0.71 0

Table 9

 9 Macroscopic approach parameters.

	Parameter	Value
	[ 2 ∕ 3 ]	1.28 × 10 17
	[ ]	7.8 × 10 8
	[ ]	25