Tomasz Byczkowski
email: tomasz.byczkowski@pwr.wroc.pl

JACEK Jakub Chorowski
email: jakub.chorowski@gmail.com

Piotr Graczyk
email: piotr.graczyk@univ-angers.fr

Jacek Malecki
email: jacek.malecki@pwr.wroc.pl

Jacek Ma

M A Lecki

Harmonic measures

Keywords: 2010 MS Classification: Primary 60J45; Secondary 60G15, 60G40 harmonic measure, Ornstein-Uhlenbeck diffusion, Girsanov theorem, hyperbolic spaces, Poisson kernel

come   L'archive ouverte pluridisciplinaire

Introduction

The purpose of the paper is to investigate some diffusions on subsets of R n . More precisely, we present a new approach, based on Girsanov's theorem, of calculating the basic object of potential theory, that is, the harmonic measure of a ball of hyperbolic Brownian motion or the Ornstein-Uhlenbeck process. The importance of the hyperbolic Brownian motion stems from the fact that it is the canonical diffusion on hyperbolic spaces; it has also some important applications in the risk theory in financial mathematics, see [START_REF] Dufresne | The distribution of a perpetuity, with application to risk theory and pension funding[END_REF] and [START_REF] Yor | Some Aspects of Brownian Motion, Part I: Some Special Functional[END_REF]. Also there is a connection with the Brownian motion in hyperbolic halfspaces (see e.g. [START_REF] Yor | Some Aspects of Brownian Motion, Part I: Some Special Functional[END_REF], [START_REF] Baldi | Stable laws arising from hitting distributions of processes on homogeneous trees and the hyperbolic half-plane[END_REF], [START_REF] Baldi | Non-symmetric hitting distributions on the hyperbolic half-plane and subordinated perpetuities[END_REF]). In a series of papers ( [START_REF] Byczkowski | Poisson kernels of half-spaces in real hyperbolic spaces[END_REF], [START_REF] Byczkowski | Poisson kernel and Green function of the ball in real hyperbolic spaces[END_REF], [START_REF] Żak | Poisson kernel and Green Function of balls for complex hyperbolic Brownian motion[END_REF]) the harmonic measure and the Green operator in hyperbolic spaces were investigated. The main objective of these efforts was to find explicit integral representation of the density functions of these objects. Representations of this kind are crucial in obtaining estimates of harmonic measures and also the Poisson integral [START_REF] Byczkowski | Poisson kernels of half-spaces in real hyperbolic spaces[END_REF]. Authors of all these papers use the Feynman-Kac formula as a main tool in describing the distribution of a stopped multiplicative functional. The present approach is based on methods related to Girsanov's theorem (see [START_REF] Cranston | Conditional transformation of drift formula and potential theory for 1 2 ∆ + b∇[END_REF]). As a result we obtain representation formulas for the harmonic measure, different from those mentioned above. The advantage of this approach is seen in Theorem 4, where we obtain the precise asymptotics of harmonic measure for large values of parameters. Another, important result is the Theorem 6, where we provide adequate representation of harmonic measure of a ball. The real hyperbolic Brownian motion is naturally related to the Ornstein-Uhlenbeck processes. The importance of the Ornstein-Uhlenbeck processes in stochastic modelization and financial mathematics is well known, see e.g. [START_REF] Lamberton | Introduction au Calcul Stochastique Appliquea la Finance[END_REF]. The potential theory of Ornstein-Uhlenbeck processes was studied in [START_REF] Cranston | Conditional transformation of drift formula and potential theory for 1 2 ∆ + b∇[END_REF], [START_REF] Alili | Representations of the first hitting time density of an Ornstein-Uhlenbeck process[END_REF] and [START_REF] Graczyk | Exit times and Poisson kernels of the Ornstein-Uhlenbeck diffusion[END_REF].

The paper is organized as follows. In Section 2 we provide the general framework for the next sections. Our approach is based on the paper [START_REF] Cranston | Conditional transformation of drift formula and potential theory for 1 2 ∆ + b∇[END_REF], where the potential theory for bounded sets D and the operator

1 2 ∆ + F • ∇ (1)
was established. In the present paper we are concentrated on some special cases of the vector fields F and on D being a ball or a half-plane. In particular, we further assume that we deal with potential vector field F on D, orthogonal to the boundary, that is F = ∇ V on D and F is orthogonal to the boundary ∂D .

(

) 2
In Section 3 we provide a closed formula for the density function P (x n , y) (Poisson kernel) of harmonic measure of half-space for the hyperbolic Brownian motion on the real hyperbolic space H n (Theorem 3) and provide an asymptotic formula for P (x n , y) (Theorem 4). In Section 4 we provide an integral representation of Poisson kernel for centered balls for hyperbolic Brownian motion on the ball model D n . We remark here that a similar representation from the paper [START_REF] Byczkowski | Poisson kernel and Green function of the ball in real hyperbolic spaces[END_REF] depends on additional conjectures on the zeros of some hypergeometric functions which so far remained unsettled. The important tool here, as in the next section, is the skew-product representation of the n-dimensional Brownian motion. In Section 5 we provide an integral representation for Poisson kernel of a ball for the classical Ornstein-Uhlenbeck process. In Appendix we collect some useful information on Bessel functions, hypergeometric and Legendre functions and on skew-product of n-dimensional Brownian motion.

2 Change of measure due to Girsanov theorem Notation: For n > 2 we denote by R n the n-dimensional Euclidean space, x, y denotes the standard inner product of x, y ∈ R n and by |x| we denote the Euclidean length of a vector x ∈ R n . The ball with center at zero and the radius r is written as B r = {x ∈ R n : |x| < r}; its boundary, which is the (n -1)-dimensional sphere, is denoted by S n-1 r = {x ∈ R n : |x| = r} and the spherical measure on S n-1 r is denoted by σ n-1 r . Furthermore, we write f (x) ∼ g(x), when lim x→b f (x)/g(x) = 1, as x → b. If for two functions f and g there exist constants c 1 , c 2 such that c 1 < f (x)/g(x) < c 2 for every x ∈ D we will write shortly f ≈ g, x ∈ D.

Throughout the paper D will stand for a domain D ⊂ R n with a smooth, connected boundary ∂D and F will be a bounded vector field which is defined on an open Lipschitz set U containing D. We assume that F is continuously differentiable up to the boundary of U and continuously vanishes on the boundary of U . We further assume that F = ∇V on D for a scalar valued function V and call the function V a potential (and F a potential vector field on D). We set F = V ≡ 0 on the complement of U . We say that the vector field F is orthogonal to the boundary ∂D of the set D if for every differentiable curve Γ : [0, 1) → ∂D we have F (Γ(s)) • Γ ′ (s) = 0 for every s ∈ [0, 1). Lemma 1. Under the above assumptions, if the potential vector field F is orthogonal to the boundary of the set D then the potential function V determined, up to a constant, by the equation

∇V (x) = F (x)
is constant on the boundary ∂D.

Proof. Set x 0 ∈ ∂D. The potential V is described by the curve integral

V (x) = γ F (r)dr + V (x 0 ),
where γ is an arbitrary continuously differentiable path beginning at x 0 and ending at x. For x ∈ ∂D we choose γ to follow the boundary of the set D i.e. γ : [0

, 1] → ∂D, γ(0) = x 0 , γ(1) = x. Then V (x) -V (x 0 ) = 1 0 F (γ(s)) • γ ′ (s)ds = 0 .
Since ∂D is connected, we obtain the conclusion.

Throughout the paper we work within the framework of the canonical representation of the process, i.e. our basic probability space is the space of all continuous R n -valued functions defined on [0, ∞) with appropriate σ-fields (see [START_REF] Ikeda | Stochastic Differential Equations and Diffusion Processes[END_REF]). The standard n-dimensional Brownian motion is denoted by

W (t) = (W 1 (t), . . . , W n (t)).
Define the process X by the SDE

dX(t) = dW (t) + F (X(t))dt , (3)
under the conditions specified above. Then X is a local diffusion on U with the generator L = 1 2 ∆+F (x)•∇. Since the field F is bounded, X can be defined as a local semimartingale (see, e.g. [START_REF] Ikeda | Stochastic Differential Equations and Diffusion Processes[END_REF]). Let τ be the first exit time of the trajectory from the set D. The harmonic measure w x on ∂D, is defined as the distribution of X(τ ) under the distribution P x of the process X starting at x ∈ D. We define a local martingale M by the formula

M (t) = t 0 F (W (s)) • dW (s). (4)
Its quadratic variation is then given by the formula

M (t) = t 0 |F (W (s))| 2 ds.
We further define the basic object of our study, namely

N (t) = exp M (t) - 1 2 M (t) . (5)
We now provide the basic formula for the harmonic measure of the process defined by (3) under some additional conditions.

Theorem 2. Under conditions stated above, suppose that X is the process defined by the system (3). Assume additionally (i) The vector field F is potential and orthogonal to the boundary of the domain D, (ii) For every t > 0

E x [exp( M (t ∧ τ ))] < ∞ , (6)
(iii) {N (t ∧ τ )} t>0 is uniformly integrable . ( 7
)
Then for x ∈ D the harmonic measure w x has the density function expressed by the formula

w x (dz) = e V (∂D)-V (x) E x exp - 1 2 τ 0 |∇V (W (s))| 2 + ∆V (W (s)) ds ; W (τ ) ∈ dz ,
where V is the potential function of the field F and V (∂D) is its value on the boundary ∂D.

Proof. According to (4) the process

F (W (t)) • dW (t) - 1 2 |F (W (t))| 2 dt = dM (t) - 1 2 d M (t)
is a local semimartingale. Writing, as in ( 5)

N (t) = exp(M (t) - 1 2 M (t))
we obtain, as an application of Itô's formula, that N (t) is a local martingale. If we define the measure Q x by dQ x dP x | Ft∧τ = N (t ∧ τ ) , then, as a consequence of Girsanov's theorem, (W, Q x ) and (X, P x ) are different descriptions of the same process, up to time τ , see [START_REF] Cranston | Conditional transformation of drift formula and potential theory for 1 2 ∆ + b∇[END_REF]. Consequently, for a continuous bounded function f defined on R n we obtain

E x f (X(t ∧ τ )) = E x (N (t ∧ τ ); f (W (t ∧ τ )))
Now, the condition [START_REF] Byczkowski | Poisson kernel and Green function of the ball in real hyperbolic spaces[END_REF] shows that the expression on the right-hand side converges to E x [N (τ ); f (W (τ ))], as t → ∞. The left-hand side converges to E x f (X(τ )), by the continuity of the process X. This indicates that indeed w x has density given by

w x (dz) = E x [N (τ ) ; W (τ ) ∈ dz] (8)
We now provide a further description of the function w x . Recall that F is the potential of the vector field B. Define

Z(t) = V (W (t)).
Applying the Itô formula we see that

Z(t) -Z(0) = t 0 ∇V (W (s)) • dW (s) + 1 2 t 0 ∆V (W (s))ds = = M (t) - 1 2 M (t) + 1 2 M (t) + 1 2 t 0 ∆V (W (s))ds.
Consequently

N (t) = exp(M (t) - 1 2 M (t)) = exp Z(t) -Z(0) - 1 2 t 0 |∇V (W (s))| 2 + ∆V (W (s)) ds . (9)
Remark that Y (τ ) = V (W (τ )) = V (∂D) since the vector field F is orthogonal to the boundary of the set D. Hence, stopping at t ∧ τ and taking expectation we get, when t → ∞

w x (dz) = e V (∂D)-V (x) E x exp - 1 2 τ 0 |∇V (W (s))| 2 + ∆V (W (s)) ds ; W (τ ) ∈ dz .
This, together with [START_REF] Byczkowski | Hitting half-spaces by Bessel-Brownian diffusions[END_REF], finishes the proof.

3 Harmonic measure of hyperbolic horocycle in H n

For every a > 0 we define H a = {x ∈ R n : x n > a}. In this section we consider the harmonic measure ω x a of the set H a for the operator

∆ LB = x 2 n n i=1 ∂ 2 ∂x 2 i -(n -2)x n ∂ ∂x n , n 2.
The motivations for studying this operator comes from hyperbolic geometry. More precisely, this operator is the Laplace-Beltrami operator associated with Riemannian metric in the half-space model H n of the real n-dimensional hyperbolic space. From geometric point of view, the set H a is an interior of the hyperbolic horocycle ∂H a = {x ∈ R n :

x n = a}. Let (B i (t)) i=1...n be an n-dimensional Brownian motion on R n with the generator d 2 dx 2 (and not 1 2

d 2
dx 2 ) i.e. the variance E 0 B 2 i (t) = 2t. Then the Brownian motion on H n , Y = (Y i ) i=1...n , can be described by the following system of stochastic differential equations

       dY 1 (t) = Y n (t)dB 1 (t) dY 2 (t) = Y n (t)dB 2 (t) . . . dY n (t) = Y n (t)dB n (t) -(n -2)Y n (t)dt. ( 10
)
By the Itô formula one verifies that the generator of the solution of this system is ∆ LB . The Laplace-Beltrami operator can be rewritten in the form ∆ LB = 2x 2 n L 1 , where

