Antoine Ayache
email: antoine.ayache@univ-lille.fr

Harmonizable Fractional Stable Motion: simultaneous estimators for the both parameters

Keywords: Estimation of Harmonizable Fractional Stable Motion AMS Subject Classification (MSC2020 database): 60G22, 60F15, 60H05 Heavy-tailed stable distributions, hamonizable fractional processes, non-ergodicity, laws of large numbers, discrete wavelet transforms

There are two classical very different extensions of the well-known Gaussian fractional Brownian motion to non-Gaussian frameworks of heavy-tailed stable distributions: the harmonizable fractional stable motion (HFSM) and the linear fractional stable motion (LFSM). As far as we know, while several articles in the literature, some of which appeared a long time ago, have proposed statistical estimators for the parameters of LFSM, no estimator has yet been proposed in the framework of HFSM. Among other things, what makes statistical estimation of parameters of HFSM to be a difficult problem is that, in contrast to LFSM, HFSM is not ergodic. The main goal of our work is to propose a new strategy for dealing with this problem and obtaining solutions of it. The keystone of our new strategy consists in the construction of new transforms of HFSM which allow to obtain, at any dyadic level, a sequence of independent stable random variables.

Introduction and statement of the two main results

A real-valued harmonizable fractional stable motion 1 (HFSM), denoted by {X(t)} t∈R , is a paradigmatic example of a continuous symmetric stable self-similar stochastic process with stationary increments. It was introduced, about 35 years ago, by Cambanis and Maejima in [START_REF] Cambanis | Two classes of self-similar stable processes with stationary increments[END_REF]. Basically, it depends on two parameters: the Hurst parameter H belonging to the open interval (0, 1), and the stability parameter α belonging to the interval (0, 2]. Among other 1 Notice that the HFSM is somtimes called harmonizable fractional stable process.

1 things, the parameter H governs roughness of sample paths of {X(t)} t∈R and its self-similarity property:

a -H X(at) t∈R d = {X(t)} t∈R , for any fixed a ∈ (0, +∞), where d = means that the equality holds in the sense of the finite-dimensional distributions. While the parameter α determines heaviness of tails of marginal distributions of {X(t)} t∈R : except in the very particular Gaussian case α = 2 in which the probability P |X(t)| ≥ z vanishes exponentially fast when z → +∞, for any other value of α, one has, for each fixed t = 0 (notice that X(0) a.s. = 0), c (t)z -α ≤ P |X(t)| ≥ z ≤ c (t)z -α , for all z ∈ [1, +∞), where 0 < c (t) < c (t) are two finite constants.

The HFSM {X(t)} t∈R is defined, for all t ∈ R, through the stable stochastic integral in the frequency domain:

X(t) := Re R e itξ -1 |ξ| H+1/α d M α (ξ) , (1.1)
where M α is a complex-valued rotationally invariant α-stable random measure with Lebesgue control measure. A detailed presentation of such a random measure and the corresponding stable stochastic integral and related topics can for instance be found in Chapter 6 of the well-known book [START_REF] Samorodnitsky | Stable Non-Gaussian Processes: Stochastic Models with Infinite Variance[END_REF]. The following remark, which provides two very important properties of this stochastic integral, will play a fundamental role in our work. α dξ.

(1.

2)

The equality (1.2) is reminiscent of the classical isometry property of Wiener integrals; in particular, it implies that Re R g n (ξ) d M α (ξ) converges to Re R g(ξ) d M α (ξ) in probability, when a sequence (g n ) n converges to g in L α (R).

(ii) Let m ∈ N be arbitrary and let f 1 , . . . , f m be arbitrary functions of L α (R) whose supports are disjoint up to Lebesgue negligible sets, then the real-valued SαS random variables Re R f 1 (ξ) d M α (ξ) , . . . , Re R f m (ξ) d M α (ξ) are independent.

In the very particular Gaussian case where the stability parameter α = 2, the HFSM represented by (1.1) reduces to the very classical Gaussian fractional Brownian motion (FBM) with Hurst parameter H, denoted by {B H (t)} t∈R . It is well-known that, up to a deterministic multiplicative constant, the Gaussian process {B H (t)} t∈R can also be represented as a moving average stochastic Wiener integral in the time domain, whose integrand is no longer the complexvalued function ξ → e itξ -1 |ξ| H+1/α but the real-valued function s → (t -s)

H-1/2 + -(-s) H-1/2 + .
One recalls, in passing, the usual convention that, for all (x, β) ∈ R 2 , (x) β + := x β when x > 0 and (x) β + := 0 else.

When the stability parameter α = 2, the HFSM in (1.1) can no longer be represented as a moving average stable stochastic integral in the time domain. Actually, it is very different from the real-valued linear fractional stable motion (LFSM) {L(t)} t∈R defined, for each t ∈ R, by

L(t) := R (t -s) H-1/α + -(-s) H-1/α + dM α (ds), (1.3)
where M α is a real-valued α-stable random measure. The large differences between the two processes {X(t)} t∈R and {L(t)} t∈R can be explained by several reasons. Let us mention two important ones of them: (1) in contrast to the process {L(t)} t∈R the process {X(t)} t∈R is not ergodic, (2) behavior of sample paths of {L(t)} t∈R and {X(t)} t∈R is far from being the same. Indeed, sample paths of {L(t)} t∈R are multifractal functions (see [START_REF] Balança | Fine regularity of Lévy processes and linear (multi)fractional stable motion[END_REF]) which become discontinuous when H ≤ 1/α and even unbounded on any interval when H < 1/α (see for example [START_REF] Samorodnitsky | Stable Non-Gaussian Processes: Stochastic Models with Infinite Variance[END_REF][START_REF] Embrechts | Selfsimilar processes[END_REF]). While those of {X(t)} t∈R are Hölder continuous of any order strictly less than H for any value of H ∈ (0, 1) (see [START_REF] Kôno | Self-similar stable processes with stationary increments[END_REF][START_REF] Kôno | Hölder continuity of sample paths of some self-similar stable processes[END_REF][START_REF] Samorodnitsky | Stable Non-Gaussian Processes: Stochastic Models with Infinite Variance[END_REF][START_REF] Embrechts | Selfsimilar processes[END_REF]), namely, for each fixed δ > 0 and T > 0, one has almost surely

sup -T ≤t <t ≤T X(t ) -X(t ) |t -t | H-δ < +∞ . (1.4)
Moreover, sample paths of {X(t)} t∈R are monofractal functions; the latter fact results from their Hölderianity property combined with Corollary 4.4 in [START_REF] Ayache | Harmonizable fractional stable fields: local nondeterminism and joint continuity of the local times[END_REF]. Also, for later purposes, one mentions that as regards their behavior at infinity, one can derive from Corollary 4.2 in [START_REF] Ayache | Stationary increments harmonizable stable fields: upper estimates on path behaviour[END_REF], that, for all fixed δ > 0, one has, almost surely, sup |t|≥1 X(t) |t| H+δ < +∞ .

(1.5)

Let us now present the main motivations behind our present work and its main contributions. Statistical estimators for the parameters H and α of the LFSM {L(t)} t∈R and related moving average stable processes have been proposed in several articles in the literature (see for instance [START_REF] Stoev | Estimation of the self-similarity parameter in linear fractional stable motion[END_REF][START_REF] Pipiras | Bounds for the covariance of functions of infinite variance stable random variables with applications to central limit theorems and wavelet-based estimation[END_REF][START_REF] Ayache | Linear mutifractional stable motion: wavelet estimation of H(•) and α parameters[END_REF][START_REF] Ayache | Uniformly and strongly consistent estimation for the Hurst function of a linear multifractional stable motion[END_REF][START_REF] Dang | Estimation of the Hurst and the stability indices of a H-selfsimilar stable process[END_REF][START_REF] Mazur | Estimation of the linear fractional stable motion[END_REF][START_REF] Ljungdahl | A minimal contrast estimator for the linear fractional stable motion[END_REF][START_REF] Ljungdahl | Multidimensional parameter estimation of heavy-tailed moving averages[END_REF]), some of which appeared a long time ago. However, as far as we know, in the framework of the HFSM {X(t)} t∈R and related harmonizable stable processes and fields no statistical estimator for any one of these two parameters has yet been proposed in the literature. Also, according to what is mentioned in Remark 1.2 (D) and on page 2 of the very recent preprint [START_REF] Basse-O'connor | Asymptotic theory for quadratic variation of harmonizable fractional stable processes[END_REF], their statistical estimation in such a framework is far from being an obvious problem due to the fact that HFSM and related harmonizable stable processes and fields are not ergodic. The main idea behind our strategy for dealing with the latter problem is to construct new transforms of HFSM which allow to obtain, at any dyadic level j ∈ N, a sequence Y j,k , k ∈ N, of independent real-valued SαS random variables whose scale parameters σ(Y j,k ), k ∈ N, are closely connected to the unknown parameters H and α of HFSM through simple formulas which are rather easy to handle. Roughly speaking, these new transforms Y j,k , (j, k) ∈ N 2 , of HFSM are at the same time inspired by discrete wavelet transforms W j,k , (j, k) ∈ N 2 , of HFSM and significantly different from them. Indeed, while W j,k is defined (see e.g. [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF][START_REF] Meyer | Wavelets and Operators[END_REF]), up to normalizing factor, as

W j,k := R ψ j,k (t)X(t) dt,
where ψ j,k (t) := ψ(2 j t -k), for all t ∈ R, with ψ being a "nice" real-valued function; we define Y j,k as

Y j,k := 2 2π R Re ψ j,k (t) X(t) dt,
where ψ j,k is the Fourier transform of ψ j,k . Also, we assume that the function ψ is given by the formula (2.4); we mention that many other choices for ψ are possible, and there is no need to impose to ψ to have any vanishing moment, while such a condition on moment(s) of ψ plays a crucial role in the case of the discrete wavelet transform W j,k . Notice that the definition of Y j,k is given in a more explicit way in (2.8).

Let us now state the main two theorems of our work.

Theorem 1.2 One assumes that α ∈ [α, 2], where the lower bound α ∈ (0, 2] is known. Let (m j ) j∈N be an arbitrary non-decreasing sequence (that is m j ≤ m j+1 , for all j ∈ N) of integers larger than 2 which satisfy the condition m j ≥ j , for all j ∈ N.

(1.6)

For each fixed γ ∈ (0, 4 -1 α), and for every j ∈ N, the statistics V m j j,γ is defined as

