Construction of a relativistic Ornstein-Uhlenbeck process

TL;DR

This paper constructs a relativistic Ornstein-Uhlenbeck process using a modified Dudley's Wiener process on the mass shell, and computes its invariant distribution and related processes.

Contribution

It introduces a novel relativistic Ornstein-Uhlenbeck process model based on Dudley's Wiener process and derives its invariant distribution.

Findings

Invariant distribution of the process computed

Model extends classical Ornstein-Uhlenbeck to relativistic context

Provides analytical tools for relativistic stochastic processes

Abstract

Based on a version of Dudley's Wiener process on the mass shell in the momentum Minkowski space of a massive point particle, a model of a relativistic Ornstein--Uhlenbeck process is constructed by addition of a specific drift term. The invariant distribution of this momentum process as well as other associated processes are computed.

Equations

$$\mathbf{p}=(p_0,p_1,p_2,p_3)\in {\mathbb{R}}^4$$

$$p_0^2-\sum _{i=1}^3{p}_i^2=m^2{c}^2,p_0\ge 0,$$

$$⟨\mathbf{p},\mathbf{q}⟩=(\mathbf{p},g_M\cdot \mathbf{q})=p_0{q}_0-\sum _{i=1}^d{p}_i{q}_i,$$

$$\mathcal{V}_{m}^{d}=\bigl{\{}\mathbf{p}\in\mathbb{R}^{1+d},\,p_{0}>0,\,\langle\mathbf{p},\mathbf{p}\rangle=m^{2}c^{2}\bigr{\}},$$

$$p_0=m\cosh (s),p=m\sinh (s)\omega ,s>0,\omega \in S^{d-1}.$$

$$\iota:(s,\theta)\mapsto\mathbf{p}(s,\theta)=\bigl{(}m\cosh(s),m\sinh(s)\,\omega(\theta)\bigr{)}.$$

$$J=m\begin{pmatrix}\sinh(s)&0\\ \cosh(s)\,\omega(\theta)&\displaystyle\sinh(s)\,\Bigl{(}\frac{\partial\omega(\theta)}{\partial\theta}\Bigr{)}\end{pmatrix},$$

$$g_d=-J^t\cdot g_M\cdot J.$$

