Paavo Salminen
email: phsalmin@abo.fi

Pierre Vallois
email: vallois@iecn.u-nancy.fr

Marc Yor

On the excursion theory for linear diffusions

Keywords: Brownian motion, last exit decomposition, local time, infinite divisibility, spectral representation, Ornstein-Uhlenbeck process AMS Classification: 60J65, 60J60. 1

We present a number of important identities related to the excursion theory of linear diffusions. In particular, excursions straddling an independent exponential time are studied in detail. Letting the parameter of the exponential time tend to zero it is seen that these results connect to the corresponding results for excursions of stationary diffusions (in stationary state). We characterize also the laws of the diffusion prior and posterior to the last zero before the exponential time. It is proved using Krein's representations that, e.g., the law of the length of the excursion straddling an exponential time is infinitely divisible. As an illustration of the results we discuss Ornstein-Uhlenbeck processes.

Introduction and preliminaries

1.1 Throughout this paper, we shall assume that X is a linear regular recurrent diffusion taking values in R + with 0 an instantaneously reflecting boundary. Let P x and E x denote, respectively, the probability measure and the expectation associated with X when started from x ≥ 0. We assume that X is defined in the canonical space C of continuous functions ω : R + → R + . Let C t := σ{ω(s) : s ≤ t} denote the smallest σ-algebra making the co-ordinate mappings up to time t measurable and take C to be the smallest σ-algebra including all σ-algebras C t , t ≥ 0. The excursion space for excursions from 0 to 0 associated with X is a subset of C, denoted by E, and given by

E := {ε ∈ C : ε(0) = 0, ∃ ζ(ε) > 0 such that ε(t) > 0 ∀ t ∈ (0, ζ(ε)) and ε(t) = 0 ∀ t ≥ ζ(ε)}.
The notation E t is used for the trace of C t on E.

As indicated in the title of the paper our aim is to gather a number of fundamental results concerning the excursion theory for the diffusion X. In Section 2 the classical descriptions, the first one due to Itô and McKean and the second one due to Williams, are presented. In Section 3 the stationary excursions are discussed and, in particular, the description due to Bismut is reviewed. After this, in Section 4, we proceed by analyzing excursions straddling an exponential time. The paper is concluded with an example on Ornstein-Uhlenbeck processes.

Our motivation for this work arose from different origins:

• First, we would like to contribute to Professor Itô's being awarded the 1st Gauss prize, by offering some discussion and illustration of K Itô's excursion theory, see [START_REF] Itô | Poisson point processes attached to Markov processes[END_REF], when specialized to linear diffusions. The present paper also illustrates Pitman and Yor's discussion (see [START_REF] Pitman | Itø's excursion theory and its applications[END_REF] in this volume) of K. Itô's general theory of excursions for a Markov process.

• In the literature there seems to be lacking a detailed discussion on the excursion theory of linear diffusions. Information available has a very scattered character, see, e.g., Williams [START_REF] Williams | Path decompositions and continuity of local time for onedimensional diffusions[END_REF], Walsh [START_REF] Walsh | Excursions and local time[END_REF], Pitman and Yor [START_REF] Pitman | A decomposition of Bessel bridges[END_REF], [START_REF] Pitman | Decomposition at the maximum for excursions and bridges of one-dimensional diffusions[END_REF], [START_REF] Pitman | On the lengths of excursions of some Markov processes[END_REF], [START_REF] Pitman | Laplace transforms related to excursions of a one-dimensional diffusion[END_REF], Rogers [START_REF] Rogers | A guided tour through excursions[END_REF], Salminen [START_REF] Salminen | On last exit decomposition of linear diffusions[END_REF]. The general theory of excursions has been developed in Itô [START_REF] Itô | Poisson point processes attached to Markov processes[END_REF], Meyer [START_REF] Meyer | Processus de Poisson ponctuels[END_REF], Getoor [START_REF] Getoor | Excursions of a Markov process[END_REF], Getoor and Sharpe [START_REF] Getoor | Last exit times and additive functionals[END_REF], [START_REF] Getoor | Last exit decompositions and distributions[END_REF], [START_REF] Getoor | Excursions of Brownian motion and Bessel processes[END_REF], [START_REF] Getoor | Excursions of dual processes[END_REF], Blumenthal [START_REF] Blumenthal | Excursions of Markov Processes[END_REF]. Although the case with Brownian motion is well studied and understood, for textbook treatments see, e.g., Revuz and Yor [START_REF] Revuz | Continuous Martingales and Brownian Motion[END_REF] and Rogers and Williams [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF], we find it important to highlight the main formulas for more general diffusions using the traditional Fellerian terminology and language.

• To generalize some recent results (see Winkel [START_REF] Winkel | Electronic foreign-exchange markets and passage events of independent subordinators[END_REF] and Bertoin, Fujita, Roynette and Yor [START_REF] Bertoin | On a particular class of self-decomposable random variables: the duration of a bessel excursion straddling an independent exponential time[END_REF]) on infinite divisibility of the distribution of the length of the excursion of a diffusion straddling an independent exponential time.

• The Ornstein-Uhlenbeck process is one of the most essential diffusions.

To present in detail formulae for its excursions is important per se. One of the key tools hereby is the distribution of the first hitting time H y of the point y from which the excursions are observed. For y = 0 this distribution can be derived via Doob's transform (see Doob [START_REF] Doob | The Brownian movement and stochastic equations[END_REF]) which connects the Ornstein-Uhlenbeck process with standard Brownian motion (see Sato [START_REF] Sato | Evaluation of the first passage time probability to a square root boundary for the Wiener process[END_REF], and Göing-Jaeschke and Yor [START_REF] Göing-Jaeschke | A clarification note about hitting times densities for Ornstein-Uhlenbeck processes[END_REF]). For arbitary y the distribution is very complicated; for explicit expressions via series expansions, see Ricciardi and Sato [START_REF] Ricciardi | First passage-time density and moments of the Ornstein-Uhlenbeck process[END_REF], Linetsky [START_REF] Linetsky | Computing hitting time densities for CIR and OU diffusions: applications to mean reverting models[END_REF] and Alili, Patie and Pedersen [START_REF] Alili | Representations of the first hitting time density of an Ornstein-Uhlenbeck process[END_REF]. We will focus on excursions from 0 to 0 and relate our work to earlier papers by Hawkes and Truman [START_REF] Hawkes | Statistics of local time and excursions for the Ornstein-Uhlenbeck process[END_REF], Pitman and Yor [START_REF] Pitman | On the lengths of excursions of some Markov processes[END_REF], and Salminen [START_REF] Salminen | On conditional Ornstein-Uhlenbeck processes[END_REF]. Due to the symmetry of the Ornstein-Uhlenbeck process around 0, it is sufficient for our purposes to consider only positive excursions -the treatment of negative ones is similar -and view the process with values in R + and 0 being a reflecting boundary.

1.2

In this subsection we introduce the basic notation and facts concerning linear diffusions needed in the sequel. A main source of information remains Itô and McKean [START_REF] Itô | Diffusion Processes and Their Sample Paths[END_REF], see also Rogers and Williams [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF], and Borodin and Salminen [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF].

(i) Speed measure m associated with X is a measure on R + which satisfies for all 0 < a < b < ∞

0 < m((a, b)) < ∞.
For simplicity, it is assumed that m does not have atoms. An important fact is that X has a jointly continuous transition density p(t; x, y) with respect to m, i.e.,

P x (X t ∈ A) = A p(t; x, y) m(dy),
where A is a Borel subset of R + . Moreover, p is symmetric in x and y, that is, p(t; x, y) = p(t; y, x). The Green or the resolvent kernel of X is defined for λ > 0 as

R λ (x, y) = ∞ 0
dt e -λt p(t; x, y).

(ii) Scale function S is an increasing and continuous function which can be defined via the identity

P x (H a < H b ) = S(b) -S(x) S(b) -S(a) , 0 ≤ a < x < b, (1)
where H . denotes the first hitting time, i.e.,