V m j j,γ := m j k=1 Y j,k γ , (1.7)
where the Y j,k 's are the real-valued SαS random variables defined in (2.8). Then, the following two results hold.

(i) Let

H j,γ := log 2 (V m j j,γ ) γ j , (1.8)
where log 2 denotes the binary logarithm (that is log 2 (2) = 1). Under the condition

lim j→+∞ log 2 (m j ) j = 0, (1.9)
H j,γ is a strongly consistent (almost surely convergent) estimator of the Hurst parameter H of HFSM.

(ii) Let

α j,γ := γ j j -log 2 (V m j j,γ )
.

(1.10)

Under the condition

lim j→+∞ log 2 (m j ) j = 1, (1.11)
α j,γ is a strongly consistent (almost surely convergent) estimator of the stability parameter α of HFSM.

Theorem 1.3 Let (m j ) j∈N be an arbitrary non-decreasing sequence of integers larger than 2 which satisfies the condition (1.6). For every j ∈ N, the statistics V m j j,log 2 is defined as

V m j j,log 2 := m j k=1 log 2 |Y j,k , (1.12)
where the Y j,k 's are the real-valued SαS random variables defined in (2.8). Then, the following two results hold.

(i) Let

H j,log 2 := V m j j,log 2 j m j . (1.13)
Under the condition (1.9), H j,log 2 is a strongly consistent (almost surely convergent) estimator of the Hurst parameter H of HFSM.

(ii) Let

α j,log 2 := j m j V m j j,log 2 . (1.14)
Under the condition (1.11), α j,log 2 is a strongly consistent (almost surely convergent) estimator of the stability parameter α of HFSM.

Remark 1.4

(i) We believe that the new strategy introduced in our present work would open the door to statistical estimation of parameters of harmonizable stable fields extending the HFSM, as for instance the harmonizable fractional stable field studied in e.g. [START_REF] Ayache | Harmonizable fractional stable fields: local nondeterminism and joint continuity of the local times[END_REF], or the harmonizable fractional stable sheet studied in e.g. [START_REF] Ayache | Wavelet series representation and geometric properties of harmonizable fractional stable sheets[END_REF].

(ii) Since, for each for each fixed j ∈ N, the random variables Y j,k , k ∈ N, have the very nice property to be independent, we believe that the centered and standardized versions of the two statistics V m j j,γ and V m j j,log 2 would be asymptotically normal. This would imply that the four estimators H j,γ , α j,γ , H j,log 2 and α j,log 2 would be asymptotically normal as well.

(iii) In our present work, the four estimators H j,γ , α j,γ , H j,log 2 and α j,log 2 are obtained from the observation of a sample path of the HFSM X in continuous time, we believe that it would be possible to extend our estimation procedures to a framework where only a discretized sample path of X is observed.

We intend to study these three issues (i), (ii) and (iii) in future works.

The remaining of our present work is organized in the following way. In Section 2, basically we show that the real-valued SαS random variables Y j,k , (j, k) ∈ N 2 , defined in (2.8), can be represented in terms of the stable stochastic integral

R • d M α (see Lemma 2.
2); a straightforward consequence of this representation is the independence property, for any fixed j ∈ N, of the random variables Y j,k , k ∈ N (see Remark 2.3). Section 3 is devoted to the proofs of Theorems 1.2 and 1.3.

The keystone

Remark 2.1 Throughout our work, we use the very classical convention that, for any function f ∈ L 1 (R), the Fourier transform f , also denoted by F(f ), is defined as

F(f )(t) = f (t) := R e -itx f (x) dx, for all t ∈ R. (2.1)
While, the inverse Fourier transform of f , denoted by F -1 (f ), is defined as

F -1 (f )(ξ) := 1 2π R e iξt f (t) dt, for all ξ ∈ R. (2.2)
It is well-known that the two maps F and F -1 can be extended to L 2 (R), and satisfy

F -1 F(g) = F F -1 (g) = g, for every g ∈ L 2 (R). (2.3)
One denotes by ψ the even piecewise linear continuous function from R into [0, 2] with defined as

ψ(ξ) := 4 1l [-4 -1 ,4 -1 ] * 1l [-4 -1 ,4 -1 ] )(ξ) = 2 1 -|2ξ| 1l [-2 -1 ,2 -1 ] (ξ), for all ξ ∈ R, (2.4)
where " * " denotes the convolution product. It is clear that ψ is compactly supported and that

supp ψ = [-2 -1 , 2 -1 ]. (2.5)
Moreover, it follows from elementary calculations and a classical property of Fourier transform that

ψ(t) = 4 -1 t -2 sin 2 (4 -1 t), for every t ∈ R, (2.6)
which in particular implies that there is a finite deterministic constant c > 0 such that

ψ(t) ≤ c 1 + |t| -2 , for every t ∈ R. (2.7)
For all (j, k) ∈ N 2 , the random variable Y j,k is defined through the pathwise Lebesgue integral

Y j,k := 2 1-j 2π R cos(2 -j kt) ψ(2 -j t)X(t) dt . (2.8)
Notice that the fact that the latter integral is well-defined and finite results from (2.7), (1.5) and continuity of sample paths of the HFSM {X(t)} t∈R .

Lemma 2.2 For all (j, k) ∈ N 2 , Y j,k (defined in (2.8
)) is a real-valued SαS random variable which can almost surely be expressed as

Y j,k = Re R ψ j,k (ξ) |ξ| H+1/α d M α (ξ) , (2.9)
where

ψ j,k (ξ) := ψ(2 j ξ + k) + ψ(2 j ξ -k), for every ξ ∈ R. (2.10)
The following remark is very fundamental.

Remark 2.3 One knows from (2.10) and (2.4) that, for all (j, k) ∈ N 2 , ψ j,k is an even piecewise linear continuous function from R into [0, 2] with compact support such that

I j,k := supp ψ j,k = -k -2 -1 2 j , -k + 2 -1 2 j k -2 -1 2 j , k + 2 -1 2 j . (2.11)
Then, in view of the fact that k ≥ 1, it turns out that the function ξ → |ξ| -H-1/α ψ j,k (ξ) belongs to the Lebegue space L α (R), with guarantees that the SαS stochastic integral in (2.9) is well-defined. More importantly, observe that, for each fixed j ∈ N, since (see (2.11)) the supports of the functions ψ j,k , k ∈ N, are disjoint (up to Lebesgue negligible sets), the corresponding (via (2.9)) SαS random variables Y j,k , k ∈ N are independent.

Proof of Lemma 2.2 Throughout the proof the positive integers j and k are arbitrary and fixed. First observe that, one can derive from (2.10) and standard calculations that the Fourier transform of ψ j,k is given by

ψ j,k (t) = 2 1-j cos(2 -j kt) ψ(2 -j t) for all t ∈ R.
(2.12) Thus, (2.8) reduces to

Y j,k = 1 2π R ψ j,k (t)X(t) dt . (2.13)
Moreover, one can derive from (2.12) and (2.7), that, for some finite deterministic constant c 1 (j) > 0, only depending on j, one has

ψ j,k (t) ≤ c 1 (j) 1 + |t| -2 , for every t ∈ R. (2.14)
For all (n, m) ∈ N × Z, the dyadic number d n,m := 2 -n m. Having introduced the latter notation, for each n ∈ N, the random variable Y n j,k is defined as the finite sum:

Y n j,k := 1 2π |m|≤4 n X(d n,m ) d n,m+1
dn,m

ψ j,k (t) dt (2.15) = Re   R |ξ| -H-1/α   1 2π |m|≤4 n d n,m+1 dn,m e idn,mξ -1 ψ j,k (t) dt   d M α (ξ)   ,
where the last equality follows from (1.1).

First step: we show that, when n goes to +∞, the random variable Y n j,k converges almost surely to the random random variable Y j,k defined through (2.13).

Let N be an arbitrary fixed positive integer. One sets

Y N j,k := 1 2π N -N ψ j,k (t)X(t) dt = 1 2π N 2 n -1 m=-N 2 n d n,m+1
dn,m

X(t) ψ j,k (t) dt (2.16)
and, for every integer n ≥ N ,

Y N,n j,k := 1 2π N 2 n -1 m=-N 2 n X(d n,m ) d n,m+1
dn,m ψ j,k (t) dt.

(2.17)

Then, using (1.4) with T = N , and the fact that

d n,m ∈ [-N, N ] when -N 2 n ≤ m ≤ N 2 n -1,
one gets, on an event of probability 1 depending only on N and denoted by Ω * N , that

Y N j,k -Y N,n j,k ≤ N 2 n -1 m=-N 2 n d n,m+1 dn,m X(t) -X(d n,m ) ψ j,k (t) dt ≤ C * N,δ N 2 n -1 m=-N 2 n d n,m+1 dn,m |t -d n,m | H-δ ψ j,k (t) dt ≤ C * N,δ R ψ j,k (t) dt 2 -n(H-δ) , (2.18)
where δ is an arbitrarily small fixed positive real number, and C * N,δ is a positive finite random variable, only depending on N and δ. Next, let Ω * be the event of probability 1 defined as

Ω * := N ∈N Ω * N . One can derive from (2.18) that lim n→+∞ Y N j,k (ω) -Y N,n j,k (ω) = 0, for all N ∈ N and ω ∈ Ω * . (2.19)
On another hand, it follows from (1.5) that there are Ω * * an event of probability 1 and C * * δ a positive finite random variable only depending on δ, such that

X(t, ω) ≤ C * * δ (ω)
Y n j,k (ω) -Y N,n j,k (ω) ≤ +∞ m=N 2 n X(d n,m , ω) d n,m+1
dn,m

ψ j,k (t) dt + -N 2 n -1 m=-∞ X(d n,m , ω) d n,m+1
dn,m

ψ j,k (t) dt ≤ C * * δ (ω) +∞ m=N 2 n d n,m H+δ d n,m+1
dn,m

ψ j,k (t) dt + -N 2 n -1 m=-∞ d n,m H+δ d n,m+1
dn,m