L 1 = 1 2 ∆ + F 1 (x) • ∇, (11)
with F 1 (x) = (0, . . . , 0, (2 -n)/(2x n )). Now, we make change of time. Namely, we write

A(u) = u 0 Y 2 n (s) ds ,
and σ(t) = inf{u > 0; A(u) > t} .

If we now write Bk (t) =

σ(t/2) 0 Y n (s) dB k (s) , k = 1, ..., n
then Bk are martingales with mutual variations < Bk , Bl > (t) = δ(k, l) t, k, l = 1, ...n so B = ( Bk ) is the standard n-dimensional Brownian motion. Substituting

Ỹk (t) = Y k (σ(t/2)) , k = 1, ..., n
we obtain that (10) transforms into the following the system of SDE

         d Ỹ1 (t) = d B1 (t) d Ỹ2 (t) = d B2 (t) . . . d Ỹn (t) = Z n (t)d Bn (t) -(n -2) dt 2Zn(t) . (12)
Again, by Itô's formula, we see that the operator L 1 is the generator of the process Ỹ = ( Ỹk ). Since the operation of change of time does not affect the exit place, the harmonic measures of the operators ∆ LB and L 1 are the same.

The potential of the vector field F 1 is given by V 1 (x) = (2 -n) ln(x n )/2. Moreover, it is easy to check that the vector field F 1 is orthogonal to the boundary of H a . Using (4) and [START_REF] Byczkowski | Hitting times of Bessel processes[END_REF], we obtain

M (t) = 2 -n 2 t 0 dW n (s) W n (s) , M (t) = n -2 2 2 t 0 ds W 2 n (s) , N (t) = W n (0) W n (t) (n-2)/2 exp - n(n -2) 8 t 0 ds W 2 n (s)
,

where

W (t) = (W 1 (t), . . . , W n (t)) denotes the standard Brownian motion in R n starting from W (0) = x. If we put τ a = inf{t > 0 : W (t) / ∈ H a } we obtain E x exp( M (t ∧ τ a )) = E x exp (n -2) 2 4 t∧τa 0 ds W 2 n (s) E x exp (n -2) 2 4 t∧τa 0 ds a 2 < ∞.
Moreover, we have

N (t ∧ τ a ) W n (0) W n (t ∧ τ a ) (n-2)/2
x n a

(n-2)/2
, t 0. Now, the results of Theorem 2 imply that

ω x a (dy) = x n a n-2 2 E x exp - n(n -2) 8 τa 0 ds W n (s) 2 ; W (τ a ) ∈ dy . (13)
The above-given formula enable us to find the density function P a (x, y), x ∈ H a , y ∈ ∂H a , of the measure ω x a (dy) with respect to the Lebesgue measure on ∂H a . The scaling property of n-dimensional Brownian motion implies the scaling property for Poisson kernels P a (x, y) = a 1-n P 1 (x/a, y/a), x ∈ H a , y ∈ ∂H a . Moreover, Brownian motion W (t) and the set H 1 are invariant under translations (x, x n ) → (x + b, x n + b), where b ∈ R n-1 . Consequently, P 1 (x, y) = P 1 ((0, x n ), (ỹ -x, 1)) for every x ∈ H 1 and y ∈ ∂H 1 . We will use these properties in further considerations to simplify the notation. We will denote by P (x n , y) = P 1 ((0, x n ), (y, 1)), where y ∈ R n-1 and τ = τ 1 .

Theorem 3. For every x n > 1 and y ∈ R n-1 we have

P (x n , y) = 1 2 ν-1 π ν+1 x ν n |y| ν-1 ∞ 0 J ν (t)Y ν (tx n ) -J ν (tx n )Y ν (t) J 2 ν (t) + Y 2 ν (t) t ν K ν-1 (t|y|) dt, (14)
where ν = (n -1)/2.

Proof. Observe that the integral appearing in (13) as well as the hitting time τ depend only on the last coordinate of the Brownian motion W (t) = ( W (t), W n (t)). Since the processes W (t) and W n (t) are independent we obtain

ω x 1 (dy) = x n-2 2 n ∞ 0 E 0 [ W (s) ∈ dy] E xn exp - n(n -2) 8 τ 0 ds W n (s) 2 ; τ ∈ ds = x n-2 2 n ∞ 0 exp(-|y| 2 /(2s)) (2πs) (n-1)/2 µ xn (ds) dy, (15)
where

µ xn (ds) = E xn exp - n(n -2) 8 τ 0 ds W n (s) 2 ; τ ∈ ds .
The Laplace transform of µ xn is given by

Lµ xn (w) = E xn exp - n(n -2) 8 τ 0 ds W n (s) 2 -w τ 0 ds = E xn e q (τ ), w 0, where q(x) = -n(n -2)/(8x 2 ) -w. The function ϕ(x n ) = E xn e q (τ
) is a gauge function for an appropriate Schrödinger operator based on generator of W n (t). Consequently, ϕ is a bounded solution of the equation

1 2 ϕ ′′ (x) - n(n -2) 8x 2 + w φ(x) = 0, x 1 such that ϕ(1) = 1. Making substitution √ xψ(x √ 2w) = ϕ(x) we reduce the above-given equation to 2x 2 wψ ′′ (x √ 2w) + x √ 2wψ ′ (x √ 2w) - (n -1) 2 4 + 2wx 2 ψ(x √ 2w) = 0,
which is the modified Bessel equation ( 33) with ν = n-1 2 . Taking into account the general form of solutions of (33), the boundary condition and boundedness of ϕ we arrive at

Lµ xn (w) = √ x n K ν (x n √ 2w) K ν ( √ 2w) , w 0. ( 16
)
Since square root can be extended to a holomorphic function on C \ (-∞, 0] and the modified Bessel function K ν has no zeros in the positive half-plane Re w 0, the Laplace transform Lµ xn (w) can also be extended to analytic function on C \ (-∞, 0]. Moreover, using the asymptotic expansion (34) we obtain that

|Lµ xn (w)| e -(xn-1) √ 2w 1 + E(x n √ 2w) 1 + E( √ 2w) 2 exp -(x n -1) 2|w| cos arg w 2 , (17)
for every w ∈ C \ (-∞, 0] such that |w| is large enough. Note that here arg w ∈ [-π, π].

In particular, Lµ xn (w) is bounded for |w| 1. These properties of Lµ xn and its analytic continuation guarantee that we can apply the inverse Laplace transform to [START_REF] Gradstein | Table of integrals, series and products[END_REF] (see [START_REF] Folland | Fourier Analysis and its applications[END_REF], Theorem 8.5, p. 267). More precisely, there exist a density function of µ xn with respect to Lebesgue measure on (0, ∞) given by the inversion formula

µ xn (s) = 1 2πi lim r→∞ 1+ir 1-ir
Lµ xn (w)e sw dw.

To compute the limit we integrate the function f s (w) = Lµ xn (w)e sw over rectangular contour surrounding the branch-cut of f s which is negative real axis. Let Γ be the positively oriented contour consisting of four horizontal segments

γ 1 = [-r + i/r, i/r], γ 2 = [-r -i/r, -i/r], γ 3 = [-r + ir, 1 + ir], γ 4 = [-r -ir, 1 -ir], three vertical segments γ 5 = [-r + i/r, -r + ir], γ 6 = [-r -i/r, -r -ir], γ 7 = [1 -ir, 1 + ir] and a semi-circle γ 8 = {|w| = 1/r, Re w > 0}.
The formula (35) implies that f s is bounded for small w such that Re w > 0. Consequently, the integral over γ 8 tends to zero when r → ∞. The boundedness of Lµ xn (w) for large w implies that for r 1 and every s > 0 we have

γ 5 + γ 6 f s (w)dw 2 sup |w| 1
Lµ xn (w) re -rs → 0, as r → ∞. Finally, using [START_REF] Ikeda | Stochastic Differential Equations and Diffusion Processes[END_REF], we obtain

γ 3 + γ 4 f s (w)dw 4 exp(-(x n -1) cos(3π/8) √ r) ∞ -1
e -su du → 0 as r → ∞. The Cauchy Theorem together with previous considerations and (36) give

µ xn (s) = 1 2πi lim r→∞ 1+ir 1-ir Lµ xn (w)e sw dw = 1 2πi lim r→∞ γ 1 + γ 2 f s (w)dw = √ x n 2πi ∞ 0 K ν (-i √ 2tx n ) K ν (-i √ 2t) - K ν (i √ 2tx n ) K ν (i √ 2t) e -st dt = - √ x n π ∞ 0 Im K ν (i √ 2tx n ) K ν (i √ 2t) e -st dt = - √ x n π ∞ 0 Im K ν (itx n ) K ν (it) te -st 2 /2 dt = √ x n π ∞ 0 J ν (t)Y ν (tx n ) -J ν (tx n )Y ν (t) J 2 ν (t) + Y 2 ν (t) te -st 2 /2 dt.
From (31) and ( 32) is easy to see that

J 2 ν (t) + Y 2 ν (t) ∼ t -1 ∨ t -2ν and J ν (t)Y ν (tx n ) -J ν (tx n )Y ν (t) J 2 ν (t) + Y 2 ν (t) C (1 + t) , t > 0 for some constant C = C(x n ) > 0.
With the use of (37) we verify

∞ 0 |J ν (t)Y ν (tx n ) -J ν (tx n )Y ν (t)| J 2 ν (t) + Y 2 ν (t) t ∞ 0 e -|y| 2 /(2s) e -st 2 /2 ds (2πs) (n-1)/2 dt C ∞ 0 (t + 1)t ν K ν-1 (t|y|)dt.
The last integral is finite by ( 34) and (35). Consequently, by the Fubini's theorem we obtain

P (x n , y) = x ν n π ∞ 0 J ν (t)Y ν (tx n ) -J ν (tx n )Y ν (t) J 2 ν (t) + Y 2 ν (t) t ∞ 0 e -|y| 2 /(2s) e -st 2 /2 ds (2πs) (n-1)/2 dt = 1 2 ν-1 π ν+1 x ν n |y| ν-1 ∞ 0 J ν (t)Y ν (tx n ) -J ν (tx n )Y ν (t) J 2 ν (t) + Y 2 ν (t) t ν K ν-1 (t|y|) dt.
This ends the proof.

The integral formula presented in Theorem 3 can be used to obtain the asymptotics of the Poisson kernel P (x n , y) as well as its sharp bounds for small x n and large |y|. Note that results given in the next theorem cover those obtained in [START_REF] Byczkowski | Poisson kernels of half-spaces in real hyperbolic spaces[END_REF] (see Theorem 4.9 and Theorem 4.10, compare also with Theorem 5.3 in [START_REF] Byczkowski | Hitting distibution of geometric Brownian motion[END_REF]). Moreover, the formula ( 14) allows to omit very laborious and sophisticated computations used in [START_REF] Byczkowski | Poisson kernels of half-spaces in real hyperbolic spaces[END_REF] to examine the behavior of P (x n , y) when |y| tends to infinity. Our approach is simpler and gives more general results.

Theorem 4. For every x 0 1 we have

P (x n , y) ∼ Γ(n/2) 2 n-2 π n/2 n-2 k=0 x k 0 • x n -1 |y| 2n-2 , x n → x 0 , |y| → ∞. (18)
Moreover, for every y 0 > 0 we have

P (x n , y) ∼ c(y 0 ) π (n+3)/2 2 (n-5)/2 • (x n -1), x n → 1, |y| → y 0 , (19)
where

c(y 0 ) = |y 0 | 1-n 2 ∞ 0 s ν K ν-1 (sy 0 )ds J 2 ν (s) + Y 2 ν (s)
.

Proof. Making a substitution t|y| = s in ( 14) we can rewrite the Poisson kernel in the following way

P (x n , y)|y| 4ν x n -1 = x ν n π ν+1 2 ν-1 ∞ 0 g ν x n , s |y| s 3ν K ν-1 (s)ds,
where

g ν (x, t) = 1 t 2ν (x -1) J ν (t)Y ν (xt) -J ν (tx)Y ν (t) J 2 ν (t) + Y 2 ν (t) , x > 1, t > 0.
Since

s 3ν K ν-1 (s) is integrable on [0, ∞)
P (x n , y)|y| 4ν x n -1 = 1 π ν 2 3ν-2 Γ 2 (ν) ∞ 0 s 3ν K ν-1 (s)ds = 2Γ(2ν)Γ(ν + 1) π ν Γ 2 (ν) and lim (x,|y|)→(x 0 ,∞) P (x n , y)|y| 4ν x n -1 = x 0 π ν+1 2 ν-1 ∞ 0 π(x 2ν 0 -1)s 3ν K ν-1 (s)ds 2 2ν (x 0 -1)Γ(ν)Γ(ν + 1)x ν 0 = (x 2ν 0 -1)Γ(2ν) π ν (x 0 -1)Γ(ν)x ν-1 0 = n-2 k=0 x k 0 Γ(2ν) π ν x ν-1 0 Γ(ν) whenever x 0 > 1.
Here we used formula (39) from Lemma 8 and relation (38). The duplication formula for gamma function gives [START_REF] Itô | Diffusion Processes and Their Sample Paths[END_REF]. In the same way, using (39), we get lim