H y := inf{t : X t = y}, y ≥ 0.
We normalize by setting S(0) = 0. Due to the recurrence assumption it holds S(+∞) = +∞. Recall that {S(X t∧H 0 ) : t ≥ 0} is a continuous local P x -martingale for every x ≥ 0 (see, e.g., Rogers and Williams [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF] p. 276). It is easily proved that S(X) = {S(X t ) : t ≥ 0} is a (recurrent) diffusion taking values in R + . The scale function associated with S(X) is the identity mapping x → x, x ≥ 0, and we say that S(X) is in natural scale. Clearly, also for S(X) the boundary point 0 is instantaneously reflecting. Using the Skorokhod reflection equation it is seen that S(X) is a P x -submartingale (cf. Meleard [START_REF] Meleard | Application du calcul stochastique à l'etude de processus de Markov réguliers sur [0, 1[END_REF] Proposition 1.4 where the semimartingale decomposition is given in case there are two reflecting boundaries).

(iii) The infinitesimal generator of X can be expressed as the generalized differential operator G = d dm d dS acting on functions f belonging to the appropriately defined domain D(G) of G (see Itô and McKean [22], Freedman [START_REF] Freedman | Brownian motion and diffusions. Holden-Day[END_REF], Borodin and Salminen [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF]). In particular, since 0 is assumed to be reflecting then

f ∈ D(G) implies that f + (0) := lim x↑0 f (x) -f (0) S(x) -S(0) = 0.
(iv) The distribution of the first hitting time of a point y > 0 has a P xdensity:

P x (H y ∈ dt) = f xy (t) dt.
This density can be connected with the derivative of a transition density of a killed diffusion obtained from X. To explain this, introduce the sample paths

X (y) t := X t , t < H y , ∂, t ≥ H y ,
where ∂ is a point isolated from R + (a "cemetary" point). Then { X (y) t

: t ≥ 0} is a diffusion with the same scale and speed as X. Let p denote the transition density of X (y) with respect to m. Then, e.g., for x > y

f xy (t) = lim z↓y p(t; x, z) S(z) -S(y) . (2)
For a fixed x and y, the mapping t → f xy (t) is continuous, as follows from the eigen-differential expansions and discussion in Itô and McKean p. 153 and 217 (see also Kent [24], [START_REF] Kent | The spectral decomposition of a diffusion hitting time[END_REF]). Recall also the following formula for the Laplace transform of H y

E x e -αHy = R α (x, y) R α (y, y) , (3)
which leads to ∞ 0 m(dx) E x e -αHy = 1 αR α (y, y) .

(v) There exists a jointly continuous family of local times

{L (y) t : t ≥ 0, y ≥ 0}
such that X satisfies the occupation time formula

t 0 ds h(X s ) = ∞ 0 h(y)L (y) t m(dy),
where h is a nonnegative measurable function (see, e.g., Rogers and Williams [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF] 49.1 Theorem p. 289). Consequently,

L (y) t = lim δ↓0 1 m((y -δ, y + δ)) t 0 1 (y-δ,y+δ) (X s ) ds.
For a fixed y introduce the inverse of L (y) via τ (y) := inf{s : L (y) s > }.

Then τ (y) = {τ (y) : ≥ 0} is an increasing Lévy process, in other words, a subordinator and its Lévy exponent is given by

E y exp(-λτ (y) ) = exp (-/R λ (y, y)) = exp(- ∞ 0 ν (y) (dv)(1 -e -λv )), (4)
where ν (y) is the Lévy measure of τ (y) . The assumption that the speed measure does not have atoms implies that τ (y) does not have a drift.

In case y = 0 we write simply L, τ and ν.

1.3 Assuming that X is started from 0 we define for t > 0

G t := sup{s ≤ t : X s = 0} and D t := inf{s ≥ t : X s = 0}. ( 5
)
The last exit decomposition at a fixed time t states that for u < t < v

P 0 (G t ∈ du, X t ∈ dy, D t ∈ dv) = p(u; 0, 0) f y0 (t -u) f y0 (v -t) du dv m(dy). (6)
In fact, this trivariate distribution is only the skeleton of a more complete body of processes:

{X u : u ≤ G t }, {X Gt+v : v ≤ t -G t }, and {X t+v : v ≤ D t -t} (7)
the distributions of which we now characterize following Salminen [START_REF] Salminen | On last exit decomposition of linear diffusions[END_REF]. For general approaches; see Getoor and Sharpe [START_REF] Getoor | Last exit times and additive functionals[END_REF], [START_REF] Getoor | Last exit decompositions and distributions[END_REF], and Maisonneuve [START_REF] Maisonneuve | Exit systems[END_REF]. Let x, y ∈ R + and u > 0 be given. Denote by (X x,u,y , P x,u,y ) the diffusion bridge from x to y of length u constructed from X, i.e., the measure P x,u,y governing X x,u,y is the conditional measure associated with X started from x and conditioned to be at y at time u. The bridge X x,u,y is a strong non-time-homogeneous Markov process defined on the time axis [0, u). For the first component in [START_REF] Chung | Excursions in Brownian motion[END_REF], we have conditionally on

G t = u {X s : 0 ≤ s < G t } d = {X 0,u,0 s : 0 ≤ s < u} (8)
For the second component in [START_REF] Chung | Excursions in Brownian motion[END_REF] consider the process X (y) as introduced in (iv) above with y = 0. We write simply X instead of X (0) . For positive x and y let X x,u,y denote the bridge from x to y of length u constructed, as above, from X. The measure P x,u,y governing X x,u,y can be extended by taking (in the weak sense) P 0,u,y := lim x↓0 P x,u,y .

We let X 0,u,y denote the process associated with P 0,u,y . Then, conditionally on G t = u and X t = y,

{X Gt+s : 0 ≤ s < t -G t } d = { X 0,t-u,y s : 0 ≤ s < t -u}. ( 9
)
For the final part in [START_REF] Chung | Excursions in Brownian motion[END_REF], by the Markov property, we have conditionally on X t = y

{X t+s : s < D t -t} d = { X s : s ≥ 0}, ( 10
)
where X 0 = y.

2 Two descriptions of the Itô measure

Description due to Itô and McKean

We discuss the description of the Itô measure n where the excursions are studied by conditioning with respect to their lifetimes. Let X be as in section 1.3 and p(t; x, y) its transition density with respect to the speed measure, in other words, P x ( X t ∈ dy) = P x (X t ∈ dy; t < H 0 ) = p(t; x, y) m(dy).

The Lévy measure ν of τ is absolutely continuous with respect to the Lebesgue measure, and the density -which we also denote by ν -is given by

ν(v) := ν(dv)/dv = lim x↓0 lim y↓0 p(v; x, y) S(x)S(y) =: p ↑ (v; 0, 0). (11)
In Section 1.3 we have introduced the bridge X x,t,y and the measure P x,t,y associated with it. The family of probability measures { P x,t,y : x > 0, y > 0} is weakly convergent as y ↓ 0 thus defining P x,t,0 for all x > 0. Intuitively, this is the process X conditioned to hit 0 at time t. Moreover, letting now x ↓ 0 we obtain a measure which we denote by P 0,t,0 which governs a non-time homogeneous Markov process X 0,t,0 starting from 0, staying positive on the time interval (0, t) and ending at 0 at time t.

Theorem 1. a. The law of the excursion life time ζ under the Itô excursion measure n is equal to the Lévy measure of the subordinator {τ } ≥0 and is given by

n(ζ ∈ dv) = ν(dv) = p ↑ (v; 0, 0) dv. (12)
b. The Itô measure can be represented as the following integral

n(dε) = ∞ 0 n(ζ ∈ dv) P 0,v,0 (dε). ( 13
)
Moreover, the finite dimensional distributions of the excursion are characterized for

0 < t 1 < t 2 < • • • < t n and x i > 0, i = 1, 2, . . . , n by n(ε t 1 ∈ dx 1 , ε t 2 ∈ dx 2 , . . . , ε tn ∈ dx n ) = m(dx 1 ) f x 1 0 (t 1 ) p(t 2 -t 1 ; x 1 , x 2 ) m(dx 2 ) (14) × . . . p(t n -t n-1 ; x n-1 , x n ) m(dx n ).
In particular, the excursion entrance law is given by

n(ε t ∈ dx) = m(dx) f x0 (t),
and it holds

n(ζ > t) = ∞ 0 n(ε t ∈ dx) = ∞ 0 m(dx) f x0 (t). (15)
Combining the formulas ( 12) and ( 13) with the last exit decomposition (6) leads to a curious relation between the transition densities p and p ↑ . Proposition 2. The functions p(t; 0, 0) and p ↑ (t; 0, 0) satisfy the identity

t 0 du p(u; 0, 0) ∞ t-u dv p ↑ (v; 0, 0) = 1. ( 16
)
Proof. From ( 12) and ( 15) we may write

∞ t dv p ↑ (v; 0, 0) = n(ζ > t) = ∞ 0 m(dx) f x0 (t).
Consequently, identity ( 16) can be rewritten as

t 0 du p(u; 0, 0) ∞ 0 m(dx) f x0 (t -u) = 1, (17)
but, in view of the last exit decomposition (6), identity [START_REF] Getoor | Excursions of dual processes[END_REF] states that the last exit from 0 when starting from 0 takes place with probability 1 before t, in other words,

P 0 (G t ≤ t) = 1,
which, of course, is trivially true.

Remark 3. For another approach to [START_REF] Getoor | Excursions of Brownian motion and Bessel processes[END_REF] notice that it follows from ( 4) and ( 11)

1 R λ (0, 0) = ∞ 0 dv p ↑ (v; 0, 0) 1 -e -λv .
Hence, from the definition of the Green kernel,

1 = ∞ 0 du e -λu p(u; 0, 0) ∞ 0 dv p ↑ (v; 0, 0)(1 -e -λv ). ( 18
)
Consequently,

1 λ = ∞ 0 du e -λu p(u; 0, 0) ∞ 0 dv e -λv ∞ v ds p ↑ (s; 0, 0) = ∞ 0 du ∞ 0 dv e -λ(u+v) p(u; 0, 0) ∞ v ds p ↑ (s; 0, 0),
from which ( 16) is easily deduced.

Description due to Williams

In the approach via the lengths of the excursions the focus is first on the time axis. In Williams' description (see Williams [START_REF] Williams | Path decompositions and continuity of local time for onedimensional diffusions[END_REF], and Rogers [START_REF] Rogers | Williams chracterization of the Brownian excursion law : proof and applications[END_REF], [START_REF] Rogers | A guided tour through excursions[END_REF]) the starting point of the analysis is on the space axis since the basic conditioning is with respect to the maximum of an excursion. To formulate the result,

let for ε ∈ E M (ε) := sup{ε t : 0 < t < ζ(ε)}.
The key element in Williams' description is the diffusion X ↑ obtained by conditioning X not to hit 0. We use the notation P ↑ and E ↑ for the measure and the expectation associated with X ↑ . To define this process rigorously set for a bounded F t ∈ C t , t > 0,

E ↑ x (F t ) := lim a↑+∞ E x (F t ; t < H a | H a < H 0 ) = lim a↑+∞ E x (F t ; t < H a < H 0 ) P x (H a < H 0 ) = lim a↑+∞ E x (F t S(X t ) ; t < H a ∧ H 0 ) S(x) ,
where the Markov property and formula (1) for the scale function are applied. The monotone convergence theorem yields

E ↑ x (F t ) = 1 S(x) P x (F t S(X t ) ; t < H 0 ) ,
in other words, the desired conditioning is realized as Doob's h-transform of X by taking h to be the scale function of X. It is easily deduced that the transition density and the speed measure associated with X ↑ are given by p ↑ (t; x, y) := p(t; x, y) S(y)S(x) , m ↑ (dy) := S(y) 2 m(dy).

We remark that the boundary point 0 is entrance-not-exit for X ↑ and, therefore, X ↑ can be started from 0 after which it immediately enters (0, ∞) and never hits 0.

Theorem 4. a. The law of the excursion maximum M under the Itô excursion measure n is given by

n(M ≥ a) = 1 S(a)
.

b. The Itô excursion measure n can be represented via

n(dε) = ∞ 0 n(M ∈ da) Q ( * ,a) (dε),
where Q ( * ,a) is the distribution of two independent X ↑ processes put back to back run from 0 until they first hit level a.

As an illustration, we give the following formula

n(1 -exp(- ζ 0 ds V (ε s ))) = ∞ 0 n(M ∈ da) 1 -E ↑ 0 (exp(- Ha 0 du V (ω u ))) 2 .
If V ≥ 0, this quantity is the Lévy exponent of the subordinator

τ 0 ds V (X s ) : ≥ 0 , that is, E exp -α τ 0 ds V (X s ) = exp -n 1 -exp -α ζ 0 ds V (ε s ) .
Comparing the descriptions of the Itô excursion measure in Theorem 1 (in particular formula [START_REF] Getoor | Last exit decompositions and distributions[END_REF]) and in Theorem 4 hints that the processes X and X ↑ have, in addition to conditioning relationship, also a time reversal relationship. This is due to Williams [START_REF] Williams | Path decompositions and continuity of local time for onedimensional diffusions[END_REF], who particularized to the case of diffusions the general time reversal result, obtained by Nagasawa [START_REF] Nagasawa | Time reversions of Markov processes[END_REF]. See also [START_REF] Revuz | Continuous Martingales and Brownian Motion[END_REF] p. 313, and [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF]

{ X s : 0 ≤ s < H 0 } d = {X ↑ Λx-s : 0 ≤ s < Λ x }, (19)
where X 0 = x and X ↑ 0 = 0.

Stationary excursions; Bismut's description

Consider the diffusion X with the time parameter t taking values in the whole of R. In the case m(R + ) < ∞ the measure governing X can be normalized to be a probability measure. Indeed, in this case the distribution of X t is for every t ∈ R defined to be

P(X t ∈ dx) = m(dx)/m(R + ) =: m(dx).
Recall from (5) the definitions of G t and D t , and introduce also ∆ t := D t -G t .

Theorem 6. Assume that m(R + ) < ∞. Then the joint distribution of t-G t and D t -t is given by

P(t -G t ∈ du, D t -t ∈ dv)/dudv = ∞ 0 m(dy) f y0 (u) f y0 (v) = ν(u + v)/m(R + ).
Consequently, for ∆ t it holds

P(∆ t ∈ du)/du = u ν(u)/m(R + ). ( 20
)
Moreover, the law of the process

{X Gt+v : v ≤ ∆ t } is given by ζ(ε) n(dε)/m(R + ), ( 21
)
where n(dε) is the Itô measure as introduced in Theorem 1 and 4 and ζ denotes the length of an excursion.

Proof. The density of (t -G t , D t -t) is derived using the time reversibility of the diffusion X, i.e.,

{X t : t ∈ R} d = {X -t : t ∈ R},
and the conditional independence given X t . The fact that the density can be expressed via the density of the Lévy measure is stated (and proved) in Proposition 12 below, see formulas [START_REF] Küchler | On spectral measures of strings and excursions of quasi diffusions[END_REF] and [START_REF] Lévy | Sur certains processus stochastiques homogénes[END_REF]. To compute the distribution of ∆ t is elementary from the joint distribution of t -G t and D t -t. For these results, we refer also to Kozlova and Salminen [START_REF] Kozlova | A note on occupation times of stationary processes[END_REF]. The statement concerning the law of {X Gt+v : v ≤ ∆ t } has been proved in Pitman [START_REF] Pitman | Stationary excursions[END_REF] (see Theorem p. 290 point (iii) and the formulation for excursions on p. 293 and 294) -all that remains for us to do is to find the right normalization constant, but this is fairly obvious, e.g., from the density of ∆ t .

If m(R + ) = ∞ the measure associated with X is still well-defined but "only" σ-finite. In this case, the distribution of X t is plainly taken to be m. From [START_REF] Itô | Poisson point processes attached to Markov processes[END_REF] it is seen that we are faced with a representation of the Itô measure via stationary excursions valid in both cases m(R + ) < ∞ and m(R + ) = ∞. We focus now on this representation as displayed in [START_REF] Itô | Diffusion Processes and Their Sample Paths[END_REF] below, and present a proof of the representation using the diffusion theory (this provides, of course, also a proof of ( 21)). We remark that in Pitman [START_REF] Pitman | Stationary excursions[END_REF] a more general case concerning homogeneous random sets is proved, and, hence, it seems worthwhile to give a "direct" proof in the diffusion case.

Theorem 7. Let F be a measurable non-negative functional defined in the excursion space E. Then up to a normalization

n(F (ε)) = E 1 ∆ t F (X Gt+s : 0 ≤ s ≤ ∆ t ) . ( 22
)
In particular, the process {X Gt+s : 0 ≤ s ≤ ∆ t } conditionally on ∆ t = v is identical in law with the excursion bridge X 0,v,0 as introduced in Section 2.1.

Proof. Without loss of generality, we take t = 0. From [START_REF] Hawkes | Statistics of local time and excursions for the Ornstein-Uhlenbeck process[END_REF] we have

n(f (ζ)) = ∞ 0 f (a) ν(a) da = E 1 ∆ 0 f (∆ 0 ) .
Therefore, it is enough (cf. Theorem 1) to prove that

n(F (ε) | ζ = u) = E (F (X G 0 +s : 0 ≤ s ≤ ∆ 0 ) | ∆ 0 = u) . ( 23
) Define for 0 ≤ s 1 < s 2 < • • • < s n ≤ u A 1,n := {X G 0 +s 1 ∈ dy 1 , . . . , X G 0 +sn ∈ dy n },
and consider

E (A 1,n | ∆ 0 = u) = ∞ y=0 u v=0 E (A 1,n , -G 0 ∈ dv, X 0 ∈ dy | ∆ 0 = u) = ∞ y=0 u v=0 E (A 1,n | ∆ 0 = u, G 0 = -v, X 0 = y) ×P (-G 0 ∈ dv, X 0 ∈ dy | ∆ 0 = u) .
From the description of the process X, the conditional independence and the equality of the laws of the past and future given X 0 , and using formula [START_REF] Hawkes | Statistics of local time and excursions for the Ornstein-Uhlenbeck process[END_REF] we obtain

P (-G 0 ∈ dv, X 0 ∈ dy | ∆ 0 = u) = 1 u ν(u) f y,0 (v) f y,0 (u -v) m(dy)dv. ( 24
)
Letting k be such that -v + s k < t < -v + s k+1 , if any, we write applying again the conditional independence

E (A 1,n | ∆ 0 = u, G 0 = -v, X 0 = y) = E (A 1,k A k+1,n | ∆ 0 = u, G 0 = -v, X 0 = y) = E (A 1,k | G 0 = -v, X 0 = y) E (A k+1,n | D 0 = u -v, X 0 = y) .
Recall from Introduction section 1.2 (iv) the notation X for the diffusion X killed when it hits 0. As in section 1.3 of Introduction we may construct the bridge X y,v,0 starting from y having the length v and ending at 0. We let P y,v,0 denote the measure associated with X y,v,0 . With these new notations,

E (A 1,k | G 0 = -v, X 0 = y) = P y,v,0 (ω v-s k ∈ dy k , . . . , ω v-s 1 ∈ dy 1 ) = 1 f y0 (v) p(v -s k ; y, y k ) m(dy k ) p(s k -s k-1 ; y k , y k-1 ) m(dy k-1 ) × . . . p(s 2 -s 1 ; y 2 , y 1 ) m(dy 1 ) f y 1 0 (s 1 )
and

E (A k+1,n | D 0 = u -v, X 0 = y) = P y,v,0 ω s k+1 -v ∈ dy k+1 , . . . , ω sn-v ∈ dy n = 1 f y0 (u -v) p(s k+1 -v; y, y k+1 ) m(dy k+1 ) × p(s k+2 -s k+1 ; y k+1 , y k+2 ) m(dy k+2 ) × . . . p(s n -s n-1 ; y n-1 , y n ) m(dy n ) f yn0 (u -s n ).
Using now [START_REF] Kent | Eigenvalue expansions for diffusion hitting times[END_REF] and formulas above we have after some rearranging and applying the symmetry of the transition density p

E (A 1,n | ∆ 0 = u) = 1 u ν(u) m(dy 1 ) f y 1 0 (s 1 ) p(s 2 -s 1 ; y 1 , y 2 ) m(dy 2 ) × . . . u 0 dv ∞ 0 m(dy) p(v -s k ; y k , y) p(s k+1 -v; y, y k+1 ) × . . . p(s n -s n-1 ; y n-1 , y n ) m(dy n ) f yn0 (u -s n ).
Performimg the integration yields

E (A 1,n | ∆ 0 = u) = 1 ν(u) m(dy 1 ) f y 1 0 (s 1 ) p(s 2 -s 1 ; y 1 , y 2 ) m(dy 2 ) × . . . p(s n -s n-1 ; y n-1 , y n ) m(dy n ) f yn0 (u -s n ),
and this means that (23) holds completing the proof.

Remark 8. The formula [START_REF] Itô | Diffusion Processes and Their Sample Paths[END_REF] was derived for Brownian motion by Bismut [START_REF] Bismut | Last exit decompositions and regularity at the boundary of transition probabilities[END_REF]. The connection with the Palm measure and stationary processes is discussed in Pitman [START_REF] Pitman | Stationary excursions[END_REF]. In fact, Bismut describes in the Brownian case the law of the process {X Gt+s : 0 ≤ s ≤ ∆ t } in terms of two independent 3-dimensional Bessel processes started from 0 and killed at the last exit time from an independent level distributed according to the Lebesgue measure (see [START_REF] Bismut | Last exit decompositions and regularity at the boundary of transition probabilities[END_REF] and [START_REF] Revuz | Continuous Martingales and Brownian Motion[END_REF] for details).

On the excursion straddling an independent exponential time

In the literature one can find several papers devoted to the properties of excursions straddling a fixed time t; first of all, Lévy's fundamental paper [START_REF] Lévy | Sur certains processus stochastiques homogénes[END_REF], which contains a lot about the zero set of Brownian motion, its (inverse) local time, excursions, and so on. See also Chung [START_REF] Chung | Excursions in Brownian motion[END_REF] starting from Lévy's paper [START_REF] Lévy | Sur certains processus stochastiques homogénes[END_REF], Durret and Iglehart [START_REF] Durrett | Functionals of Brownian meander and Brownian excursion[END_REF], and Getoor and Sharpe [START_REF] Getoor | Excursions of Brownian motion and Bessel processes[END_REF], [START_REF] Getoor | Excursions of dual processes[END_REF]. In fact, the last exit decomposition (6) lies in the heart of these studies (see Getoor and Sharpe [START_REF] Getoor | Last exit times and additive functionals[END_REF], [START_REF] Getoor | Last exit decompositions and distributions[END_REF]). However, it seems to us that excursions straddling an exponential time are not so much analyzed. Here we make some remarks on this subject.

Let T be an exponentially distributed random variable with parameter α > 0, independent of X, and define (iii) The law of the process {X G T +v : v ≤ ∆ T } is given by

1 Φ(α) 1 -e -αζ(ε) n(dε). ( 26
)
where n(dε) is the Itô measure associated with the excursions away from 0 for X and ζ denotes the length of an excursion.

Proof. Let F 1 and F 2 be two nonnegative functionals of continuous functions and consider

E 0 (F 1 (X u : u ≤ G T ) F 2 (X G T +v : v ≤ ∆ T )) = α ∞ 0 dt e -αt E 0 (F 1 (X u : u ≤ G t ) F 2 (X Gt+v : v ≤ ∆ t )) = αE 0 τ τ - dt e -αt F 1 (X u : u ≤ τ -) F 2 (X τ -+v : v ≤ τ -τ -) = E 0 ∞ 0 d e -α τ F 1 (X u : u ≤ τ ) × n(dε) 1 -e -α ζ(ε) F 2 (ε s : s ≤ ζ(ε)) = Φ(α)E 0 ∞ 0 d e -α τ F 1 (X u : u ≤ τ ) × n(dε) 1 -e -α ζ(ε) Φ(α) F 2 (ε s : s ≤ ζ(ε)).
where the third equality is based on the properties of the Poisson random measure associated with the excursions (see Revuz and Yor [START_REF] Revuz | Continuous Martingales and Brownian Motion[END_REF] Master Formula p. 475).

Remark 10. Notice that letting α → 0 in (26) and using

lim α→0 α Φ(α) = lim α→0 α R α (0, 0) = 1/m(R + ) ( 27
)
yield the probability law of the excursion straddling a fixed time in the stationary setting, cf. ( 21) in Theorem 6.

As a corollary of Theorem 9 we have the following results which show that after conditioning the quantities do not depend on α, and in this context α is entirely "contained" in G T and ∆ T . The formulas should be compared with (8), [START_REF] Durrett | Functionals of Brownian meander and Brownian excursion[END_REF], and [START_REF] Dym | Gaussian processes, function theory, and the inverse spectral problem[END_REF]. The distributions of G T and ∆ T are given, respectively, in ( 40) and ( 35) below.

Corollary 11. For any nonnegative functionals F 1 and F 2 of continuous functions it holds

E 0 (F 1 (X u : u ≤ G T ) | G T = g) = E 0 (F 1 (X u : u ≤ g) | X g = 0) = E 0,g,0 (F 1 (ω u : u ≤ g)) (28)
and

E 0 (F 2 (X G T +v : v ≤ ∆ T ) | ∆ T = h) = E ↑ 0 (F 2 (X v : v ≤ h) | X h = 0) . = E 0,h,0 (F 2 (ω s : s ≤ h)) (29)
Proof. The statement ( 28) can be obtained from the corresponding result for fixed time as presented in [START_REF] Doob | The Brownian movement and stochastic equations[END_REF]. Also ( 29) can be derived from the fixed time result but we prefer to present here a proof based on the Master Formula. For this consider for δ > 0

E 0 F 2 (X G T +v : v ≤ ∆ T ) 1 {h≤∆ T <h+δ} = n(dε) 1 -e -α ζ(ε) Φ(α) F 2 (ε s : s ≤ ζ(ε)) 1 {h≤ζ(ε)<h+δ} .
Using the description ( 13) of the Itô excursion law we obtain

E 0 F 2 (X G T +v : v ≤ ∆ T ) 1 {h≤∆ T <h+δ} = h+δ h du ν(u) 1 -e -α u Φ(α) E 0,u,0 (F 2 (ω s : s ≤ u)).
Applying the explicit form of the distribution of ∆ T as given in [START_REF] Meyer | Processus de Poisson ponctuels[END_REF] and letting δ ↓ 0 leads to (29).

P 0 (T -G T ∈ du, D T -T ∈ dv) = dudv α R α (0, 0) e -αu ∞ 0 m(dy) f y0 (u) f y0 (v). ( 30
) = α e -αu ν(u + v) Φ(α) dudv. ( 31
)
In particular,

ν(u + v) = ∞ 0 m(dy) f y0 (u) f y0 (v). ( 32
)
Proof. From ( 6),

P 0 (t -G t ∈ du, D t -t ∈ dv) = dudv p(t -u; 0, 0) 1 {u≤t,v≥0} ∞ 0 m(dy) f y0 (u) f y0 (v),
and, hence,

P 0 (T -G T ∈ du, D T -T ∈ dv) = dudv α ∞ u dt e -αt p(t -u; 0, 0) ∞ 0 m(dy) f y0 (u) f y0 (v),
from which (30) follows. To derive [START_REF] Lévy | Sur certains processus stochastiques homogénes[END_REF], we apply again the Master Formula (see Revuz and Yor [START_REF] Revuz | Continuous Martingales and Brownian Motion[END_REF] p. 475). For this, let (u, v) → ϕ(u, v) be a nonnegative and Borel measurable function and define

Q(ϕ) := E 0 (ϕ(T -G T , D T -T )).
Letting τ denote the inverse of the local time L at 0 we have

Q(ϕ) = E 0   ≥0 ϕ(T -τ -, τ -T ) 1 {τ -<T <τ }   = E 0 R 2 + ϕ(T -τ , z + τ -T ) 1 {τ <T <τ +z} ν(z) dzd ,
since {( , τ ) : ≥ 0} is a Poisson point process with Lévy measure d dν, and T is independent of {τ : ≥ 0}. Apply next that T is exponentially distributed to obtain

Q(ϕ) = E 0 R 2 + ν(z)dzd τ +z τ dt α e -αt ϕ(t -τ , z + τ -t) = αE 0 R 3 + ν(z)e -α(x+τ ) ϕ(x, z -x) 1 {x≤z} dxdzd ,
where we have substituted x = t -τ . Furthermore, setting y = z -x yields

Q(ϕ) = αE 0 R 3 + ν(y + x)e -α(x+τ ) ϕ(x, y) dxdyd = αE 0 ∞ 0 e -α τ d R 2 + ϕ(x, y) e -α x ν(y + x) dxdy,
and (31) follows now easily from (4). The equality ( 32) is an immediate consequence of ( 30) and [START_REF] Lévy | Sur certains processus stochastiques homogénes[END_REF].

Corollary 13. 1. The densities for T -G T , D T -T, and ∆ T are given, respectively, by

P 0 (T -G T ∈ du)/du = α Φ(α) e -αu ∞ u ν(z)dz, (33)
P 0 (D T -T ∈ dv)/dv = α Φ(α) e αv ∞ v e -αz ν(z)dz, (34)
and

P 0 (∆ T ∈ da)/da = (1 -e -αa )ν(a) Φ(α) . ( 35
)
2. The joint density of T -G T and ∆ T is

P 0 (T -G T ∈ du, ∆ T ∈ da)/du da = α Φ(α) e -αu ν(a), u ≤ a. ( 36
)
3. The density of T -G T conditionally on ∆ T = a is

P 0 (T -G T ∈ du | ∆ T = a)/du = α 1 -e -αa e -αu , u ≤ a. ( 37
)
Proposition 14. The joint Laplace transform of G T and D T is given by

E 0 (exp (-γ 1 G T -γ 2 D T )) = Φ(γ 2 + α) -Φ(γ 2 ) Φ(γ 1 + γ 2 + α) . (38)
In particular,

E 0 e -γ∆ T = Φ(γ + α) -Φ(γ) Φ(α) , ( 39
)
and the random variables G T and ∆ T are independent. The density of G T is given by P 0 (G T ∈ du)/du = Φ(α) e -αu p(u; 0, 0). ( 40)

Proof. The formula [START_REF] Pitman | On the lengths of excursions of some Markov processes[END_REF] for the density of G T is obtained from ( 6) by integrating. The independence of G T and ∆ T follows immediately from [START_REF] Pitman | A decomposition of Bessel bridges[END_REF].

To derive the joint Laplace transform of G T and D T , consider

E 0 (exp (-γ 1 G T -γ 2 D T )) = v u dt α e -αt e -γ 1 u-γ 2 v P 0 (G t ∈ du, D t ∈ dv).
Applying the last exit decomposition formula ( 6) yields

E 0 (exp (-γ 1 G T -γ 2 D T )) = ∞ 0 du e -γ 1 u p(u; 0, 0) ∞ u dv e -γ 2 v v u dt α e -αt × ∞ 0 m(dy) f y0 (t -u) f y0 (v -t) = ∞ 0 du e -γ 1 u p(u; 0, 0) ∞ 0 da e -γ 2 (a+u) a+u u
dt α e -αt p(u; 0, 0)

× ∞ 0 m(dy) f y0 (t -u) f y0 (a + u -t) = ∞ 0 du e -(γ 1 +γ 2 )u p(u; 0, 0) ∞ 0 da e -γ 2 a a 0
db α e -α(b+u) p(u; 0, 0)

× ∞ 0 m(dy) f y0 (b) f y0 (a -b) = α ∞ 0 du e -(γ 1 +γ 2 +α)u p(u; 0, 0) ∞ 0 m(dy) ∞ 0 da e -γ 2 a × a 0 db e -αb f y0 (b) f y0 (a -b) = αR γ 1 +γ 2 +α (0, 0) ∞ 0 m(dy) E y e -(γ 2 +α)H 0 E y e -γ 2 H 0 .
To proceed, we have

∞ 0 m(dy) E y e -(γ 2 +α)H 0 E y e -γ 2 H 0 = 1 R γ 2 +α (0, 0)R γ 2 (0, 0) ∞ 0 m(dy) R γ 2 +α (y, 0)R γ 2 (y, 0).
The integral term in this expression can be evaluated:

∞ 0 m(dy) R α+γ 2 (y, 0)R γ 2 (y, 0) = ∞ 0 m(dy) ∞ 0 dt e -(α+γ 2 )t p(t; y, 0) ∞ 0 ds e -γ 2 s p(s; y, 0) = ∞ 0 dt e -(α+γ 2 )t ∞ 0 ds e -γ 2 s p(t + s; 0, 0) = ∞ 0 dt e -(α+γ 2 )t ∞ t du e -γ 2 (u-t) p(u; 0, 0) = ∞ 0 du e -γ 2 u 1 -e -αu α p(u; 0, 0), = 1 α (R γ 2 (0, 0) -R γ 2 +α (0, 0)) .
where the Chapman-Kolmogorov equation and the symmetry of the transition density p is applied, and by [START_REF] Kent | The spectral decomposition of a diffusion hitting time[END_REF] this completes the proof.

Remark 15. 1. From Proposition 12 it is seen that the density of ∆ T can also be written in the form

P 0 (∆ T ∈ da)/da = α Φ(α) ∞ 0 m(dy) a 0 db e -αb f y0 (b) f y0 (a -b),
which taking into account [START_REF] Meyer | Processus de Poisson ponctuels[END_REF] leads to the identity

(1 -e -αa ) α ν(a) = ∞ 0 m(dy) a 0 db e -αb f y0 (b) f y0 (a -b).
Let here α → 0 to obtain

ν(a) = ∞ 0 m(dy) a 0 db a f y0 (b) f y0 (a -b). (41)
It is interesting to compare this expression with the following one obtained from ( 32)

ν(a) = ∞ 0 m(dy) f y0 (b) f y0 (a -b). (42)
The fact that the right hand sides of ( 41) and ( 42) do not depend on b can also be explained via the Chapman-Kolmogorov equation.

2. We may study distributions associated with G t , D t and ∆ t in the stationary case, i.e., if m(R + ) < ∞, by letting α → 0, as observed in Remark 10.

From Proposition 12 and Corollary 13 we deduce the following results take t = 0:

P(-G 0 ∈ du, D 0 ∈ dv)/dudv = 1 m(R + ) ν(u + v). P(-G 0 ∈ du)/du = P(D 0 ∈ du)/du = 1 m(R + ) ∞ u ν(v) dv, P(∆ 0 ∈ da)/da = 1 m(R + ) a ν(a).
Moreover, letting

Z T := (T -G T )/∆ T then (Z T , ∆ T ) converges in distribu- tion as α → 0 to (U, ∆)
, where U and ∆ are independent with U uniformly distributed on (0, 1) and ∆ is distributed as ∆ 0 (cf. Theorem 6).

Infinite divisibility

In the paper by Bertoin et al. [START_REF] Bertoin | On a particular class of self-decomposable random variables: the duration of a bessel excursion straddling an independent exponential time[END_REF] it is proved that the distribution of ∆ T for a Bessel process with dimension d = 2(1 -α), 0 < α < 1, is infinitely divisible (in fact, self-decomposable) and the Lévy measure associated with this distribution is computed. In this section we show that the distribution of ∆ T is infinitely divisible in general, i.e., for all regular and recurrent diffusions. Moreover, we also prove that the distributions of T -G T and D T -T have this property. The key to these results is the Krein representation of the density of the Lévy measure ν (see Knight [START_REF] Knight | Characterization of the Lévy measures of inverse local times of gap diffusion[END_REF], Kent [START_REF] Kent | The spectral decomposition of a diffusion hitting time[END_REF], Küchler and Salminen [START_REF] Küchler | On spectral measures of strings and excursions of quasi diffusions[END_REF], and in general on Krein's theory of strings Kotani and Watanabe [START_REF] Kotani | Krein's spectral theory of strings and generalized diffusion processes[END_REF], Dym and McKean [START_REF] Dym | Gaussian processes, function theory, and the inverse spectral problem[END_REF]) according to which

ν(a) = ∞ 0 e -az M (dz), (43)
where the measure M has the properties

∞ 0 M (dz) z(z + 1) < ∞ and ∞ 0 M (dz) z = ∞.
. Proof. As seen from ( 33), [START_REF] Meleard | Application du calcul stochastique à l'etude de processus de Markov réguliers sur [0, 1[END_REF], and ( 35 Applying the Krein representation ( 43) in [START_REF] Maisonneuve | Exit systems[END_REF] yields

P 0 (T -G T ∈ du)/du = α R α (0, 0) e -αu ∞ u da ∞ 0 M (dz) e -az = α Φ(α) e -αu ∞ 0 M (dz) z e -uz = α Φ(α) ∞ 0 M (dz) z e -(α+z)u = ∞ 0 (α + z) e -(α+z)u M (dz), with M α (dz) = α Φ(α) M (dz) z(α + z) . ( 44
)
The claim of the theorem follows now from the fact that mixtures of exponential distributions are infinitely divisible. (see Bondesson [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF]). For D T -T we compute similarly from (34) via the Krein representation

P 0 (D T -T ∈ dv)/dv = α Φ(α) e αv ∞ v e -αa ν(a)da. = ∞ 0 z e -zv M (dz).
To analyze the distribution of ∆ T we use the Krein representation in [START_REF] Meyer | Processus de Poisson ponctuels[END_REF] to obtain

P 0 (∆ T ∈ da)/da = 1 Φ(α) ∞ 0 e -za -e -(α+z)a M (dz). ( 45
)
Notice that for a ≥ 0

f (a; z, α) = z(α + z) α e -za -e -(α+z)a
is a probability density as a function of a. In fact, letting T 1 and T 2 be two independent exponentially distributed random variables, with respective parameters z and α + z, then the sum T 1 + T 2 has the density f (a; z, α). In other words, the distribution of T 1 + T 2 is a so called gamma convolution.

Next we notice that letting Π z,α (dx) := z(α + z) α x -2 dx, z < x < α + z we may represent the distribution of T 1 + T 2 as a mixture of Gamma(2)distributions as follows

f (a; z, α) = ∞ 0 x 2 a e -xa Π z,α (dx). ( 46
)
Combining the representation ( 46) with (45) yields

P 0 (∆ T ∈ da)/da = α Φ(α) ∞ 0 f (a; z, α) z(α + z) M (dz) = ∞ 0 x 2 a e -xa Π α (dx), (47)
where Π α is a probability measure on R + given for any Borel set A in R + by

Π α (A) = ∞ 0 M α (dz)Π z,α (A). ( 48
)
The claim that the distribution of ∆ T is infinitely divisible follows now from ( 47) by evoking the result that mixtures of Gamma(2)-distributions are infinitely divisible (see Kristiansen [START_REF] Kristiansen | A proof of Steutel's conjecture[END_REF]).

Remark 17. 1. Recall from Bondesson [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF] that a probability distribution F on R + is called a generalized gamma convolution (GGC) if its Laplace transform can be written as

∞ 0 e -sa F (da) = exp -µs + ∞ 0 log t t + s U (dt) , (49)
where µ ≥ 0 and U is a measure on (0, ∞) satisfying

(0,1] | log t|U (dt) < ∞ and (1,∞) U (dt) t < ∞.
It is known see [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF] Theorem 4.1.1 p. 49 that if β is the total mass of U then the distribution F in ( 49) is a mixture of Gamma(β)-distributions 2. The distribution of the length ∆ t of an excursion straddling a fixed time t for a stationary diffusion (with stationary probability distribution) is given in Theorem 6 [START_REF] Hawkes | Statistics of local time and excursions for the Ornstein-Uhlenbeck process[END_REF] as

P(∆ t ∈ da) = a ν(a) m(R + ) da.
Also in this case the distribution of ∆ t is a mixture of Gamma(2)-distributions and, hence, it is infinitely divisible. In fact,

P(∆ t ∈ da)/da = ∞ 0 z 2 a e -za M (dz).
where the probability measure M is given in terms of the Krein measure M via M (dz) = M (dz)/(m(R + )z 2 ).

5 Case study: Ornstein-Uhlenbeck processes

In this section we give some explicit formulas for excursions from 0 to 0 associated with Ornstein-Uhlenbeck processes. It is possible to obtain such formulas due to the symmetry of the Ornstein-Uhlenbeck process around 0. Analogous results for excursions from an arbitrary point x to x are less tractable.

Basics

Let U denote the Ornstein-Uhlenbeck diffusion with parameter γ > 0, i.e., U is the solution of the SDE

dU t = dB t -γU t dt with U 0 = u,
and most of the time, but not always, we take u = 0. Recall that the speed measure and the scale function of U can be taken to be m(dx) := 2 e -γx 2 dx and S(x) :=

x 0 e γy 2 dy, respectively. Moreover, see [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF] p. 137, the Green kernel of Ornstein-Uhlenbeck process with respect to the speed measure is given for x ≥ y by

R λ (x, y) = Γ(λ/γ) 2 √ γπ exp γx 2 2 D -λ/γ x 2γ × exp γy 2 2 D -λ/γ -y 2γ ,
where D denotes the parabolic cylinder function. In particular, since

D -λ/γ (0) = √ π 2 λ/(2γ) Γ ((λ + γ)/(2γ)) -1
, we have, after some manipulations,

R λ (0, 0) = √ π Γ(λ/γ) 2 √ γ 2 λ/(2γ) Γ ((λ + γ)/(2γ)) -2
.

Consequently, using the formula

Γ(x) = 2 x-1 √ π Γ((x + 1)/2) Γ(x/2) we obtain R λ (0, 0) = 1 Φ(λ) = Γ(λ/(2γ)) 4 Γ((λ + γ)/(2γ)) . (50)
We remind also that U can be represented as the deterministic time change (Doob's transformation) of Brownian motion via

U t = e -γt (u + β at ) ,
where β is a standard Brownian motion and a t := (e 2γt -1)/2γ (see Doob [START_REF] Doob | The Brownian movement and stochastic equations[END_REF]).

Killed Ornstein-Uhlenbeck processes

We consider now the Ornstein-Uhlenbeck process killed at the first hitting time of 0, and denote this process by U . Let Y be the diffusion on R + satisfying the SDE

dY t = dB t + 1 Y t -γY t dt, Y 0 = y > 0.
Recall that Y may be described as the radial part of the three-dimensional Ornstein-Uhlenbeck process. In [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF] p. 138 the basic properties of such processes are presented. In particular, we record that 0 is an entrance-notexit boundary point and the process is positively recurrent its stationary distribution being the Maxwell distribution, i.e., the distribution with the density proportional to the speed measure of Y, that is, m Y (dx) := 2x 2 e -γx 2 dx, x > 0.

We remark that there is a misprint in [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF] and can be computed from the transition density of a Bessel process using Doob's transform (for an approach via inverting the Laplace transform see Giorno et al. [START_REF] Giorno | Some remarks on the Rayleigh process[END_REF]). In Salminen [START_REF] Salminen | On conditional Ornstein-Uhlenbeck processes[END_REF] it is proved that

P x ( U t ∈ dy) = P x (U t ∈ dy, t < H 0 ) = e -γt p Y (t; x, y) ϕ γ (y) ϕ γ (x) m Y (dy), (51)
where ϕ γ (x) = 1/x is the unique (up to multiplicative constants) decreasing positive solution of the ODE associated with Y killed at rate γ:

1 2 u (x) + 1 x -γx u (x) = γu(x).
From ( 51) we obtain Proposition 18. The transition density (with respect to its speed measure m) of the Ornstein-Uhlenbeck killed at the first hitting time of 0 is given by

p(t; x, y) = √ γ e γt/2 2πsinh(γt) exp - γe -γt (x 2 + y 2 ) 2sinh(γt) sinh γxy sinh(γt) . (52)
Combining the expression of the transition density in [START_REF] Williams | Path decompositions and continuity of local time for onedimensional diffusions[END_REF] with formula (2) yields the distribution of H 0 (see also Sato [START_REF] Sato | Evaluation of the first passage time probability to a square root boundary for the Wiener process[END_REF] and Going-Jaeschke and Yor [START_REF] Göing-Jaeschke | A clarification note about hitting times densities for Ornstein-Uhlenbeck processes[END_REF]).

Proposition 19. The density of the first hitting time of 0 for the Ornstein-Uhlenbeck process {U t } is given by The density of the Lévy measure of the inverse local time at 0 is obtained by applying formula [START_REF] Erdélyi | Higher transcendental functions[END_REF] (see also Hawkes and Truman [20]). Moreover, using [START_REF] Sato | Evaluation of the first passage time probability to a square root boundary for the Wiener process[END_REF] in formula (4) leads to an explicit expression for the Bernstein function associated with the inverse local time at 0. Next we display the distributions of ∆ T , T -G T , and D T -T. Recall that these distributions are infinitely divisible and the densities are expressable via the density of the Lévy measure, as stated in Corollary 13 formulae [START_REF] Maisonneuve | Exit systems[END_REF] and [START_REF] Meleard | Application du calcul stochastique à l'etude de processus de Markov réguliers sur [0, 1[END_REF], and in Theorem 16. To simply the notation, we take γ = 1. Therefore, it seems motivated to compute the measure M (cf. [START_REF] Revuz | Continuous Martingales and Brownian Motion[END_REF]) in this representation for Ornstein-Uhlenbeck processes.

f x0 (t) = γ 3/2 x e γt/2 √ 2π(sinh(γt)) 3/2 exp - γe -γt x 2 2sinh(γt) . ( 53
To start with, we give the spectral representation of the transition density of p of the Ornstein-Uhlenbeck process killed at the first hitting time of 0. Instead of computing from scratch, we exploit the spectral representation for p Y (with γ = 1) as presented in Karlin and Taylor [START_REF] Karlin | A second course in stochastic processes[END_REF] p. 333:

p Y (t; x, y) = ∞ n=0 w -1 n,1/2 e -2nt L (1/2) n (x 2 ) L (1/2) n (y 2 ), (58)
where {L

The normalization (59) coincides with the normalization in Erdelyi et al. [START_REF] Erdélyi | Higher transcendental functions[END_REF] (see formula (2) p. 188 where the notation for the norm is h n ). Therefore, from [START_REF] Erdélyi | Higher transcendental functions[END_REF] formula [START_REF] Getoor | Excursions of a Markov process[END_REF] p. 189 we have

L (1/2) n (0) = n + 1 2 n (61)
and, consequently (cf. ( 53)), we obtain the spectral representation for the density of the first hitting time of 0

f x0 (t) = ∞ n=0 w -1 n,1/2 e -(2n+1)t xL (1/2) n (x 2 ) L (1/2) n (0). = 2 √ π ∞ n=0 e -(2n+1)t xL (1/2) n (x 2 ). ( 62
)
To find the spectral representation for the density of the Lévy measure we apply formula (54) which yields

ν(t) = 2 √ π ∞ n=0 n + 1 2 n e -(2n+1)t . (63)
In view of (43), we have Proposition 22. The measure M in the Krein representation of ν for the Ornstein-Uhlenbeck process is given by

M (dz) = 2 √ π ∞ n=0 n + 1 2 n δ {2n+1} (dz),
where δ is the Dirac measure.

Notice that ν(t) = 2 √ π e -t 1 -e -2t -3/2 , and, hence, (63) is obtained also from the MacLaurin expansion of x → (1 -x) -3/2 evaluated at x = e -2t .

p. 35 . 5 .

 355 Proposition Let for a given x > 0 Λ x := sup{t : ω(t) = x} denote the last exit time from x. Then

G

  T := sup{s ≤ T : X s = 0}, D T := inf{s ≥ T : X s = 0}, and ∆ T := D T -G T . The Lévy exponent of the inverse local time at 0 is denoted by Φ, in other words, E 0 (exp(-λτ )) = exp (-Φ(λ))Recall the relation (cf. (4) with y = 0)Φ(λ) R λ (0, 0) = 1. (25)4.1 Last exit decomposition at TWe begin by discussing the last exit decomposition at the exponential time T.Theorem 9. (i) The processes{X u : u ≤ G T } and {X G T +v : v ≤ ∆ T } are independent.(ii) The law of {X u : u ≤ G T } may be described as follows: (a) L T := L G T is exponentially distributed with mean 1/Φ(α), (b) The process {X u : u ≤ G T } conditionally on L T = is distributed as {X u : u ≤ τ } under the probability exp (-α τ + Φ(α)) P 0 .

4. 2 Proposition 12 .

 212 On the distribution of (G T , D T ) In this section the distributions of T -G T , D T -T and ∆ T := D T -G T are studied in detail. The joint distribution of T -G T and D T -T is given by

Theorem 16 .

 16 The distributions of T -G T , D T -T and ∆ T are infinitely divisible.

  ), the intrinsic term in the densities of T -G T , D T -T and ∆ T is the density ν(a) of the Lévy measure of the inverse local time at 0. We consider first the distribution of T -G T .

) 5 . 3

 53 Lévy measure of inverse local time and densities of ∆ T , T -G T and D T -T

Proposition 20 .

 20 The density of the Lévy measure of the inverse local time at Let {τ : ≥ 0} be the inverse local time at 0. Then E 0 (exp(-λτ )) = exp -4 Γ ((λ + γ)/2γ) Γ(λ/2γ) .

Proposition 21 .√ π e 2a e 2a - 1 - 3 / 2 ,du e -αu 2 √ π e 2u e 2u - 1 - 3 / 2 . (57) 5 . 4

 2113213254 With Φ(α) as in[START_REF] Sato | Evaluation of the first passage time probability to a square root boundary for the Wiener process[END_REF], the distributions of ∆ T , T -G T and D T -T are given, respectively, byP 0 (∆ T ∈ da)/da = 1 -e -αa Φ(α)2 (55)P 0 (T -G T ∈ da)/da = α e -αa Φ(α) 2 √ π e 2a -1 -1/2 , (56)andP 0 (D T -T ∈ da)/da = α e αa Φ(α)∞ a The Krein measure As seen in Section 4.3, the Krein representation plays a central rôle in the proof of infinite divisibility of the distributions of T -G T , D T -T, and ∆ T .

:

  n = 0, 1, 2, . . . } is the famly of Laguerre polynomials with parameter 1consider the symmetric density with respect to the speed measure m Y . From (51) and (58) the spectral representation of p is now obtained immediately and is given by p(t; x, y) =

  p. 139; the stationary distribution in the general case is not a χ 2 -distribution. The transition density of Y with respect to its speed measure m Y is

	p Y (t; x, y) =	√ γ e 3γt/2 2π sinh(γt) xy	exp -	γe -γt (x 2 + y 2 ) 2 sinh(γt)	sinh	γxy sinh(γt)

Acknowledgement. We thank Lennart Bondesson for co-operation concerning gamma convolutions.


G Vermot Des Roches

O Chiello

E Balmes

X Lorang

Benchmarking Signorini and exponential contact laws for an industrial train brake squeal application

Contact and friction play a key role in many industrial applications and their modeling is becoming accessible in vibration applications. Two strategies are commonly considered: idealized or perfect contact which assumes no interpenetration of structures at contact points and functional representations that use pressurepenetration relationships. The two strategies imply widely different solvers so that the choice between the two is almost always based on a priori decisions. For an application to the prediction of unstable modes in railway brake squeal, the objective of this paper is to objectively compare results obtained with the two models and to show how this comparison can be used to determine parameters of functional laws giving the same results as idealized contact.

Introduction

Simulating interactions of multiple components in mechanical assemblies is critical in many industrial applications. This requires contact and possibly friction models which have become accessible to simulation due to advances in solvers and computational power. Train brake squeal during station parking operations generates noise disturbances that can raise issues in the vicinity of stations potentially limiting train exploitation. The AcouFren project managed by SNCF (the french railway company) aims at providing, through simulation, noise indicators at early brake design stages.

Contact modeling strategies are widely documented in the literature, from which two trends are clearly found. The oldest representation is due to Signorini from the 30's, it idealizes a perfect contact between two ideally smooth bodies, and constrains the relative displacement of contact points. The second representation strategy is functional and authorizes a level of interpenetration between bodies. First measurements were performed in the 60's, from which Greenwood and Williamson [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF] developed models where contact stiffness is function of the body interpenetration. State-of-the-art measurements as performed by Nogueira et al. [START_REF] Nogueira | Experimental validations of elastic to plastic asperitybased models using normal indentations of rough surfaces[END_REF] still show that this model can be relevant, mostly for hard surfaces. Direct ultrasonic measurement of contact stiffnesses can also be performed. Biwa et al. [START_REF] Biwa | Stiffness evaluation of contacting surfaces by bulk and interface waves[END_REF] thus showed that contact stiffness is function of the contact pressure, mostly for low interfacial loading.

Physically speaking, contact between two bodies can be interpreted as the balance between their respective stiffness and asperity compression, as illustrated in figure 1. When one material is much softer than the other, like for brake systems, local compression of asperities becomes non-negligible at the macroscopic scale, generating local interpenetration of the nominal surfaces (with sometimes plastic deformation of asperity). It can then be relevant to identify a law p(g) linking penetration of nominal surfaces and contact pressure.

Figure 1: Representation of contact between two bodies, body stiffness and asperities Numerically speaking, implementation of functional contact laws using continuous non linear functions of displacement is easier since classical non-linear algorithms can then be used. The ideal contact law does not feature any mechanical parameter to be identified, but its numerical implementation requires specific algorithms. Their convergence is often very sensitive and fine tuning of numerical parameters is needed. Besides, numerical convergence is often obtained with a tolerance, so that interpenetration can still occur in the results, with levels that are difficult to control.

Choosing a contact strategy is thus difficult and the general feeling is that these different contact models yield different results. The aim of this paper is to present an objective comparison of Signorini and exponential functional contact laws in application to an industrial train brake squeal simulation. Section 2 starts by presenting numerical implementation details of functional contact, based on the work of Vermot des Roches [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF]. Signorini implementation is then presented, based on the work of Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF]. Section 3 then demonstrates that an equivalence between both contact laws is possible, and a functional law of the exponential type is proposed. Section 4 eventually presents a validation applied to industrial railway brake models.

Contact formulations for squeal simulations

Squeal simulations in the frequency domain is nowadays widely deployed in industry, and is based on the sliding perturbation method. This was used by Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF], Lorang [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant: Application au crissement des freins de TGV[END_REF] for train brake squeal applications, or by Vola et al. [START_REF] Vola | Friction and instability of steady sliding: squeal of a rubber/glass contact[END_REF] to study rubber/glass instabilities in sliding steady states, and by Vermot des Roches in [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF] for automotive brake squeal application. Contact formulations and adapted numerical schemes for squeal simulations are here presented.

Contact handling for finite element models

The representation of contact between two solids as formalized by Signorini is defined by a contact direction N along the outward normal at a contact point of one of the bodies (the chosen one, S1, is commonly called master, as contact forces will be evaluated on its surface). Figure 2 illustrates the contact configuration. The distance between the two bodies along the contact normal is called the gap, here noted g. The convention is that the gap is positive when contact is open, or not effective. When contact is closed, or effective, a contact pressure p exists between the bodies, applied at contact points, equal in amplitude on both solids. The contact law proposed by Signorini, illustrated in figure 2, defines the condition of non-interpenetration as

     g ≥ 0 p ≥ 0 (g).(p) = 0 (1)
Quite differently, functional contact laws define a pressure-gap relationship. These can represent physical phenomena due to interaction between asperities, but have also been considered as a numerical regularization trick. Both approaches have their limitations, so that in the end the most convenient to reproduce the behavior of interest should be retained. Classical laws are illustrated in figure 3. Numerically the gap is a linear observation of relative displacement between two surfaces along a normal, and can thus be written using an observation matrix [C N OR ]. Noting N the number of system DOF, and N c the number of contact points, it is written with a possible gap offset {g 0 } as

{g} Nc×1 = [C N OR ] Nc×N {q(u)} N ×1 -{g 0 } Nc×1 (2)
The force resulting from the gap-pressure relationship is then defined as an integral over the contact surface

{q} T {f N } = Γ û(q)N p(q)dS e j {û(q)} T {N } p(x j , q)ω (e) j J (e) (x j ) (3)
where f N is the global contact force, p the contact pressure, q a virtual displacement, q the displacement, x (e) j are the integration points of current element e, J (e) (x j ) the Jacobian of the shape transformation (surface associated to each integration point) and ω (e) j the weighting associated with the integration rule for element e.

In practice generalized contact forces can be recovered at DOF using the transpose of the gap observation (2):

{f N } N ×1 = [C N OR ] T N ×Nc ω (e) j J (e) (x j ) p(x j , q) Nc×1 (4)
The Coulomb law is expressed at each contact point j as

|f T j | < µ|p(x j )| ⇔ w j = 0 |f T j | = µ|p(x j )| ⇔ f T j = µ|p(x j )| w j |w j | (5)
where f T j is the friction constraint, µ the friction coefficient, and w j the sliding velocity. Computation of the sliding velocity requires computation of differential velocities between slave and master surfaces in the plane orthogonal to the contact normal. A tangential displacement observation matrix [C T AN ] can be used

{w} 2Nc×1 = [C T AN ] (2Nc×N ) { q} N ×1 = [C T AN ] 1 [C T AN ] 2 { q} N ×1 (6)
The tangential displacement is observed in the two directions of the plane defined by the contact normal. It can be useful to distinguish observations in both directions when a natural orientation of contact forces exist. This is for example the case for this paper application as a steady sliding state exists.

Friction generalized forces can be recovered, from the local friction constraint f T j by

{f T } N ×1 = [C T AN ] T N ×2Nc ω (e) j J (e) (x j ) f T j (p(x j , q)) 2Nc×1 (7)
Given a multibody system with contact interactions, the system must then verify

[M el ] {q} + [C el ] { q} + [K el ] {q} = {f ext } + {f N L } (8)
where [M el ], [C el ], [K el ] are respectively the system elastic mass, damping and stiffness matrices, {q} the displacement vector. The non linear force vector can be decomposed as the sum of contact {f N } and friction forces {f T }, {f N L } = {f N } + {f T }.

Brake squeal simulations are based on a steady sliding state, considering that all points in effective contact are sliding, thus verifying the second equation of the Coulomb law [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF]. A tangent model can then be derived from the steady sliding state to obtain the system stability from complex mode analysis. The same hypothesis of constant sliding for points in effective contact is used, hence the name of sliding perturbations.

Due to the priviledged sliding direction of the steady sliding state, observation matrices [C T AN ] i are constructed for modal analysis so that [C T AN ] 1 is along the direction of the sliding velocity of the steady sliding state at each contact point, and [C T AN ] 2 is along the orthogonal direction.

For 3D solids, Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF] noted that the sliding friction tangent state also features a damping term due to the possibility for the friction force to change its orientation in vibration. A planar friction damping matrix