$$g_d($$

$$\displaystyle=m^{2}\begin{pmatrix}\cosh(s)^{2}\,\omega(\theta)^{t}\cdot\omega(\theta)-\sinh(s)^{2}&\sinh(s)\cosh(s)\,\omega(\theta)^{t}\cdot\bigl{(}\frac{\partial\omega(\theta)}{\partial\theta}\bigr{)}\\[8.61108pt] \sinh(s)\cosh(s)\,\omega(\theta)\cdot\bigl{(}\frac{\partial\omega(\theta)}{\partial\theta}\bigr{)}^{t}&\sinh(s)^{2}\bigl{(}\frac{\partial\omega(\theta)}{\partial\theta}\bigr{)}^{t}\cdot\bigl{(}\frac{\partial\omega(\theta)}{\partial\theta}\bigr{)}\end{pmatrix}$$

$$=m^2\left(\begin{array}{cc}1& 0\\ 0& \sinh (s{)}^2{g}_{S^{d-1}}(\theta )\end{array}\right),$$

$$\begin{split}d\text{vol}_{d}(s,\theta)&=m\,\sinh(s)^{d-1}\,ds\,d\sigma_{S^{d-1}}(\theta)\\ &=m\,\sinh(s)^{d-1}\bigl{(}\det g_{S^{d-1}}(\theta)\bigr{)}^{1/2}\,ds\,d\theta,\end{split}$$

$$\begin{array}{rl}{\Delta }_d& =\frac{1}{m^2\sinh (s{)}^{d-1}}\frac{\partial }{\partial s}\sinh (s{)}^{d-1}\frac{\partial }{\partial s}+\frac{1}{m^2\sinh (s{)}^2}{\Lambda }_{S^{d-1}}\\ & =\frac{1}{m^2}\frac{{\partial }^2}{\partial s^2}+\frac{d-1}{m^2}\coth (s)\frac{\partial }{\partial s}+\frac{1}{m^2\sinh (s{)}^2}{\Lambda }_{S^{d-1}},\end{array}$$

$${\Delta }_2=\frac{1}{m^2}\frac{{\partial }^2}{\partial s^2}+\frac{1}{m^2}\coth (s)\frac{\partial }{\partial s}+\frac{1}{m^2\sinh (s{)}^2}{\Lambda }_{S^{d-1}}\frac{{\partial }^2}{\partial {\phi }^2},$$

$$\begin{split}\Delta_{3}=\frac{1}{m^{2}}\,\frac{\partial^{2}}{\partial s^{2}}&+\frac{2}{m^{2}}\,\coth(s)\,\frac{\partial}{\partial s}\\ &+\frac{1}{m^{2}\sinh(s)^{2}}\,\Bigl{(}\frac{1}{\sin(\theta)}\,\frac{\partial}{\partial\theta}\,\sin(\theta)\,\frac{\partial}{\partial\theta}+\frac{1}{\sin(\theta)^{2}}\,\frac{\partial^{2}}{\partial\varphi^{2}}\Bigr{)}.\end{split}$$

$$A_{d}f(s)=\frac{1}{2m^2}{f}^{\prime \prime }(s)+\frac{d-1}{2m^2}\coth (s)f^{\prime }(s),s\in (0,+\infty ),$$

$$d{S}_t=\frac{d-1}{2m^2}\coth (S_t)dt+\frac{1}{m}d{W}_t,$$

$$\tau(t)=\int_{0}^{t}\bigl{(}m\sinh(S_{r})\bigr{)}^{-2}\,dr$$

$$B_{t}=\bigl{(}S_{t},\Theta_{\tau(t)}\bigr{)},\qquad t\in\mathbb{R}_{+}$$

$$b_0(s)=\frac{d-1}{2m^2}\coth (s),s>0,$$

$$b_0(s)=\frac{d-1}{2}\coth (s),s>0,$$

$$\rho(s)=\exp\Bigl{(}2\int_{s}^{a}b_{0}(u)\,du\Bigr{)}=\frac{\sinh(a)^{d-1}}{\sinh(s)^{d-1}}.$$

$${\kappa }_a(s)$$

$${\kappa }_{\infty }(s)$$

$${\int }_0^a\rho (s)ds$$

$${\int }_a^{\infty }\rho (s)ds$$

$$I_a$$

$$I_{\infty }$$

$$\int_{0}^{a}\Bigl{(}1+\frac{(d-1)\cosh(s)}{2\sinh(s)}\Bigr{)}\sinh(s)^{d-1}\Bigl{(}\int_{s}^{a}\sinh(u)^{-(d-1)}\,du\Bigr{)}\,ds.$$

$$s\le \sinh (s)\le \sinh (a)s.$$

$$I_{a}\leq\sinh(a)^{d-1}\int_{0}^{a}\Bigl{(}1+\frac{d-1}{2}\,\coth(x)\Bigr{)}\,s^{d-1}\Bigl{(}\int_{s}^{a}u^{-(d-1)}\,du\Bigr{)}\,ds.$$

$$I_{a}\leq\sinh(a)^{d-1}\int_{0}^{a}\Bigl{(}1+\frac{d-1}{2}\,\coth(s)\Bigr{)}\,s\ln\Bigl{(}\frac{a}{s}\Bigr{)}\,ds<+\infty,$$

$$\begin{split}I_{a}\leq\frac{\sinh(a)^{d-1}}{d-2}\int_{0}^{a}\Bigl{(}1+\frac{d-1}{2}\,&\coth(s)\Bigr{)}\,s^{d-1}\\[2.15277pt] &\times\,\bigl{(}s^{-(d-2)}-a^{-(d-2)}\bigr{)}\,ds<+\infty.\end{split}$$

Full text

††copyright: ©2017: J. Potthoff, R. Schrader

1. Introduction

In 1930 Ornstein and Uhlenbeck [32] introduced the stochastic process which afterwards carried their name in order to treat the case where the particle undergoing a motion of Brownian type has a surrounding medium which is a rarefied gas instead of a liquid. They argued that in this case one has to take into account the friction that the particle experiences by hitting the gas molecules, which they called Doppler friction. As a consequence, they proposed an equation of Langevin type for the velocity of the particle instead of for its position. They proved that the velocity of the so defined motion admits a stationary state which is described by a centered normal density. For a discussion of the Ornstein–Uhlenbeck process from the physics point of view, especially in comparison to the Einstein–Smoluchowski theory of Brownian motion, we refer the interested reader to [28]. For example in [22] one can find a treatment within the context of Itô’s theory of stochastic differential equations.

The ground breaking paper [4] by Dudley in 1965 was the first in which a relativistic Wiener process has been constructed. Since then a large amount of literature on diffusion processes in the frameworks of special and general relativity has been published. We refer the interested reader especially to the overview papers [7, 11], to the literature quoted there, and also to [1, 5, 6, 10, 12, 13, 14, 15, 16, 19, 20].

A construction of a relativistic Ornstein–Uhlenbeck process is provided by [3], based on a relativistic formulation of the Langevin equation. In [15, 16] Haba has studied a variety of relativistic diffusion processes, and among them also processes of Ornstein–Uhlenbeck type. The relation of these papers to the model we construct in the present manuscript has still to be worked out.

We consider the momentum of the particle as the basic dynamical quantity. Therefore we consider the relativistic four momentum

[TABLE]

of a massive point particle, and the special theory of relativity demands that $\mathbf{p}$ has to be a point of the mass shell, that is the condition

[TABLE]

has always to be fulfilled, where $m$ is the mass of the particle, and $c$ is the speed of light in vacuum. It turns out that the mass shell is a $3$ dimensional Riemannian manifold, and which therefore is equipped with a canonical, positive definite Laplace–Beltrami operator. Therefore we first construct a stochastic process of Wiener type on the mass shell, which has this Laplace–Beltrami operator as its generator, and the resulting process is a version of Dudley’s Wiener process [4]. Next, via an Itô stochastic differential equation, we add a drift of a specific form (cf. section 4) to this Wiener process in order to imitate the drift which has been introduced by Ornstein and Uhlenbeck to model the Doppler friction. As consequence, we obtain a stochastic process which we call relativistic Ornstein–Uhlenbeck momentum process, which moves on the mass shell and admits a stationary state. The relativistic Ornstein–Uhlenbeck velocity process is then defined as prescribed by special relativity, namely as the space components of the momentum process divided by the energy process (times $c^{2}$).

The plan of the article is as follows. In section 2 we discuss the mass shell as a Riemannian manifold and calculate the associated Laplace–Beltrami operator. The Wiener process on the mass shell is constructed in section 3, while the relativistic Ornstein–Uhlenbeck process is treated in section 4. In section 5 various invariant measures (or stationary states) are computed for the momentum and the velocity process. In appendix A we describe our simulation method for the stochastic differential equations.

Acknowledgement. JP gratefully acknowledges helpful discussions with Leif Döring and Andreas Neuenkirch. The authors are very much indebted to Zbiegniew Haba for pointing out the work by R.M. Dudley, J. Franchi, Y. Le Jan, and the references [15, 16].

2. The Mass Shell and its Laplacian

Throughout this article we consider a space dimension $d$, which is greater or equal to two. For points $\mathbf{p}$ in $\mathbb{R}^{1+d}$ we write their cartesian coordinates as $\mathbf{p}=(p_{0},p)$ with $p_{0}\in\mathbb{R}$, $p=(p_{1},\dotsc,p_{d})\in\mathbb{R}^{d}\!$. $\mathbb{R}^{1+d}$ is equipped with the Minkowski metric tensor $g_{M}=\text{diag}(1,-1,\dotsc,-1)$, and inner product

[TABLE]

where “$\,\cdot\,$” denotes matrix multiplication, and $(\,\cdot\,,\,\cdot\,)$ stands for the euclidean product in $\mathbb{R}^{1+d}\!$. For $m>0$ define the mass shell

[TABLE]

where $c$ is the speed of light in vacuum. From now on we shall use physical units such that $c=1$.

Clearly, if $\mathbf{p}\in\mathcal{V}_{m}^{d}$ then $p_{0}\geq m$ holds true, and for given $p_{0}>m$, $p$ belongs to the $d-1$ dimensional sphere $S^{d-1}_{\rho}$ of radius $\rho=(p_{0}^{2}-m^{2})^{1/2}$.

It is convenient to coordinatize $\mathcal{V}_{m}^{d}$ by hyperbolic coordinates $(s,\omega)\in\mathbb{R}_{+}\times S^{d-1}$, where $S^{d-1}$ denotes the $(d-1)$–dimensional unit sphere. Namely — with the exception of the apex $\mathbf{p}=(m,0,\dotsc,0)$ — every point $\mathbf{p}=(p_{0},p)$ in $\mathcal{V}_{m}^{d}$ can uniquely be written as

[TABLE]

At the apex, i.e., for $s=0$, we simply leave $\omega$ undefined.

To make this more concrete, we let $S^{d-1}$ be parametrized in the usual way by angles $\theta_{1}$, …, $\theta_{d-1}$ with $\theta_{k}\in[0,\pi)$, $k=1$, …, $d-2$, and $\theta_{d-1}\in[0,2\pi]$. Set $\theta=(\theta_{1},\dotsc,\theta_{d-1})$, and consider the mapping

[TABLE]

We may consider this mapping as an immersion of $\mathcal{V}_{m}^{d}$ into $\mathbb{R}^{1+d}$. In block form its Jacobian reads

[TABLE]

and \bigl{(}\partial\omega(\theta)/\partial\theta\bigr{)} is the Jacobian of the embedding of the unit sphere $S^{d-1}$, coordinatized by the angles $\theta$, into $\mathbb{R}^{d}$. With the immersion $\iota$ we pull the Minkowski metric $g_{M}$ back on $\mathcal{V}_{m}^{d}$ yielding a metric $g_{d}$, which in matrix form is given by

[TABLE]

The minus sign is chosen for later convenience, and the superscript “$t$” stands for transposition. Hence

[TABLE]

because $\omega(\theta)^{t}\cdot\omega(\theta)=1$, which is also the reason that the off-diagonal terms vanish. Moreover, $g_{S^{d-1}}(\theta)$ is the usual Riemannian metric tensor of the unit sphere $S^{d-1}$ in $d$ dimensions, written as a matrix parametrized by the angle variables $\theta=(\theta_{1},\dotsc,\theta_{d-1})$.

Observe that $g_{d}$ is positive definite, providing a Riemannian metric on $\mathcal{V}_{m}^{d}$. The associated volume element is

[TABLE]

where $d\sigma_{S^{d-1}}$ is the Riemannian surface element of the sphere $S^{d-1}$. The usual well-known formula (e.g., [18], [25]) for the Laplace–Beltrami operator $\Delta_{d}$ on $\mathcal{V}_{m}^{d}$ relative to $g_{d}$ yields

[TABLE]

where $\Lambda_{S^{d-1}}$ is the standard Laplace–Beltrami operator on the sphere $S^{d-1}$. For $d=2$ we find the explicit form

[TABLE]

while for the case $d=3$ of the physical Minkowksi space it reads

[TABLE]

3. Wiener Process on the Mass Shell

A natural way to define a Wiener process on a Riemannian manifold is as a stochastic process whose generator is one half times the canonical Laplace–Beltrami operator on the manifold. The history of Wiener and — more generally — diffusion processes on manifolds is quite long, and probably the first papers where those by Yosida [34, 35], and by Itô [23]. A turning point in this development was the construction by Eells, Elworthy [8], and Malliavin [27] (see also [9]), based on the rolling map of Cartan. The resulting construction of a diffusion on the orthonormal frame bundle is since then one of the basic methods, and it can be found in many textbooks such as [9, 22, 17, 21].

Dudley [4] has been the first who constructed a Wiener process on $\mathcal{V}_{m}^{d}$. His construction uses the theory of convolution semigroups on homogeneous spaces. Here we employ a different method which is as follows. The specific form (3) of the Laplace–Beltrami operator $\Delta_{d}$, namely the fact the differential operator in the $s$–variable involves no dependency on the angle variables, and the $s$–dependence of the last term only appears in form of a factor, suggests another possible construction via stochastic differential equations and a skew product, e.g., [24, Sect. 7.15], which we carry out now.

In fact, suppose that $\Theta=(\Theta_{t},\,t\in\mathbb{R}_{+})$ is a Wiener process on the sphere $S^{d-1}$ (i.e., with generator $1/2\,\Lambda_{S^{d-1}}$), and $S=(S_{t},\,t\in\mathbb{R}_{+})$ is a stochastic process on $(0,+\infty)$ with continuous paths, and generator $A_{d}$ defined by

[TABLE]

for $f\in C^{2}((0,+\infty))$. The Itô stochastic differential equation (SDE) associated with the generator $A_{d}$ is

[TABLE]

where $W$ is a standard Wiener process on the real line. Below we shall prove the existence and uniqueness of solutions of this equation. Consider the stochastic time scale

[TABLE]

then the skew product

[TABLE]

of $S$ and $\Theta$ defines a path-continuous stochastic process on $\mathcal{V}_{m}^{d}$ whose generator is $1/2\,\Delta_{d}$, i.e., a Wiener process on $\mathcal{V}_{m}^{d}$. The proof is the same as for the spherical Wiener process in Section 7.15 of [24]. In fact there one can also find a method for the construction of the Wiener process $\Theta$ on $S^{d-1}$ by successive applications of skew products. Another possibility for the construction of the Wiener process $\Theta$ on $S^{d-1}$ is Stroock’s method ([31], cf. also [21]): In this case one takes a standard Wiener process in the euclidean space $\mathbb{R}^{d}$, starts it on the embedded sphere $S^{d-1}$, and projects the infinitesimal increments of the euclidean Wiener process onto the sphere by the usual orthogonal projections. The resulting stochastic differential equation can be solved, and yields another version of the process $\Theta$.

Therefore, in order to complete our first task, viz. to construct a Wiener process on the mass shell $\mathcal{V}_{m}^{d}$, it remains to construct the process $S$ on $(0,+\infty)$ with generator $A_{d}$ as in (6). In other words, we want to prove that the Itô SDE (7) has unique solutions (in which precise sense will be clarified further below). Due to the singularity of the drift term of the SDE (7)

[TABLE]

at $s=0$, one cannot employ the standard theorems on the existence and uniqueness of SDE’s, as they can be found in, e.g., [22], [26], [30]. However, by using results in the book [2] by Cherny and Engelbert we can prove following

Theorem 3.1**.**

For every starting point $s_{0}\in(0,+\infty)$ the SDE (7) has a pathwise unique strong solution $S=(S_{t},\,t\in\mathbb{R}_{+})$ with paths which are $P$–a.s. strictly positive. Moreover, the solutions are transient in the following sense: For every $a>0$ and every initial condition $s_{0}>a$ the event $\{T_{a}=+\infty\}$, where $T_{a}$ is the first hitting time of $a$ by $S$, has strictly positive probability, and on this set $\lim_{t\to+\infty}S_{t}=+\infty$ $P$–a.s.

Remark 3.2*.*

We quickly (and somewhat roughly) recall the definition of strong existence and pathwise uniqueness — for an in depth overview of the various notions of existence and uniqueness of solutions of stochastic differential equations and their interrelations we refer the interested reader to, e.g., [2, Sect. 1.1]. The existence of a strong solution of (7) means that for any given Wiener process $W$ on the real line, there exists a solution $S$ of (7) which is adapted to the filtration generated by $W$. Thus the paths of $S$ can be considered as adapted functionals of the paths of $W$. Pathwise uniqueness means that if $S$ and $S^{\prime}$ are two solutions defined on the same probability space with the same initial condition, and with the same driving Wiener process $W$, then $P(S_{t}=S^{\prime}_{t},\,t\in\mathbb{R}_{+})=1$.

For the proof of theorem 3.1 show first two lemmas. In order to simplify our notation for the following discussion we shall temporarily set $m=1$. The first step is to prove the analogue statement as in theorem 3.1 for weak existence and uniqueness:

Lemma 3.3**.**

For every starting point $s_{0}\in(0,+\infty)$ the SDE (7) has a weak solution $S=(S_{t},\,t\in\mathbb{R}_{+})$ which is unique in law, and with paths which are $P$–a.s. strictly positive. Moreover, the solutions are transient in the following sense: For every $a>0$ and every initial condition $s_{0}>a$ the event $\{T_{a}=+\infty\}$, where $T_{a}$ is the first hitting time of $a$ by $S$, has strictly positive probability, and on this set $\lim_{t\to\infty}S_{t}=+\infty$ $P$–a.s.

Remark 3.4*.*

Also here we first want to quickly recall the meaning of the existence and uniqueness statement. That the stochastic differential equation (7) has a weak solution roughly speaking means that on some filtered probability space there exists a pair $(S,W)$ of adapted processes so that the integrated version of (7) holds true. Uniqueness in law of the solution means that if $(S,W)$ and $(S^{\prime},W^{\prime})$ are two such pairs (possibly defined on different probability spaces) with the same initial condition, then the laws of $S$ and $S^{\prime}$ coincide. Furthermore we remark in passing that the existence of a weak solution is equivalent to the existence of the associated martingale problem (e.g., [2, Theorem 1.27] or [22, Proposition IV.2.1]).

Proof of lemma 3.4.

We show that the conditions of theorems 2.16, 4.2, and part (viii) of theorem 4.6 in [2] hold true. First we remark that for every $a>0$ the drift

[TABLE]

obviously belongs to $L^{1}_{{\rm loc}}([a,+\infty))$ so that the origin $s=0$ is indeed an isolated singularity in the sense of [2, Sect. 2.1]. We fix some $a>0$ for the remainder of this proof. Next we compute the density $\rho$ of the so-called scale functions. Since in our case the diffusion coefficient is equal to $1$, $\rho$ is given by

[TABLE]

Furthermore we define the scale functions

[TABLE]

Since $d\geq 2$, we clearly have from (12) that

[TABLE]

Moreover we claim that the following are true:

[TABLE]

The integral $I_{a}$ is equal to

[TABLE]

Since $s\mapsto\sinh(s)$ is convex, we have for all $s\in(0,a]$ the inequalities

[TABLE]

Therefore we can estimate as follows

[TABLE]

For $d=2$ we get

[TABLE]

while for $d\geq 3$ we find

[TABLE]

Hence (15) is proved, and together with (13) this shows that the conditions of theorem 2.16 in [2] are fulfilled. As a consequence we get the statement that the SDE (7) has for every starting point $s_{0}>0$ a unique weak solution up to the first hitting time $T_{a}$, and the solution is a.s. strictly positive. Moreover — and this will be more important below — the singularity at $s=0$ is of type 3 in the nomenclature of [2, p. 37].

Next we show (16). We have

[TABLE]

With the inequalities

[TABLE]

where $c_{a}=(1-\exp(-2a))/2$, we obtain (recall that $d\geq 2$)

[TABLE]

and (16) is proved. Together with (13) this result shows that the hypotheses of theorem 4.2 in [2] are satisfied. This entails that for every start point $s_{0}\geq a$ there exists a unique weak solution of (7) up to the first hitting time $T_{a}$, and the solution is transient in the sense stated in the lemma. Moreover, we obtain that the behavior of the SDE (7) at infinity is of type B as defined in [2, p. 82].

Finally, with the result that the SDE (7) has a (right) singularity of type 3 and the behavior of type B at infinity, we can apply theorem 4.6.(viii) in [2], which implies the statement of the lemma. (We remark that the statement of the above quoted theorem in [2] is formulated there for a two-sided singularity, but actually the properties of the SDE on the negative half axis do not enter the statement nor its proof at all. In order to bring our situation precisely into the one discussed in chapter 4 of [2], we simply could interpret the SDE (7) as one formulated on all of $\mathbb{R}\setminus\{0\}$, and we would get the same result.) ∎

Next we show

Lemma 3.5**.**

The solutions of the SDE (7) are pathwise unique.

Proof.

This statement is a direct consequence of the fact that the drift $x\mapsto b_{0}(x)=\coth(x)$ is decreasing on $(0,+\infty)$, cf. Example 5.2.4 in [26]. ∎

Lemma 3.5 allows the application of the Yamada–Watanabe theorem [33] (cf. also [22, Theorem 1.1, Chap. IV] or [30, Theorem 1.7, Chap. IX]) which entails that we even have strong solutions, and thereby concludes the proof of theorem 3.1.

Having established the existence of Wiener processes on the mass shells $\mathcal{V}_{m}^{d}$, $d\geq 2$, we now turn to the special cases $d=2$, $d=3$, and provide more explicit expressions descriptions of theses processes in terms of stochastic differential equations in hyperbolic as well as cartesian coordinates.

3.1. The Case $d=2$

Consider formula (4) for the Laplacian on $\mathcal{V}_{m}^{2}$. Thus the associated Itô stochastic differential equations for stochastic processes $S$, and $\Phi$ in the $s$, resp. $\varphi$ coordinates are

[TABLE]

where $W^{1}$ and $W^{2}$ are independent standard one dimensional Wiener processes. Of course, the solutions of the SDE for $\Phi$ have to be taken modulo $2\pi$. We transform these equations into three dimensional cartesian coordinates using Itô calculus. A straightforward computation with Itô’s formula yields the following stochastic differential equations for the cartesian components $P=(P_{0},P_{1},P_{2})$

[TABLE]

where we have set $r(t)=\sqrt{\mathstrut P_{1}(t)^{2}+P_{2}(t)^{2}}$. An application of Itô calculus yields the associated generator in cartesian coordinates acting on smooth functions on $\mathbb{R}^{3}$:

[TABLE]

where $\partial_{i}$, $i=0$, $1$, $2$, denotes the usual partial derivative in the $i$–th coordinate direction. We want to point out the appearance of a first order term with the linear “drift” coefficient function $\mathbf{p}\mapsto 1/m^{2}\mathbf{p}$ in the generator $L$.

3.2. The Case $d=3$

From the form (5) of $\Delta_{3}$ we deduce the following system of stochastic differential equations for coordinate processes $S$, $\Theta$, $\Phi$:

[TABLE]

where $W^{1}$, $W^{2}$, and $W^{3}$ are independent a standard one dimensional Wiener processes. It is clear that also here the solutions of the equation for $\Phi$ have to taken modulo $2\pi$. A straightforward — even though somewhat lengthy — calculation with Itô’s formula gives the following stochastic differential equations in cartesian coordinates of $\mathbb{R}^{4}$:

[TABLE]

In the last equations we have set $R(t)=\sqrt{\mathstrut P_{1}(t)^{2}+P_{2}(t)^{2}+P_{3}(t)^{2}}$, and $r(t)$ is as above. The generator has in cartesian coordinates the following form

[TABLE]

and also in this case we remark the linear drift term with a linear coefficient function $\mathbf{p}\mapsto 3/2m^{2}\mathbf{p}$.

4. Relativistic Ornstein–Uhlenbeck Process

Based on the Wiener process constructed on the mass shell $\mathcal{V}_{m}^{d}$ in the previous section, we shall construct here stochastic processes on $\mathcal{V}^{d}_{m}$ which resemble the standard Ornstein–Uhlenbeck process. As we have recalled in section 1, in the usual euclidean setting the Ornstein–Uhlenbeck process is constructed by adding (via a stochastic differential equation) a linear drift term to a standard Wiener process, which pushes the Wiener process back towards the origin. As a consequence, the classical Ornstein–Uhlenbeck process has an invariant distribution which is given by a centered normal law.

Consider first the special cases $d=2$, $3$, and the SDE’s  (20), (23), for the Wiener processes on the mass shell. Clearly, one cannot simply add linear drift terms to these SDE’s, since there is no guarantee that the resulting process would continue to live on the mass shells — actually, as our computations below show, this will definitely not be the case.

Instead we introduce — for general space dimension $d$ — an additional drift term into the SDE (7) in hyperbolic coordinates. Then we take the skew product of this new process in the $s$-coordinate with a standard Wiener process on the unit sphere $S^{d-1}$ as in (8), (9). Transforming this process via (1) into a stochastic process with values in $\mathbb{R}^{1+d}$ we obtain a process which lives on the mass shell $\mathcal{V}^{d}_{m}$, if started thereon.

It turns out that a simple, natural choice for the additional drift term is given by $s\mapsto-\gamma/2m^{2}\,\tanh(s)$, where $\gamma$ is some non-negative constant. So we consider now the SDE

[TABLE]

with

[TABLE]

Hence for $\gamma>d-1$ we have a backward drift which is asymptotically constant with value $(d-\gamma-1)/2m^{2}$. The choice of this additional drift term has two advantages: For one, it turns out that in the special cases $d=2$, $d=3$ the SDE’s in cartesian coordinates will be supplemented with almost linear drifts, which are directed towards the origin and compensate the linear outward drifts which we had observed for Wiener processes in the SDE’s (20), (23). Therefore this shows some similarity with the construction of the classical Ornstein–Uhlenbeck process. Moreover, for $\gamma$ large enough this additional drift yields the existence of an invariant state for the resulting process, which can be computed explicitly (as well as some other invariant states, see section 5).

For the question of existence and uniqueness of solutions of (25) we have the following

Theorem 4.1**.**

For every initial condition $S_{0}=s_{0}\in(0,+\infty)$, the stochastic differential equation (25) has a pathwise unique, strong solution which is a.s. strictly positive for all times. For $\gamma\in[0,d-1)$ the solution is transient in the same sense as in theorem 3.1. For $\gamma\geq d-1$ the solution is recurrent in the sense that if $a>0$ and $s_{0}>a$ then $P$–a.s. $T_{a}<+\infty$.

Proof.

The proof is quite similar to the one of theorem 3.1, so we only sketch it. Again we temporarily put $m^{2}=1$ for notational simplicity. In this case the scale density $\rho$ becomes

[TABLE]

Therefore the estimations leading to the inequalities (13) and (15) are completely unaffected by the additional smooth, bounded function $s\mapsto\cosh(s)$, and we find again that the singularity at $s=0$ is of type 3. However, instead of (14) we this time get the following for any fixed $a>0$:

[TABLE]

This shows that for $\gamma\geq d-1$ we get now type A for the behavior at $+\infty$, as defined on p. 82 in [2]. For $\gamma\in[0,d-1)$ we have to estimate the present analogue of $I_{\infty}$, see (16). To this end, we use in addition to (17) the trivial bounds $1/2\,\exp(s)\leq\cosh(s)\leq\exp(s)$. The result is $I_{\infty}=+\infty$. Hence for $\gamma\in[0,d-1)$ the behavior at infinity is again of type B. Now we apply once more theorem 4.6.(viii) in [2] to conclude that for every initial condition $S_{0}=s_{0}>0$ we have the existence of a strictly positive weak solution which is unique in law.

Observe that the drift $b_{\gamma}$ is monotone decreasing on $(0,+\infty)$, so that by the same argument as in the proof of lemma 3.5 pathwise uniqueness of the solutions holds true. Another application of the Yamada–Watanabe theorem provides us with the existence of a strong solution for every initial condition $S_{0}=s_{0}\in(0,+\infty)$.

Finally we remark that theorem 4.1, [2, p. 81], states that the behavior of type A at infinity of the SDE entails that the solutions are recurrent in the sense of the theorem. ∎

As in section 3, let $\Theta=(\Theta_{t},\,t\in\mathbb{R}_{+})$ be a standard Wiener process on the $d-1$ dimensional unit sphere $S^{d-1}$, and define the stochastic time scale $\tau$ as in (8) where this time we choose for $S$ the process defined by the SDE (25). Consider the skew product

[TABLE]

We transform this process with the equations (1) into a stochastic process \mathbf{P}=\bigl{(}\mathbf{P}(t),\,t\in\mathbb{R}_{+}\bigr{)}, \mathbf{P}(t)=\bigl{(}P_{0}(t),P(t)\bigr{)}, on the mass shell $\mathcal{V}^{d}_{m}$ written in cartesian coordinates:

[TABLE]

We call the process $\mathbf{P}$ the relativistic Ornstein–Uhlenbeck momentum process in $1+d$ dimensions. The relativistic Ornstein–Uhlenbeck velocity process V=\bigl{(}V(t),\,t\in\mathbb{R}_{+}\bigr{)} in $1+d$ dimensions is then defined as

[TABLE]

(Recall that we work with physical units so that the speed of light $c$ in the vacuum is equal to $1$. In other units, we have an additional factor $c$ on the right hand side.)

Similarly as for the Wiener process, which we treated in section 3, for the cases of dimensions $d=2$ and $d=3$, we give an alternative, more explicit description of the relativistic Ornstein–Uhlenbeck processes in cartesian coordinates in terms of stochastic differential equations instead of using the skew product.

For $d=2$ we replace the first equation in (19) by (25) and transform them into cartesian coordinates with a straightforward computation using Itô’s formula. This yields the following SDE’s for the components of $\mathbf{P}$

[TABLE]

where we have set $r(t)=\sqrt{\mathstrut P_{1}(t)^{2}+P_{2}(t)^{2}}$. Thus, for $\gamma\geq 2$ the original outward drift of the Wiener process is compensated, while for $\gamma>2$ we have an effective drift towards the origin, and except for the term $\gamma/2P_{0}(t)\,dt$ this drift acts in a linear way as for the classical Ornstein–Uhlenbeck process. The additional non-linear term in the equation for $P_{0}$ takes care that the process stays on the mass shell. Note however, that this term is bounded from above by $\gamma/2m$ since $P_{0}(t)\geq m$ on the mass shell.

For $d=3$ we obtain

[TABLE]

where $R(t)=\sqrt{\mathstrut P_{1}(t)^{2}+P_{2}(t)^{2}+P_{3}(t)^{2}}$, and $r(t)$ is as in the case $d=2$ above. So in this case we have to have $\gamma\geq 3$ in order to compensate the outward drift of the Wiener process, and for $\gamma>3$ we have as above an almost linear drift term pushing the motion towards the origin.

The generators of these processes are those obtained for the Wiener process plus the additional drift terms derived above, namely for $d=2$

[TABLE]

and for $d=3$

[TABLE]

where $L_{2}$ and $L_{3}$ are as in (21),  (24) respectively. Note that $p_{0}\geq m$ so that the non-linear drift coefficients $\gamma/2p_{0}$ in the time direction are bounded from above by $\gamma/2m$, and they asymptotically vanish as $p_{0}$ tends to $+\infty$.

5. Invariant Measures

Define a measure $\mu_{d}$ on $(\mathbb{R}_{+},\mathcal{B}(\mathbb{R}_{+}))$ by

[TABLE]

Lemma 5.1**.**

For every $\gamma>d-1$ and every initial condition $S_{0}$, the solution $S=(S_{t},\,t\in\mathbb{R}_{+})$ of the stochastic differential equation (25) has the following invariant measure

[TABLE]

where

[TABLE]

Proof.

Consider the generator of $S=(S_{t},\,t\in\mathbb{R}_{+})$:

[TABLE]

with

[TABLE]

Up to the factor $1/2$, $L_{d,0}$ is the part of the Laplace–Beltrami operator (3) involving the $s$–derivatives. Therefore, by the construction of the Laplace–Beltrami operator is symmetric with respect to the measure $\mu_{d}$ on $(\mathbb{R}_{+},\mathcal{B}(\mathbb{R}_{+}))$. Hence the adjoint $L_{d,\gamma}^{*}$ of $L_{d,\gamma}$ with respect to $\mu_{d}$ acts on smooth functions as the differential operator given by

[TABLE]

Rewrite $L_{d,\gamma}^{*}$ as follows

[TABLE]

An elementary computation shows that

[TABLE]

finishing the proof. ∎

As long as $\gamma>d-1$, we may equivalently consider the function

[TABLE]

as the Lebesgue density of the invariant measure of the stochastic process $S=(S_{t},\,t\in\mathbb{R}_{+})$. This is in particular useful, when we want to compare the theoretical result of lemma 5.1 with simulations of the process. Figure 1 shows some of the results of simulation experiments we carried out, and which are described in more technical detail in appendix A. In each of these experiments we have put $m^{2}=1$, and simulated $5\times 10^{3}$ (numerical approximations of) the paths of the process $S$ for a relatively long time (see appendix A), and plotted the resulting histograms of the final positions (in blue) versus the Lebesgue density (in red) derived above. The plots show a very reasonable agreement of the theoretical density with the histograms, as could be expected.

Theorem 5.2**.**

For every $\gamma>d-1$, the stochastic process given by (28) admits the invariant measure given by

[TABLE]

Proof.

This follows directly from lemma 5.1, together with the observation that the invariant measure of the Wiener process on the unit sphere $S^{d-1}$ is the uniform law on $S^{d-1}$:

[TABLE]

where the coefficient in front of the surface element $d\sigma_{S^{d-1}}$ is the inverse of the total area of $S^{d-1}$. ∎

In a slightly informal manner the invariant measure (38) of the relativistic Ornstein–Uhlenbeck momentum process can be written in cartesian coordinates $\mathbf{p}=(p_{0},p)\in\mathbb{R}^{1+d}$ as

[TABLE]

where the constant is the same as in (38), and $\delta$ is the Dirac delta function.

For simplicity let us put $m^{2}=1$ in the sequel. From lemma 5.1 we directly get the following

Corollary 5.3**.**

For $d\geq 2$ and $\gamma>d-1$ the energy process $P_{0}$ has an invariant density $\varphi_{P_{0}}$ with respect to Lebesgue measure on $(\mathbb{R},\mathcal{B}(\mathbb{R}))$ given by

[TABLE]

The Ornstein–Uhenbeck speed process $|V|=\sqrt{\mathstrut V_{1}^{2}+\dotsb V_{d}^{2}}$ has an invariant density $\varphi_{|V|}$ with respect to Lebesgue measure on $([0,1],\mathcal{B}([0,1]))$ given by

[TABLE]

Remark 5.4*.*

As formula (41) shows, the parameter $\gamma$ must actually be chosen to be strictly larger than $d+1$ in order that the particle undergoing this process cannot attain the speed of light with strictly positive probability (cf. also figure 3a).

Figures 2 and 3 show the long term histograms of the final values of $5\times 10^{3}$ simulated paths of $P_{0}$, $|V|$ respectively, in comparison with the marginal densities (40), (41) respectively.

For the remainder of this section we assume in addition that $d=3$, i.e., that we are in the physical Minkowski space. Then it is straightforward to compute also the marginal invariant densities of the momentum processes $P_{i}$, $i=1$, $2$, $3$, explicitly:

Corollary 5.5**.**

For $d=3$ and $\gamma>2$ the cartesian components $P_{i}$, $i=1$, $2$, $3$, of the Ornstein–Uhlenbeck momentum process have marginal invariant densities $\varphi_{P_{i}}$ with respect to Lebesgue measure on $(\mathbb{R},\mathcal{B}(\mathbb{R}))$ given by

[TABLE]

where $n_{\gamma}$ is the normalization constant.

Figure 4 illustrates the result of corollary 5.5 with simulations of the value of $P_{1}$ for large times, $d=3$ and various values of $\gamma$.

Also for components $V_{i}$, $i=1$, $2$, $3$, of the Ornstein–Uhlenbeck velocity process $V$ (see (30)) it is straightforward to calculate their marginal invariant densities:

Corollary 5.6**.**

For $d=3$ and $\gamma>2$ the cartesian components $V_{i}$, $i=1$, $2$, $3$, of the Ornstein–Uhlenbeck velocity process have marginal invariant densities $\varphi_{V_{i}}$ relative to Lebesgue measure given by a symmetric Beta law on $[-1,1]$ with parameter $\gamma/2$. Explicitly:

[TABLE]

Figure 5 shows histograms of $V_{2}$ resulting from the simulation of the momentum process, where the red line is the graph of the Beta density (43).

Remark that for small velocities $v$, the density in (43) is (up to normalization) close to a centered normal distribution:

[TABLE]

so that at least for the cartesian components of the velocity process we obtain a certain compatibility with the invariant density of the classical, non-relativistic Ornstein–Uhlenbeck process.

Appendix A Simulations

Consider the stochastic differential equations (7) and (25). As argued in sections 3 and 4, the solutions do not leave the interval $(0,+\infty)$ when started there. However, if one tries to simulate paths with a naive scheme, such as the Euler–Maruyama scheme, it is not possible to prevent all paths from crossing the singularity of the drift at $s=0$ into the region $(-\infty,0)$. The reason is of course, that one actually simulates a random walk, and the discrete increments do have the possibility to cross the singularity at $s=0$. A simulation scheme, called backward Euler–Maruyama (BEM) scheme, which does prevent this crossing has been provided by Neuenkirch and Szpruch in [29]. The conditions formulated in [29] for their results to hold are fulfilled by the SDE’s considered in the present article. The scheme is of the form

[TABLE]

where $b$ is the drift, and the increments $\Delta W_{t+\Delta t}$ of the Wiener process are — as usual — independent centered normal variates with variance $\Delta t$. Observe that in order to compute an increment of $s$ from one time step to the next, one has to numerically solve an implicit problem. We implemented this scheme in Scilab, 111http://www.scilab.org and had to observe that sometimes we still obtained paths which crossed the singularity of the drift at the origin. A careful analysis showed that this is due to the fact that Scilab’s fsolve routine does not in all cases find the correct solution of the implicit problem. We believe that this is so because probably Scilab’s fsolve is based on Newton’s method, which is well-known to fail under certain circumstances. In order to get a functioning scheme for the SDE’s (7), (25), we therefore supplemented Scilab’s fsolve with a bisection method for those cases, where a jump across the singularity had occured. Actually, a similar consideration had to be done for the SDE of $\Theta$ in (22), in which case the drift has singularities at $\theta=0$ and $\theta=\pi$.

For each of the histograms in figures 1 and 2, we generated with the above described method samples of $5\times 10^{3}$ paths, with $m^{2}=1$, $\Delta t=2^{-6}$, and let the paths develop until time $\tau=100$, i.e., altogether over $2^{6}\times 100$ time steps.

Afterword by JP. During the work on this manuscript my mentor and coauthor Robert Schrader passed away. Robert was a truly outstanding scientist, a charismatic teacher, and a wonderful colleague and friend — I miss him very much.