Asymptotic expansion of stationary distribution for reflected Brownian motion in the quarter plane via analytic approach

Brownian motion in R 2 + with covariance matrix $\Sigma$ and drift $\mu$ in the interior and reflection matrix R from the axes is considered. The asymptotic expansion of the stationary distribution density along all paths in R 2 + is found and its main term is identified depending on parameters ($\Sigma$, $\mu$, R). For this purpose the analytic approach of Fayolle, Iasnogorodski and Malyshev in [12] and [36], restricted essentially up to now to discrete random walks in Z 2 + with jumps to the nearest-neighbors in the interior is developed in this article for diffusion processes on R 2 + with reflections on the axes.


Introduction and main results
1.1.Context.Two-dimensional semimartingale reflecting Brownian motion (SRBM) in the quarter plane received a lot of attention from the mathematical community.Problems such as SRBM existence [39,40], stationary distribution conditions [19,22], explicit forms of stationary distribution in special cases [7,8,19,23,30], large deviations [1,7,33,34] construction of Lyapunov functions [10], and queueing networks approximations [19,21,31,32,43] have been intensively studied in the literature.References cited above are non-exhaustive, see also [42] for a survey of some of these topics.Many results on two-dimensional SRBM have been fully or partially generalized to higher dimensions.
In this article we consider stationary SRBMs in the quarter plane and focus on the asymptotics of their stationary distribution along any path in R 2 + .Let Z(∞) = (Z 1 (∞), Z 2 (∞)) be a random vector that has the stationary distribution of the SRBM.In [6], Dai et Myazawa obtain the following asymptotic result: for a given directional vector c ∈ R 2 + they find the function f c (x) such that where • | • is the inner product.In [7] they compute the exact asymptotics of two boundary stationary measures on the axes associated with Z(∞).In this article we solve a harder problem arisen in [6, §8 p.196], the one to compute the asymptotics of where c ∈ R 2 + is any directional vector and B ⊂ R 2 + is a compact subset.Furthermore, our objective is to find the full asymptotic expansion of the density π(x 1 , x 2 ) of Z(∞) as x 1 , x 2 → ∞ and x 2 /x 1 → tan(α) for any given angle α ∈]0, π/2[.Our main tool is the analytic method developed by V. Malyshev in [36] to compute the asymptotics of stationary probabilities for discrete random walks in Z 2  + with jumps to the nearestneighbors in the interior and reflections on the axes.This method proved to be fruitful for the analysis of Green functions and Martin boundary [26,28], and also useful for studying some joining the shortest queue models [29].The article [36] has been a part of Malyshev's global analytic approach to study discrete-time random walks in Z 2 + with four domains of spatial homogeneity (the interior of Z 2 + , the axes and the origin).Namely, in the book [35] he made explicit their stationary probability generating functions as solutions of boundary problems on the universal covering of the associated Riemann surface and studied the nature of these functions depending on parameters.G. Fayolle and R. Iasnogorodski [11] determined these generating functions as solutions of boundary problems of Riemann-Hilbert-Carleman type on the complex plane.Fayolle, Iansogorodski and Malyshev merged together and deepened their methods in the book [12].The latter is entirely devoted to the explicit form of stationary probabilities generating functions for discrete random walks in Z 2  + with nearest-neighbor jumps in the interior.The analytic approach of this book has been further applied to the analysis of random walks absorbed on the axes in [26].It has been also especially efficient in combinatorics, where it allowed to study all models of walks in Z 2  + with small steps by making explicit the generating functions of the numbers of paths and clarifying their nature, see [38] and [27].
However, the methods of [12] and [36] seem to be essentially restricted to discrete-time models of walks in the quarter plane with jumps in the interior only to the nearest-neighbors.They can hardly be extended to discrete models with bigger jumps, even at distance 2, nevertheless some attempts in this direction have been done in [13].In fact, while for jumps at distance 1 the Riemann surface associated with the random walk is the torus, bigger jumps lead to Riemann surfaces of higher genus, where the analytic procedures of [12] seem much more difficult to carry out.Up to now, as far as we know, neither the analytic approach of [12], nor the asymptotic results [36] have been translated to the continuous analogs of random walks in Z 2 + , such as SRBMs in R 2 + , except for some special cases in [2] and in [16].This article is the first one in this direction.Namely, the asymptotic expansion of the stationary distribution density for SRBMs is obtained by methods strongly inspired by [36].The aim of this work goes beyond the solution of this particular problem.It provides the basis for the development of the analytic approach of [12] for diffusion processes in cones of R 2 + which is continued in the next articles [17] and [18].In [18] the first author and K. Raschel make explicit Laplace transform of the invariant measure for SRBMs in the quarter plane with general parameters of the drift, covariance and reflection matrices.Following [12], they express it in an integral form as a solution of a boundary value problem and then discuss possible simplifications of this integral formula for some particular sets of parameters.The special case of orthogonal reflections from the axes is the subject of [17].Let us note that the analytic approach for SRBMs in R 2 + which is developed in the present paper and continued by the next ones [17] and [18], looks more transparent than the one for discrete models and deprived of many second order details.Last but not the least, contrary to random walks in Z 2  + with jumps at distance 1, it can be easily extended to diffusions in any cones of R 2 via linear transformations, as we observe in the concluding remarks, see Section 5.3.1.2.Reflected Brownian motion in the quarter plane.We now define properly the twodimensional SRBM and present our results.Let r 11 r 12 r 21 r 22 ∈ R 2×2 be a reflection matrix.
Columns R 1 and R 2 represent the directions where the Brownian motion is pushed when it reaches the axes, see Figure 1.Proposition 2. The reflected Brownian motion Z(t) associated with (Σ, µ, R) is well defined, and its stationary distribution Π exists and is unique if and only if the data satisfy the following conditions: r 11 > 0, r 22 > 0, r 11 r 22 − r 12 r 21 > 0, The proof and some more detailed statements can be found in [24,41,20].From now on we assume that conditions (1) and (2) are satisfied.The stationary distribution Π is absolutely continuous with respect to Lebesgue measure as it is shown in [22] and [4].We denote its density by π(x 1 , x 2 ).
1.3.Functional equation for the stationary distribution.Let A be the generator of (W t + µt) t 0 .For each f ∈ C 2 b (R 2 + ) (the set of twice continuously differentiable functions f on R 2 + such that f and its first and second order derivatives are bounded) one has Let us define for i = 1, 2, D i f (x) = R i |∇f that may be interpreted as generators on the axes.We define now ν 1 and ν 2 two finite boundary measures with their support on the axes: for any Borel set B ⊂ R 2 + , By definition of stationary distribution, for all t 0, E Π [f (Z(t))] = R 2 + f (z)Π(dz).A similar formula holds true for ν i : Therefore ν 1 and ν 2 may be viewed as a kind of boundary invariant measures.The basic adjoint relationship takes the following form: The proof can be found in [22] in some particular cases and then has been extended to a general case, for example in [5].We now define ϕ(θ) the two-dimensional Laplace transform of Π also called its moment generating function.Let for all θ = (θ 1 , θ 2 ) ∈ C 2 such that the integral converges.It does of course for any θ with θ 1 0, θ 2 0. We have set θ|Z = θ 1 Z 1 + θ 2 Z 2 .Likewise we define the moment generating functions for ν 1 (θ 2 ) and ν 2 (θ 1 ) on C: Function ϕ 2 (θ 1 ) exists a priori for any θ 1 with θ 1 0. It is proved in [6] that it also does for θ 1 with θ 1 ∈ [0, 1 ], up to its first singularity 1 > 0, the same is true for ϕ 1 (θ 2 ).The following key functional equation (proven in [6]) results from the basic adjoint relationship (3).

Notation. We write the asymptotic expansion
It will be more convenient to expand π(r cos α, r sin α) as r → ∞ and α → α 0 .We give our final results in Section 5, Theorems 22-25: we find the expansion of π(r cos α, r sin α) as r → ∞ and prove it uniform for α fixed in a small neighborhood O(α In this section, Theorem 4 below announces the main term of the expansion depending on parameters (µ, Σ, R) and a given direction α 0 .Next, in Section 1.5 we sketch our analytic approach following the main lines of this paper in order to get the full asymptotic expansion of π.We present at the same time the organization of the article.Now we need to introduce some notations.The quadratic form γ(θ) is defined in (4) via the covariance matrix Σ and the drift µ of the process in the interior of R 2 + .Let us restrict ourselves on θ ∈ R 2 .The equation γ(θ) = 0 determines an ellipse E on R 2 passing through the origin, its tangent in it is orthogonal to vector µ, see Figure 2. Stability conditions (1) and (2) imply the negativity of at least one of coordinates of µ, see [6,Lemma 2.1].In this article, in order to shorten the number of pictures and cases of parameters to consider, we restrict ourselves to the case µ 1 < 0 and µ 2 < 0, (6) although our methods can be applied without any difficulty to other cases, we briefly sketch some different details at the end of Section 2.4.
Let us note that the exponents in Theorem 4 are the same as in the large deviation rate function found in [7,Thm 3.2].The same phenomenon is observed for discrete random walks, cf.[36] and [25].
1.5.Sketch of the analytic approach.Organization of the paper.The starting point of our approach is the main functional equation (4) valid for any θ = (θ 1 , θ 2 ) ∈ C 2 with θ 1 0, θ 2 0. The function γ(θ 1 , θ 2 ) in the left-hand side is a polynomial of the second order of θ 1 and θ 2 .The algebraic function Θ 1 (θ 2 ) defined by γ(Θ 1 (θ 2 ), θ 2 ) ≡ 0 is 2-valued and its Riemann surface S θ 2 is of genus 0. The same is true about the 2-valued algebraic function Θ 2 (θ 1 ) defined by γ(θ 1 , Θ 2 (θ 1 )) = 0 and its Riemann surface S θ 1 .The surfaces S θ 1 and S θ 2 being equivalent, we will consider just one surface S defined by the equation γ(θ 1 , θ 2 ) = 0 with two different  2),(3), (4) coverings.Each point s ∈ S has two "coordinates" (θ 1 (s), θ 2 (s)), both of them are complex or infinite and satisfy γ(θ 1 (s), θ 2 (s)) = 0.For any point s = (θ 1 , θ 2 ) ∈ S, there exits a unique point s = (θ 1 , θ 2 ) ∈ S with the same first coordinate and there exists a unique point s = (θ 1 , θ 2 ) ∈ S with the same second coordinate.We say that s = ζs, i.e. s and s are related by Galois automorphism ζ of S that leaves untouched the first coordinate, and that s = ηs, i.e. s and s are related by Galois automorphism η of S that leaves untouched the second coordinate.Clearly ζ 2 = Id, η 2 = Id and the branch points of Θ 1 (θ 2 ) and of Θ 2 (θ 1 ) are fixed points of ζ and η respectively.The ellipse E is the set of points of S where both "coordinates" are real.The construction of S and definition of Galois automorphisms are carried out in Section 2.
Next, unknown functions ϕ 1 (θ 2 ) and ϕ 2 (θ 1 ) are lifted in the domains of S where {s ∈ S : θ 2 (s) 0} and {s ∈ S : θ 1 (s) 0} respectively.The intersection of these domains on S is non-empty, both ϕ 2 and ϕ 1 are well defined in it.Since for any s = (θ 1 (s), θ 2 (s)) ∈ S we have γ(θ 1 (s), θ 2 (s)) = 0, the main functional equation (4) implies: Using this relation, Galois automorphisms and the facts that ϕ 1 and ϕ 2 depend just on one "coordinate" (ϕ 1 depends on θ 2 and ϕ 2 on θ 1 only), we continue ϕ 1 and ϕ 2 explicitly as meromorphic on the whole of S. This meromorphic continuation procedure is the crucial step of our approach, it is the subject of Section 3.1.It allows to recover ϕ 1 and ϕ 2 on the complex plane as multivalued functions and determines all poles of all its branches.Namely, it shows that poles of ϕ 1 and ϕ 2 may be only at images of zeros of γ 1 and γ 2 by automorphisms η and ζ applied several times.We are in particular interested in the poles of their first (main) branch, and more precisely in the most "important" pole (from the asymptotic point of view, to be explained below), that turns out to be at one of points ζθ * * or ηθ * defined above.The detailed analysis of the "main" poles of ϕ 1 and ϕ 2 is furnished in Section 3.2.
Let us now turn to the asymptotic expansion of the density π(x 1 , x 2 ).Its Laplace transform comes from the right-hand side of the main equation (4) divided by the kernel γ(θ 1 , θ 2 ).By inversion formula the density π(x 1 , x 2 ) is then represented as a double integral on {θ : θ 1 = θ 2 = 0}.In Section 4.1, using the residues of the function 1 γ(θ 1 ,•) or 1 γ(•,θ 2 ) we transform this double integral into a sum of two single integrals along two cycles on S, those where θ 1 (s) = 0 or θ 2 (s) = 0. Putting (x 1 , x 2 ) = re α we get the representation of the density as a sum of two single integrals along some contours on S: We would like to compute their asymptotic expansion as r → ∞ and prove it to be uniform for α fixed in a small neighborhood O(α 0 ), α 0 ∈]0, π/2[.These two integrals are typical to apply the saddle-point method, see [15,37].The point θ(α) ∈ E defined above is the saddle-point for both of them, this is the subject of Section 4.2.The integration contours on S are then shifted to new ones Γ θ 1 ,α and Γ θ 2 ,α which are constructed in such a way that they pass through the saddle-point θ(α), follow the steepest-descent curve in its neighborhood O(θ(α)) and are "higher" than the saddle-point w.r.t. the level curves of the function θ(s) | e α outside O(θ(α)), see Section 4.3.Applying Cauchy Theorem, the density is finally represented as a sum of integrals along these new contours and the sum of residues at poles of the integrands we encounter deforming the initial ones: Here P α (resp.P α ) is the set of poles of the first order of ϕ 1 (resp.ϕ 2 ) that are found when shifting the initial contour I + θ 1 to the new one Γ θ 1 ,α (resp.I + θ 2 to Γ θ 2 ,α ), all of them are on the arc {s 0 , θ(α)} (resp.{θ(α), s 0 }) of ellipse E.
The asymptotic expansion of integrals along Γ θ 1 ,α and Γ θ 2 ,α is made explicit by the standard saddle-point method in Section 4.4.The set of poles P α ∪P α is analyzed in Section 4.5 .In Case (1) of Theorem 4 this set is empty, thus the asymptotic expansion of the density is determined by the saddle-point, its first term is given in Theorem 4. In Cases (2), (3) and (4) this set is not empty.The residues at poles over P α ∪ P α in (13) bring all more important contribution to the asymptotic expansion of π(re α ) than the integrals along Γ θ 1 ,α and Γ θ 2 ,α .Taking into account the monotonicity of function θ | e α on the arcs {s 0 , θ(α)} and on {θ(α), s 0 }, they can be ranked in order of their importance: clearly, the term associated with a pole p is more important than the one with p if p | e α < p | e α .In Case (2) (resp.(3)) the most important pole is ζθ * * (resp.ηθ * ), as announced in Theorem 4. In Case (4) the most important of them is among ζθ * * and ηθ * , as stated in Theorem 4 as well.The expansion of integrals in (13) along Γ θ 1 ,α and Γ θ 2 ,α via the saddle-point method provides all smaller asymptotic terms than those coming from the poles.Section 5 is devoted to the results: they are stated from two points of view in Sections 5.1 and 5.2 respectively.First, given an angle α 0 , we find the uniform asymptotic expansion of the density π(r cos(α), r sin(α)) as r → ∞ and α ∈ O(α 0 ) depending on parameters (Σ, µ, R): Theorems 22 -25 of Section 5.1 state it in all cases of parameters (1)-(4).Second, in Section 5.2, given a set parameters (Σ, µ, R), we compute the asymptotics of the density for all angles α 0 ∈]0, π/2[, see Theorems 26-28.Remark.The constants mentioned in Theorem 4 and all those in asymptotic expansions of Theorems 22-28 are specified in terms of functions ϕ 1 and ϕ 2 .In the present paper we leave unknown these functions in their initial domains of definition although we carry out explicitly their meromorphic continuation procedure and find all their poles.In [18] the first author and K. Raschel make explicit these functions solving some boundary value problems.This determines the constants in asymptotic expansions in Theorems 4, 22-28 .
Future works.The case of parameters such that ζθ * * = θ(α) and ηθ * ∈ {s 0 , θ(α){ or the case such that ηθ * = θ(α) and ζθ * * ∈ {s 0 , θ(α){ are not treated in Theorem 4. Theorem 25 gives a partial result but not at all as satisfactory as in all other cases.In fact, in these cases the saddle-point θ(α) coincides with the "main" pole of ϕ 1 or ϕ 2 .The analysis is then reduced to a technical problem of computing the asymptotics of an integral whenever the saddle-point coincides with a pole of the integrand or approaches to it.We leave it for the future work.
In the cases α = 0 and α = π/2, the tail asymptotics of the boundary measures ν 1 and ν 2 has been found in [7] and the constants have been specified in [18].It would be also possible to find the asymptotics of π(r cos α, r sin α) where r → ∞ and α → 0 or α → π/2.This problem is reduced to obtaining the asymptotics of an integral when the saddle-point θ(0) or θ(π/2) coincides with a branch point of the integrand ϕ 1 or ϕ 2 .It can be solved by the same methods as in [26] for discrete random walks.

Riemann surface S
2.1.Kernel γ(θ 1 , θ 2 ).The kernel of the main functional equation can be written as The equation γ(θ 1 , θ 2 ) ≡ 0 defines a two-valued algebraic function Θ 1 (θ 2 ) such that γ(Θ 1 (θ 2 ), θ 2 ) ≡ 0 and a two-valued algebraic function Θ 2 (θ 1 ) such that γ(θ 1 , Θ 2 (θ 1 ) ≡ 0. These functions have two branches: , and that are both real and of opposite signs: . We can compute: 2.2.Construction of the Riemann surface S. We now construct the Riemann surface S of the algebraic function Θ 2 (θ 1 ).For this purpose we take two Riemann spheres ), and we glue them together along the borders of these cuts, joining the lower border of the cut on S 1 θ 1 to the upper border of the same cut on S 2 θ 1 and vice versa.This procedure can be viewed as gluing together two half-spheres, see Figure 6.The resulting surface S is homeomorphic to a sphere (i.e., a compact Riemann surface of genus 0) and is projected on the Riemann sphere C ∪ {∞} by a canonical covering map h θ 1 : S → C ∪ {∞}.In a standard way, we can lift the function Θ 2 (θ 1 ) to S, by setting In a similar way one constructs the Riemann surface of the function Θ 1 (θ 2 ), by gluing together two copies S 1 θ 2 and S 2 θ 2 of the Riemann sphere S cut along ).We obtain again a surface homeomorphic to a sphere where we lift function Θ 1 (θ 2 ).
Since the Riemann surfaces of Θ 1 (θ 2 ) and Θ 2 (θ 1 ) are equivalent, we can and will work on a single Riemann surface S, with two different covering maps h θ 1 , h θ 2 : S → C ∪ {∞}.Then, for s ∈ S, we set θ 1 (s) = h θ 1 (s) and θ 2 (s) = h θ 2 (s).We will often represent a point s ∈ S by the pair of its coordinates (θ 1 (s), θ 2 (s)).These coordinates are of course not independent, because the equation γ(θ 1 (s), θ 2 (s)) = 0 is valid for any s ∈ S. One can see S with points This contour goes from s ∞ to s ∞ and back as well, but via s − 2 and s + 2 .Let E be the set of points of S where both coordinates θ 1 (s) and θ 2 (s) are real.Then 5 and 7, it contains all branch points s ± 1 and s ± 2 .2.3.Galois automorphisms η and ζ.Now we need to introduce Galois automorphisms on S. For any s ∈ S \ s ± 1 there is a unique θ 1 and vice versa.On the other hand, whenever . With the previous notations we now define the mappings ζ : S → S and η : S → S by Following [35] we call them Galois automorphisms of S. Then ζ 2 = η 2 = Id, and 1 and s + 1 (resp.s − 2 and s + 2 ) are fixed points for ζ (resp.η).It is known that conformal automorphisms of a sphere (that can be identified to C ∪ ∞) are transformations of type z → az+b cz+d where a, b, c, d are any complex numbers satisfying ad − bc = 0.The automorphisms ζ and η, which are conformal automorphisms of S, have each two fixed points and are involutions (because ζ 2 = η 2 = Id).We can deduce from it that ζ (resp.η) is a symmetry w.r.t. the axis A 1 (resp.A 2 ) that passes through fixed points s − 1 and s + 1 (resp.s − 2 and s + 2 ).In other words ζ (resp.η) is a rotation of angle π, around D 1 (resp.A 2 ), see Figure 8.Let us draw the axis A orthogonal to the plane generated by the axes A 1 and A 2 and passing through the intersection point of A 1 and A 2 .We denote by β the angle between the axes A 1 and A 2 .Automorphisms ηζ and ζη are then rotations of angle 2β and −2β around the axis A. This axis goes through points s ∞ and s ∞ which are fixed points for ηζ and ζη, see Figure 8.
In the particular case Σ = Id, we have ηζ = ζη, the axes A 1 and A 2 are orthogonal.We deduce that β = π 2 and that ηζ and ζη are symmetries w.r.t. the axis A.
Assume now that the interior drift has both coordinates negative, i.e. (6).From what said above, The intersection ∆ 1 ∩ ∆ 2 consists of two connected components, both bounded by I − θ 1 and I − θ 2 .The union ∆ 1 ∪ ∆ 2 is a connected domain, but not simply connected because of the point s 0 .The domain ∆ 1 ∪ ∆ 2 ∪ s 0 is open, simply connected and bounded by I + θ 1 and I + θ 2 , see Figure 9.We set ∆ = ∆ 1 ∪ ∆ 2 .Note that in the cases of stationary SRBM with drift µ having one of coordinates non-negative, the location of contours For example, assume that µ 2 > 0. Then RΘ − 2 (θ 1 ) < 0 and RΘ + 2 (θ 1 ) 0, the contour where the second coordinate is negative, while I + θ 1 goes from s ∞ to s ∞ crossing E at s 0 = (0, 0).In order to shorten the number of cases and pictures, we restrict ourselves in this paper to the case (6) of both coordinates of µ negative, although all our methods work in these other cases as well.

Parametrization of S.
It is difficult to visualize on three-dimensional sphere different points, contours, automorphisms and domains introduced above that will be used in future steps.For this reason we propose here an explicit and practical parametrisation of S. Namely we identify S to C ∪ {∞} and in the next proposition we explicitly define h θ 1 and h θ 2 two recoveries introduced in Section 2.2.Such a parametrisation allows to visualize better in two dimensions the sphere S ≡ C ∪ {∞} and all sets we are interested in, as we can see in Figure 10.Proposition 5. We set the following covering maps The equation γ(θ 1 (s), θ 2 (s)) = 0 is valid for any s ∈ S. Galois automorphisms can be written s , and ηζ (resp.ζη) is a rotation around s ∞ ≡ 0 of angle 2β (resp.−2β) according to counterclockwise direction.
Proof.We set This parametrization is practical because it leads to a similar rational recovery h θ 2 .In order to make the equation γ(θ 1 (s), θ 2 (s)) = 0 valid for any s ∈ S we naturally set and we are going to show that θ 2 (s) = ( s e iβ + e iβ s ) where β = arccos − σ 12 √ σ 11 σ 22 .We note that d(θ 1 (s)) is the opposite of the square of a rational fraction Then we have Furthermore this parametrization leads to simple expressions for Galois automorphisms η and ζ.We derive immediately that θ Next we search η as an automorphism of the form ηs = K s .Since θ 2 (s) is of the form θ 2 (s) = us + v s + w with constants u, v, w defined by ( 14), then After setting .
Then we obtain (the last equality follows from elementary geometric properties of an ellipse).It implies concluding the proof.
Figure 10 shows different sets we are interested in according to the parametrization we have just introduced.We have We can determine the equation of the analytic curves of pure imaginary points of θ i .We have cos(a) cosh(b).It follows that We can easily notice that 3. Meromorphic continuation of ϕ 1 and ϕ 2 on S 3.1.Lifting of ϕ 1 and ϕ 2 on S and their meromomorphic continuation.
Lemma 6. Functions ϕ 1 and ϕ 2 (defined on ∆2 and ∆1 respectively) can be meromorphically continued on ∆ ∪ {s 0 } by setting Proof.The open set ∆ 1 ∩ ∆ 2 is non-empty and bounded by the curve 16) is valid for s ∈ ∆ 1 ∩∆ 2 .It allows us to continue functions ϕ 1 and ϕ 2 as meromorphic on ∆ as stated in this lemma.The functional equation ( 16) is then valid on the whole of ∆, as well as the invariance formulas (18).
Functions ϕ 1 and ϕ 2 can be then of course continued to ∆. Moreover we have the following lemma.
We may now in the same way, using formulas ( 21) and ( 22), continue function ϕ 1 (s) (resp.ϕ 2 (s)) as meromorphic on (ηζ) 2 ∆ , (ηζ) 3 ∆ (resp. Function ϕ 2 (s) can be continued as meromorphic on ζη∆ , (ζη) 2 ∆ , • • • (ζη) n ∆ by the formulas : Proof.We proceed by induction on k = 1, 2, . . .n.For k = 1, this is the subject of the previous lemma.For any k = 2, . . ., n, assume the formula (29) for any s ∈ (ηζ Proceeding as in Lemma 10 by rotations, we will continue ϕ 1 soon on the first half of S, that is S 1 θ 2 , then the whole of S and go further, turning around S for the second time, for the third, etc up to infinity.In fact, each time we complete this procedure on one of two halves of S, we recover a new branch of the function ϕ 1 as function of θ 2 ∈ C. So, going back to the complex plane, we continue this function as multivalued and determine all its branches.The same is true for ϕ 2 if we proceed by rotations in the opposite direction.This procedure could be presented better on the universal covering of S, but for the purpose of the present paper it is enough to complete it only on one-half of S, that is to study just the first (main) branch of ϕ 1 and ϕ 2 .We summarize this result in the following theorem.We recall that S = S 1 θ 1 ∪ S 2 θ 1 and we denote by S 1 θ 1 the half that contains s 0 (and not s 0 , as ζs 0 = s 0 ).In the same way S = S 1 θ 2 ∪ S 2 θ 2 and we denote by S 1 θ 2 the half that contains s 0 (and not s 0 , as ηs 0 = s 0 ), see Figure 10.
Proof.It is a direct corollary of Lemma 7 and Lemma 10.

3.2.
Poles of functions ϕ 1 and ϕ 2 on S. It follows from meromorphic continuation procedure that all poles of ϕ 1 (s) and ϕ 2 (s) on S are located on the ellipse E, they are images of zeros of γ 1 and γ 2 by automorphisms η and ζ applied several times.Then all poles of all branches of ϕ 2 (s) ). Notations of arcs on E. Let us remind that we denote by {s 1 , s 2 } an arc of the ellipse E with ends at s 1 and s 2 not passing through the origin, see Theorem 4. From now on, we will denote in square brackets ]s 1 , s 2 [ or [s 1 , s 2 ] an arc of E going in the anticlockwise direction from s 1 to s 2 .
In order to compute the asymptotic expansion of stationary distribution density, we are interested in poles of ϕ 1 on the arc ]s 0 , s + 2 [ and in those of ϕ 2 on the arc ]s + 1 , s 0 [.To determine the main asymptotic term, we are particularly interested in the pole of ϕ 1 (θ 2 (s)) on ]s 0 , s + 2 [ closest to s 0 and in the one of ϕ 2 (θ 1 (s)) on ]s + 1 , s 0 [ closest to s 0 .We identify them in this section.
We remind that θ * is a zero of γ 1 (s) on E different from s 0 and that θ * * is a zero of γ 2 (s) on E different from s 0 .Their coordinates are Lemma 12. ( Let us check that the numerator in ( 35) is non zero, this will prove the statement (1) the lemma.
The reasoning for θ * is the same.
If parameters (Σ, µ, R) are such that θ 2 (s + 1 ) > 0 and ηζθ p ∈]ηs + 1 , s 0 [, then either ηζθ p ∈ [s 0 , s 0 ] or ηζθ p ∈]s 0 , θ p [.In the first case we have (39) as previously from where γ 2 (ζθ p ) = 0 and the pole ζθ p is of the first order.Let us turn to the second case ηζθ p ∈]θ p , s 0 [ for which we will use the formula (37).The pole θ p being the closest to s 0 , then ηζθ p can not be a pole of ϕ 2 on ]θ p , s 0 [.It can neither be a pole of ϕ 2 on [s 0 , s 0 ], since this function is initially well defined on this segment.Hence in formula (37) ϕ 2 (ηζθ p 1 ) = ∞ for ηζθ p ∈]θ p , s 0 [.It follows from (37) that either γ 2 (ζθ p ) = 0 or γ 1 (ηζθ p ) = 0 and if these two equalities do not hold simultaneously, then pole θ p must be of the first order.
The proof in the case (ii) is symmetric.
Figure 13 gives two illustrations of Lemmas 12 and 13.On the left figure the parameters are such that θ 1 (s + 2 ) > 0 and θ 2 (s + 1 ) > 0. Let us look at zeros θ * of γ 1 and θ * * of γ 2 different from s 0 .We see θ * ∈]s + 2 , s 0 [, then ηθ * is the first candidate for the closest pole of ϕ 1 to s 0 on ]s 0 , s + 2 [.We also see θ We will also need the following two lemmas.

Contribution of the saddle-point and of the poles to the asymptotic expansion
4.1.Stationary distribution density as a sum of integrals on S. By the functional equation (4) and the inversion formula of Laplace transform (we refer to [9] and [3]), the density π(x 1 , x 2 ) can be represented as a double integral We now reduce it to a sum of single integrals.
Lemma 16.For any where and Proof.By inversion formula ( 40) Now it suffices to show the following formulas Let us prove (43).For any )) has two zeros Θ + 2 (θ 1 ) and Θ − 2 (θ 1 ).Their real parts are of opposite signs: (Θ − 2 (θ 1 )) < 0 and (Θ + 2 (θ 1 )) > 0. Thus for any fixed of the argument θ 2 has two poles on the complex plane C θ 2 , one at Θ − 2 (θ 1 ) with negative real part and another one at the purely imaginary segment [−iR, iR] and the half of the circle with radius R and center 0 on C θ 2 , see Figure 14.For R large enough Θ + 2 (θ 1 ) is inside the contour.The integral of γ 2 (θ)e −x 2 θ 2 γ(θ) over this contour taken in the counter-clockwise direction equals the residue at the unique pole of the integrand: for all large enough R 0. (45) Let us take the limit of this integral as R → ∞: The last term equals We note that sup [. Then by dominated convergence theorem the limit (47) equals 0 as R → ∞.Hence, due to (45) and ( 46) that proves (43) for any θ 1 ∈ iR \ {0}.The proof of (44) is analogous.Note also that the integral is absolutely convergent.In fact sup Then the integral is absolutely convergent.This concludes the proof of formula (41).The proof of ( 42) is completely analogous.
Remark.These integrals are equal to those on the Riemann surface S along properly oriented contours I + θ 1 and I + θ 2 respectively.Thanks to the parametrization of Section 2.5 we have Then we can write for x = (x 1 , x 2 ) ∈ R 2 + the density π(x 1 , x 2 ) as a sum of two integrals on S : where α ∈]0, π/2[, Our aim now is to find the asymptotic expansion of π(r cos(α), r sin(α)), that is the one of the sum (49) as r → ∞ and to prove that for any α 0 ∈]0, π/2[ this asymptotic expansion is uniform in a small neighborhood O(α 0 ) ∈]0, π/2[.These integrals are typical to apply the saddle-point method, see [15] or [37].Let us study the function θ(s) | e α on S and its critical points.

4.3.
Shifting the integration contours.Our aim now is to shift the integration contours I + θ 1 and I + θ 2 in (49) up to new contours Γ θ 1 ,α and Γ θ 2 ,α respectively which coincide with γ α in a neighborhood of θ(α) on S and are "higher" than θ(α) in the sense of level curves of the function When shifting the contours we should of course take into account the poles of the integrands and the residues at them.
Let us construct Γ θ 1 ,α and Γ θ 2 ,α .We set where V (α) > 0 will be defined later.Then the end points of Γ 1,+ θ 1 ,α are z α,+ and Z α,+ where θ 1 ,α coincides with I + θ 1 from Z 0 α,+ up to infinity : We define in the same way One can visualize this contour on Figure 16 : in the left picture it is drawn on parametrized S, in the right picture it is projected on the complex plane C θ 1 .The contour Γ θ 2 ,α is constructed analogously with respect to Let us recall that poles of ϕ 1 (s) and ϕ 2 (s) on S may occur only at E. Let us also recall the convention that an arc }a, b{ on E is the one with ends a and b which does not include s 0 = (0, 0).Notation of the sets of poles P α and P α .Let P α be the set of poles of the first order of the function ϕ 2 (θ 1 (s)) on the arc }θ(α), s 0 {.Let P α be the set of poles of the first order of the function ϕ 1 (θ 2 (s)) on the arc }θ(α), s 0 {.
Then the following lemma holds true.
Proof.It follows from the assumption of the lemma that θ(α) is not a pole of ϕ 1 (θ 2 (s)) neither of ϕ 2 (θ 1 (s)) for any α in a small enough neighborhood O(α 0 ).Then we use the representation of the density (49) and apply Cauchy theorem shifting the contours to Γ θ 1 ,α and Γ θ 2 ,α .
In order to find the asymptotic expansion of the density π(r cos(α), r sin(α)), we have to evaluate now the contribution of the residues at poles in (51) and the one of integrals along shifted contours Γ θ 1 ,α and Γ θ 2 ,α .This is the subject of the next two sections.4.4.Asymptotics of integrals along shifted contours Γ θ 1 ,α and Γ θ 2 ,α .To finish the construction of Γ θ 1 ,α and Γ θ 2 ,α , it remains to specify V (α).For that purpose we consider closer the function f α (s) = θ(s) | e α = θ 1 (s) cos α + θ 2 (s) sin α.Let us define the projection of this function on C θ 1 : (ii) There exist constants d 1 0, d 2 > 0 and V > 0 such that: Proof.We compute : with the discriminant d(u where ω − (u + iv) et ω + (u + iv) are defined as ω − (u + iv) = arg(θ + 1 − u − iv) and ω + (u + iv) = arg(u + iv − θ − 1 ), see Figure 17.We have Thus Both statements of the lemma follow directly from this representation.We may now choose V (α) and such that in accordance with notations of Lemma 19.This concludes the construction of Γ θ 1 ,α and Γ θ 2 ,α .The asymptotic expansion of integrals along these contours is given in the following lemma.The main contribution comes from the integrals along γ α , while all other parts of integrals are proved to be exponentially negligible by construction.
In the first theorem parameters (Σ, µ, R) are such that the asymptotic expansion is determined by the saddle-point.
In Theorem 23 θ(α 0 ) is assumed not to be a pole of ϕ 1 and neither of ϕ 2 , that is why Lemma 18 applies.Nevertheless, it may happen (for a very few angles and under some sets of parameters) that θ(α 0 ) is a pole of one of these functions.In this case the following theorem holds true.
Then for any δ > 0 there exists a small enough neighborhood O(α 0 ) such that π(r cos(α), r sin(α)) ∼ Furthermore, the main term in this expansion is the same as in Theorem 23, cases (i), (ii) and (iii).

Concluding remarks.
Let us remark that the approach of this article applies to the SRBM in any cone of R 2 .Thanks to a linear transformation T ∈ R 2×2 , it is easy to transform Z(t), a reflected Brownian motion of parameters (Σ, µ, R) in a cone into T Z(t) a reflected Brownian motion of parameters (T ΣT t , T µ, T R) in the quarter plane.For example if the initial cone is the set {(x, y)|x 0 and y ax} for some a > 0, we may just take T = 1 − 1 a 0 1 .The process T Z(t) lives in a quarter plane.Then the approach of this article applies and its results can be converted to the initial cone by the inverse linear transformation.The analytic approach for discrete random walks is essentially restricted to those with jumps to the nearest neighbors in the interior of the quarter plane.Since a linear transformation can not generally keep the length of jumps, this procedure does not work in the discrete case.That is why the analytic approach in R 2 has a more general scope of applications.
The first term in (81) is just the Laplace transform of the function h 1 (θ 1 ) = g 1 (θ 1 ) 1 , its asymptotics is determined by the smallest real singularity of h 1 (θ 1 ), see e.g.[9].This may be either the branch point θ + 1 of ϕ 2 (θ 1 ), or the smallest pole of h 1 (θ 1 ) on ]0, θ + 1 [ whenever it exists, the natural candidates are ζθ * * , ζηθ * due to Lemmas 12-15 or a point θ c = (θ c 1 , θ c 2 ) such that θ c 2 = Θ + 2 (θ c 1 ) = θ c 1 c 2 c 1 .To determine the asymptotics of the second integral in (81), we shift the integration contour to the new one passing through the saddle-point Θ 1 (θ + 2 ) and take into account the poles of the integrand we encounter, the most important of these poles are those listed above.The asymptotics of two terms in (82) is determined in the same way.Combining all these results together we derive the main asymptotic term depending on the parameters that can be either e This analysis leads to the results of [6].

Figure 1 .
Figure 1.Drift µ and reflection vectors R 1 and R 2

Figure 6 .Figure 7 .
Figure 6.Construction of the Riemann surface S

Figure 9 .
Figure 9. Pure imaginary points of S s and then ηζ(s) = e 2iβ s, ζη(s) = e −2iβ s.It follows that ηζ and ζη are just rotations for angles 2β et −2β respectively.By symmetry considerations we can now rewrite
R −1/2 or R −3/2 with some constant, or e * ) i , i = 1, 2 preceding by some constant and the factor R in some critical cases.