[C f 2 ]
derived from the expression of friction forces thus appears, written

δf T = [C f 2 ] δ q = [C T AN ] T 2 \ µ |w j | ω (e) j J (e) (x j )p(x j , q) \ [C T AN ] 2 δ q (9)
Physically, the higher the contact pressure, the more difficult the sliding direction variation. This is due to the priviledged sliding direction of the steady sliding state.

Functional contact

For functional contact formulations, classical non-linear methods can be applied, since contact-friction forces only depend on the system state. Equation ( 8) can then be writen using {f N L } = {f N L (q)}.

Functional contact statics

The Newton method, presented in figure 4, is the most classical method to resolve regular non linear statics. Difficulties in Newton resolution algorithms are mostly on the numerical side. Although mathematical convergence is obtained theoretically when the system is regular enough, hard laws like contact implementations are challenging.

Initial state : q 0

Residue computation The system Jacobian [J] derived from [START_REF] Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF] in statics can be written

{r k } = -[K el ] {q} k + {f ext } + f N L ((q) k ) Correction [J] {∆q} = {r k } {q} k+1 = {q} k -{∆q} Convergence r k fext ≤ no
[J] = [K el ] + [K c (g)] (10)
Considering a steady sliding state, the tangent friction perturbation is neglected. Noting ω j J(x j ) the surface associated with each contact point. the tangent contact state defining [K c (g)] writes