ψ j,k (t) dt ≤ C * * δ (ω) +∞ m=N 2 n d n,m+1 dn,m |t| H+δ ψ j,k (t) dt + -N 2 n -1 m=-∞ d n,m+1 dn,m 1 + |t| H+δ ψ j,k (t) dt ≤ 2c 1 (j)C * * δ (ω) +∞ N 1 + t H+δ-2 dt ≤ 2(1 -H -δ) -1 c 1 (j)C * * δ (ω)N -(1-H-δ) . ( 2
Y j,k (ω) -Y N j,k (ω) ≤ {|t|≥N } ψ j,k (t) X(t) dt ≤ 2c 1 (j)C * * δ (ω) +∞ N 1 + t H+δ-2 dt ≤ 2(1 -H -δ) -1 c 1 (j)C * * δ (ω)N -(1-H-δ) . (2.22)
Finally, observe that, for all positive integers N and n ≥ N , and for each ω ∈ Ω * ∩ Ω * * (the event Ω * ∩ Ω * * is clearly of probability 1), using the triangle inequality, (2.21) and (2.22), one gets that

Y j,k (ω) -Y n j,k (ω) ≤ Y j,k (ω) -Y N j,k (ω) + Y N j,k (ω) -Y N,n j,k (ω) + Y N,n j,k (ω) -Y n j,k (ω) ≤ 4(1 -H -δ) -1 c 1 (j)C * * δ (ω)N -(1-H-δ) + Y N j,k (ω) -Y N,n j,k (
ω) . Thus, one can derive from (2.19) that, for all positive integers N and for each ω

∈ Ω * ∩ Ω * * , lim sup n→+∞ Y j,k (ω) -Y n j,k (ω) ≤ 4(1 -H -δ) -1 c 1 (j)C * * δ (ω)N -(1-H-δ) . ( 2
Y j,k (ω) -Y n j,k (ω) = 0, for all ω ∈ Ω * ∩ Ω * * ,
which shows that Y n j,k converges almost surely to Y j,k , when n tends to +∞.

Second step: we show that for each fixed γ ∈ (0, α), when n goes to +∞, the SαS random variable Y n j,k (see (2.15)) converges in L γ (Ω) (Ω being the underlying probability space) to the SαS random variable Yj,k defined as

Yj,k := Re R ψ j,k (ξ) |ξ| H+1/α d M α (ξ) .
(2.24)

Observe that, since ψ j,k is a continuous function on R (see (2.10) and (2.4)) belonging to

L 1 (R) ∩ L 2 (R) whose Fourier transform ψ j,k belongs to L 1 (R) ∩ L 2 (R)
R n j,k = 1 2π Re R |ξ| -H-1/α |m|≤4 n d n,m+1
dn,m e idn,mξ -e itξ ψ j,k (t) dt

+ {t / ∈[-2 n ,2 n +2 -n )} e itξ -1 ψ j,k (t) dt d M α (ξ) . ( 2
σ(R n j,k ) α ≤ 2 α U n j,k + V n j,k , (2.30)
where

U n j,k := R |ξ| -αH-1 |m|≤4 n d n,m+1
dn,m e i(dn,m-t)ξ -1 ψ j,k (t) dt α dξ (2.31) and

V n j,k := R |ξ| -αH-1 {|t|≥2 n } e itξ -1 ψ j,k (t) dt α dξ.
(2.32)

For deriving appropriate upper bounds for U n j,k and V n j,k , we will make use of the inequality e iθ -1 β ≤ min |θ| β , 2 β , for all θ ∈ R and β ∈ (0, +∞).

(2.33)

Let B 0 := [-1, 1] and B 1 := R \ B 0 . For l ∈ {0, 1}, one sets U n,l j,k := B l |ξ| -αH-1 |m|≤4 n d n,m+1 dn,m e i(dn,m-t)ξ -1 ψ j,k (t) dt α dξ (2.34) and V n,l j,k := B l |ξ| -αH-1 {|t|≥2 n } e itξ -1 ψ j,k (t) dt α dξ. (2.35)
In view of (2.31) and (2.32), one clearly has that

U n j,k = U n,0 j,k + U n,1 j,k and V n j,k = V n,0 j,k + V n,1 j,k . (2.36)
Let us now bound U n,0 j,k and V n,0 j,k . To this end, one introduces the positive finite constant c 2 defined as

c 2 := B 0 |ξ| α2 -1 (1-H)-1 dξ. ( 2
U n,0 j,k ≤ c 2 2 2 -1 (1-H) R ψ j,k (t) dt α 2 -nα2 -1 (1+H) (2.39) and V n,0 j,k ≤ c 2 2 2 -1 (1-H) {|t|≥2 n } t 2 -1 (1+H) ψ j,k (t) dt α .
(2.40)

Moreover, one can derive from (2.14) that

c 3 (j, k) := c 2 2 2 -1 (1-H) R ψ j,k (t) dt α < +∞ (2.41)
and that

{|t|≥2 n } t 2 -1 (1+H) ψ j,k (t) dt (2.42) ≤ 2c 1 (j) +∞ 2 n 1 + t 2 -1 (H-3) dt ≤ 4(1 -H) -1 c 1 (j)2 -n2 -1 (1-H) . Setting c 4 (j) := c 2 2 2 -1 (1-H) 4(1 -H) -1 c 1 (j) α
, and putting together (2.39), (2.40), (2.41) and ( 2.42), one obtains, for all n ∈ N, that

U n,0 j,k ≤ c 3 (j, k)2 -nα2 -1 (1+H) and V n,0 j,k ≤ c 4 (j)2 -nα2 -1 (1-H) .
(2.43)

Let us now bound U n,1 j,k and V n,1 j,k . To this end, one introduces the two positive finite constants c 5 and c 6 defined as

c 5 := B 1 |ξ| -α2 -1 H-1 dξ and c 6 := B 1 |ξ| -αH-1 dξ.
(2.44)

One can derive from (2.44), (2.33) with β = 2 -1 H (for U n,1 j,k ) and β = 1 (for V n,1 j,k ), (2.34) and (2.35) with l = 1, and the inequality (2.38), that one has, for every n ∈ N, Conclusion: combining the main result obtained in the first step with the main one derived in the second step, it turns out that the almost sure equality (2.9) is valid.

U n,1 j,k ≤ c 5 2 (1-2 -1 H) R ψ j,k (t) dt α 2 -nα2 -1 H = c 5 (j, k)2 -nα2 -1 H (2.45) and V n,1 j,k ≤ c 6 2 {|t|≥2 n } ψ j,k (t) dt α ≤ c 6 (j)2 -nα , ( 2

Proofs of the two main results

Lemma 3.1 Under the sole condition (1.6) on the sequence (m j ) j∈N , one has

lim j→+∞ V m j j,γ E(V m j j,γ ) = 1, for all fixed γ ∈ (0, 4 -1 α), (3.1)
where the convergences holds almost surely.

Lemma 3.2 When the sequence (m j ) j∈N satisfies the two conditions (1.6) and (1.9), or when it satisfies the two conditions (1.6) and (1.11), one has

lim j→+∞ V m j j,log 2 E(V m j j,log 2 ) = 1, (3.2)
where the convergences holds almost surely.

For proving the two Lemmas 3.1 and 3.2, one needs several preliminary results.

Remark 3.3 For any fixed α ∈ (0, 2], let Z be an arbitrary real-valued SαS random variable with scale parameter σ(Z) ≥ 0. Then, for all γ ∈ (0, α), one has

E |Z| γ = c 0 (γ, α)σ(Z) γ , (3.3)
where the positive finite constant c 0 (γ, α) := E |W | γ , which only depends on γ and α, is the absolute moment of order γ of an arbitrary real-valued SαS random variable W with scale parameter equals to 1. Moreover, when σ(Z) > 0, letting c 1 (α) be the finite constant only depending α defined as c 1 (α) := E log 2 |W | , one has that

E log 2 (|Z|)) = log 2 (σ(Z)) + c 1 (α), (3.4)
where log 2 is the binary logarithm (that is log 2 (2) = 1).

Proof Using the fact that Φ Z and Φ W , the characteristic functions of the random variables Z and W , satisfy (see e.g. [START_REF] Samorodnitsky | Stable Non-Gaussian Processes: Stochastic Models with Infinite Variance[END_REF]), for all λ ∈ R,

Φ Z (λ) = exp -σ(Z) α |λ| α and Φ W (λ) = exp -|λ| α ,
it turns out that the random variable Z has the same distribution as the random variable σ(Z)W . The two equalities (3.3) and (3.4) are straightforward consequences of the latter fact.

Lemma 3.4 Let ν(α) the positive finite constant defined as

ν(α) := 2 R ψ(η) α dη 1/α = 2 2 -1 -2 -1 ψ(η) α dη 1/α , (3.5)
where the last equality results from the fact that

[-2 -1 , 2 -1 ] is the support of ψ (see (2.4)).
The scale parameter of Y j,k satisfies, for all (j, k)

∈ N 2 , ν(α)2 jH k + 2 -1 -(H+1/α) ≤ σ(Y j,k ) ≤ ν(α)2 jH k -2 -1 -(H+1/α) . (3.6)
Proof One knows from (2.9), (2.10), Remark 1.1 and (2.4) that, for all (j, k) ∈ N 2 , one has

σ(Y j,k ) α = 2 -j (-k+2 -1 ) 2 -j (-k-2 -1 ) ψ(2 j ξ + k) α |ξ| αH+1 dξ + 2 -j (k+2 -1 ) 2 -j (k-2 -1 ) ψ(2 j ξ -k) α |ξ| αH+1 dξ
Then, using the fact that ψ is an even function and the change of variable η = 2 j ξ -k i.e ξ = 2 -j (η + k), one gets that

σ(Y j,k ) α = 2 2 -j (k+2 -1 ) 2 -j (k-2 -1 ) ψ(2 j ξ -k) α (ξ) αH+1 dξ = 2 jαH+1 2 -1 -2 -1 ψ(η) α (η + k) αH+1 dη,
which implies that

2 jαH+1 k+2 -1 -(αH+1) 2 -1 -2 -1 ψ(η) α dη ≤ σ(Y j,k ) α ≤ 2 jαH+1 k-2 -1 -(αH+1) 2 -1 -2 -1 ψ(η) α dη.
Thus, in view of (3.5), it turns out that (3.6) is valid.

Lemma 3.5 For each fixed γ ∈ (0, 4 -1 α), the two positive finite constants µ (γ, α) and µ (γ, α) are defined as

µ (γ, α) := 1 -(3/5) 1/4 ν(α) γ c 0 (γ, α) and µ (γ, α) := 8ν(α) γ c 0 (γ, α), (3.7)
where ν(α) and c 0 (γ, α) are the same positive constants as in (3.5) and (3.3). Then, one has, for all j ∈ N,

µ (γ, α) 2 jγH m 1-γ(H+1/α) j ≤ E V m j j,γ ≤ µ (γ, α) 2 jγH m 1-γ(H+1/α) j . (3.8)
Notice that, the condition that γ ∈ (0, 4 -1 α) implies that

1 > 1 -γ(H + 1/α) > 4 -1 . (3.9)
Remark 3.6 Let γ ∈ (0, 4 -1 α) be arbitrary and fixed. Two straightforward consequences of Lemma 3.5 are the following.

(i) When the conditions (1.6) and (1.9) hold, one has

lim j→+∞ log 2 E(V m j j,γ ) γ j = H. (3.10)
(ii) When the conditions (1.6) and (1.11) hold, one has

lim j→+∞ log 2 E(V m j j,γ ) γ j = 1 γ - 1 α . (3.11)
Proof of Lemma 3.5 First observe that, one can derive from (1.7), (3.3) and (3.6) that, for all j ∈ N,

ν(α) γ c 0 (γ, α)2 jγH m j k=1 k+2 -1 -γ(H+1/α) ≤ E V m j j,γ ≤ ν(α) γ c 0 (γ, α)2 jγH m j k=1 k-2 -1 -γ(H+1/α) .
(3.12) Moreover, using (3.9) and the inequality m j ≥ 2, one has that

m j k=1 k + 2 -1 -γ(H+1/α) ≥ m j k=1 k+1 k x + 2 -1 -γ(H+1/α) dx = m j +1 1 x + 2 -1 -γ(H+1/α) dx ≥ m j + 3 2 1-γ(H+1/α) - 3 2
1-γ(H+1/α) ≥ m j + 3 2 1-γ(H+1/α) - 3 5
1-γ(H+1/α)

m j + 3 2 1-γ(H+1/α) ≥ 1 -(3/5) 1/4 m 1-γ(H+1/α) j (3.13)
and

m j k=1 k -2 -1 -γ(H+1/α) = 2 m j k=1 k k-2 -1 k -2 -1 -γ(H+1/α) dx ≤ 2 m j k=1 k k-2 -1 x -2 -1 -γ(H+1/α) dx = 2 m j k=1 k-2 -1 k-1 x -γ(H+1/α) dx ≤ 2 m j k=1 k k-1
x -γ(H+1/α) dx = 2

m j 0 x -γ(H+1/α) dx ≤ 8 m 1-γ(H+1/α) j . (3.14)
Putting together (3.12), (3.13) and (3.14), one obtains (3.8).

Lemma 3.7

The following two results hold (i) Under the conditions (1.6) and (1.9), one has that

lim j→+∞ 1 j m j E(V m j j,log 2 ) = H. (3.15)
(ii) Under the conditions (1.6) and (1.11), one has that

lim j→+∞ 1 j m j E(V m j j,log 2 ) = 1 α . (3.16)

Proof

First notice that, in view of (3.4) and (1.12), the expectation E(V m j j,log 2

) can be expressed, for all j ∈ N, as

E(V m j j,log 2 ) = c 1 (α)m j + m j k=1 log 2 σ(Y j,k ) .
(3.17) Also notice that (3.6) implies that, for every j ∈ N, From now on, one focuses on the proof of (3.19). Using the inequality log 2 (1 + x) ≤ 4|x|, for all x ∈ [-2 -1 , +∞), one obtains that S p,m j -

(H + 1/α) m j k=1 log 2 (k -2 -1 ) ≤ H(j m j ) + m j log 2 (ν(α)) - m j k=1 log 2 σ(Y j,k ) (3.18) ≤ (H + 1/α) m j k=1 log 2 (k + 2 -1
m j k=2 log 2 (k) ≤ m j k=1 log 2 1 + p2 -1 k ≤ 2 m j k=1 1 k ≤ 2 + 2 log(m j ).
Thus, for proving (3.19), it turns out that it is enough to show that

lim j→+∞ 1 m j log 2 (m j ) m j k=2 log 2 (k) = 1.
(3.21)

One clearly has, for all j ∈ N, We are now in position to prove the Lemmas 3.1 and 3.2 as well as Theorems 1.2 and 1.3 which are the two main results of our article.

m j 1 log 2 (x) dx ≤ m j k=2 k k-1 log 2 (k)dx = m j k=2 log 2 (k) = m j k=2 k+1 k log 2 (k)dx ≤ m j +1 2 log 2 (x) dx.
Proof of Lemma 3.1 First notice that using Markov inequality, for each j ∈ N, one has

P V m j j,γ E(V m j j,γ ) -1 ≥ m -ρ j = P V m j j,γ -E(V m j j,γ ) ≥ m -ρ j E(V m j j,γ ) ≤ m 4ρ j × E V m j j,γ -E(V m j j,γ ) 4 E(V m j j,γ ) 4 , (3.25)
where ρ is a fixed positive constant small enough, which will be chosen more precisely later. Let now provide an appropriate upper bound for the expectation E V

m j j,γ -E(V m j j,γ ) 4 .
Noticing that, for each (j, k) ∈ N 2 , the SαS random variable Y j,k is equal in distribution to σ(Y j,k )W , where W is an arbitrary real-valued SαS random variable with scale parameter equals to 1, one gets, for q ∈ {1, 2}, that

E |Y j,k | γ -E(|Y j,k | γ ) 2q = c q σ(Y j,k ) 2qγ , (3.26)
where the positive finite constant c

q := E |W | γ -E(|W | γ ) 2q
does not depend on (j, k). Then, one can derive from (1.7), (3.26), and the fact that, for any fixed j ∈ N, the centered random variables |Y j,k | γ -E(|Y j,k | γ ), k ∈ N, are independent (see the very fundamental Remark 2.3), and from Lemma 3.4 that

E V m j j,γ -E(V m j j,γ ) 4 = m j k 1 ,..., k 4 =1 E 4 l=1 |Y j,k | γ -E(|Y j,k | γ ) ≤ m j k=1 E |Y j,k | γ -E(|Y j,k | γ ) 4 + 6 m j k=1 E |Y j,k | γ -E(|Y j,k | γ ) 2 2 = c 2 m j k=1 σ(Y j,k ) 4γ + 6c 2 1 m j k=1 σ(Y j,k ) 2γ 2 ≤ c 3 m j k=1 σ(Y j,k ) 2γ 2 ≤ c 4 2 j4γH m j k=1 k -2 -1 -2γ(H+1/α) 2 ≤ c 4 2 j4γH 4 γ(H+1/α) + m j 1 x -2 -1 -2γ(H+1/α) 2 (3.27) ≤ c 5 2 j4γH 1 + 1l {1} 2γ(H + 1/α) log m j + m 1-2γ(H+1/α) j 2 ,
where c 3 , c 4 and c 5 are three positive finite constants not depending on j. Next, combining (3.27) with the first inequality in (3.8), one gets, for all j ∈ N,

m 4ρ j × E V m j j,γ -E(V m j j,γ ) 4 E(V m j j,γ ) 4 ≤ c 6 m -4(1-ρ+γ(H+1/α)) j 1 + 1l {1} 2γ(H + 1/α) log m j + m 1-2γ(H+1/α) j 2 (3.28) ≤ c 6 m -2(1-ρ+γ(H+1/α)) j + 1l {1} 2γ(H + 1/α) m -2(1-ρ+γ(H+1/α)) j log m j + m -(1-2ρ) j 2 ≤ 3c 6 m -4(1-ρ+γ(H+1/α)) j + 1l {1} 2γ(H + 1/α) m -4(1-ρ+γ(H+1/α)) j log 2 m j + m -2(1-2ρ) j .
Next, notice that, in view of (3.9), the positive constant ρ can be chosen small enough so that one has 4 1 -ρ -γ(H + 1/α) > 1 and 2(1 -2ρ) > 1 .

(3.29)

Then, it follows from (3.28), (3.29), (1.6) and (3.25) that

+∞ j=1 P V m j j,γ E(V m j j,γ ) -1 ≥ m -ρ j < +∞.
Therefore, (3.1) results from Borel-Cantelli Lemma.

Proof of Lemma 3.2 First notice that similarly (3.25), it can be shown that, for each j ∈ N, Noticing that, for each (j, k) ∈ N 2 , the SαS random variable Y j,k is equal in distribution to σ(Y j,k )W , where W is an arbitrary real-valued SαS random variable with scale parameter equals to 1, one gets, for q ∈ {1, 2}, that

P V m j j,log 2 E(V m j j,log 2 ) -1 ≥ m -1/8 j ≤ m 1/2 j × E V m j j,log 2 -E(V m j j,log 2 ) 4 E(V m j j,log
E log 2 |Y j,k | -E(log 2 |Y j,k |) 2q = c 0,q , (3.31)
where the positive finite constant c 0,q := E log

Remark 1. 1 (

 1 i) The stable stochastic integral R • d M α is a linear map on the Lebesgue space L α (R) such that, for any deterministic function g ∈ L α (R), the real part Re R N g(ξ) d M α (ξ) is a real-valued Symmetric α-Stable (SαS) random variable with a scale parameter satisfying σ Re R g(ξ) d M α (ξ)

m j 1 log 2 2 log 2

 1222 x → x log(x)-x / log(2) is a primitive of x → log 2 (x), elementary calculations show that (x) dx = m j log 2 (m j ) + 1 -(x) dx = (m j + 1) log 2 (m j + 1(3.23) and (3.24) with (3.22), one obtains (3.21).

2

 2

) 4 .j j,log 2 -

 42 (3.30) Let now provide an appropriate upper bound for the expectation E V m

  It follows from Remark 1.1, (2.25) and the inequality |a + b| α ≤ 2 α |a| α + |b| α , for all comlex numbers a and b, that, for every n ∈ N,

	Thus, in view of (3.3) in Remark 3.3 in the next section, for proving that lim n→+∞	E |R n j,k | γ = 0,
	for each fixed γ ∈ (0, α), it is enough to show that the scale parameter of the SαS random
	variable R n j,k satisfies	lim n→+∞	σ(R n j,k ) = 0.	(2.29)
				.28)

  Notice

	that the last
	inequality in (2.46) is obtained similarly to (2.42),
	Finally, putting together (2.30), (2.36), (2.43), (2.45) and (2.46), it follows that (2.29) is
	satisfied.

.46) where c 5 (j, k) and c 6 (j) are two finite constants not depending on n.

  ).

	Since	lim j→+∞	m j j m j	= 0,
	combining (3.17) and (3.18), it turns out that for proving the lemma it is enough to show that,
	for p ∈ {-1, 1}, one has
	lim j→+∞	S p,m j m j log 2 (m j )	= 1, for p ∈ {-1, 1},	(3.19)
	where	m j
		S p,m j :=	log 2 (k + p2 -1 ).	(3.20)
		k=1
	Observe that when the conditions imposed to the sequence (m j ) j∈N in Part (i) of the lemma
	hold, then (3.19) implies that
		lim j→+∞	S p,m j j m j	= 0.
	Also, observe that when the conditions imposed to the sequence (m j ) j∈N in Part (ii) of the
	lemma hold, then (3.19) entails that
		lim j→+∞	S p,2 j j m j	= 1.

  2 |W | -E(log 2 |W |)2q does not depend on (j, k). Then, one can derive from (1.12),(3.31), and the fact that, for any fixed j ∈ N, the centered random variables log2 |Y j,k | -E(log 2 |Y j,k | γ ), k ∈ N, are independent (see the very fundamental Remark 2.3), that |Y j,k | -E(log 2 |Y j,k |) |Y j,k | -E(log 2 |Y j,k |) |Y j,k | -E(log 2 |Y j,k |)Next, combining (3.32) with Lemma 3.7, and (3.30), under the assumptions of Lemma 3.2, it follows that (almost sure convergence), which means that H j,γ , defined in (1.8), is a strongly consistent estimator of the Hurst parameter H. On another hand, when the conditions (1.6) and (1.11) are satisfied, it follows from (3.33), (3.1) and (3.11) that Proof of Theorem 1.3 One clearly has, for all j large enough, (almost sure convergence), which means that H j,log 2 , defined in (1.13), is a strongly consistent estimator of the Hurst parameter H. On another hand, when the conditions (1.6) and (1.11) are satisfied, it follows from (3.34), (3.2) and (3.16) that which means that α j,log 2 , defined in(1.14), is a strongly consistent estimator of the stability parameter α.

			V j m j m j j,log 2	=	E(V j,log 2 m j j m j	)	×	m j j,log 2 E(V V m j j,log 2	)	.	(3.34)
	On one hand, when the conditions (1.6) and (1.9) hold, one can derive from (3.34), (3.2) and
	E V log 2 ≤ m j j,log 2 -E(V m j j,log 2 ) 4 = m j k 1 ,..., k 4 =1 E 4 l=1 E log 2 4 + 6 m j m j (3.15) that lim j→+∞ m j V j,log 2 j m j = H , lim j→+∞ V m j j,log 2 j m j = 1 α , (almost sure convergence), E log 2 2	2
	k=1									k=1
	= c 0,2 m j + 6c 2 0,1 m 2 j .								(3.32)
		+∞ j=1	P	m j j,log 2 V E(V m j j,log 2	)	-1 ≥ m	-1/8 j	< +∞.
		m j j,γ ) γ j	=	1 γ j	log 2	V E(V m j j,γ m j j,γ )	+	log 2 (E(V j,γ )) m j γ j	.	(3.33)
	On one hand, when the conditions (1.6) and (1.9) hold, one can derive from (3.33), (3.1) and
	(3.10) that	lim j→+∞ = H , lim log 2 (V m j j,γ ) γ j j→+∞ log 2 (V m j j,γ ) γ j = 1 γ -1 α	, (almost sure convergence),
	which means that α j,γ , defined in (1.10), is a strongly consistent estimator of the stability
	parameter α.

Therefore, (3.2) results from Borel-Cantelli Lemma.

Proof of Theorem 1.2 One clearly has, for all γ ∈ (0, 4 -1 α) and for every j ∈ N, log 2 (V

Acknowledgements

The author acknowledges partial support from the Labex CEMPI (ANR-11-LABX-0007-01), the GDR 3475 (Analyse Multifractale et Autosimilarité), and the Australian Research Council's Discovery Projects funding scheme (project number DP220101680).


Keywords: Asymptotic analysis, Finite Element Method, Laplace equations

1 Introduction

In the domains with a rough thin layer, numerical difficulties appear due to the complex geometry of the rough layer when computing the steady-state potentials. We present here how these difficulties may be avoided by replacing this rough layer by appropriate transmission conditions. Particularly, we show that considering only the mean effect of the roughness is not sufficient to obtain the potential with a good accuracy.

Statement of the problem

Ω O 0 ε O 1 O m ε Γ Γ ε Figure 1: Geometry of the problem.
Let Ω be a smooth bounded domain of R 2 with connected boundary ∂Ω. For ε > 0, we split Ω into three subdomains: O1 , O m ε and O 0 ε . O 1 is a smooth domain strictly embedded in Ω (see Fig. 1). We denote by Γ its connected boundary. The domain O m ε is a thin oscillating layer surrounding O 1 . We denote by Γ ε the oscillating boundary of

O m ε : Γ ε = ∂O m ε \ Γ. The domain O 0 ε is defined by: O 0 ε = Ω \ (O 1 ∪ O m ε ).
We also denote by O 0 = Ω \ O 1 . Two piecewise-constant conductivities on the domain Ω have to be defined:

σ(z) =      σ 1 , if z ∈ O 1 , σ m , if z ∈ O m ε , σ 0 , if z ∈ O 0 ε . σ(z) = σ 1 , if z ∈ O 1 , σ 0 , if z ∈ Ω \ O 1 .
where σ 1 , σ m and σ 0 are given positive constants 1 .

Let u ε and u 0 be defined by:

∇. (σ∇u ε ) = 0, in Ω, u ε | ∂Ω = g, , ∇. σ∇u 0 = 0, in Ω, u 0 | ∂Ω = g, (1)
where g is a sufficiently smooth boundary data. We present how to define the potential u 1 such that u ε is approached by

u ε = u 0 + εu 1 + o(ε 3/2 ) for ε tending to zero 2 .

Heuristics of the derivation of the conditions

Suppose Γ is a smooth closed curve of R2 of length 1 and parameterize it by the curvilinear coordinate Γ = {Ψ(θ), θ ∈ [0, 1]} . Let n be the (outward) normal to ∂O 1 . Γ ε is described by

Γ ε = {Ψ(θ) + εf (θ/ε)n(θ), θ ∈ [0, 1]},
where f is a smooth 1-periodic and positive function, which describes the roughness of the layer.

Boundary layer corrector in the infinite strip

The key-point of the derivation of the equivalent transmission conditions consists in taking advantage of the periodicity of the roughness. This is performed by unfolding and upscaling the rough thin layer into the infinite strip R × [0, 1]. Define the closed curves C 1 and C 0 , which are trigonometrically oriented by

C 0 = {0} × [0, 1], C 1 = {(f (y), y), ∀y ∈ [0, 1]} .
The outward normals to C 0 and C 1 equal

n C0 = 1 0 , n C1 = 1 1 + (f ′ (y)) 2 1 -f ′ (y) . (2)
According to [START_REF] Ciuperca | Approximate transmission conditions through a rough thin layer. The case of the periodic roughness[END_REF] there exists a unique couple (A 0 , a 0 ) where A 0 is a continuous vector field and a 0 is constant such that

A 0 is 1-periodic in y, ∆A 0 = 0, in R × [0, 1], (3a)
σ 0 ∂ n A 0 | C + 1 -σ m ∂ n A 0 | C - 1 = (σ m -σ 0 )n C1 , (3b)
σ m ∂ n A 0 | C + 0 -σ 1 ∂ n A 0 | C - 0 = -(σ m -σ 0 )n C0 , (3c)
A 0 → x→-∞ 0, A 0 -a 0 → x→+∞ 0, (3d)
where the convergences at infinity are exponential. We emphasize that a 0 is not imposed but is a floating potential.

Approximate transmission conditions

Our transmission conditions are then obtained with the help of the constant vectors D 1 and D 2 defined by:

D 1 = (σ 0 -σ m ) 1 0 f (y) dy n C0 + 1 0
A 0 (f (y), y) dy

+ (σ m -σ 1 ) 1 0 A 0 (0, y) dy -σ 0 a 0 , D 2 = (σ m -σ 0 ) 1 0 A 0 (f (y), y)f ′ (y) dy - 1 0 f (y) dy 0 1 .
The potential u 1 is then defined by3 :

               ∆u 1 = 0, in O 0 ∪ O 1 , u 1 | ∂Ω = 0, σ∂ n u 1 Γ = -κD 1 • ∂ n u 0 | Γ+ ∂ t u 0 | Γ + + D 2 • ∂ t ∂ n u 0 | Γ+ ∂ t u 0 | Γ + , u 1 Γ = a 0 • ∂ n u 0 | Γ+ ∂ t u 0 | Γ + ,
where ∂ t and ∂ n denote the tangential and the normal derivatives along Γ and κ is the curvature of Γ. We emphasize that our conditions are different than if we would only consider the mean effect of the roughness. In this case, denoting by f the mean of f , the conditions would be (see [START_REF] Ciuperca | Approximate transmission conditions through a rough thin layer. The case of the periodic roughness[END_REF][START_REF] Poignard | Approximate conditions replacing thin layer[END_REF][START_REF] Poignard | Approximate transmission conditions through a weakly oscillating thin layer[END_REF]):

σ∂ n u 1 Γ = (σ 0 -σ m )f ∂ 2 t u 0 | Γ , u 1 Γ = σ 0 -σ m σ m f ∂ n u 0 | Γ + .

Numerical simulations

In order to verify the convergence rate stated in Section 1, we consider a problem where the geometry and the boundary conditions are ε-periodic. The computational domain Ω is delimited by the circles of radius 2 and of radius 0.2 centered in 0, while O 1 is the intersection of Ω with the concentric disk of radius 1. The rough layer is then described by f (y) = 1+1/2 sin(y). One period of the domain is shown Fig 2(a). The Dirichlet boundary data is identically 1 on the outer circle and 0 on the inner circle. The mesh generator Gmsh [START_REF] Geuzaine | Gmsh mesh generator[END_REF] and the finite element library Getfem++ [START_REF] Renard | avenue de l'Europe -38334 Montbonnot Saint-Ismier Centre de recherche INRIA Lille -Nord Europe : Parc Scientifique de la Haute Borne -40, avenue Halley -59650 Villeneuve d'Ascq Centre de recherche INRIA Nancy -Grand Est : LORIA, Technopôle de Nancy-Brabois -Campus scientifique 615, rue du Jardin Botanique -BP 101 -54602 Villers-lès[END_REF] enables us to compute the four potentials u ε , u 0 , u 1 and u 1 . The rough thin layer is supposed slighty insulating. The conductivities σ 0 , σ 1 and σ m respectively equal to 3, 1 and 0.1. The computed coefficients4 issued from Problem (3) are given in Table 1. The numerical convergence rates for the H 1 -norm in O 1 of the three following errors u ε -u 0 , u ε -u 0 -εu 1 and u ε -u 0 -ε u 1 as ε goes to zero are given Figure 2. As predicted by the theory, the rates are close to 1 for the order 0 and for the order 1 with the mean effect, whereas it is close to 2 for the "real" order 1 equal to u εu 0εu 1 .

Figure 3: H 1 -Error in the cytoplasm vs ε for three approximate solutions.

1 4 2 2 . 1 2 . 2

 142122 Statement of the problem . . . . . . . . . . . . . . . . . . . . . . Heuristics of the derivation of the conditions 5 Boundary layer corrector in the infinite strip . . . . . . . . . . . 5 Approximate transmission conditions . . . . . . . . . . . . . . . . 6

Figure 2 :

 2 Figure 2: Representation of one period of the domain and the corresponding errors with approximate solutions u 0 and u 0 + εu 1 . ε = 2π/30. Do not consider the error in the rough layer because a proper reconstruction of the solution in it is not currently implemented.

Table 1 :

 1 Coefficients issued from the solution to problem (3). 3 significant digits are kept.

RR n°6935

The same following results hold if σ 0 , σ 1 , and σm are given complex numbers with imaginary parts (and respectively real parts) with the same sign.INRIA

The notation o(ε

3/2 ) means that ' ' u ε -(u 0 + εu 1 ) ' ' goes to zero faster than ε 3/2 as ε goes to zero. We refer to Theorem 1.1 of[START_REF] Ciuperca | Approximate transmission conditions through a rough thin layer. The case of the periodic roughness[END_REF] for a precise description of the involved norms and the accuracy of the convergence.RR n°6935

We denote by [w] Γ the jump of a function w on Γ.

The convergences at the infinity in Problem (3) are exponential hence we just have to compute problem (3) for |x| ≤ M , with M large enough to obtain a 0 with a good accuracy.


Connections, disconnections and reconnections in socio-hydrosystems:

Towards a better integration of the vertical dimension?

Anne-Lise Boyer, Postdoctoral Fellow, CNRS Labex DRIIHM, IRL iGlobes at the University of Arizona -2023 AAG Annual Meeting, Denver, CO

Abstract:

The relationship between societies and their hydrosystems can be understood as a series of connections and disconnections that occur in all domains (e.g., hydrological, geomorphological, ecological, etc.) and dimensions (e.g, longitudinal, vertical, lateral). The use of watercourses creates material and symbolic connections which leads to also consider the sociological, political, aesthetic, and spiritual dimensions of the socio-hydrosystem.

Technical progress and the multiplication of infrastructures have contributed to the invisibilization of water resources. The issue of visibility is to be understood as lying at the intersection of relations of perceptions and relations of power which shapes decision-making.

As a result, recent water and environmental policies emphasize the importance of connectivity between the various compartments of the socio-hydrosystems, and focus especially on reconnections at all levels. Successfully restoring the socioecological connectivity is tied to (in)visibility dynamics that must be taken into account. The objective is to analyze the impacts of the (in)visibility of challenges such as emerging pollutants and hydrological changes in the context of climate change.

A network of Human -Environment Observatories:

To respond to the complexity of the socioecosystemic crises that increasingly affect an anthropized and globalized planet, the Human-Environment Observatories (OHM) are a research infrastructure created by the Institute for Ecology and Environment of the CNRS (France). The first one was created in 2007. Today there are 13 OHM, spread across metropolitan France (6), overseas (2) and abroad (5). All OHMs are built upon:

• a founding event triggering of human origin that disrupts the entire functioning of an ecosystem;

• the convergence of Earth, Life, Human and Social Sciences in an interdisciplinary approach to address the complexity of the socio-ecosystems under study;

• the development of scientific research to understand the mechanisms at play and to better inform political and economic decision-makers and citizens in their choices for the resolution of ongoing environmental crises. Reconnection : recovery of lost linkages over time, e.g. restoration projects and the responses they trigger in the system.

Attachment and stronger emotional reactions, e.g. mourning and grieving in the face of environmental changes

The social connectivity of invisible subsurface flows: towards a better integration of the vertical dimension in socio-hydrosystems studies o The importance of vertical connectivity for the sustainability of groundwater-dependent ecosystems -two-way exchanges between the surface and water located in the underground, e.g., surface water infiltration into the alluvial aquifer or exfiltration of phreatic water from the hillslopes aquifer to supply the stream and other floodplain water bodies.

-Lagoons but also rivers and their adjacent areas (e.g., alluvial plain) are considered groundwater-dependent ecosystems : about 30% of streamflow comes from groundwater o Groundwater, shallow groundwater, and the hyporheic zone -The existing hydrologic exchanges between the different water bodies and the vertical linkages are not always fully appreciated: difficulties to assess groundwater volumes, recharge rates and sources + low recognition of the linkages between groundwater and many surface water ecosystems.

-Groundwaters are heterogenous and ubiquitous because of their extensive distribution , "both on top of and within the earth with a varying degrees of connectivity that depends on geology [e.g., shallow aquifer, deep aquifer, fossil aquifer, etc.] and human action [e.g., withdrawals]" (Walsh, 2022)

-The hyporheic zone: the region of sediment and porous space beneath and alongside a stream bed, where shallow groundwater and surface water mix -Groundwater-dependent socio-hydrosystems face a particularly volumetric challenge: groundwater depletion and contamination, and their consequences on the surface (e.g., subsidence, drinking water supply availability, irrigated agriculture and food production, social conflicts, and land use transformations, etc.) -In the Anthropocene, as interests in conserving and restoring the subsurface flows and their interactions with the surface strengthen, verticality comes to matter in the process of making sense of this continuum between surface and subsurface by articulating with existing geomorphic, hydrologic, and ecologic analyses but also with current modes of governing and policing resources and people.

Different factors for a social recognition of the hyporheic zone and the importance of shallow groundwater:

-the presence/ absence of a healthy riparian forest: wildlife refuges, important recreation areas, or just shade in desert cities parched by the urban heat island effect.

-dry and wet participatory surveys, efforts are made to recover the loss of socio-cultural connectivity to "the rhythm of rivers" (Wantzen, 2022).

-materializing waterscapes with a focus on the interface between groundwater and surface water to disrupt the loose regulation of extractive exploitation of streams and groundwater (Wölfle-Hazard, 2022).

-integration of unconventional research methods and traditional knowledge, environmental justice issues and research related to health and livelihood.

  HYDECO: an interdisciplinary research project 2022 -2024 Monitoring the evolutions of 8 socio-hydrosystems in the face of climate change DISCONNECTIONS RECONNECTIONS Materials and Methods: Lake Arbu monitoring, Pyrenees, November 2022 Hydroelectric dam on the Rhine River Salinity control gate, Etang de l'Or, FR Rice fields, Aveiro, Portugal Waterfront in Lyon, Rhône River, FR o Hydrological connectivity : "the water-mediated fluxes of material, energy, organisms within and amongst components (Kondolf et al., 2006) o In 5 dimensions: longitudinal, lateral, and vertical dimensions -+ over time -+ social connectivity : communication and movement of people, goods, ideas, and cultures along and across rivers Disconnections = lost linkages Kondolf, G. M., et al. "Process-based ecological river restoration: visualizing three-dimensional connectivity and dynamic vectors to recover lost linkages." Ecology and society 11.2 (2006).

o

  Reclaiming the hyporheic zone: the social connectivity of subsurface flows o Verticality, the subsurface, and why it matters in the Anthropocene Cienega Creek, Empire Ranch, Las Cienegas Conservation Area (Southern Arizona, April 2022) Wellfield of the Crépieux Charmy water production facility (Lyon, FR, September 2019)-The "vertical turn" and "volumetric thinking" in social sciences: the close articulation of horizontality and verticality + integration of variables of depth, height, and distance.


Amine Ammar
email: amine.ammar@ujf-grenoble.fr

Francisco Chinesta
email: francisco.chinesta@paris.ensam.fr

A Particle Strategy for Solving the Fokker-Planck Equation Modelling the Fiber Orientation Distribution in Steady Recirculating Flows Involving Short Fiber Suspensions

In this work we have analyzed the application of a deterministic approximation of the diffusion term in the Fokker-Planck equation using smooth particles for computing its steady solution in a steady recirculating flow. The main idea of this approach lies in the introduction of the Fokker-Planck diffusion term into the advection one, which allows to proceed in a Lagrangian deterministic manner without a mesh support requirement.

Introduction

As indicated by Chaubal et al. [START_REF] Chaubal | Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part 1: Method[END_REF], "complex fluid" is the term commonly used to describe a wide class of liquid-like materials, in which the relaxation time towards the equilibrium state occurs sufficiently slowly that significant changes in the microstructural configuration, and thus in the macroscopic properties, can be induced by the flow. Viscoelastic fluids or short fiber suspensions may be considered as examples of complex fluids. The Fokker-Planck formalism is a commonly used description of kinetic theory problems, for describing the evolution of the configuration distribution function. This function represents the probability of finding the microstructural element in a particular configuration.

In the case of a short fiber suspension, the configuration distribution function (also known as orientation distribution function) gives the probability of finding the fiber in a given direction. Obviously, this function depends on the physical coordinates (space and time) as well as on the configuration coordinates, that taking into account the rigid character of the fibers, are defined on the surface of the unit sphere. Thus, we can write Ψ (x,t,p), where x defines the position of the fiber center of mass, t the time and p the unit vector defining the fiber orientation. The evolution of the distribution function is given by the Fokker-Planck equation

dΨ dt = - ∂ ∂p (Ψ ṗ)+ ∂ ∂p D r ∂Ψ ∂p (1.1)
where d/dt represents the material derivative, D r is a diffusion coefficient and ṗ is the fiber rotation velocity. When the fibers are assumed with an ellipsoidal shape and the suspension is dilute enough, the rotation velocity can be obtained from the Jeffery's equation

ṗ = Ω p + kDp -k(p T D p) p (1.2)
where Ω and D are the vorticity and the strain rate tensors respectively, associated with the fluid flow undisturbed by the presence of the fiber, and k is a scalar which depends on the fiber aspect ratio λ (ratio between the fiber length and the fiber diameter)

k = λ 2 -1 λ 2 +1 (1.3)
Although this work focuses in the flow of short fiber suspensions, the numerical procedures here developed could be applied successfully to other fluids whose microstructure is described by similar kinetic theory models. This is the case for example of some viscoelastic fluids. In the FENE model [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF], the probability distribution function depends on the physical coordinates and on the conformation coordinates that in this case are defined by the end-to-end molecule vector q. The equation governing its evolution results:

dΨ dt = - ∂ ∂q (Ψ q)+ 1 2 ∂ ∂q ∂Ψ ∂q (1.4)
where the evolution of q can be evaluated from

q = Gradv q - 1 2 F (q) (1.5)
with the spring force F given by

F = 1 1 - q 2 b 2 q q (1.6)
being b the maximum molecule length. The FENE-P model is obtained by averaging the square norm of the molecule length in the expression of the spring force [START_REF] Bird | Polymer solution rheology based on a finitely extensible bead-spring chain model[END_REF].

In these cases the Fokker-Planck equations remain linear, even if the advection fields depend in a non-linear manner on the configuration coordinates (see for instance Eqs. (1.2) and (1.5)).

Another more complex example consists in the model describing the behavior of Liquid Crystals Polymers (LCP) due to Doi. This model has been treated using a particle technique in [START_REF] Chaubal | Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part 1: Method[END_REF]. In this case the microstructural element is a rigid axis-symmetric rod of infinite aspect ratio. Under some simplifying assumptions the orientation distribution (also defined on the surface of the unit sphere) evolves according to

dΨ dt = - ∂ ∂q (Ψ ṗ)+ ∂ ∂p D r ∂Ψ ∂p + Ψ ∂U ∂p (1.7)
where the advection field ṗ is given by

ṗ = Gradv p -(p T D p) p (1.8)
As the potential U depends on the second moment of the distribution function p ⊗ p defined by

p ⊗ p = p ⊗ p Ψ (p) dp (1.9)
and the Fokker-Planck equation results in this case non-linear. Many of the experimental and industrial flows show recirculating areas or recirculate themselves. For example, many rheometric devices involve this type of flows, whose steady and recirculating character introduce some additional difficulties in their numerical simulation. Actually, the Fokker-Planck equation which defines an advection problem in physical coordinates, is supposed to have a steady solution in these steady recirculating flows but neither boundary conditions nor initial conditions are known in such flows.

In a former paper [START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF] the discretisation of the advection dominated Fokker-Planck equation, governing the fiber orientation in short fiber suspension flows, was carried out using a particle technique, where the diffusion term was modelled from random motions. It was pointed out that the number of fibers required in this stochastic simulation to describe the fiber distribution increases significantly with the diffusion coefficient D r . Thus, it was argued that for practical applications the use of the particle method in the framework of a stochastic simulation, is restricted to very slight diffusion effects. When the diffusion becomes dominant, continuous approximations using a fixed mesh seem to be suitable, but in this case accurate stabilizations are required for dealing with small diffusion effects, and a lack of accuracy is noticed in the treatment of the advection dominated case.

Chaubal et al. [START_REF] Chaubal | Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part 1: Method[END_REF] propose the use of the SPH (Smooth Particle Hydrodynamics) to solve the dynamics of a liquid crystalline polymer (LCP) using the Doi's model. Thus, the diffusion term is treated in a determinist manner, and high accurate results were obtained due to the meshless and Lagrangian character of the SPH technique considered by the authors.

In this work we examine the application of this kind of techniques for solving the steady Fokker-Planck equation in steady recirculating flows, but we will limit our attention to the 2D linear Fokker-Planck equations.

A Particle Discretisation

In this section we will consider the simplest linear form of the Fokker-Planck equation (1.1), that in the 2D case results

dΨ dt = - ∂ ∂ϕ (Ψ φ)+ ∂ ∂ϕ D r ∂Ψ ∂ϕ (2.1)
where Ψ (x,t,ϕ), φ =φ(x,t,ϕ) and D r is assumed constant. Due to the steady character of the flow kinematics we can remove the temporal variable in the expression of φ. A general approximation of the steady probability distribution function at point x 0 can be written as

Ψ (x 0 ,ϕ) ≈ i=N i=1 α i (x 0 ) F i (ϕ) (2.2)
Due to the linearity and homogeneity of Eq. (2.1) we can compute its solution for each function F i (ϕ) along the pathline related to the point x 0 . These solutions are denoted by Ψ i (x(x 0 ,t 0 ; t),ϕ), where the notation x(x 0 ,t 0 ; t) refers to the position at time t of a particle located at point x 0 at time t 0 . The trajectory is defined by the equation

x(x 0 ,t 0 ; t)= x 0 + t t0 v(x(x 0 ,t 0 ; t ′ )) dt ′ (2.3)
Thus, the general solution of Eq. (2.1) results

Ψ (x(x 0 ,t 0 ; t),ϕ) ≈ i=N i=1 α i (x 0 ) Ψ i (x(x 0 ,t 0 ; t),ϕ) (2.4)
The particle method lies in taking the Dirac masses in Eq. (2.2). Thus, if we consider the N directions defined by

ϕ i =(i -1) h ϕ , i ∈ [1, ••• ,N], with h ϕ =2π/N , Eq. (2.2) becomes Ψ (x 0 ,ϕ) ≈ i=N i=1 α i (x 0 ) δ(ϕ -ϕ i ) (2.5)
where α i (x 0 ) represents the "mass" of the particle aligned in the ϕ i direction.

The Advection Equation

In the case of neglecting all the diffusion effects, the solution Ψ i (x(x 0 ,t 0 ; t),ϕ) results

Ψ i (x(x 0 ,t 0 ; t),ϕ)=δ(ϕ -ϕ(x 0 ,t 0 ,ϕ i ; t)) (2.6)
where ϕ(x 0 ,t 0 ,ϕ i ; t) denotes the orientation at time t of a fiber located at time t 0 in x 0 and whose orientation was defined by ϕ i . Obviously, the spatial location of a fiber refers to the position of its center of mass, that at time t is given by Eq. (2.3).

In the pure advection case ϕ(x 0 ,t 0 ,ϕ i ; t) becomes

ϕ(x 0 ,t 0 ,ϕ i ; t)= ϕ i + t t0 φ(x(x 0 ,t 0 ; t ′ ),ϕ(x 0 ,t 0 ,ϕ i ; t ′ )) dt ′ (2.7)
Thus, Eqs. (2.3) and (2.7) describe the position and orientation of each fiber along the flow trajectory. Several discretisations of these equations exist, being the simplest one the backward Newton method, that consider the time interval [t 0 ,t] divided into M intervals [t m ,t m+1 ] of length Δt such that M Δt = t -t 0 , and the fibers updating given by

⎧ ⎨ ⎩ x(x 0 ,t 0 ; t m+1 )=x(x 0 ,t 0 ; t m )+v(x(x 0 ,t 0 ; t m )) Δt ϕ(x 0 ,t 0 ,ϕ i ; t m+1 )=ϕ(x 0 ,t 0 ,ϕ i ; t m )+ +φ(x(x 0 ,t 0 ; t m ),ϕ(x 0 ,t 0 ,ϕ i ; t m )) Δt, ∀i ∈ [1, ••• ,N]
(2.8) where

x(x 0 ,t 0 ; t 0 )=x 0 ϕ(x 0 ,t 0 ,ϕ i ; t 0 )=ϕ i , ∀i ∈ [1, ••• ,N]
(2.9)

Steady Recirculating Flows

In the case of a steady recirculating flow the particles come back to the departure position after a time T , that corresponds with the period of the considered trajectory. Thus, we can write

x(x 0 ,t 0 ; T )=x(x 0 ,t 0 ; t 0 )=x 0 (2.10) However, the final orientation will be, in general, different to the initial one, i.e.

ϕ(x 0 ,t 0 ,ϕ i ; T ) = ϕ i , ∀i (2.11)
On the other hand, the steady solution of the probability distribution at point x 0 , Ψ (x 0 ,ϕ), requires its periodicity along the closed streamline

i=N i=1 α i (x 0 ) δ(ϕ -ϕ(x 0 ,t 0 ,ϕ i ; T )) = i=N i=1 α i (x 0 ) δ(ϕ -ϕ i )
(2.12) but this expression, (Eq. (2.12)), cannot be used in its present form because the Dirac masses are concentrated in different angles at t 0 and T .

To use this expression we need to transfer the masses concentrated in the directions ϕ(x 0 ,t 0 ,ϕ i ; T ) to the initial ones ϕ i . In [START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF] this transfer is performed from each angle ϕ(x 0 ,t 0 ,ϕ i ; T ) towards the neighbor directions, being the mass transferred to each neighbor direction proportional to the distance between them. Thus, we can finally write (see [START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF] for more details)

i=N i=1 α i (x 0 ) δ(ϕ -ϕ(x 0 ,t 0 ,ϕ i ; T )) = i=N i=1 j=N j=1 β ij α j (x 0 ) δ(ϕ -ϕ i ) (2.13)
where β ij depends in the considered probability transfer from ϕ(x 0 ,t 0 ,ϕ i ; T ) to ϕ j . If a linear transfer is used, then it results

β ij = 1-|ϕ(x 0 ,t0,ϕi;T )-ϕj | hϕ if |ϕ(x 0 ,t 0 ,ϕ i ; T ) -ϕ j |≤h ϕ 0 if |ϕ(x 0 ,t 0 ,ϕ i ; T ) -ϕ j | >h ϕ
The periodicity condition becomes in this case

i=N i=1 j=N j=1 β ij α j (x 0 ) δ(ϕ -ϕ i )= = i=N i=1 α i (x 0 ) δ(ϕ -ϕ i )= i=N i=1 j=N j=1 δ ij α j (x 0 ) δ(ϕ -ϕ i ) (2.14)
where δ ij is the unit tensor. The previous equation implies

j=N j=1 (β ij -δ ij )α j (x 0 )=0 , ∀i ∈ [1, ••• ,N] (2.15)
that with the normality condition

1= 2π 0 Ψ (x 0 ,ϕ) dϕ = i=N i=1 α i (x 0 ) (2.16)
allows to compute the N coefficients α i (x 0 ), and then, the steady solution of the probability distribution at point x 0

Ψ (x 0 ,ϕ) ≈ i=N i=1 α i (x 0 ) δ(ϕ -ϕ i ) (2.17)

Introducing Diffusion Effects Using a Deterministic Approximation

When diffusion effects are retained in the kinetic model the application of a particle technique becomes delicate as we will describe in this section. The first possibility for introducing of the diffusion effects lies in the use of a stochastic technique [START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF]. Thus, a great amount of particles are considered at point x 0 aligned in each direction ϕ i . Then, each fiber is subjected to three actions during its movement along its closed trajectory: (i) the advection of its center of mass, (ii) the rotation induced by the term ṗ and (iii) the rotation associated with the diffusion that is modelled from a random motion. In spite of the simplicity of its computational implementation and the very accurate results obtained for slight diffusion coefficients, the great number of particles required when the diffusion coefficient increases makes this strategy useless for treating the problems encountered in practical applications.

To circumvent these computational drawbacks, we transform the linear Fokker-Planck advection-diffusion equation in a pure advection problem. For this purpose we define the new advection field ṗ * ṗ * =ṗ -D r ∂Ψ ∂p

Ψ

(2.18) that introduced in the pure advection Fokker-Planck equation

dΨ dt = - ∂ ∂p (Ψ ṗ * ) (2.19)
leads to the standard linear advection-diffusion Fokker-Planck equation (1.1).

In the 2D case developed in the previous section, both equations result in the following purely-advection problem φ * =φ -D r ∂Ψ ∂ϕ

Ψ

(2.20)

dΨ dt = - ∂ ∂ϕ (Ψ φ * ) (2.21)
Now, Eq. (2.21) could be used for computing the different solutions associated with the Dirac's distributions δ(ϕ -ϕ i ), ∀i, but two difficulties appear: (i) the derivative and the ratio of linear combinations of Dirac masses cannot be defined in a proper way, and (ii) due to the distribution function derivative required to compute the modified advection field φ * in Eq. (2.20), a coupling between the different fibers takes place, and consequently, in this case, the evolution of each fiber cannot be computed independently of the other ones.

In order to circumvent the first difficulty we introduce a smoothed approximation of the Dirac mass given by

ζ(x)= e -x 2 √ π (2.22) which verifies the normality condition ∞ -∞ ζ(x) dx =1 (2.23)
To adapt its sharpness we modify this function in the following manner

ζ ǫ (x)= ζ x ǫ ǫ (2.24)
which verifies also the normality condition. Fig. 2.1 depicts this function for three values of the ǫ parameter. In this form, the gradients and ratios can be computed in a proper way. On the other hand to avoid the second limitation two simple possibilities exist: (i) a direct procedure that compute N times, the evolution history of N particles along the closed trajectory, and (ii) an iteration algorithm that computes the evolution of each particle along the closed trajectory independently of the other ones due to the fact that the coupling term is evaluated from the solution at the previous iteration. In the following sections we describe both procedures.

A Direct Procedure.

We write the searched solution at point x 0 in a matrix form

Ψ (x 0 ,ϕ) ≈ i=N i=1 α i (x 0 ) ζ ǫ (ϕ -ϕ i )=α(x 0 ) • ζ ǫ (ϕ) (2.25)
where the unknown coefficients α i (x 0 ), contained in the vector α(x 0 ), must be determined in order to satisfy the periodicity condition imposed by both the steady character of the searched solution and the steady and recirculating character of the flow. We define the following N alpha-vectors:

α k (x 0 )= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ α k 1 (x 0 ) . . . α k k-1 (x 0 ) α k k (x 0 ) α k k+1 (x 0 ) . . . α k N (x 0 ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ν . . . ν μ ν . . . ν ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , ∀k ∈ [1, ••• ,N] (2.26)
where ν and μ can be chosen arbitrarily with the only restriction of verifying the normality condition, i.e. (N -1)ν+μ = 1. We can notice that the uncoupled procedure, applied when diffusion effects are neglected, uses ν = 0 and μ = 1. Introducing this notation, Eq. (2.25) can be rewritten as

Ψ (x 0 ,ϕ) ≈ k=N k=1 β k α k (x 0 ) • ζ = k=N k=1 β k Ψ k (x 0 ,ϕ) (2.27)
Now, we can compute the evolution along a closed trajectory of each solution

Ψ k (x 0 ,ϕ)=α k (x 0 ) • ζ ǫ (ϕ) (2.28)
which in fact requires the tracking of N particles. For this purpose, we need to integrate the evolution of the orientation for each fiber (represented by the subscript i) in each problem (noted by the superscript k). If we consider, for a sake of simplicity, the backward Newton method again, it results

ϕ(x 0 ,t 0 ,ϕ k i ; t m+1 )=ϕ(x 0 ,t 0 ,ϕ k i ; t m )+ φ(x(x 0 ,t 0 ; t m ),ϕ(x 0 ,t 0 ,ϕ k i ; t m )) Δt - -D r ∂ ∂ϕ j=N j=1 α k j (x 0 )ζ ǫ (ϕ -ϕ(x 0 ,t 0 ,ϕ k j ; t m )) ϕ=ϕ(x 0 ,t0,ϕ k i ;tm) j=N j=1 α k j (x 0 )ζ ǫ (ϕ(x 0 ,t 0 ,ϕ k i ; t m ) -ϕ(x 0 ,t 0 ,ϕ k j ; t m )) Δt (2.29) ∀k ∈ [1, ••• ,N], ∀i ∈ [1, ••• ,N], ∀m ∈ [0, ••• ,M -1]
After a complete turn of period T we obtain

ϕ(x 0 ,t 0 ,ϕ k i ; T ) ∀i ∈ [1, ••• ,N], ∀k ∈ [1, ••• ,N] (2.30)
Now the weights associated to each orientation ϕ(x 0 ,t 0 ,ϕ k i ; T ), β k α k i (x 0 ), must be transferred to the angles ϕ i to enforce the periodicity. This operation can be accurately carried out by using the algorithm proposed in [START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF] as previously described.

This algorithm offers the steady solution of the problem from the resolution of the N problems, requiring each one of them, the integration of the orientation evolution of N fibers along the closed trajectory. In this form the algorithm results of order N 2 .

A Fixed Point Iteration Algorithm. In this case we start from the solution obtained assuming D r = 0. This solution can be written as

Ψ 0 (x(x 0 ,t 0 ; t),ϕ) ≈ i=N i=1 α 0 i (x 0 ) ζ ǫ (ϕ -ϕ(x 0 ,t 0 ,ϕ i ; t)), ∀t ∈ [t 0 ,T] (2.31)
Now, the iteration k results

ϕ k (x 0 ,t 0 ,ϕ i ; t m+1 )=ϕ k (x 0 ,t 0 ,ϕ i ; t m )+ φ(x(x 0 ,t 0 ; t m ),ϕ k (x 0 ,t 0 ,ϕ i ; t m ))Δt -D r ∂ ∂ϕ Ψ k-1 (x(x 0 ,t 0 ; t),ϕ) ϕ=ϕ k (x 0 ,t0,ϕi;tm) Ψ k-1 (x(x 0 ,t 0 ; t),ϕ k (x 0 ,t 0 ,ϕ i ; t m )) Δt (2.32) for all i ∈ [1, ••• ,N], and all m ∈ [0, ••• ,M -1],
where in this case the superscript k ≥ 1 refers to the iteration.

This algorithm results of order N but the convergence is not ensured. Some improvements can be introduced, as for example the introduction of the diffusion coefficient step by step.

Numerical Examples

In this section two numerical examples involving steady recirculating flows will be considered. The first one is defined by the following kinematics:

v = u v = -y x 2 + y 2 x x 2 + y 2 (3.1)
The steady solution of the probability distribution is searched at point (0, 1). For this problem the exact solution can be computed due to the solution symmetry as described in [START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF]. Thus, for a diffusion coefficient D r =0 .2 and fibers with an aspect ratio related to k =0 .6, the fixed point algorithm proposed in the section 2.3 with ǫ =0 .5 and N = 72 converges in around 10 iterations to the reference solution as depicted in Fig. 3.2. Now, the direct method described in the section 2.3 is applied. Obviously, the numerical solution accuracy will depends on the μ coefficient (for μ ∈ [3, 5] the extact and the computed solutions are in very good agreement). Figs. 3.3 and 3.4 depicts the fiber distribution related to μ =1 .2/N and μ = 10/N respectively as well as the exact solution (dashed line). The case of μ =0 (where N -1 fibers with equal weights are considered) is depicted in Fig. 3.5. The solution accuracy can be improved by increasing the number of particles N involved in the simulation. A first order convergence is noticed when the number of particles N is increased. In these examples we don't take μ in the interval [3, 5] because in that case the numerical and reference solutions are completely superposed. Fig. 3.6 depicts the fiber distribution in some points on the closed streamline related to point (0, 1) computed using the following model parameters: k =0 .8, ǫ =0 .5, D r =0 .2 and N = 72. The fiber orientation is then defined on the unit circle, where the orientation of each fiber is depicted by a small circle. Thus, for each one of the eight points considered, the orientation of the 72 fibers involved in the simulation are depicted. The second example concerns the flow kinematics induced by two-eccentric rotating cylinders. The velocity field was computed assuming a Newtonian fluid behavior and using the well known lubrification hypothesis (justified by the fact that the flow gap is much more lower than the cylinder radius). The inner cylinder is rotating with a velocity of 0.1 rad/s. The cylinder radius are 0.01 m. and 0.015 m., being the distance between the both cylinder axes of 0.0031 m. The model parameters are k =0.8, D r =0.2, ǫ =0.5 and N =0.72. Figs. 3.7 and 3.8 depict the fiber orientation distribution along a streamline located in the secondary vortex and the main recirculation respectively.

Conclusions

In this work we have analyzed the application of a deterministic approximation of the diffusion term in the Fokker-Planck equation using smooth particles for computing its steady solution in a steady recirculating flow. The main idea of this approach lies in the introduction of the Fokker-Planck diffusion term into the advection one, which allows to proceed in a Lagrangian deterministic manner without a mesh support requirement. In spite of its simplicity, and from the first numerical results here shown, this technique seems to be very accurate.

Figure 2 . 1 .

 21 Figure 2.1. Representation of the cut-off function for different values of the ǫ parameter

Figure 3 . 2 .Figure 3 . 3 .Figure 3 . 4 .ΨFigure 3 . 5 .

 32333435 Figure 3.2. Numerical solution computed at point (0, 1) using the fixed point algorithm

Figure 3 . 6 .

 36 Figure 3.6. Fiber orientation distribution along the streamline related to the point (0, 1)

Figure 3 . 7 .Figure 3 . 8 .

 3738 Figure 3.7. Fiber orientation distribution a streamline located in the secondary vortex of a two-eccentric rotating cylinders device