(xn,|y|)→(1,y 0 ) P (x n , y)|y| 4ν x -1 = 1 π ν+1 2 ν-1 ∞ 0 2s ν |y 0 | 2ν K ν-1 (s)ds π J 2 ν s |y 0 | + Y 2 ν s |y 0 |
. which, by substitution s = ty 0 , proves [START_REF] Lamberton | Introduction au Calcul Stochastique Appliquea la Finance[END_REF].

As a consequence of Theorem 4 we obtain the following sharp bounds for the Poisson kernel for small x n and large |y|. Similar results have been obtained recently (see Theorem 11 in [START_REF] Byczkowski | Hitting times of Bessel processes[END_REF]). The results are more general (there is no restriction for x n and |y|) however the methods of proof are much more complicated.

Corollary 5. We have

P (x n , y) ≈ x n -1 |y| 2n-2 , 1 < x n 2, |y| 1.
Proof. Existence and positivity of the limits proved in Theorem 4 imply that for every x ∈ [1, 2] there exist ε x > 0 and Y x > 0 and strictly positive constants c 1 (x), c 2 (x) such that

c 1 (x) x ′ -1 |y| 2n-2 P (x ′ , y) c 2 (x) x ′ -1 |y| 2n-2 for every x ′ ∈ [1, 2] and y ∈ R n-1 satisfying |x -x ′ | < ε x and |y| > Y x . Since the family of intervals (x -ε x , x + ε x
) is an open cover of [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Alili | Representations of the first hitting time density of an Ornstein-Uhlenbeck process[END_REF], we can choose the finite subcover {(x i -ε x i , x i + ε x i ), i = 1, . . . m}. Putting Y = max{Y x i : i = 1, . . . , m}, c 1 = min{c 1 (x i ) : i = 1, . . . , m} and c 2 = max{c 2 (x i ) : i = 1, . . . , m} we get

c 1 x n -1 |y| 2n-2 P (x n , y) c 2 x n -1 |y| 2n-2 (20)
for every x n ∈ [1, 2] and |y| > Y . Observe that formula (15) implies positivity of P (x n , y). Moreover, by [START_REF] Lamberton | Introduction au Calcul Stochastique Appliquea la Finance[END_REF], we get that the function P (x n , y)|y| 2n-2 /(x n -1) can be continuously extended to the strictly positive function on the compact set [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Alili | Representations of the first hitting time density of an Ornstein-Uhlenbeck process[END_REF]×[1, Y ] and consequently [START_REF] Lebedev | Special functions and their applications[END_REF] is true also for |y| ∈ [1, Y ] (with possibly different constants c 1 and c 2 ). This ends the proof.

Harmonic measure of hyperbolic balls in D n

In this section we consider the harmonic measure of the balls associated with the operator

∆ LB = (1 -|x| 2 ) 2 4 n i=1 ∂ 2 ∂x 2 i + (n -2) 1 -|x| 2 2 n i=1 x i ∂ ∂x i . ( 21
)
This operator appears naturally as the Laplace-Beltrami operator on the ball model of the real hyperbolic space D n (see [START_REF] Byczkowski | Poisson kernel and Green function of the ball in real hyperbolic spaces[END_REF] for more details).

In particular, if B = (B k ) is the standard n-dimensional Brownian motion then the following system of SDE:

dY k (t) 1 -|Y (t)| 2 = dB k (t) + 2(n -2)Y k (t)dt, k = 1, . . . , n.
describes a diffusion with values in the real hyperbolic space D n , with the generator 2∆ LB . As in the case of half-space model, we perform a change of time defined by:

A(u) = u 0 (1 -|X(s)| 2
) ds , and σ(t) = inf{u; A(u) > t} .

Then the process defined by

Ỹk (t) = Y k (σ(t))
is the diffusion with values in D n with the generator L 2 .

As in the case of H n in Section 4, the harmonic measures of the operators ∆ LB and L 2 are the same.

We consider now the harmonic measure ω x r of a ball B r = {x ∈ R n : |x| < r}, r < 1, supported on the boundary of B r which is the sphere S n-1 r of radius r. We denote by P r (x, y) the Poisson kernel of B r , i.e. the density of the measure ω x r with respect to the (n -1)-dimensional spherical measure σ n-1 r . As in the previous section we can write the Laplace-Beltrami operator in the form (

1-|x| 2 ) 2 2 L 2 , where L 2 = 1 2 ∆ + F 2 (x) • ∇, with F 2 (x) = n-2 1-|x| 2 (x 1 , . . . , x n ).
The positivity of the factor (1 -|x| 2 ) 2 /2 implies that the harmonic functions on the ball B r for the operators (21) and L 2 are exactly the same and consequently the harmonic measures coincide.

Moreover, the vector field F 2 is orthogonal to the sphere S n-1 r and its potential function is V 2 (x) = 2-n 2 ln(1 -|x| 2 ). We denote by τ r = inf{t > 0; W (t) / ∈ B r } the first exit time of Brownian motion W (t) from a ball B r . Note that τ r depends only on the Euclidean norm of W . The martingale M , related to the vector field F 2 , and its quadratic variation are

M (t) = (n -2) t 0 W i (s)dW i (s) 1 -|W (s)| 2 , M (t) = (n -2) 2 t 0 |W (s)| 2 ds (1 -|W (s)| 2 ) 2 .
Observe that the condition ( 6) is fulfilled in this case, since

E x exp (n -2) 2 t∧τr 0 |W (s)| 2 ds (1 -|W (s)| 2 ) 2 E x exp (n -2) 2 r 2 (t ∧ τ r ) (1 -r 2 ) 2 < ∞ .
By [START_REF] Byczkowski | Hitting times of Bessel processes[END_REF], the kernel N (t) is of the form

N (t) = 1 -|W (0)| 2 1 -|W (t)| 2 n-2 2 exp - n(n -2) 2 t 0 ds (1 -|W (s)| 2 ) 2 ,
and it is now evident that {N (t ∧ τ r )} t>0 is uniformly bounded in t so the condition ( 7) holds. Applying Theorem 2 we obtain

w x r (dy) = 1 -|x| 2 1 -r 2 n-2 2 E x exp - n(n -2) 2 τr 0 ds (1 -|W (s)| 2 ) 2 ; W (τ r ) ∈ dy , (22)
From now on we assume that x = 0. For x = 0, from the rotational invariance of the Laplace-Beltrami operator, we easily obtain that ω x r is just

σ n-1 r /σ n-1 r (S n-1 r
). Recall the skew-product representation of the Brownian motion

W (t) = R (ν) (t)Θ(A (ν) (t)),
where R (ν) is the Bessel process with index ν = n 2 -1 staring from |x| and Θ is spherical Brownian motion on S n-1 1 independent from R (ν) (see Appendix). Using the fact that τ r depends only on R (ν) we get that W (τ r ) = R (ν) (τ r )Θ(A (ν) (τ r )), where Θ is independent from R (ν) (τ r ) and A (ν) (τ r ). Applying this decomposition to formula [START_REF] Żak | Poisson kernel and Green Function of balls for complex hyperbolic Brownian motion[END_REF] we get

ω x r (dy) = 1 -|x| 2 1 -r 2 n-2 2 ∞ 0 P x |x| (Θ t ∈ dy)E |x| exp τr 0 q(R (ν) s )ds ; A (ν) (τ r ) ∈ dt ,
where q(y) = -n(n-2) 2(1-y 2 ) 2 . Rotational invariance of spherical Brownian motion implies that the harmonic measure ω x r is axially symmetric with axis x. As a consequence we get that its density P r (x, y) depends only on the cosine of angle between starting point x and the point y, i.e. P r (x, y) = Pr x,

x, y |x||y| .

If we consider the sets of the form

A = {η ∈ S n-1 1 : x,η |x| ∈ (a, b)},
where -1 < a < b < 1 we get (for a definition of the process S see Appendix)

P x |x| (Θ t ∈ A) = P 1 (S t ∈ (a, b)) = b a p S t (1, z)m(dz) = 2 b a p S t (1, z)(1 -z 2 ) (n-3)/2 dz,
where p S t is defined in (51). From the other side, using the spherical coordinates we obtain

ω x r (rA) = rA P r (x, y)dσ n-1 r (y) = r n-1 σ n-2 r (S n-2 1 ) cos φ∈(a,b) Pr (x, cos φ) sin n-2 φ dφ = nπ n-1 2 r n-1 Γ n+1 2 b a Pr (x, z)(1 -z 2 ) (n-3)/2 dz.
Comparing both sides we get the following formula for the Poisson kernel P r (x, y)

P r (x, y) = Γ 2 n+1 2 π n-1 2 nr n-1 1 -|x| 2 1 -r 2 n-2 2 ∞ 0 p S t 1, x, y |x||y| µ |x| (dt), (23)
where

µ y (dt) = E y exp τr 0 q(R (ν) (s))ds ; A (ν) (τ r ) ∈ dt , y ∈ (0, r].
The formula for p S t can be computed from the appropriate formula for the transition density function for Θ, which is given in terms of spherical harmonics and that approach leads to the series representation for P r (x, y) presented in [START_REF] Byczkowski | Poisson kernel and Green function of the ball in real hyperbolic spaces[END_REF]. However, we want to compute the Laplace transform of p S t which is so called λ-Green function of the process S

G λ (x, 1) = ∞ 0 e -λt p S t (1, x)dt, x ∈ (-1, 1),
and we do it directly. From the general theory (see for example [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF] Chapter II for short resume) a function G λ is described by solutions of the second-order differential equation

1 -x 2 2 u ′′ (x) - n -1 2 x u ′ (x) = λu(x), x ∈ (-1, 1). ( 24
)
Note that the expression on the left-hand side is just Gu(x), where G is just the generator of S described in (50). More precisely, we have

G λ (x, 1) = ϕ λ (1)ψ λ (x) w λ , x ∈ (-1, 1),
where ϕ λ is a decreasing and ψ λ is an increasing solution of (24) such that ϕ - λ (1) = 0 and ψ + λ (-1) = 0. The boundary conditions for the derivatives follows from the fact that non-singular points -1 and 1 are reflecting. Here f + and f -denote the right and left derivative with respect to the speed function s(x). The Wronskian w λ is given by

w λ = ψ - λ (x)ϕ λ (x) -ψ λ (x)ϕ - λ (x)
and it does not depend on x. Putting x = 1 in the above-given formula and using the boundary conditions we obtain that

G λ (x, 1) = ϕ λ (1)ψ λ (x) ψ - λ (1)ϕ λ (1) -ψ λ (1)ϕ - λ (1) = ψ λ (x) ψ - λ (1) , x ∈ (-1, 1).
This implies that G λ (x, 1) is uniquely described as a solution of (24) such that u + (-1) = 0 and u -(1) = 1. Making a substitution u(x) = f (z) with z = 1+x 2 in the equation ( 24) we reduce it to the following hypergeometric equation

z(1 -z)f ′′ (z) + n -1 2 -(n -1)z f ′ (z) -2λf (z) = 0. with α = n-2 2 -A(λ), β = n-2 2 +A(λ), γ = n-1 2 , where A(λ) = 1 2 (n -2) 2 -8λ. Consider the function h λ (x) = 2 F 1 n -2 2 -A(λ), n -2 2 + A(λ); n -1 2 ; 1 + x 2 .
The above-given computation implies that the function h λ is a solution of (24). Using ( 44) and (45) we compute the derivative of this function with respect to the scale function s(x) in the following way

(1 -x 2 ) n-1 2 d dx h λ (x) = (1 -x 2 ) n-1 2 2λ n -1 2 F 1 n 2 -A(λ), n 2 + A(λ); n + 1 2 ; 1 + x 2 = 2 n+1 2 λ n -1 (1 + x) n-1 2 2 F 1 1 2 + A(λ), 1 2 -A(λ); n + 1 2 ; 1 + x 2 .
The first equality and the fact that the hypergeometric function 2 F 1 is equal to 1 at zero implies h + λ (-1) = 0. Using the second equality and ( 43) we obtain

h - λ (1) = 2 n λ n -1 2 F 1 1 2 + A(λ), 1 2 -A(λ); n + 1 2 ; 1 = 2 n-1 λΓ 2 n-1 2 Γ n-2 2 -A(λ) Γ n-2 2 + A(λ)
.

Moreover, using once again (45) and the definition (47) we can express function h λ in terms of the Legendre function of the first kind

h λ (x) = 1 -x 2 3-n 2 2 F 1 1 2 -A(λ), 1 2 + A(λ); n -1 2 ; 1 + x 2 = (1 -x 2 ) 3-n 4 2 3-n 2 Γ n -1 2 P 3-n 2 A(λ)-1 2 (-x).
Finally we have just obtained that

G λ (x, 1) = B n (λ)(1 -x 2 ) 3-n 4 P 3-n 2 A(λ)-1 2 (-x),
where

B n (λ) = Γ n-2 2 -A(λ) Γ n-2 2 + A(λ) 2 n+1 2 λΓ n-1 2 .
The second part of the formula (23) relates to the measure µ y . Observe that µ y depends only on the Bessel process R (ν) , which is one-dimensional diffusion. For every w 0 the Laplace transform Lµ y (w) is

E y exp - n(n -2) 2 τr 0 ds (1 -(R (ν) (s)) 2 ) 2 -w τr 0 ds (R (ν) (s)) 2 = E y e g (τ r ),
where r ) is by definition the gauge function for the Schrödinger operator based on the generator of the process R (ν) and the non-positive potential g. From the Feynman-Kac formula ϕ is a solution of the Schrödinger equation. Using (49) we obtain that ϕ is a bounded solution to the following second-order differential equation

g(y) = -n(n-2) 2(1-y 2 ) 2 -w y 2 . The function ϕ(y) = E y e g (τ
1 2 ϕ ′′ (y) + n -1 2y ϕ ′ (y) - n(n -2) 2(1 -y 2 ) 2 + w y 2 ϕ(y) = 0, y ∈ [0, r), w 0, ( 25
)
with the boundary condition ϕ(r) = 1. Substituting ϕ(y) = y 1-n/2 ψ 1+y 2 1-y 2 we obtain

ϕ ′ (y) = 1 - n 2 y -n/2 ψ 1 + y 2 1 -y 2 + 4y 2-n/2 (1 -y 2 ) 2 ψ ′ 1 + y 2 1 -y 2 , ϕ ′′ (y) = n 2 n 2 -1 y -1-n/2 ψ 1 + y 2 1 -y 2 + 16y 3-n/2 (1 -y 2 ) 4 ψ ′′ 1 + y 2 1 -y 2 + + 4y 1-n/2 (1 -y 2 ) 3 (3 -n + (n + 1)y 2 )ψ ′ 1 + y 2 1 -y 2 .
Putting the above given formulas to the differential equation ( 25) and dividing both sides by the factor

-2y 1-n/2 (1-y 2 ) 2 give 0 = -4y 2 (1 -y 2 ) 2 ψ ′′ 1 + y 2 1 -y 2 - 2(1 -y 4 ) (1 -y 2 ) 2 ψ ′ 1 + y 2 1 -y 2 + + (1 -y 2 ) 2 4y 2 (n -2) 2 4 + 2w + n(n -2) 4 ψ 1 + y 2 1 -y 2 .
Moreover, putting z = 1+y 2 1-y 2 and using the equality 1

-1+y 2 1-y 2 2 = -4y 2 (1-y 2 ) 2 lead to the following differential equation for ψ (1 -z 2 )ψ ′′ (z) -2zψ ′ (z) + ν(ν + 1) - A(-w) 2 1 -z 2 ψ(z) = 0, z 1. with ν = n 2 -1 and A(-w) = 1 2 (n -2) 2 + 8w
. This is the Legendre's differential equation ( 46). Thus, the general solution of (25) is given by

ϕ(y) = c 1 • y 1-n/2 P -A(-w) ν 1 + y 2 1 -y 2 + c 2 • y 1-n/2 Q -A(-w) ν 1 + y 2 1 -y 2 , y ∈ [0, r],
where c 1 and c 2 are absolute constants. Using (47) and (48) one can easily check that the function y

-ν P -A(-w) ν 1+y 2 1-y 2 is bounded on the interval 1, 1+r 2 1-r 2 in contrast to the function y -ν Q -A(-w) ν 1+y 2
1-y 2 , which is unbounded in the neighborhood of 1. Thus c 2 = 0 and the boundary condition ϕ(r) = 1 gives

Lµ |x| (w) = r |x| n/2-1 P -A(-w) ν 1+|x| 2 1-|x| 2 P -A(-w) ν 1+r 2 1-r 2 , |x| r, w 0. (26)
Now observe that for every complex number w such that Re(w)

> -ν 2 2 = -(n-2) 2 8 |Lµ |x| (w)| E |x| exp - n(n -2) 2 τr 0 ds (1 -(R (ν) s ) 2 ) 2 -Re(w) τr 0 ds (R (ν) s ) 2 E |x| exp (n -2) 2 8 τr 0 ds (R (ν) s ) 2 = E |x| exp ν 2 2 τr 0 ds (R (ν) s ) 2 = r x ν .
The last equality follows from (see [START_REF] Borodin | Handbook of Brownian Motion -Facts and Formulae[END_REF] 2.20.4 p.407)

P (ν) x τr 0 ds (R (ν) (s)) 2 ∈ dy = r x ν ln(r/x) √ 2πy 3/2 exp - ν 2 y 2 - ln 2 (r/x) 2y dy, y > 0.
In particular Lµ |x| (-ν 2 /2) is finite. This implies that the formula

Lµ |x| (w) = ∞ 0 e -wt µ |x| (dt)
defines a holomorphic function in the complex half-plane Re(w) > -v 2 /2. Moreover, for |z| < 1 the function 2 F 1 (α, β; γ; z)/Γ(α) as a function of α is analytic function on C. Using this fact and the representation of P µ ν in terms of hypergeometric function 2 F 1 we deduce that the function on the right-hand side of ( 26) is a meromorphic function in the half-plane Re(w) > -ν 2 /2. In fact, the equality (26) implies that the ratio of Legendre functions is analytic for Re(w) > -ν 2 /2 and consequently the function in the denominator has no zeros in this region. Compare this result with Conjecture 5.2 in [START_REF] Byczkowski | Poisson kernel and Green function of the ball in real hyperbolic spaces[END_REF]. Moreover, we have just proved that (26) holds whenever Re(w) > -ν 2 /2.

Now let c = -(n-2) 2 16 . We have ∞ 0 p S t 1, x, y |x||y| µ |x| (dt) = ∞ 0 p S t 1, x, y |x||y| 1 2πi c+i∞ c-i∞ e zt Lµ |x| (z)dz dt = 1 2πi c+i∞ c-i∞ Lµ |x| (z) ∞ 0 e zt p S t 1, x, y |x||y| dt dz = 1 2πi c+i∞ c-i∞ Lµ |x| (z)G z 1, x, y |x||y| dtdz.
Taking into account the previously found formulas for the Laplace transform Lµ |x| (z) and the Green function G z (1, x) we finally obtain Theorem 6. For every x ∈ B r , x = 0 and y ∈ S n-1 r the Poisson kernel P r (x, y) is given by the following formula

2Γ n+1 2 π n-1 2 nr n-1 1 -|x| 2 1 -r 2 r |x| ν sin 3-n 2 ϕ 2πi c+i∞ c-i∞ P -A(-z) ν 1+|x| 2 1-|x| 2 P -A(-z) ν 1+r 2 1-r 2 B n (z)P 3-n 2 A(z)-1 2 (-cos ϕ)dz,
where

A(z) = 1 2 (n -2) 2 -8z, B n (z) = Γ( n-2 2 -A(z))Γ( n-2 2 +A(z)) 2 n+1 2 zΓ( n-1 2 )
, c = -ν 2 /4 and ϕ is an angle between x and y,

Harmonic measure of Ornstein-Uhlenbeck process

As in the previous chapter for fixed r > 0 we denote B r = {x ∈ R n : |x| < r}. Consider a vector field F 3 (x) = λx, where λ > 0. F 3 is a potential vector field with the potential function V 3 (x) = λ|x| 2 /2 and as in the previous cases it is orthogonal to the boundary of B r . The corresponding martingale M and its quadratic variation are

M (t) = λ n i=1 t 0 W i (s)dW i (s), M (t) = t 0 λ 2 |W (s)| 2 ds.
The validity of [START_REF] Byczkowski | Poisson kernels of half-spaces in real hyperbolic spaces[END_REF] in this case follows from

E x [exp M (t ∧ τ r )] = E x exp t∧τr 0 λ 2 |W (s)| 2 ds E x [exp[λ 2 r 2 (t ∧ τ r )]] < ∞,
where τ r = inf{t > 0 :

W (t) / ∈ B r }. Since N (t) = exp λ(|W (t)| 2 -|W (0)| 2 ) 2 exp - 1 2 τ 0 (λ 2 |W (s)| 2 + 2nλ)ds ,
all the assumptions of Theorem 2 are satisfied and consequently we obtain that the harmonic measure w x r (dy) of B r for the operator

L 3 = 1 2 ∆ + λx • ∇
is given by

w x r (dy) = exp λ(r 2 -|x| 2 ) 2 E x exp - 1 2 τr 0 (λ 2 |W (s)| 2 + 2nλ)ds ; W (τ r ) ∈ dy .
Computations in this case mimic those introduced in the previous section so we omit some details and present only a sketch of the argumentation. For x = 0 the skew-product representation of the Brownian motion allows us to write

ω x r (dy) = e λ 2 (r 2 -|x| 2 ) ∞ 0 P x |x| (Θ t ∈ dy)E |x| exp τr 0 q(R (ν) (s))ds ; τr 0 ds (R (ν) (s)) 2 ∈ dt ,
where q(y) = -λ 2 2 |x| 2 -nλ and consequently the Poisson kernel P r (x, y) is given by As previously, the Laplace transform Lµ y (w) given by

P r (x, y) = 2Γ n+1 2 π (n-1)/2 nr n-1 e λ 2 (r 2 -|x| 2 ) ∞ 0 p S t 1, x,
E y exp τr 0 - λ 2 2 (R (ν) s ) 2 -nλ - w (R (ν) s ) 2 ds = E y e g (τ r ), (27)
where g(y) = -λ 2 2 y 2 -nλ -w y 2 , can be identified (by applying the Feynman-Kac formula) as a bounded solution of the following Schrödinger equation

1 2 ϕ ′′ (y) + n -1 2y ϕ ′ (y) - λ 2 2 y 2 + nλ + w y 2 ϕ(y) = 0, y ∈ (0, r), w 0,
with the boundary condition ϕ(r) = 1. Setting ϕ(y) = y -n/2 f (λy 2 ) we reduce this equation to the Whittaker equation

f ′′ (x) + f (x) - 1 4 - n 2x - n(n -4) 16 + w 2 1 x 2 = 0, with parameters k = -n 2 and µ = √ (n-2) 2 +8w 4
. Consequently

ϕ(y) = y -n 2 c 1 M (k, µ, λy 2 ) + c 2 W (k, µ, λy 2 ) ,
where M and W are Whittaker functions (see [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] 13.1.32, 13.1.33 p.505). The boundedness of ϕ implies that c 2 = 0 and the boundary condition ϕ(r) = 1 gives

Lµ y (w) = r y n 2 M (k, µ, λy 2 ) M (k, µ, λr 2 ) . (28)
If we look at (27), the probabilistic definition of Lµ y (w), by the same argument as previously gives that Lµ y (w) can be extended to holomorphic function in the complex halfplane Re(w) > -(n-2) 2

8

. Since the Whittaker functions are well defined in this region the formula (28) is also satisfied in this region. As before we use the Laplace inverse formula and for c = - Theorem 7. For every x ∈ B r , x = 0 and y ∈ S n-1 r the Poisson kernel P r (x, y) is given by the following formula

2Γ n+1 2 π n-1 2 nr n-1 1 -|x| 2 1 -r 2 n-2 2 r |x| n 2 sin 3-n 2 ϕ 2πi c+i∞ c-i∞ M (k, µ, λy 2 ) M (k, µ, λr 2 ) B n (z)P 3-n 2 A(z)-1 2 (-cos ϕ)dz,
where ϕ is an angle between x and y, A(z) = 1 2 (n -2) 2 -8z and

B n (z) = Γ n-2 2 -A(z) Γ n-2 2 + A(z) 2 n+1 2 zΓ n-1 2 .

Appendix

For convenience of the Reader we collect here basic information about Bessel functions, hypergeometric functions and other special functions appearing throughout the paper. Mainly we follow the exposition given in [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] and [START_REF] Erdelyi | Higher Transcendental Functions[END_REF], where we refer the Reader for more details (see also [START_REF] Lebedev | Special functions and their applications[END_REF] and [START_REF] Gradstein | Table of integrals, series and products[END_REF]).

Bessel functions

The Bessel functions J ν (z) and Y ν (z) are independent solutions of the Bessel equation

z 2 y ′′ (z) + zy ′ (z) + (z 2 -ν 2 )y(z) = 0, ν ∈ R.
The Wronskian of the pair (J ν (z), Y ν (z)) is equal to 2/(πz) (see [START_REF] Lebedev | Special functions and their applications[END_REF], p.113). The derivatives of Bessel functions can be expressed by itself in the following way

J ′ ν (x) = J ν-1 (x) - ν x J ν (x), x > 0, (29)
Y ′ ν (x) = Y ν-1 (x) - ν x Y ν (x), x > 0. ( 30
)
For every ν > 0 we have (see [START_REF] Lebedev | Special functions and their applications[END_REF] 5.16 p.134-135).

J ν (x) ∼ x ν 2 ν Γ(ν + 1) , Y ν (x) ∼ - 2 ν Γ(ν) π 1 x ν , x → 0 + . (31) J ν (x) ∼ 2 xπ cos(x - νπ 2 - π 4 ) , Y α (x) ∼ 2 xπ sin(x - νπ 2 - π 4 ), x → ∞.(32)
The modified Bessel functions are independent solutions to the modified Bessel equation

z 2 y ′′ (z) + zy ′ (z) -(ν 2 + z 2 )y(z) = 0. ( 33
)
The following asymptotic expansion for the function K ν holds (see [START_REF] Gradstein | Table of integrals, series and products[END_REF] 8.451 ( 6))

K ν (z) = π 2z e -z (1 + E(z)), |E(z)| = O(|z| -1 ) as |z| → ∞. (34)
whenever | arg z| 3π/2 . The behavior of K ν near zero is described by (see [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF], 9.6.9)

K ν (z) ≈ 2 ν-1 Γ(ν) z ν , Re z > 0. ( 35
)
The connection between modified Bessel function of purely imaginary argument and Bessel functions is given by

K ν (ix) = - iπ 2 e -iνπ/2 (J ν (x) -iY ν (x)), x > 0. (36)
Finally, we recall the integral representation of K ν (see [START_REF] Gradstein | Table of integrals, series and products[END_REF], 8.432 ( 7))

K ϑ (z) = z ϑ 2 ∞ 0 exp - t + z 2 /t 2 t -ϑ-1 dt, z > 0, ϑ ∈ R, (37)
as well as the formula ([16] 6.561 formula 16, p. 676)

∞ 0 x µ K ν (x)dx = 2 µ-1 Γ 1 + µ + ν 2 Γ 1 + µ -ν 2 , µ + 1 > ν > 0. ( 38
)
For every ν > 1/2 we introduce the following function of two variables

g ν (x, t) = 1 t 2ν (x -1) J ν (t)Y ν (tx) -J ν (tx)Y ν (t) J 2 ν (t) + Y 2 ν (t) , x > 1, t > 0.
It is obvious that g ν is a continuous function on (1, ∞) × R + . However, the most crucial for the considerations given in Chapter 3 are the following asymptotic properties of g ν .

Lemma 8. Set x 0 1, t 0 > 0. Then lim (x,t)→(x 0 ,0) g ν (x, t) = π n-2 k=0 x k 0 2 2ν Γ(ν)Γ(ν + 1)x ν 0 , (39) lim
(x,t)→(1,t 0 ) g ν (x, t) = 2 πt 2ν 0 (J 2 ν (t 0 ) + Y 2 ν (t 0 )) . ( 40
)
Proof. From the Lagrange theorem, there exist θ 1 , θ 2 ∈ (1, x), depending on x and such that

J ν (t)Y ν (tx) -J ν (tx)Y ν (t) tx -t = J ν (t) Y ν (tx) -Y ν (t) tx -t -Y ν (t) J ν (tx) -J ν (t) tx -t = J ν (t)Y ′ ν (tθ 1 ) -J ′ ν (tθ 2 )Y ν (t).
Obviously, when x → 1 then θ 1 , θ 2 also tend to 1. Furthermore, since the Wronskian of (J ν (z), Y ν (z)) is 2/πz we get that lim

(x,t)→(1,t 0 ) J ν (t)Y ν (tx) -J ν (tx)Y ν (t) tx -t = J ν (t 0 )Y ′ ν (t 0 ) -J ′ ν (t 0 )Y ν (t 0 ) = 2 πt 0 ,
which proves (40). If we use the recurrent formulas for the Bessel function derivatives (29) and (30) we obtain that J

ν (t)Y ′ ν (tθ 1 ) -J ′ ν (tθ 2 )Y ν (t) is equal to J ν (t)Y ν-1 (tθ 1 ) - ν t J ν (t)Y ν (tθ 1 ) θ 1 - J ν (tθ 1 )Y ν (t) θ 2 -J ν-1 (tθ 2 )Y ν (t)
Multiplying the last expression by t, taking (x, t) → (1, 0) and using (31) it is easy to see that first two summands tends to zero and the last one tends to 2/π. Since, by (31), we have lim

t→0 + t 2ν (J 2 ν (t) + Y 2 ν (t)) = 2 2ν Γ 2 (ν)/π 2 and consequently lim (x,t)→(1,0) g ν (x, t) = π 2 2ν-1 Γ 2 (ν) ,
which is (39) for x 0 = 1. For x 0 > 1 relation (39) follows directly from (31) lim (x,t)→(x 0 ,0)

g ν (x, t) = π(x ν 0 -x -ν 0 ) ν2 2ν Γ 2 (ν) = π n-2 k=0 x k 0 2 2ν Γ(ν)Γ(ν + 1)x ν 0 .
Note that Lemma 8 together with (31) and (32) imply that the function g ν can be extended to the continuous function on [1, ∞) × [0, ∞) which is bounded whenever x is bounded, i.e. for every R > 1 there exists C(R) > 0 such that

|g ν (x, t)| C(R), (x, t) ∈ [1, R] × [0, ∞).
(41)

Hypergeometric Function and Legendre Functions

For γ = -1, -2, . . . the hypergeometric function is defined by

2 F 1 (α, β; γ; z) = ∞ k=0 (α) k (β k ) (γ k )k! z k , |z| < 1,
Here (a) k = Γ(a + k)/Γ(a).

The derivative of 2 F 1 is given by (see [START_REF] Erdelyi | Higher Transcendental Functions[END_REF]

and the following elementary relation holds (see [START_REF] Erdelyi | Higher Transcendental Functions[END_REF] vol.1 p.105 2.9(2))

2 F 1 (α, β; γ; z) = (1 -z) γ-α-β 2 F 1 (γ -α, γ -β; γ; z).

The Legendre functions are solutions of Legendre's differential equation ; a + 3/2; 1 x 2 ), x > 1, respectively. The functions P b a and Q b a are independent solutions of (46).

Skew-product representation of Brownian motion

We now introduce the spherical Brownian motion on the unit sphere S n-1 1 ⊆ R n as a diffusion on S n-1 1 with the generator being the one half of the Laplace-Beltrami operator ∆ S n-1 1 of the manifold S n-1

1

. It is well-known that the following formula holds

∆ S n-1 1 = (sin φ) 2-n ∂ ∂ φ [(sin φ) n-2 ∂ ∂ φ ] + (sin φ) -2 ∆ S n-2 1 ,
where the angle φ is the angle between the pole and the given point on the sphere and ∆ S 1 1 = ∂ 2 ∂ φ 2 . Now, if we consider the action of ∆ S n-1

1
on function depending only on φ, this reduces to the generator of the Legendre process LEG(d):

1 2 ∆ S n-1 1 = 1 2 (sin φ) 2-n ∂ ∂ φ [(sin φ) n-2 ∂ ∂ φ ] = 1 2 ∂ 2 ∂ φ 2 + n -2 2 cot φ ∂ ∂ φ
Changing variable cos φ = t we obtain

1 -t 2 2 ∂ 2 ∂ t 2 - n -1 2 t ∂ ∂ t .
We now invoque the classical skew-product representation of the n-dimensional Brownian motion (see, e.g. [START_REF] Itô | Diffusion Processes and Their Sample Paths[END_REF] 7.15) stating that it can be represented as the product of R (ν) = {R

∂ 2 ∂ r 2 + n -1 2r ∂ ∂ r , (49)
and independent spherical Brownian motion Θ = {Θ(t); t 0} on S n-1

1
with time changed according to the formula

A (ν) (t) = t 0 ds (R (ν) (s)) 2 .
Moreover, for x = 0, we introduce the process S = {S(t); t 0} defined by S(t) =

x,Θ(t)

|x| . The process S describes the cosine of the angle between the starting point x and the spherical Brownian motion Θ. Consequently, the cosine between the starting point x and W (t) is just S(A(t)). The skew-product representation and the previous given considerations imply that S is independent from the Bessel process R (ν) and the generator of S is given by

G = 1 -t 2 2 ∂ 2 ∂ t 2 - n -1 2 t ∂ ∂ t (50)
with a domain D G = {u ∈ C 2 [-1, 1]; u ′ (-1) = u ′ (1) = 0}. Three basic characteristics of the diffusion: the speed measure, the scale function and the killing measure are described by the following relations (see also [START_REF] Byczkowski | Hitting half-spaces by Bessel-Brownian diffusions[END_REF]) m(dx) = 2(1-x 2 ) (n-3)/2 dx, s ′ (x) = (1-x 2 ) (1-n)/2 , k(dx) = 0. Moreover, the points -1 and 1 are non-singular reflecting points. We denote by p S t (x, y) the transition density function with respect to the speed measure, i.e.

2 ∈

 2 dt , y ∈ (0, r].

  Lµ |x| (z)G z 1,x, y |x||y| dtdz.

Function 2 F 2 F 1

 221 1 is a solution of the hypergeometric equationz(1 -z)u ′′ (z) + [γ -(α + β + 1)z]u ′ (z) -αβu(z) = 0 (42)regular at z = 0. Whenever Re (γ -α -β) > 0 we have (see[START_REF] Erdelyi | Higher Transcendental Functions[END_REF] vol.1 p.104 2.8(46)) (α, β; γ; 1) = Γ(γ)Γ(γ -α -β) Γ(γ -α)Γ(γ -β) .

( 1 - 2 F 1 (

 121 z 2 )u ′′ (z) -2zu ′ (z) + [a(a + 1) -b 2 (1 -z 2 ) -1 ]u(z) = 0. (46)Making appropriate substitution it can be reduced to the hypergeometric equation (42) and consequently its solutions are given in terms of hypergeometric function. More precisely, the Legendre function of the first and second kind are defined by (see[START_REF] Erdelyi | Higher Transcendental Functions[END_REF] vol.1 p.122 3.2(3) and p.143 3.4(6)) -a, a + 1; 1 -b;

t

  ; t 0} -the Bessel process BES(n), ν = n/2 -1 with the generator 1 2

P

  x (S(t) ∈ A) = A p S t (x, y)m(dx), A ∈ B[-1, 1].(51)

  and (41) gives boundedness of g ν x n , s

	x n < R, we can apply Lebesgue dominated convergence theorem to get	|y|	for
	lim (xn,|y|)→(1,∞)

  vol.1 p.102 2.8(20))

	d dz	2 F 1 (α, β; γ; z) =	αβ γ	2 F 1 (α + 1, β + 1; γ + 1; z).

Research supported by Polish Ministry of Science and Higher Eduction grant N N201 3731 36 and l'Agence Nationale de la Recherche grant no. ANR-09-BLAN-0084-01.


Oštir Krištof

Nuninger Laure

Paleorelief detection and modellinga case of study in eastern Languedoc (France)

Introduction

The main objective of the paper is paleo digital elevation model (paleoDEM) detection and modelling, a multidisciplinary approach involving geodesy, archaeology, and paleoenvironmental studies. The digital elevation model (DEM) is an important data layer to understand settlement pattern and socioenvironmental contexts and evolution that archaeologists are using together with GIS analysis. It can, for example help to analyze viewsheds, to define travel time between places, to simulate ancient paths, etc. Considering a period of several millennia, some topographical areas can be considered without fundamental change. Within littoral zones or fluvial areas, however, the situation is different. Sedimentation and erosion change the relief permanently, due to alluvial and sea level dynamic. A nowadays-flat area could be hilly in prehistoric times and the application of present DEM for the analysis is therefore unsuitable, since it would distort the expected result completely. Our case of study area is situated in Languedoc (France), on the littoral plain of Mauguio between the cities of Montpellier to the west and Nîmes to the north (Fig. 1). This area covers around 250 km 2 . From an archaeological point of view, the zone is very well explored with almost 350 sites inhabited from Bronze Age to Middle Ages, showing also land division evidences obtained from aerial photo interpretation. In addition, in 1990s several geoarchaeological surveys were done in the main deltas of Lez and Vidourle rivers. Unfortunately the geo-archaeological coring, performed to understand paleochannel profiles in the Vidourle delta, is concentrated in three zones only. Considering a regional scale, only three average points could be used in further analyses. The geoarchaeological approach in itself brings indices about paleorelief but since the measurements are represented as points they are not adequate to produce a proper historical DEM. Additionally the approach is slow and relatively expensive. We have therefore used remote sensing, optical image processing and radar inter-ferometry, in combination with spatial analysis to detect the paleorelief.

The first step was to create an accurate present DEM with a resolution suitable for analyses. Afterwards, the DEM, as well as SPOT and Landsat satellite images, was processed with several filtering techniques to detect "paleofeatures" -evidences of paleorelief. Features were digitalised and their spatial distribution was statistically tested according to archaeological data distribution. At last, we compared and tried to qualify the paleofeatures with the results obtained by photo-interpretation, classification and simulation.

Paleorelief interpretation

Existing technology does not allow a direct detection of paleorelief. Nevertheless, satellite imagery and current digital elevation models can be used to detect the "paleofeatures" and to relate them to past relief.

The digital elevation model of the study area has been made as a combination (weighted sum) of all available sources [START_REF] Podobnikar | Data integration for the DTM production[END_REF]. Firstly, a DEM has been produced with radar interferometry (InSAR, [START_REF] Oštir | Producing digital elevation models with radar interferometry[END_REF]) from two ERS radar images, acquired on 1996-04-04 and 1996-04-05. The images were acquired in tandem mode by ERS-1 and ERS-2 with a time difference of only one day. The coherence of the image pair is very high in the whole study area, with exception of water bodies (sea, rivers, channels, etc.). Interferometric processing produced a digital elevation model with a resolution of 25 m. The vertical accuracy of the relief is approximately 8 m and the positional accuracy is almost 5 m (one fifth of pixel size). The model has been only moderately filtered and prepared for combination with other DEMs.

Additionally the following DEMs were used: IGN DEM (with the resolution of 50 m, provided by IGN), Aster DEM (30 m, made from Aster stereo satellite images), SRTM DEM (90 m, produced with radar interferometry from the Shuttle Radar Topography Mission data). These models have different resolution and different coverage of the area: from ~60% for Aster to almost 100% for InSAR and SRTM (both are produced with interferometry and contain some areas with missing data) and complete for IGN. Standard deviations of all DEMs for the study area are less than 10 m (Table 1); Aster and InSAR are systematically higher because they include vegetation cover.

The final (combined) model is a weighted sum of partial models [START_REF] Podobnikar | Data integration for the DTM production[END_REF]. It takes the rough shape of IGN DEM and adds details of others, particularly InSAR. Its resolution is 25 m and it is very close to the control points (-0.2 m, standard deviation 3.7 m, Table 1, Fig. 1). The model is rather homogenous and could be used for further analysis and paleorelief interpretation. Optical satellite images were an additional important interpretation source. A panchromatic image acquired by SPOT 4 in 1999-04-08 and a Landsat 7 ETM+ image from 2001-08-13 were used. Optical images (Landsat and SPOT) were enhanced prior interpretation. We have been trying to observe paleofeatures indirectly through detection of edges, particularly related to humidity and vegetation anomalies. Landsat with its seven bands gives much better spectral information than SPOT while SPOT has a better spatial resolution (SPOT 4 pan with 10 m resolution has been used only).

Simple and advanced edge detection has been performed, including Sobel and high-pass filtering supplemented with edge thresholding. Special attention has been given to filtering in the direction of the Sun to expose small variations in the direction of incidence radiation. The position of the Sun has been computed from the image acquisition time and the filter has been oriented in the selected direction.

To increase the spatial and preserve spectral resolution of Landsat resolution merge of panchromatic and selected multispectral channels has been performed [START_REF] Švab | High-resolution image fusion: methods to preserve spectral and spatial resolution[END_REF]. A similar procedure has been repeated with panchromatic SPOT and multispectral Landsat data. It has to be mentioned the huge difference between 15 m resolution of panchromatic Landsat and 10 m of SPOT. The latter is much sharper, what is probably the influence of different sensor types, simplifying the observation of anomalies.

The pre-processed DEM and satellite imagery enabled paleofeature delineation. Manual on screen digitalization has been performed and all possible features have been included. The following data has been mainly used to digitalize potential features (Fig. 2):

• panchromatic and edge filtered SPOT,

• multispectral and pan-sharpened Landsat, edge filtered panchromatic Landsat, and • shaded and edge filtered combined DEM. All possible anomalies have been considered and the features have been over-digitalized and later cleaned (Fig. 3). In the iterative procedure different images were used to digitize and recognize features. Maps in the scale 1:25.000 were used to verify all the vectors by eliminating existent relief characteristics (roads, buildings, etc.) and water courses (rivers, channels, etc.). SPOT and shaded combined DEM were most useful in feature delineation and all other layers were used only to update or check the results. Much better results are expected with a digital elevation model produced from lidar imaging [START_REF] Challis | Airborne laser altimetry in alluviated landscapes[END_REF]) that is scheduled for December 2006.

Thematical identification

To test if the detected paleofeatures can be linked to real paleorelief, an analysis of archaeological sites in the area has been performed (Fig. 3). The distribution of sites in relation to the features (buffer zones, site proximity analysis) has been observed for Prehistory, Roman period and medieval period. The sites were compared to random point distribution. For all known sites the distance from digitalized paleofeatures has been determined. Only sites that are less than certain distance (2.5 km) from the features have been considered. It has been statistically proven that the average distance of sites for all periods, especially for Prehistory, is smaller than for random points (Fig. 4a). For Prehistory the average distance is only approximately 450 m, while for random points it is more than 800 m, almost twice as large. One can also observe that the majority of the sites lie less then 500 m from the detected features (Fig. 4b). Despite of the connection between evidence of paleofeatures and archaeological sites, we do not know their exact nature. Three different approaches were performed to define types of paleofeatures. First, the digitized layer was overlaid with geology, pedology, cadastre with contour lines at 1 m interval, and historical aerial photography (from allied missions in 1944). By observation of small areas chosen as test sites, several paleofeatures were identified as paleochannels, edges of natural paleodepressions, ancient roads, terraces, dried valleys etc. Nevertheless, it was noticed that some of the features are still related to the artefacts observed on the IGN DEM (stages caused by contour line interpolation) and some features remain unresolved.

The second approach was related to paleocoastline or sea border (Fig. 3). To ascertain if any of the features can be interpreted as coastline, we used simulations taking into account geoarchaeological information. A small number of cores with an acceptable spatial distribution permitted the generation of a very basic model with an approximated subsiding of the alluvial floor [START_REF] Nuninger | Contribution à la modélisation des paléo-reliefs de la plaine littorale de l'étang de Mauguio (Languedoc, France): premières approches par télédétection. In Temps et espaces de l'Homme en société. Analyses et modèles spatiaux en archéologie[END_REF]. The result shows a rather good general connection and by adding archaeological sites distribution, we were able to notice a relation between chronology and local variation. In fact, reconstitution of paleo-shore during antiquity suggests that the sea-level was approximately 1 m above the present level. However, the archaeological site locations suggest less elevated paleo-shore which is better connected to detected paleofeatures. Even if the scale of layers is different, we have to consider the local variation of sea level, due to hydro sedimentary dynamics. Thus, features detected (Fig. 5) show two steps of shore evolution from Bronze Age to Antiquity. In this case, a shore seems to have moved back until Antiquity and Bronze Age settlements were probably under water during that time.

Fig. 5. Evolution of paleo sea border from Bronze Age to Antiquity

At last, in order to automatically classify the detected evidence of paleorelief, a multivariate analysis (factor analysis with classification 1 ) was processed [START_REF] Coulon | Contribution à la recherche des paléoreliefs dans la plaine littorale de Mauguio : analyse spatiale et statistique des données pour l'élaboration d'une typologie des indices de paléoreliefs[END_REF]. First, the detected features were rasterised and filtered to delete each cells corresponding to DEM artificial stages and thus to avoid noise due to IGN DEM artefacts. Remaining cells were grouped as new entities (1020) according to neighbourhood and their thematic homogeneity (mostly soils). Then, each entity was described by several criteria:

• Archaeological density was calculated as kernel density taking into account duration of occupation (settlements). E.g. 1 settlement occupied during period 100 to 199 AD and period 200 to 299 AD counts as 2 occupations for density calculation. The factor analysis results grouped the detected paleofeatures in 7 classes for 1020 entities. We will focus to the most significant only. In spite of containing dip 6% of all entities, class 1 is very interesting. It groups features with a strong archaeological density, a very strong compactness, a principal aspect to the east and a concentration on the oldest soils (61% "fersiallitiques").

Features from this class can be interpreted as micro-knoll with a high probability. Class 7 and 8 (14% and 6%) can be interpreted as class 1 with lower probability. On the contrary within class 6 (12%), one can observe a weak archaeological density (78% from null to mean value), a weak compactness (88%), principal aspect to north or north-east (67%) and most recent soils (70%). Features from this class are probably in connection with alluvial terrace or dried valley. At last, the class 5 (10%) is represented mostly by no archaeological occupation (82%), a very strong compactness (70%) and most recent soil (85%), especially soils with hydromorphological characteristics (41%). This class should bring together small depressions and paleomeanders.

The multivariate analysis is a first stage to sort automatically the detected paleofeatures. Nevertheless, such assumptions should be verified by control on the original data (IGN map, cadastre, aerial photograph) and in the field. Then, if the classification can be validated, criteria should be integrated directly within the process of detection.

Conclusions

With the case study we have proven that remote sensing can help in the detection of paleorelief. While it is not possible to detect past environment directly, several data layers -in our case most successfully the enhanced panchromatic SPOT 4 and the combined DEM -can be used to observe paleofeatures and to relate them to relief. The described approach will not replace precise work of archaeological and paleoenvironmental investigation, but can simplify the interpretation.

The results obtained present a first step and should be improved with further processing and additional data. The project will continue with advanced processing of high resolution satellite imagery and aerial photography. Furthermore, lidar scanning is scheduled to provide better elevation data in the flat areas.

Fig. 1 .

 1 Fig. 1. Location of the study area.

Fig. 2 .

 2 Fig. 2. Enhanced SPOT (a) and Landsat (b) satellite images were used together with DEM (c) to detect paleofeatures.

Fig. 3 .

 3 Fig. 3. Digitized and cleaned paleofetures are shown in relation to archaeological sites.

Fig. 4 .

 4 Fig. 4. Average distance to of archaeological sites to paleofeatures for different periods (a) and percentage of the sites within certain distance (b) compared with random points.

Table 1 .

 1 Statistical comparison of available DEMs and the final combined DEM

		Control	IGN	Aster	SRTM	InSAR	Combined
		points	DEM 50	DEM	DEM	DEM	DEM
	Number	110	110	67	108	105	110
	Percentage	100%	100%	61%	98%	95%	100%
	Average	81.4	80.9	82.4	81.6	90.0	81.2
	Min Diff		-13	-13	-25	-15	-15
	M a x D i f f		8	2 4	1 1	2 9	8
	Average Diff		-0.5	9.1	-1.2	4.7	-0.2
	StdDev Diff		3.7	8.4	5.0	7.8	3.7

The factor analysis and classification used are an « Analyse Factorielle des Correspondances -AFC » with a « Classification Ascendante Hiérarchique -CAH »

Acknowledgments

The project was supported by CNRS within the framework of a four months associate researcher position for Krištof Oštir (Laboratory ThéMA, UMR 6049, Besançon, France). We would like to thank Jean-François Berger, Cécile Jung and M. Coulon for their collaboration and access to their data. Ziga Kokalj helped with figures and text corrections.


Julien Toulouse

Wuming Zhu

Andreas Savin

Georg Jansen

János G Ángyán

Closed-shell ring coupled cluster doubles theory with range separation applied on weak intermolecular interactions

We explore different variants of the random phase approximation (RPA) to the correlation energy derived from closed-shell ring-diagram approximations to coupled cluster doubles theory. We implement these variants in range-separated density-functional theory, i.e. by combining the long-range random phase approximations with short-range density-functional approximations. We perform tests on the rare-gas dimers He 2 , Ne 2 , and Ar 2 , and on the weakly interacting molecular complexes of the S22 set of Jurečka et al. [Phys. Chem. Chem. Phys. 8, 1985 (2006)]. The two best variants correspond to the ones originally proposed by Szabo and Ostlund [J. Chem. Phys. 67, 4351 (1977)]. With range separation, they reach mean absolute errors on the equilibrium interaction energies of the S22 set of about 0.4 kcal/mol, corresponding to mean absolute percentage errors of about 4%, with the aug-cc-pVDZ basis set.

I. INTRODUCTION

In the last decade, there has been a revived interest in the random phase approximation (RPA) and other related approximations for calculating the electron correlation energy of atomic, molecular and solid-state systems . One particularly appealing feature of RPA is its correct description of dispersion forces at large separation [42][43][44]. However, RPA is a poor approximation to short-range correlations [1], and, in a Gaussian localized basis, RPA calculations have a slow convergence with respect to the basis size [2]. A promising strategy is thus to combine a long-range RPA-type approximation with a short-range density-functional approximation [15-17, 28, 31, 33], hence avoiding the inaccurate description and slow basis-set convergence of short-range correlations in RPA.

Among the different formulations of RPA, the one based on a ring-diagram approximation to coupled cluster doubles (CCD) theory [14,[45][46][47] is particularly attractive since it avoids the numerical integration over the adiabatic connection, and in principle is amenable to a fast algorithm [14]. However, due to the fact that the ring approximation breaks the antisymmetry property of the coupled-cluster amplitudes, several non-equivalent variants of ring CCD can be constructed, especially when the exchange terms are included. In this paper, we explore these various ring CCD variants for closed-shell systems, and show that some of them correspond to the RPA correlation energy expressions originally proposed by Szabo and Ostlund [48,49]. We apply these closed-shell ring CCD variants in the context of range-separated densityfunctional theory, and test them on rare-gas dimers and on the weakly interacting molecular complexes of the S22 set of Jurečka et al. [50]. * Electronic address: julien.toulouse@upmc.fr

II. THEORY

We first show how to rigorously combine a long-range CCD calculation with a short-range density functional (for details on range-separated density-functional theory, see e.g. Refs. 33,51,52). We start from a self-consistent rangeseparated hybrid (RSH) calculation [52]

E RSH = min Φ { Φ| T + Vext + Ŵlr ee |Φ + E sr Hxc [n Φ ]}, (1)
where T is the kinetic energy operator, Vext is the external potential operator (e.g., nuclei-electron interaction), Ŵlr ee is a long-range electron-electron interaction operator, E sr Hxc [n] is the associated short-range Hartree-exchange-correlation density functional, and Φ is a single-determinant wave function with density n Φ . The long-range interaction is constructed with the error function, w lr ee (r) = erf(µr)/r, where µ is a parameter whose inverse gives the range of the separation. The minimizing RSH single-determinant wave function is denoted by Φ 0 and its associated (approximate) density by n 0 . In principle, the exact ground-state energy can be obtained from the RSH energy by adding the long-range correlation energy

E lr c E = E RSH + E lr c . (2)
Several formally exact expressions can be derived for E lr c . The one that is most convenient for applying coupled-cluster theory is 3), the last two terms are the variation of the energy functional, ∆E sr Hxc = E sr Hxc [n] -E sr Hxc [n 0 ], and the variation of the associated potential expectation value due to the variation of the density from the RSH one to the exact one, ∆n = nn 0 . The contribution of these last two terms is expected to be small since it is of second order in ∆n

E lr c = Ψ lr | Ĥlr [n]|Ψ lr -Φ 0 | Ĥlr [n]|Φ 0 +∆E sr Hxc -v sr Hxc [n](r) ∆n(r)dr, (3)
∆E sr Hxc -v sr Hxc [n](r) ∆n(r)dr = - 1 2 δE sr Hxc [n 0 ] δn(r)δn(r ′ ) ∆n(r)∆n(r ′ )drdr ′ + O(∆n 3 ). ( 4
)
Using a spin-unrestricted CCD ansatz (see Appendix A for a review of CCD theory) for the long-range wave function, |Ψ lr CCD = exp T2 |Φ 0 , where T2 = (1/4) i jab (t ab i j ) lr â † a âi â † b â j is the cluster operator for double excitations written in terms of the long-range amplitudes (t ab i j ) lr , and occupied (i, j) and virtual (a, b) RSH spin-orbital creation and annihilation operators, we approximate the long-range correlation energy as

E lr c,CCD = Φ 0 | Ĥlr [n 0 ]|Ψ lr CCD -Φ 0 | Ĥlr [n 0 ]|Φ 0 . (5)
In Eq. ( 5), the variation of the density has been neglected, n ≈ n 0 (and thus the contribution of Eq. ( 4) vanishes), which seems appropriate if we define the coupled-cluster density as the projected one, Φ 0 |n(r)|Ψ lr CCD = Φ 0 |n(r)|Φ 0 = n 0 (r), which does not vary since the CCD wave function does not contain single excitations. The long-range correlation energy can be calculated as, for real spin orbitals,

E lr c,CCD = 1 4 tr B lr T lr = 1 2 tr K lr T lr , (6)
where B lr ia, jb = ab||i j lr and K lr ia, jb = ab|i j lr are the matrices of antisymmetrized and non-antisymmetrized two-electron integrals with long-range interaction w lr ee (r), respectively, and T lr ia, jb = (t ab i j ) lr is the amplitude matrix. The second equality in Eq. ( 6) is due to the antisymmetry property of the coupledcluster amplitudes T lr ia, jb with respect to the exchange of the indices i and j. These amplitudes can be determined by the usual coupled-cluster equations, replacing the normal Hamiltonian by the long-range one Ĥlr [n 0 ], which amounts to using the RSH orbital eigenvalues and the long-range two-electron integrals. The present range-separated CCD method can be seen a special case of the more general range-separated coupledcluster approach of Goll et al. [53] which also includes single excitations and possibly perturbative triples.

We now consider the ring-diagram approximation for closed-shell systems. A number of closed-shell ring CCD variants can be defined. In the ring approximation without exchange terms, the direct RPA (dRPA, also sometimes referred to as RPA or time-dependent Hartree) amplitudes for spin-singlet excitations, 1 T lr dRPA , are obtained by the following Riccati equation [14] 1 K lr + 1 L lr 1 T lr dRPA + 1 T lr dRPA 1 L lr + 1 T lr dRPA 1 K lr 1 T lr dRPA = 0, (7) with the spin-adapted matrices 1 K lr ia, jb = 2 ab|i j lr and 1 L lr ia, jb = ∆ǫ ia, jb + 1 K lr ia, jb , where ∆ǫ ia, jb = (ǫ aǫ i )δ i j δ ab is the matrix of the RSH orbital eigenvalue differences (i, j and a, b refer now to occupied and virtual spatial orbitals, respectively). Contracting the dRPA amplitudes with the nonantisymmetrized two-electron integrals 1 K lr gives the dRPA long-range correlation energy (also referred to as dRPA-I in Ref. 54)

E lr c,dRPA = 1 2 tr 1 K lr 1 T lr dRPA . (8)
Contracting the dRPA amplitudes with the spin-singletadapted antisymmetrized two-electron integrals 1 B lr ia, jb = 2 ab|i j lr -ab| ji lr gives the dRPA+SOSEX (or just SOSEX for short) long-range correlation energy [23,28]

E lr c,SOSEX = 1 2 tr 1 B lr 1 T lr dRPA . (9)
Similarly, in the ring approximation with exchange terms, that we will refer to as RPAx (also sometimes referred to as RPA or time-dependent Hartree-Fock), the singlet and triplet amplitudes 1 T lr RPAx and 3 T lr RPAx are obtained by the equations

1 B lr + 1 A lr 1 T lr RPAx + 1 T lr RPAx 1 A lr + 1 T lr RPAx 1 B lr 1 T lr RPAx = 0, (10)
and

3 B lr + 3 A lr 3 T lr RPAx + 3 T lr RPAx 3 A lr + 3 T lr RPAx 3 B lr 3 T lr RPAx = 0, (11)
where 1 A lr ia, jb = ∆ǫ ia, jb + 2 ib|a j lr -ib| ja lr , 3 A lr ia, jb = ∆ǫ ia, jb -ib| ja lr , and 3 B lr ia, jb = -ab| ji lr . Using these amplitudes in the CCD correlation energy expression of Eq. ( A9) gives what we call the RPAx-II long-range correlation energy (see, also, Refs. 48, 49, 54, 56)

E lr c,RPAx-II = 1 4 tr 1 B lr 1 T lr RPAx + 3 3 B lr 3 T lr RPAx , (12)
which is equivalent to the plasmon formula expression of McLachlan and Ball [57].

E lr c,RPAx-SO2 = 1 2 tr 1 K lr 1 T lr RPAx . (13)
Equations ( 12) and ( 13 This last variant is the one preferred by Szabo and Ostlund because in a supermolecule approach it consistently gives a dispersion coefficient C 6 identical to the one given by the Casimir-Polder formula applied with the RPAx polarizabilities of the fragments, which is not the case for the other variants RPAx-II and RPAx-SO2. On the other hand, among the three RPAx methods proposed here, RPAx-SO2 has the advantage of involving only singlet excitations and thus is not subject to triplet instabilities. The RPAx method of Refs. 15, 33, that we will rename RPAx-I here, is yet another alternative correlation energy expression that involves only singlet excitations, but for which, as far as we know, the numerical integration over the adiabatic connection cannot be avoided (although in practice a single-point quadrature works well [31]). It can be shown that the SOSEX, RPAx-I, RPAx-II, RPAx-SO1 and RPAx-SO2 correlation energies all correctly reduce to the MP2 correlation energy at second order in the electron-electron interaction, but dRPA does not. Finally, we note that an another RPAx correlation energy variant first proposed by Fukuda et al. [61] and defined as 2E c,RPAx-II -E c,MP2 has also been discussed in the literature [41,[47][48][49]. It obviously correctly reduces to the MP2 correlation energy at second order, but numerical experience [41] shows that this variant gives very inaccurate correlation energies. [START_REF] Werner | Molpro, version 2008.2, a package of ab initio programs[END_REF], implementing equations ( 7)- (14). We first perform a self-consistent RSH calculation with the shortrange Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional of Ref. 64 and add the long-range RPA correlation energies calculated with RSH orbitals. The range separation parameter is taken at µ = 0.5 bohr -1 , according to previous studies [65], without trying to readjust it. For the rare-gas dimers, we also carry out full-range RPA calculations using PBE orbitals [66] for comparison. The Riccati equation ( 7) is solved by decomposing the matrix 1 L lr into diagonal and off-diagonal parts and iteratively extracting 1 T lr dRPA from its product with the diagonal part and updating it in the other terms, and similarly for Eqs. ( 10) and (11). For RPAx-I calculations, the adiabatic-connection integration is performed by a 8-point Gauss-Legendre quadrature for the rare-gas dimers, and by a single-point quadrature [Eq. ( 14) of Ref. 31] for the S22 set. We use the correlation-consistent basis sets of Dunning [67,68]. Core electrons are kept frozen (i.e. only excitations of valence electrons are considered). Basis set superposition error (BSSE) is removed by the counterpoise method. The geometries of the complexes of the S22 set are taken from Ref. 50. The geometries of the isolated monomers are fixed to those in the dimers, thus the so-called monomer deformation energy is not included in the interaction energy. For each method, mean error (ME), mean absolute error (MAE) and mean absolute percentage error (MA%E) are given using as a reference the CCSD(T) values extrapolated to the complete basis set (CBS) limit of Takatani et al. [62].

E lr c,RPAx-SO1 = 1 2 tr 1 B lr 1 T lr RPAx -3 T lr RPAx . (14)
In our present, most basic implementation, the computational cost of all the RPA methods used here formally scales as N 3 v N 3 o for large basis sets, where N v and N o are the numbers of virtual and occupied orbitals, respectively. The computational cost of the CCD (or CCSD) method without the ring approximation is higher and it scales as N 4 v N 2 o for large basis sets [69]. Of course, far better scalings should be obtained by using integral-direct methods and resolution-of-identity/Choleskydecomposition techniques [14].

IV. RESULTS

The interaction energy curves of He 2 , Ne 2 , and Ar 2 calculated by the full-range and range-separated RPA methods are compared in Fig. 1. We use the large aug-cc-pV6Z basis set to ensure that the full-range calculations are converged. Full-range dRPA and SOSEX strongly underestimate the interaction energies, while full-range RPAx-II and RPAx-SO1 strongly overestimate them. The best full-range methods are RPAx-I and RPAx-SO2, which is in agreement with the recent study of Heßelmann [41,[START_REF]Ref. 41, AC-RPA, NRPA1, NRPA3, and NRPA4 refer to what we call here RPAx-I, RPAx-II, RPA-SO2, and RPA-SO1, respectively[END_REF]. In passing, we note that the full-range RPAx-I method better performs for Ne 2 and Ar 2 when using PBE orbitals than when using HF orbitals, as done in Ref. 33. Range separation greatly improves the accuracy of all the RPA variants. However, range-separated dRPA and SOSEX still underestimate the interaction energies, and range-separated RPAx-II significantly overestimates the interaction energy of Ar 2 . Range-separated RPAx-I, RPAx-SO1, and RPAx-SO2 give the most reasonable interaction energy curves.

The interaction energies for the complexes of the S22 set calculated with the range-separated RPA methods with the aug-cc-pVDZ basis set are given in Table I. For comparison, range-separated CCD results (without the ring approximation) are also reported. Although the aug-cc-pVDZ basis set may appear small, range-separated RPA methods are weakly dependent on the basis size [15,33], and indeed it was estimated in Ref. 31 that when going from the aug-cc-pVDZ to the augcc-pVTZ basis set the range-separated RPAx-I interaction energies of the S22 set are lower by at most 7%, and the corresponding total MA%E decreases by less than 2%. Therefore, we believe that the aug-cc-pVDZ basis set is sufficient to compare the different range-separated RPA methods.

The S22 set includes seven hydrogen-bonded complexes (HB7 subset), eight weakly interacting complexes with predominant dispersion contributions (WI8 subset), and seven mixed complexes featuring also multipole interactions (MI7 subset). The trends are quite different for the HB7 subset on the one hand, and the WI8 and MI7 subsets on the other hand. It was previously argued that the general overestimation of the interaction energies of hydrogen-bonded complexes is due to the approximate short-range density functional [31,[START_REF] Goll | [END_REF]. The fact that dRPA and SOSEX give the smallest MAEs for the HB7 subset is thus not believed to be significant but rather due to a compensation of errors between an underestimated longrange contribution and an overestimated short-range contribution. This is corroborated by the relatively large overestimation of the interaction energies of this subset by rangeseparated CCD which should most accurately describe the long-range correlation energies. We will thus focus our analysis on the WI8 and MI7 subsets.

For the WI8 and MI7 subsets, dRPA gives largely underestimated interaction energies, with MA%Es of 34.3% and 15.8%, respectively. SOSEX barely improves dRPA with MA%Es of 31.2% and 13.5%, respectively. This may not be surprising since, in the limit of large separation, SOSEX only adds exponentially decaying exchange interactions between the monomers, but does not change the coupled-cluster amplitudes and thus does not change the polarizabilities of the monomers. The RPAx-I method of Refs. 15, 33, which incorporates exchange effects in the monomers, greatly reduces the underestimation of the interaction energies, with MA%Es of 13.9% and 5.0%, respectively. The RPAx-II variant, which may be seen as the most straightforward way of defining a closed-shell ring CCD with exchange terms, is disappointingly inaccurate. It overestimates the interaction energies by about the same amount that dRPA underestimates them. Finally, the two variants RPAx-SO1 and RPAx-SO2 give remarkably accurate interaction energies, with MA%Es of 3.1% and 2.7% for RPAx-SO1, and 2.5% and 2.2% for RPAx-SO2. They are globally as accurate as range-separated CCD without the ring approximation. However, it must be noted that RPAx-SO1 and RPAx-SO2 tend to overestimate dispersion energies, while RPAx-I underestimates them. Therefore, increasing the basis size will likely increase the MA%Es of RPAx-SO1 and RPAx-SO2, while it will decrease the MA%E of RPAx-I.

V. CONCLUSION

We have studied various RPA variants that can be cast in the form of closed-shell ring CCD approximations. We have tested these variants with range separation, i.e. by combining a long-range RPA-type approximation with a shortrange density-functional approximation, on rare-gas dimers and on the weakly interacting complexes of the S22 set. Among all these variants, the ones first proposed by Szabo and Ostlund [48,49], called here RPAx-SO1 [Eq. ( 14)] and RPAx-SO2 [Eq. ( 13)], give the most accurate dispersion energies. The other variants tend to either strongly underestimate (dRPA and SOSEX) or strongly overestimate (RPAx-II) the interaction energies. For comparison, we have also reported results from the RPAx-I method of Refs. 15, 33, which is not based on a ring CCD approximation but on the adiabatic connection formula, and which gives reasonable interaction energies as well. From a practical point of view, RPAx-SO2 appears to be the most convenient variant since, contrary to RPAx-I, it does not use any numerical adiabatic-connection integration and, contrary to RPAx-SO1, it involves only singlet excitations and is thus not subject to triplet instabilities.

CCD correlation energy in spin-orbital basis

The spin-unrestricted CCD wave function ansatz is

|Ψ CCD = exp T2 |Φ , ( A1
)
where |Φ is a single-determinant reference wave function, and T2 is the cluster operator for double excitations which is written in a spin-orbital basis as

T2 = 1 4 i jab t ab i j â † a âi â † b â j , (A2)
where i, j and a, b refer to occupied and virtual spin-orbitals, respectively, and the amplitudes t ab i j must be antisymmetric with respect to any exchange of two indices: t ab i j = -t ab ji = -t ba i j = t ba ji . The CCD correlation energy is obtained by the transition formula

E CCD c = Φ| Ĥ|Ψ CCD -Φ| Ĥ|Φ = Φ| Ĥ T2 |Φ = 1 4 i jab ab||i j t ab i j = 1 4 tr [BT] , (A3)
where B ia, jb = ab||i j are the antisymmetrized two-electron integrals over real spin orbitals and T ia, jb = t ab i j is the amplitude matrix. Using the antisymmetry of the amplitudes, the CCD correlation energy can also be written as

E CCD c = 1 2 i jab ab|i j t ab i j = 1 2 tr [KT] , (A4)
where K ia, jb = ab|i j are the two-electron integrals.

CCD correlation energy in spatial-orbital basis for closed-shell systems a. Expression in terms of the singlet and triplet amplitudes

For spin-restricted closed-shell calculations, all the matrices in the spin-orbital excitation basis encountered so far (e.g., A,B,K,T) have the following spin block structure

C =               C ↑↑,↑↑ C ↑↑,↓↓ 0 0 C ↓↓,↑↑ C ↓↓,↓↓ 0 0 0 0 C ↑↓,↑↓ C ↑↓,↓↑ 0 0 C ↓↑,↑↓ C ↓↑,↓↑               , (A5)
and can be brought to a block-diagonal spin-adapted matrix C = U T C U by the orthogonal transformation

U = 1 √ 2               1 1 0 0 1 -1 0 0 0 0 1 1 0 0 1 -1              
In practice, the CCD correlation energy is normally calculated starting from the spin-restricted closed-shell CCD wave function ansatz

|Ψ CCD = exp R T2 |Φ , (A12)
where the restricted cluster operator R T2 is written in a spatialorbital basis as

R T2 = 1 2 i jab R t ab i j Êai Êb j , (A13)
where Êai = â † a↑ âi↑ + â † a↓ âi↓ is the singlet excitation operator and R t ab i j are the restricted amplitudes which must be symmetric with respect to the exchange of both i, j and a, b, i.e., R t ab i j = R t ba ji , but not antisymmetric with respect to the exchange of only two indices. The CCD correlation energy is obtained by the transition formula

E CCD c = Φ| Ĥ|Ψ CCD -Φ| Ĥ|Φ = Φ| Ĥ R T2 |Φ = i jab (2 ab|i j -ab| ji ) R t ab i j = tr 1 B R T , ( A14
)
where R T ia, jb = R t ab i j .

c. Expression in terms of the singlet and triplet restricted amplitudes

Another equivalent correlation energy expression can be obtained by decomposing the restricted amplitudes into spinsinglet and spin-triplet components. Indeed, the restricted cluster operator can be decomposed as (see, e.g., Ref. 73, 74) R T2 = 1 2 i jab 1,R t ab i j Ŝ aib j + 3,R t ab i j Taib j , (A15) where 1,R t ab i j are singlet restricted amplitudes 1,R t ab i j = R t ab ji + R t ab i j , (A16) which are totally symmetric (i.e., 1,R t ab i j = 1,R t ab ji = 1,R t ba i j = 1,R t ba ji ), and 3,R t ab i j are the triplet restricted amplitudes 3,R t ab i j = R t ab ji -R t ab i j , (A17) which are totally antisymmetric (i.e., 3,R t ab i j = -3,R t ab ji = -3,R t ba i j = 3,R t ba ji ). In Eq. (A15), Ŝ aib j is the singlet doubleexcitation operator Using the symmetry properties of 1,R t ab i j and 3,R t ab i j , it is easy to check that Eqs. (A13) and (A15) are equivalent. Combining Eqs. (A16) and (A17) leads to the decomposition of the restricted amplitudes into spin components R t ab i j = 1 2 1,R t ab i j -3,R t ab i j , (A22) and the CCD correlation energy [Eq. (A14)] can thus be written as

E CCD c = 1

  where Ψ lr is the ground-state wave function of the long-range interacting Hamiltonian Ĥlr [n] = T + Ŵlr ee + Vext + Vsr Hxc [n] with the short-range Hartree-exchange-correlation potential operator Vsr Hxc [n] = v sr Hxc [n](r) n(r)dr written with the density operator n(r) and v sr Hxc [n](r) = δE sr Hxc [n]/δn(r). The long-range wave function Ψ lr is associated with the exact density n. In Eq. (

  ) are not equivalent because the ring approximation does not preserve the antisymmetry of the amplitudes with respect to the exchange of two spin-orbital indices. Using the same amplitudes in place of the singlet and triplet restricted amplitudes in the CCD correlation energy expression of Eq. (A23) gives another alternative RPAx correlation energy corresponding to the first approximation proposed bySzabo and Ostlund [Eq. (3.20) of Ref. 49, or Eq. (17) of Ref. 48] which is an alternative zeroth iteration of the selfconsistent RPA scheme

2 FIG. 1 :

 21 FIG. 1: (Color online) Interaction energy curves of He 2 , Ne 2 , and Ar 2 calculated by the full-range (left) and range-separated (right) RPA methods with the aug-cc-pV6Z basis set. Cubic splines are used to interpolate between the calculated points. The accurate curves are from Ref. 55.

2 tr 1

 1 B 1,R T -3,R T . (A23)This corresponds to the definition of singlet and triplet contributions to the correlation energy,E S c = (1/2)tr[ 1 B 1,R T] and E T c = -(1/2)tr[ 1 B 3,R T]. By using the symmetry properties of 1,R T and 3,R T, one can show that they are equivalent to the more usual expressions in terms of the restricted amplitudes (see, e.g., Ref. 75): E S c = (1/4)tr[( 1 B -3 3 B) R T] and E T c = (3/4)tr[( 1 B + 3 B) R T].

TABLE I :

 I Interaction energies (in kcal/mol) for the complexes of the S22 set from the range-separated RPA methods with the aug-cc-pVDZ basis set. For comparison, range-separated CCD results (without the ring approximation) are also reported. The geometries of complexes are taken from Ref. 50 and the reference interaction energies in the rightmost column are taken as the CCSD(T)/CBS estimates of Ref.62. Mean errors (ME), mean absolute errors (MAE) and mean absolute percentage errors (MA%E) are given.

	No. Complex	dRPA SOSEX RPAx-I RPAx-II RPAx-SO1 RPAx-SO2 CCD Reference
	Hydrogen-bonded complexes (HB7)
	1 (NH 3 ) 2	-2.87	-2.92 -3.07	-3.26	-3.20	-3.18 -3.20	-3.17
	2 (H 2 O) 2	-5.16	-5.23 -5.33	-5.42	-5.40	-5.39 -5.41	-5.02
	3 Formic acid dimer	-20.30 -20.55 -20.81 -20.98	-20.94	-20.98 -20.94	-18.80
	4 Formamide dimer	-16.51 -16.68 -17.03 -17.48	-17.32	-17.27 -17.25	-16.12
	5 Uracil dimer C 2h	-21.03 -21.36 -21.80 -22.58	-22.04	-22.15 -22.00	-20.69
	6 2-pyridoxine/2-aminopyridine -17.07 -17.28 -17.81 -18.89	-18.25	-18.20 -18.08	-17.00
	7 Adenine/thymine WC	-16.53 -16.73 -17.29 -18.25	-17.81	-17.69 -17.62	-16.74
	ME	-0.28	-0.46 -0.80	-1.33	-1.06	-1.05 -0.99	0.00
	MAE	0.42	0.53	0.83	1.33	1.06	1.05 0.99	0.00
	MA%E	3.7%	4.3% 5.6%	8.6%	6.8%	6.6% 6.4%	0.0%
	Complexes with predominant dispersion contribution (WI8)
	8 (CH 4 ) 2	-0.30	-0.31 -0.42	-0.56	-0.53	-0.51 -0.51	-0.53
	9 (C 2 H 4 ) 2	-0.97	-1.02 -1.28	-1.66	-1.52	-1.47 -1.45	-1.50
	10 Benzene/CH 4	-0.92	-0.98 -1.23	-1.75	-1.47	-1.43 -1.40	-1.45
	11 Benzene dimer C 2h	-1.27	-1.38 -2.05	-4.28	-2.72	-2.61 -2.40	-2.62
	12 Pyrazine dimer	-2.99	-3.10 -3.78	-6.12	-4.49	-4.34 -4.14	-4.20
	13 Uracil dimer C 2	-8.22	-8.46 -9.38 -11.93	-10.25	-10.13 -9.94	-9.74
	14 Indole/benzene	-2.58	-2.75 -3.70	-7.12	-4.64	-4.48 -4.17	-4.59
	15 Adenine/thymine stack	-9.38	-9.68 -10.97 -15.14	-12.23	-12.02 -11.72	-11.66
	ME	1.21	1.08	0.43	-1.53	-0.20	-0.09 0.07	0.00
	MAE	1.21	1.08	0.43	1.53	0.20	0.13 0.13	0.00
	MA%E	34.3% 31.2% 13.9% 31.7%	3.1%	2.5% 3.8%	0.0%
	Mixed complexes (MI7)
	16 Ethene/ethyne	-1.31	-1.36 -1.48	-1.67	-1.58	-1.57 -1.55	-1.51
	17 Benzene/H 2 O	-2.90	-2.96 -3.16	-3.52	-3.34	-3.30 -3.29	-3.29
	18 Benzene/NH 3	-1.83	-1.88 -2.11	-2.57	-2.33	-2.29 -2.27	-2.32
	19 Benzene/HCN	-4.20	-4.31 -4.54	-4.98	-4.72	-4.71 -4.65	-4.55
	20 Benzene dimer C 2v	-1.92	-2.00 -2.39	-3.40	-2.77	-2.70 -2.61	-2.71
	21 Indole/benzene T-shape	-4.54	-4.65 -5.17	-6.57	-5.66	-5.57 -5.44	-5.62
	22 Phenol dimer	-6.48	-6.62 -7.07	-8.16	-7.49	-7.43 -7.35	-7.09
	ME	0.56	0.47	0.17	-0.54	-0.11	-0.07 -0.01	0.00
	MAE	0.56	0.47	0.17	0.54	0.11	0.09 0.11	0.00
	MA%E	15.8% 13.5% 5.0% 13.6%	2.7%	2.2% 2.6%	0.0%
	total ME	0.53	0.40 -0.04	-1.15	-0.44	-0.39 -0.30	0.00
	total MAE	0.75	0.71	0.47	1.15	0.44	0.41 0.40	0.00
	total MA%E	18.7% 17.0% 8.4% 18.6%	4.1%	3.7% 4.3%	0.0%
	III. COMPUTATIONAL DETAILS
	All calculations have been done with a development version
	of MOLPRO 2008

Acknowledgments

This work was supported by ANR (French national research agency) via contract number ANR-07-BLAN-0272 (Wademecom). We thank D. Mukherjee (Kolkata, India) for discussions.

Appendix A: CCD correlation energy

In this appendix, we review several equivalent CCD correlation energy expressions, in view of justifying the different ring CCD variants.