∂f N (q + δq) ∂δq = [C N OR ] T \ ω j J(x j ) \ ∂p(g + δg) ∂δg [C N OR ] = [C N OR ] T \ k cj (g) \ [C N OR ] (11)
For hard laws, tangent state variation can be very brutal, which is a drawback of such algorithms. Indeed, when contact opening occurs in a loaded area at iteration k, tangent contact state becomes locally null, and correction can become very large (with possibly contact penetrations). Divergence can then occur quickly as the non-linear forces returned for significant penetrations lead to numerically infinite forces.

A full set of adaptations can be applied to the original Newton scheme, yielding numerous variations of quasi-Newton methods. For hard contact, cut-backs on corrections based on the norm of the mechanical residue, and controlled release on the contact part of the Jacobian for opening points are usually working and have been used here. Implementation details is key to performance, but are outside the scope of this paper.

Functional contact sliding perturbation formulation

Functional contact formulations are defined as function of the gap, such that a tangent state can be derived for each contact point by deriving the gap/pressure relation, illustrated in figure 3. A key aspect to note is that for similar contact pressure levels, different functional contact laws can yield very different contact stiffnesses. Contact law calibration in statics and dynamics can thus be very different.

Tangent contact stiffness has been defined in [START_REF] Balmes | Vermot des Roches[END_REF]. Defining tangent friction states is less direct as it couples normal and tangential directions, which is typically non-symmetric. Friction expression ( 5) indeed shows that a variation of contact force has a direct effect on friction force, whereas friction force variation may occur without effect on contact force.

For sliding states, friction force is explicitly defined and only depends on contact force at the same point.

The tangent friction stiffness is then the tangent contact stiffness [START_REF] Balmes | Vermot des Roches[END_REF] scaled by the friction coefficient µ, and transfers normal displacement to planar force in the sliding direction. Since planar displacements are free for a fixed contact state, this effect is not considered for real mode computation. The tangent sliding friction coupling stiffness K slide nlf is thus defined as

K slide nlf = 0 for real mode computation [C T AN ] T 1 \ µk cj (g) \ [C N OR ] for complex mode computation (12)
Equation ( 11) illustrates one of the main interest of functional contact laws, which is the representaton of a contact stiffness distribution depending on contact pressure. Such variation can be non negligible for low pressures [START_REF] Biwa | Stiffness evaluation of contacting surfaces by bulk and interface waves[END_REF]. This is impossible with Signorini contact, where this pattern can only be assessed by evaluating contact opening thresholds as function of complex mode amplitudes.

Signorini, or idealized contact

Idealized contact verifying the exact Signorini condition ( 1) is non regular, hence the common name of nonsmooth dynamics. {f N L } has no specific link to the displacement, and the pair ({q} , {f N L }) has to be found using appropriate solvers.

Signorini statics

Solution strategies are usually based on Lagrange resolutions. The non-symmetric nature of friction requires a methodology refinement, leading to Augmented-Lagrange methods. De Saxcé and Feng developed the bipotential method [START_REF] Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF]. Alart et al. presented a generalization of the Newton method adapted to Lagrange formulations introducing Gauss-Siedel like algorithms [START_REF] Jourdan | A gauss-siedel like algorithm to solve frictional contact problems[END_REF].

The solution chosen here was developped by Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF]. It is based on the status assumption of each contact points (open/closed gap). Closed contact points generate a displacement constraint under which resolution is linear. Noting with superscript (e) the observation restricted to conctact points in effectively closed contact, one resolves the problem in the orthogonal subspace T (e) that verifies the constraint

{q} = T (e) {q e } with [C N OR ] (e) T (e) = 0 (13)
The subspace T (e) is unique, any {q e } is defined by the basis chosen for T (e) which may be arbitrary. Equation ( 8) can thus be rewritten in statics as

T (e) T [K el ] T (e) {q e } = T (e) T ({f ext } + {f N } + {f T }) (14)
Working with pressures is mandatory, one thus introduces an operator passing from generalized forces at DOF to local pressures at contact points [C N P ] (based on the pseudo inversion of relation ( 4)), and to friction local constraints [C T P ] (pseudo inverse operator of ( 7)). Contact forces verify the equilibrium and are opposed in direction and of equal norm at each observation point. Contact pressures are then expressed as {p(

x j )} = [C N P ] ([K el ] {q} -{f ext }) (15)
and are null in T (e) . Sliding friction forces [START_REF] Vola | Friction and instability of steady sliding: squeal of a rubber/glass contact[END_REF] can thus be written (with -1 denoting a pseudo inverse)

{f T } = µ [C T P ] -1 [C N P ] ([K el ] {q} -{f ext }) (16)
Equation ( 14) projected in T (e) becomes after simplification

T (e) T [I] -µ [C T P ] -1 [C N P ] [K el ] T (e) {q} = T (e) T [I] -µ [C T P ] -1 [C N P ] {f ext } (17)
Assuming contact points at nodes and first order contact elements between compatible meshes, one can directly work using contact resultants. Matrices [C N P ] and [C T P ] are then equivalent to observation matrices [C N OR ] and [C T AN ] of equations ( 2) and ( 6), generated using integration rules at nodes. This simplification of implementation was used in this paper application.

It must be noted that the factor [I] -µ [C T P ] -1 [C N P ] is not symmetric. Direct resolution of (17) can then be challenging for large problems, like the million-DOF system presented in section 4. Use of the latest PARDISO librairies [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF], implemented in SDT [START_REF] Balmes | Vermot des Roches[END_REF] was required here.

Contact status must be predicted, and can be updated depending on the result of (17). This is formalized by considering the group of contact nodes in effective contact in (17) at resolution iteration k, noted C k .

At iteration k, contact points violating Signorini conditions in the solution of (17) are classed in group D k such that

     x j ∈ C k , p(x j ) k < 0 or x j / ∈ C k , g(x j ) k < 0 (18)
Status updating allows forming a new [C N OR ] (ek) and consists in switching the status of points in D k ,

C k+1 = (C k -C k ∩ D k ) ∪ (D k -D k ∩ C k ) (19)
Iterative status correction can be resolved by an Uzawa algorithm, presented in figure 5. This method convergence has been proved without friction and by switching contact statuses one by one. Its performance is however very satisfying in the application presented, with very few iterations even for complex systems and high friction coefficients.

Prediction : Initial contact status C 0 Mechanical resolution with C k [K el ] {q} k = {f ext } -{f N L } [C N OR ] (ek) {q} k = 0 Check solution      {p(x j )} k = [C N P ] [K el ] {q} k -{f ext } {g(x j )} k = [C N OR ] {q} k GenerateD k Convergence D k = ∅? no, C k+1

Signorini modal analysis

Following the exact contact condition of Signorini, the system tangent state only depends on the contact surface, materialised by C ∞ , as presented in 2.3.1. The mechanical problem is here formulated as

λ 2 [M el ] + λ ([C el ] + [C f 2 ]) + [K el ] {q} = {f N L } [C N OR ] (e∞) {q} = 0 (20)
As for static resolution, projection in the subspace T (e∞) verifying the contact constraints yields

λ 2 M (e) µ + λ C (e) µ + [C f 2 ] + K (e) µ {q e } = 0 (21)
where matrices X (e) µ are the so-called sliding matrices, defined as

X (e) µ = I -µ [C T P ] -1 [C N P ] T (e∞) T [X el ] T (e∞) (22)

Definition of an exponential contact law

The two contact models presented in section 2 are physically limited, so that the choice of one or the other should be based on the ability to reproduce specific behavior and on the performance of algorithms used to resolve equations. The objective of this section is to demonstrate how the ability to reproduce specific modal behavior can be used to calibrate the parameters of an exponential pressure/gap law that reproduces the idealized behavior.

The clamped/sliding block

Railway industry pads can be of different shapes and materials, the sample pad shown in figure 6 is used as an example here. A single friction block (made of friction material), solidar to a metallic ring making the link to the backplate, is extracted for the illustration. The friction material has a non-negligible loss-factor. Here, this damping is modelled by a material Rayleigh damping for easy compatibility for transient simulations, even though better models [START_REF] Bianchi | Using modal damping for full model transient analysis[END_REF] are probably needed. When computing steady sliding states and sliding perturbations, the friction blocks effective contact area is supposed to be fully sliding (although some contacts can be opened). A sample demonstration of the blocks behavior is thus to consider them in clamped/sliding condition. Ideal contact condition is materialized in this case by clamping normal displacements of the friction block bottom side. Figure 6 shows some modes associated to the pad block in such conditions.

To demonstrate the effect of contact stiffness (possibly derived from a functional contact law), normal spring elements with variable stiffness are added to the blocks bottoms instead of clamping. The previous case would then correspond to a spring with infinite stiffness.

Figure 7 presents complex mode frequency and damping results of the clamped sliding block as function of contact stiffness. An S-shaped frequency evolution can be observed. Saturation occurs around 10 4 M P a, after which frequencies are equivalent to an ideal contact implementation. Mode damping evolution as function of contact stiffness can be observed in figure 7. Damping rates increase with frequency as a consequence of the stiffness proportionality coefficient used to form the material Rayleigh damping matrix.

For low contact stiffnesses, ground coupling is small and the friction block behaves as free, that-is-to-say the whole strain energy is located in the block. Damping ratio of the clamped/sliding block is thus linked to the damping ratio of the block itself. The same observation can be done for high contact stiffnesses, when coupling is strong and the friction block behaves as clamped. The expected damping ratio increase is obtained regarding the frequency increase, materialized by the dotted line in figure 7right.

In the curve inflection area, a non negligible strain energy is located in the coupling springs, so that a non negligible part of the total system strain energy is outside the friction block. The damping ratio thus varies non linearly, and can even decrease when coupling becomes non negligible while still being low. This observation is very general regarding impact of components strain energy distribution in damping of assemblies. More detailed illustrations of this effect can be found in [START_REF] Bianchi | Using modal damping for full model transient analysis[END_REF].

Exponential contact law definition

As a large variety of functional contact laws could be calibrated, one decided here to focus on the exponential type. This choice is based on the expected transient behavior, where contact opening transitions and penetration saturation effects also need to be properly represented for a potentially large range of pressures. Numerical results obtained in [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF] with such law were satisfying.

Given {g} the gap vector (opposite of the penetration), obtained from an observation of the relative displacement of both surfaces along the local contact normal, pressure {p} is defined as

{p} = p 0 e -λ{g} (23)
The exponential law has two parameters p 0 and λ, whose effect on contact behavior is illustrated in figure 8.

To interpret these curves, values have to be observed as function of a reference force level, conceptually named F . Horizontal intersections are thus iso-values regarding system loading outside the contact area.

Parameter p 0 only drives the curve offset, thus the level of gap obtained for a given force. Identification of this parameter can be obtained if measurements of surface rougthness are accessible, providing acceptable values of interpenetrations. In this paper, coefficients are only numerically calibrated so p 0 can be set arbitrarily. Choosing a p 0 such that significant contact pressures are obtained for closed (or negative) gaps, simplifies implementation and allows easier physical interpretations. Parameter λ impacts contact stiffness at a given force level, the higher, the stiffer. This parameter itself can thus be updated to choose a contact stiffness that complies with the fact that the frequencies in clamped/sliding conditions are converged towards the perfect contact behavior. Looking at the diagrams of figure 7, target contact stiffnesses values can be pin-pointed in figure 9. For calibration, the level of contact pressure generated by the system must be evaluated, either by using a status implementation or by choosing a linear contact law for a preliminary computation, or by exploiting the system command (this was done here). Noting p 1 the reference pressure, contact pressure and contact stiffness k c (g) (first order derivative of contact pressure) expressions yield

p 1 = p 0 e -λg k c (g) = λp 0 e -λg (24)
The product λg can be substituted in the first equation of (24) so that

λ = k c (g) p 1 (25)
Choosing a target value of k = 10 7 M P a in figure 9, and knowing that the contact force command on the TGV brake system is here of 5kN over a surface of 21.4 10 3 mm 2 , one obtains λ = 4.28 10 7 mm -1 . Higher values could also be used for security margin regarding system variability. It must however be kept in mind that values with too high stiffness may alter numerical conditioning of the stability problem.

For comparison purposes, a converged value kc = 10 5.3 M P a at the limit, and a non converged value kc = 10 4 M P a will also be tested. Correspondances between λ and target stiffnesses are reported in table 3.2.

Application to industrial railway brakes

The french high speed train, TGV, features 4 disc brake systems on each axle, as presented in figure 10. A single brake system is here modelled by finite elements on which subassemblies (the disc, rig and pads) have been updated as part of the AcouFren project. The full model is free-meshed using second order ten nodes tetraedrons, yielding between 500,000 and 1,000,000 DOF depending on the fitted pad model.

The mesh is kinematically positioned depending on the pad thickness to have zero gap at rest. Disc sections underlying the pads are remeshed to obtain compatible contact interfaces. Braking force is applied by an actuator at the system rear, which translates the pad holder towards the disc. The translation is achieved kinematically by the rotations of levers and rods, constituting the brake rig. The illustration will focus on complex mode evaluation for station parking operations. Comparisons between calibrated exponential laws, defined in section 3.2 and the idealized contact implementation are presented for steady sliding solutions in section 4.1 and complex modes in section 4.2.

Steady sliding simulations

A contact force command of 5kN is applied by the actuator. Static results are compared between Signorini contact and the family of exponential contact laws defined in table 3.2. In the following, contact nodes define nodes of the mesh in the contact area, while Gauss contact points define the Gauss contact integration points used for functional contact implementation.

Figure 12 presents contact results in terms of pressure and gaps. Points with gaps opened at contact nodes are presented in figure 12a. It can be seen that no global difference rises for these displacements. Although contact pressure from an exponential law is strictly positive, exponential decay as function of opening is sufficient not to alter the behavior of unloaded contact points. Figure 12c presents contact pressure at contact nodes. Contact pressures at the interface are here computed a posteriori from displacements to allow comparison between implementations, using a nodal resultant to nodal pressure operator ([C N P ] in (15), integrated at Gauss points). Opened contact points will thus show negative pressures at the interface, which would correspond to depressions on a membrane. Such observation shows that pressure levels seen by the structure outside the interface are very similar. Forces transmitted by the interface are naturally recovered in all cases, which is physically natural.

Global results in figure 12c can be put into perspective with the local results of figure 12d. Contact pressures are here plotted at the Gauss contact points for the functional laws only. Clear differences occur for low pressure areas, with a convergence pattern as function of λ. Such fluctuations for very low pressures are however not percieved by the structure.

Figure 13 presents absolute differences in displacement between the Signorini response and functional contact laws. It can be seen that for all cases, maximum differences scale in nanometers are a thousand times smaller than displacement scale.

Figure 13: Absolute displacement differences in nanometers between Signorini and functional static results, from top to bottom along the x, y and z axis, from left to right in increasing λ. Strain energy densities of differential displacements are presented for elements with highest levels only Orientation axis references are provided in figure 10. Convergence in the z axis is direclty linked to gap observations of figure 12b, and is clearly established as function of λ. Displacements in the x, y plane are linked to the rig deformation, depending on the pad holder displacement. Larger differences can here occur, with very limited strain energy errors for converged cases. Although static displacement seems better for the middle law, contact stiffness distribution is not as good, as will be discussed in section 4.2.

Sliding perturbation results

Complex modes are directly computed based on the formulations of section 2. Stability diagrams are presented in figure 14 for 500 modes, providing a frequency bandwidth of 0 -10kHz. The global stability diagram of figure 14b shows the effect of material Rayleigh damping, with very large damping ratios on the first modes (mass factor) and an average damping ratio linear increase with frequency (stiffness factor).

The first four most unstable modes are plotted in figure 14a and14c. Modes 27 and 69 are pad/disc coupling modes, with effective displacement of the pad holder. Modes 81 and 59 are planar modes, instability coming from friction block modes (possibly coupled with in-plane disc modes). No unstable modes of the global structure were found for high frequencies in this configuration -Rayleigh damping coefficients still needed experimental fitting in this study's model.

Globally, stability diagrams fit relatively well for all exponential calibrations. Larger differences are obvious for the soft value, highlighting a convergence as function of λ. Zooms in are provided in figures 14d, 14e and 14f. Structure modes (with clear displacement of the brake rig) are in the 2kHz range, and show correct fittings, although the soft law presents some peculiarities.

Modes in the 5.5kHz and 9kHz ranges mostly feature friction block modes. Slight differences between converged λ values and Signorini are mostly due to contact stiffness variation of low contact pressure areas.

A near perfect match stability diagram is here only obtained for the hardest value of λ.

Comparison between unstable mode shapes with MAC is finally presented in figures 14g,h,i. Due to the large number of modes computed and model size, each modal basis weights over 3 GBytes, making direct comparisons difficult. Unstable modes are of higher interest due to their propensity to generate squeal, hence the choice to compare only these shapes, with a detection threshold set under a damping ratio of -10 -3 %.

Convergence as function of λ is well observed regarding the number of detected unstable modes and their shape correlation, figure 14i displays a squared matrix with unit diagonal. More unstable modes with no equivalent in the Signorini results are found for the two softest λ values.

Conclusion

Simulation of structure-structure interaction is important in many industrial applications, requirering performant implementation of relevant contact models. Physically, transmission of contact forces between structures should naturally show correct contact pressures independently of the contact law.

Idealized contact avoids identification of physical parameters, but also lacks representation of contact stiffness variation as function of contact pressure. Functional laws require parameters, which can be very difficult to obtain experimentally. An exponential law was retained for its advantage of a quick decay of contact forces for opened gaps, while keeping relevant saturation patterns. For this law it was shown here that numerical calibration can be performed by considering the apparent stiffness at interfaces.

For the computation of complex modes of squealing brakes using the sliding perturbation framework, ideal and exponential contact were then shown, for properly calibrated parameters, to give nearly identical results in pressures and predicted modes. With two models leading to equivalent results, choices can then focus on implementation performance. An ongoing extension of the work presented here is the use of this benchmark to analyze transient simulations using methods presented in [START_REF] Vermot | Frequency and time simulation of squeal instabilities[END_REF].

Figure 2 :

 2 Figure 2: Contact normal (N) and gap (g) definition between solids, Signorini and Coulomb laws

Figure 3 :

 3 Figure 3: Classical contact laws, and their derivatives (for contact stiffness)

Figure 4 :

 4 Figure 4: Implementation of a static Newton resolution scheme

Figure 5 :

 5 Figure 5: Implementation of an Uzawa resolution algorithm for the sliding Signorini-Coulomb problem using status

Figure 6 :

 6 Figure 6: Left: A sample railway braking pad. The friction blocks in orange are linked to the backplate in green, using rings in blue. A dovetail joint (in red) fixes the pad in the brake rig. Right: Sample block modes in clamped/sliding conditions with ideal contact for a ring fitted block

Figure 7 :

 7 Figure 7: Frequency (left) and damping (middle) evolution as function of the contact stiffness applied to the lining block of the brake pad. Right: zoom in on the first modes.

  (a) p0 = 10 4 , λ ∈ (10 2 ; 10 10 )m -1 (b) λ = 1.3.10 9 , p0 ∈ (10 2 ; 10 7 )P a

Figure 8 :

 8 Figure 8: Exponential contact law evolution as function of its parameters

  Table 3.2 : k, λ correspondances

Figure 9 :

 9 Figure 9: Identification of target contact stiffnesses of the clamped/sliding block, and identified λ

  Figure 10: The TGV brake system

Figure 11

 11 presents global steady sliding solution with ideal contact implementation. Levels of displacement and contact force resultants are presented. Largest displacements are seen by the lever with values around 140µm. Contact force distribution greatly varies over the contact surface. Friction blocks of the front end are the most loaded, unloading of the trailing end seems complicated and is function of the pad holder fixations.

Figure 11 :

 11 Figure 11: Steady sliding state absolute displacement of a TGV brake using Signorini contact implementation, contact resultant at nodes of the disc top side

Figure 12b presents gapFigure 12 :

 12 Figure12bpresents gap results for closed contact. Great differences can here be observed, the softer the law, the deeper the penetration. Penetration levels seem reasonable compared to global displacements for converged λ values only.

(a) Most unstable modes 1 and 2 (

 2 b) Global stability diagram (c) Most unstable modes 3 and 4 (d) Zoom around 2kHz (e) Zoom around 5.5kHz (f) Zoom around 9kHz (g) λ = 4.83.10 4 (h) λ = 9.64.10 5 (i) λ = 4.28.10 7

Figure 14 :

 14 Figure 14: Stability diagrams of Signorini solutions and functional law solutions. Most unstable TGV modes, zooms around frequency ranges of interest and unstable modes (ζ ≤ -10 -3 %) MAC between Signorini and exponential laws solutions

Acknowledgements

The authors wish to thank the members of the AcouFren project, SNCF, IFSTTAR-LTE, SDTOOLS, VI-BRATEC, ENPC, ECL-LTDS, ALSTOM TRANSPORT, BOMBARDIER, FAIVELEY, and funder ADEME under the convention 0966C0281.