Skew Brownian motion with dry friction: The Pugachev-Sveshnikov approach

TL;DR

This paper studies a generalized skew Brownian motion model incorporating dry friction, analyzing its probability distribution and occupation time using the Pugachev-Sveshnikov approach, with applications across physics and mathematics.

Contribution

It introduces a comprehensive analysis of skew Brownian motion with dry friction using the Pugachev-Sveshnikov method, extending existing models in stochastic processes.

Findings

Derived the probability distribution of the process.

Analyzed the occupation time on the positive half line.

Provided insights into the process's behavior under dry friction.

Abstract

The Caughey-Dieness process, also known as the Brownian motion with two valued drift, is used in theoretical physics as an advanced model of the Brownian particle velocity if the resistant force is assumed to be dry friction. This process also appears in many other fields, such as applied physics, mechanics, astrophysics, and pure mathematics. In the present paper we are concerned with a more general process, skew Brownian motion with dry friction. The probability distribution of the process itself and of its occupation time on the positive half line are studied. The approach based on the Pugachev-Sveshnikov equation is used.

Equations

$$dX(t)=-2\mu \mathrm{sign}\left(X(t)\right)dt+(2\alpha -1)d{L}_X^0(t)+\sqrt{2}dW(t),t>0,X(0)=0,$$

$$L_X^0(t)=\underset{\epsilon \to +0}{\lim }\frac{1}{2\epsilon }{\int }_0^t\text{mathds}{1}_{(-\epsilon ,+\epsilon )}(X(s))d[X{]}_s,$$

$$\mathcal{I}(t)={\int }_0^t\text{mathds}{1}_{(0,+\infty )}(X(s))ds.$$

$$\{\begin{array}{rl}& dX(t)=-2\mu \mathrm{sign}\left(X(t)\right)dt+(2\alpha -1)d{L}_X^0(t)+\sqrt{2}dW(t),\\ & d\mathcal{I}(t)=\text{mathds}{1}_{(0,+\infty )}(X(s))ds,X(0)=\mathcal{I}(0)=0.\end{array}$$

$$\frac{\partial E}{\partial t}+(z_1^2-i{z}_2/2)E+(2\mu z_1-z_2/2)\hat{E}-2i(2\alpha -1)z_1{\Psi }_0=0,E(z_1,z_2;0)=1,$$

$$\hat{E}=\frac{1}{\pi }v.p.\underset{-\infty }{\overset{+\infty }{\int }}\frac{E{|}_{z_1=s}}{s-z_1}ds,{\Psi }_0=\frac{1}{2\pi }v.p.\underset{-\infty }{\overset{+\infty }{\int }}E{|}_{z_1=s}ds.$$

$$\Phi (\zeta ,z_2;t)=\frac{1}{2\pi i}\underset{-\infty }{\overset{+\infty }{\int }}\frac{E{|}_{z_1=s}}{s-\zeta }ds,{\Phi }^{\pm }(z,z_2;t)=\underset{\zeta \to z\pm i0}{\lim }\Phi (\zeta ,z_2;t),\mathrm{Im}z=0.$$

$${\Phi }^{+}-{\Phi }^{-}=E,{\Phi }^{+}+{\Phi }^{-}=-i\hat{E}.$$

$$(z_1^2+2\mu i{z}_1+p-i{z}_2){\tilde{\Phi }}^{+}-i(2\alpha -1)z_1{\tilde{\Psi }}_0-\frac{1}{2}=(z_1^2-2\mu i{z}_1+p){\tilde{\Phi }}^{-}+i(2\alpha -1)z_1{\tilde{\Psi }}_0+\frac{1}{2}.$$

$${\tilde{\Phi }}^{\pm }=\frac{G_0+z_1{G}_1\pm i(2\alpha -1)z_1{\tilde{\Psi }}_0\pm 1/2}{z_1^2\pm 2i\mu z_1+p-(1\pm 1)i{z}_2/2}.$$

$$G_0+2i{\nu }^{+}\alpha G_1+1/2=0,G_0+2i{\varkappa }^{-}(1-\alpha )G_1-1/2=0.$$

$${\tilde{E}}_X(z_1,p)=\tilde{E}(z_1,0;p)=\frac{1}{i{\varkappa }^{-}}\left(\frac{\alpha }{z_1+i{\varkappa }^{+}}-\frac{1-\alpha }{z_1-i{\varkappa }^{+}}\right),$$

$${\tilde{E}}_{\mathcal{I}}(z_2,p)=\tilde{E}(0,z_2;p)=\frac{\alpha {\varkappa }^{+}-(1-\alpha ){\nu }^{-}}{{\nu }^{-}{\varkappa }^{+}((1-\alpha ){\varkappa }^{-}-\alpha {\nu }^{+})}.$$

$$f_X(x,t)=\left(\frac{1}{\sqrt{\pi t}}{e}^{-\frac{(|x|+2\mu t{)}^2}{4t}}+\mu e^{-2\mu |x|}\mathrm{Erfc}\frac{|x|-2\mu t}{2\sqrt{t}}\right)\cdot \{\begin{array}{rlrl}& \alpha ,& & x>0,\\ & 1-\alpha ,& & x<0.\end{array}$$

$$f_X^{\infty }(x)=f_X(x,+\infty )=2\mu e^{-2\mu |x|}\{\begin{array}{rlrl}& \alpha ,& & x>0,\\ & 1-\alpha ,& & x<0.\end{array}$$

$$S(w_1,w_2)=\frac{\alpha (w_2+\mu )+(1-\alpha )(w_1+\mu )}{(w_1+\mu )(w_2+\mu )(\alpha (w_1-\mu )+(1-\alpha )(w_2-\mu ))},$$

$${\tilde{E}}_{\mathcal{I}}(z_2,p)=S(\sqrt{{\mu }^2+p-i{z}_2},\sqrt{{\mu }^2+p}).$$

$$f_{\mathcal{I}}(y,t)={\mathfrak{L}}_{p,t}^{-1}\left[{\mathfrak{L}}_{q,y}^{-1}\left[S(\sqrt{{\mu }^2+p+q},\sqrt{{\mu }^2+p})\right]\right],$$

$$f_{\mathcal{I}}(y,t)={\mathfrak{L}}_{p,t}^{-1}\left[{\mathfrak{L}}_{q,y}^{-1}\left[S(\sqrt{{\mu }^2+p+q},\sqrt{{\mu }^2+p})\right]\right]=e^{-{\mu }^{2}t}{\mathfrak{L}}_{p,t}^{-1}\left[{\mathfrak{L}}_{q,y}^{-1}\left[S(\sqrt{p+q},\sqrt{p})\right]\right]$$

$$=e^{-{\mu }^{2}t}{\mathfrak{L}}_{p,t}^{-1}\left[e^{-py}{\mathfrak{L}}_{q,y}^{-1}\left[S(\sqrt{q},\sqrt{p})\right]\right]=\chi (t-y)e^{-{\mu }^{2}t}{\mathfrak{L}}_{p,t-y}^{-1}\left[{\mathfrak{L}}_{q,y}^{-1}\left[S(\sqrt{q},\sqrt{p})\right]\right]$$

$$=\frac{\chi (t-y)e^{-{\mu }^{2}t}}{4\pi (y(t-y){)}^{3/2}}\underset{0}{\overset{+\infty }{\int }}\underset{0}{\overset{+\infty }{\int }}{\mathfrak{L}}_{p,s_2}^{-1}\left[{\mathfrak{L}}_{q,s_1}^{-1}\left[S(q,p)\right]\right]s_1{s}_2{e}^{-\frac{s_1^2}{4y}-\frac{s_2^2}{4(t-y)}}d{s}_{1}d{s}_2$$

$$=\frac{4\chi (t-y)e^{-{\mu }^{2}t}}{\pi \sqrt{y(t-y)}}\underset{0}{\overset{+\infty }{\int }}\underset{0}{\overset{+\infty }{\int }}{\mathfrak{L}}_{p,2\sqrt{t-y}{s}_2}^{-1}\left[{\mathfrak{L}}_{q,2\sqrt{y}{s}_1}^{-1}\left[S(q,p)\right]\right]s_1{s}_2{e}^{-s_1^2-s_2^2}d{s}_{1}d{s}_2.$$

$$f_{\mathcal{I}}(y,t)=\frac{4e^{-{\mu }^{2}t}}{\pi \sqrt{y(t-y)}}\underset{0}{\overset{+\infty }{\int }}\underset{0}{\overset{+\infty }{\int }}\chi (2\sqrt{y}{s}_1,2\sqrt{t-y}{s}_2)s_1{s}_2{e}^{-s_1^2-s_2^2}d{s}_{1}d{s}_2,0<y<t,$$

$$\chi (s_1,s_2)=\frac{1-\alpha }{\alpha }{e}^{-\mu \left(s_1-\frac{2-\alpha }{\alpha }{s}_2\right)}{\chi }^{+}(s_1,s_2)+\frac{\alpha }{1-\alpha }{e}^{\mu \left(\frac{1+\alpha }{1-\alpha }{s}_1-s_2\right)}{\chi }^{-}(s_1,s_2),$$

$${\chi }^{+}(s_1,s_2)=\text{mathds}{1}_{(0,+\infty )}(\alpha s_1-(1-\alpha )s_2),{\chi }^{-}(s_1,s_2)=1-{\chi }^{+}(s_1,s_2).$$

$$f_X(x,t)=\left(\frac{1}{\sqrt{\pi t}}{e}^{-\frac{(|x|+2\mu t{)}^2}{4t}}+\mu e^{-2\mu |x|}\mathrm{Erfc}\frac{|x|-2\mu t}{2\sqrt{t}}\right)\cdot \{\begin{array}{rlrl}& \alpha ,& & x>0,\\ & 1-\alpha ,& & x<0,\end{array}$$

$$f_X^{\infty }(x)=f_X(x,+\infty )=2\mu e^{-2\mu |x|}(\alpha \text{mathds}{1}_{(0,+\infty )}(x)+(1-\alpha )\text{mathds}{1}_{(-\infty ,0]}(x)),$$

$$f_{\mathcal{I}}(y,t)=\frac{4e^{-{\mu }^{2}t}}{\pi \sqrt{y(t-y)}}\underset{0}{\overset{+\infty }{\int }}\underset{0}{\overset{+\infty }{\int }}\chi (2\sqrt{y}{s}_1,2\sqrt{t-y}{s}_2)s_1{s}_2{e}^{-s_1^2-s_2^2}d{s}_{1}d{s}_2,0<y<t$$

$$\chi (s_1,s_2)=\frac{1-\alpha }{\alpha }{e}^{-\mu \left(s_1-\frac{2-\alpha }{\alpha }{s}_2\right)}{\chi }^{+}(s_1,s_2)+\frac{\alpha }{1-\alpha }{e}^{\mu \left(\frac{1+\alpha }{1-\alpha }{s}_1-s_2\right)}{\chi }^{-}(s_1,s_2),$$

$${\chi }^{+}(s_1,s_2)=\text{mathds}{1}_{(0,+\infty )}(\alpha s_1-(1-\alpha )s_2),{\chi }^{-}(s_1,s_2)=1-{\chi }^{+}(s_1,s_2).$$

Full text

1 Introduction

Brownian motion plays an important role in statistical physics and other applied areas of science. The classical physical theory of Brownian motion was developed by Einstein and Smoluchovwski[1] in the beginning of 20th century. In this theory, a weightless Brownian particle is driven by the random force resulting from the collision with molecules. Also, the friction between this particle and the molecules, in other words the resistant force due to collision, is assumed to be viscous, that is, proportional to the particle’s velocity. Mathematically, these assumptions lead to the Wiener process. The major drawback of this model is that almost all the trajectories of the particle are continuous but nowhere differentiable, so the velocity cannot be defined properly. Later on, this problem was solved by Ornstein and Uhlenbeck [2], who developed a refined physical theory that accounts for the particle’s inertia.

A more evolved model of Brownian motion uses an assumption that along with viscous, dry friction is present, which is independent of the particle’s speed, but depends on the motion direction. Such a model in its simplest form, when the viscous friction is absent, is called Brownian motion with dry friction and was first studied by Caughey and Dienes [3] in 1960s. They considered a process similar to Ornstein–Uhlenbeck’s but the linear term was replaced by the sign function. This process turned out to be useful in analysis of different phenomena in control theory [4], seismic mechanics [5], communication systems theory [6], radio physics [7], nonlinear stochastic dynamics [8], and also in pure mathematics [11, 10, 9]. In 2000s, other applications of the Caughey–Dienes process emerged. The major interest was in nanofrictional systems [12], particles separation [13], ratchets [14, 15], granular motors [16], dynamics of granular media [18, 17], and in dynamics of droplets on moving surfaces [22, 19, 20, 21]. Some additional publications on the subject can be found in the author’s work [23].

The present paper deals with the model of skew Brownian motion with dry friction, or more precisely with the so-called skew Caughey–Dienes process, which can be understood as the velocity of a particle in this model. In what follows, for simplicity, we limit ourselves to studying skewing at zero only.

In order to understand how skewing works we consider an excursion of the process $V(t)$; an excursion is the part of the process’s trajectory, located between two time points at which $V(t)$ vanishes, that is to say, between two consequent stops of a particle. An excursion of the Caughey–Dienes process without skewing is known to be a sufficiently good approximation of the real Brownian particle behavior. However, note that for such a process, due to symmetry, positive and negative excursions appear with the same probability $0.5$, which is not always the case. Generally speaking, in reality the probability $\alpha$ to have a positive excursion is different from $0.5$, that is, some skewing takes place. In this way we come to the model of the skew Caughey–Dieness process, where the probability $\alpha$ is an additional parameter of the model, which has to be estimated from the experimental data. Note that due to the above interpretation of the skewing, it is not that difficult to do.

Skew processes are tightly connected with diffusion in discontinuous media and has numerous applications in several fields of physics [24]. For example, they describe shock acceleration of charged particles in a magnetic field [25] and diffusion in geophysical problems associated with the study of inhomogeneous porous media [26]. It should also be noted that such processes arise in some special randomly perturbed Hamiltonian systems [27]. Besides, the skewing procedure was studied from purely mathematical standpoint for a number of typical stochastic processes. The most detailed results are available for the process of skew Brownian motion [24], for the Ornstein–Uhlenbeck process [28], and for the Cox–Ingersoll–Ross process [29]. Unfortunately, the skew Caughey–Dienes process is not that well known, and this, in particular, motivates our study.

Generally, calculating the probabilistic characteristics of a skew diffusion is considered to be a difficult mathematical problem [10]. The standard tools to solve it are the Fokker–Planck–Kolmogorov equation, Feynman-Kac equation, or the random walk approximation. We propose an alternative way, based on the Pugachev–Sveshnikov equation, which was used in authors’ previous works [23, 30]. From our perspective, this approach allows one to get to the result faster, and in a more algorithmic way.

2 Skew Caughey–Dienes process

Skewing is intimately connected with the notion of local time, and it is possible to proof that the skew Caughey–Dienes process described in the introduction can be represented as a unique strong solution [24] of the following stochastic differential equation

[TABLE]

where $\alpha\in(0,1)$ is the skewing parameter, the probability of an excursion to be positive. By $W(t)$ we denote a standard Wiener process starting at zero, and $L_{X}^{0}(t)$ is the symmetric local time of $X(t)$ at the level zero:

[TABLE]

where $[X]_{s}=2s$ is the so-called quadratic variation of $X(s)$. Further, we will be interested in the positive half-line occupation time of $X(t)$

[TABLE]

One can think of $X(t)$ and $\mathcal{I}(t)$ as of the components of the vector diffusion process $(X(t),\mathcal{I}(t))$ governed by the system of SDEs

[TABLE]

In what follows, we derive explicit formulas for the probability density function of (4), following ideas from [30].

It can be shown [30] that the characteristic function $E(z_{1},z_{2};t)=\mathbb{E}\!\left[e^{i(z_{1}X(t)+z_{2}\mathcal{I}(t))}\right]$ of the process $(X(t),\mathcal{I}(t))$ satisfies the equation

[TABLE]

where we adopt the short notation $\hat{E}(z_{1},z_{2};t)$ and $\Psi_{0}(z_{2},t)$:

[TABLE]

For $\mathrm{Im}\ \zeta\neq 0$ let us introduce the Cauchy-type integral $\Phi(\zeta,z_{2};t)$ and its limit values $\Phi^{\pm}(z,z_{2};t)$ on the real axis from upper and lower half-planes (with respect to the first argument):

[TABLE]

It is well known that $\Phi(\cdot,z_{2};t)$ is analytic when $\mathrm{Im}\ \zeta\neq 0$, and that $\Phi^{\pm}$ satisfy the Sokhotski–Plemelj formulas when $\mathrm{Im}\ z=0$:

[TABLE]

Clearly, one can rewrite (5) in terms of $\Phi^{\pm}$, that gives the condition of the Riemann boundary value problem, and now we need to recover analytic functions $\Phi^{\pm}$, from this condition. Applying the Laplace transform with respect to $t$, and denoting its argument by $p$, we get to the formula

[TABLE]

The Laplace transforms are labeled with tildes above the functions.

Note that the left-hand side of (9) can be analytically continued for all $z_{1}\in\mathbb{C}$ such that $\mathrm{Im}\ z_{1}>0$, also the right-hand side can be analytically continued for all $z_{1}\in\mathbb{C}$ such that $\mathrm{Im}\ z_{1}<0$. Since they match when $\mathrm{Im}\ z_{1}=0$, they turn out to be the elements of the same entire function of argument $z_{1}\in\mathbb{C}$ . Assuming that $\Phi^{\pm}(z,z_{2};t)=O(\frac{1}{|z|})$ when $z\to\infty$ for $\mathrm{Im}\ z\gtrless 0$, by generalized Liouville’s theorem one can realize that this entire function is actually linear: $G_{0}(z_{2},t)+z_{1}G_{1}(z_{2},t)$. This leads to the equality

[TABLE]

Note that the denominator in (10) has zeros $i\nu^{\pm}=i(-\mu\pm\sqrt{\mu^{2}+p-iz_{2}})$ and $i\varkappa^{\pm}=i(\mu\pm\sqrt{\mu^{2}+p})$ such that $\mathrm{Im}\ \!\!(i\nu^{\pm})\gtrless 0$ and $\mathrm{Im}\ \!\!(i\varkappa^{\pm})\gtrless 0$. At the same time, $\tilde{\Phi}^{\pm}$ should be analytic in upper and lower half-planes, therefore, the singularities at $i\nu^{+}$ and $i\varkappa^{-}$ are removable. This gives a system of linear equations to determine $G_{0}$ and $G_{1}$. Also, taking into account the definition of $\Psi_{0}$ in (6) and performing integration in (10) one can find that $\tilde{\Psi}_{0}=-iG_{1}$. After that, the system of linear equations for $G_{0}$ and $G_{1}$ can be written in the following form

[TABLE]

Finally, substituting $G_{0}$ and $G_{1}$ from (11) into (10), one can get to the final expression for $\tilde{E}$, using the first of the Sokhotski–Plemelj formulas (8). Particularly, after necessary simplifications one can get the Laplace transform of the characteristic function of $X(t)$ and $\mathcal{I}(t)$:

[TABLE]

Straightforward computations using the table of Laplace and Fourier transforms for the first expression in (12) gives us the probability density function of $X(t)$:

[TABLE]

Letting $t$ tend to $+\infty$, one can obtain the steady-state probability density function

[TABLE]

The plot of the probability density function of $X(t)$ for different $\alpha$ is given in Fig. 1. Note that for $\alpha\neq 0.5$ the curve has a jump at zero, so that the probability that $X(t)>0$ is equal to $\alpha$.

Handling the second line in (12) takes somewhat more effort. First, we introduce the function:

[TABLE]

then the expression for $\tilde{E}_{\mathcal{I}}(z_{2},p)$ can be written as

[TABLE]

After that, we note that since $\mathcal{I}$ is non-negative, instead of using the inverse Fourier transform with respect to $z_{2}$, one can use the inverse Laplace transform with respect to $q=-iz_{2}$:

[TABLE]

where $\mathfrak{L}^{-1}_{p,t}\left[\cdot\right]$ is the inverse Laplace transform with respect to the parameter $p$, and the argument of the original is $t$.

Now we write the following chain of equalities:

[TABLE]

In the first and second lines we used well-known properties of the Laplace transform, in the line three we applied Efros’s theorem [31], and in the last line we made a change of variables in the double integral.

As the function $S$ in (15) is rational, it is easy to calculate $\mathfrak{L}^{-1}_{p,2\sqrt{t-y}s_{1}}\left[\mathfrak{L}^{-1}_{q,2\sqrt{y}s_{2}}\left[S(p,q)\right]\right]$. Skipping cumbersom but trivial in nature computations, we get to the final formula for the probability density function of $\mathcal{I}$

[TABLE]

where $\mathrm{Erfc}\ \!(\cdot)$ is the complementary error function, and

[TABLE]

Not performing any other simplifications of this formula, we illustrate the generic shape of the scaled occupation time $\mathcal{T}(t)=\mathcal{I}(t)/t$ in Fig. 2. The probability density function of $\mathcal{T}(t)$ is given by $f_{\mathcal{T}}(y,t)=tf_{\mathcal{I}}(ty,t)$.

It is worth noticing that unlike the well-known arcsine law, the distribution in Fig. 2 is unimodal, and its mode can be controlled by the parameter $\alpha$. This completely agrees with the physical interpretation of the skew Caughey–Dienes process given in the introduction.

At the end of the section, we formulate our result in the form of the following theorem.

Theorem**.**

The PDF of $X(t)$, the steady-state PDF of $X(t)$, and the PDF of the positive half-line occupation time $\mathcal{I}(t)$ are given by:

[TABLE]

[TABLE]

[TABLE]

where $\mathrm{Erfc}\ \!(\cdot)$ is the complementary error function, and

[TABLE]

3 Conclusions

We derived explicit formulas for the probability density function of the skew Caughey–Dienes process and its occupation time on the positive half-line, which generalizes known results for the regular Caughey–Dienes process. In fact, more general result was obtained for the Laplace transform of the joined characteristic function. Essentially, our approach is based on the reduction to a Riemann boundary value problem, and clearly it can be used to find the characteristics of more general SDEs with piecewise linear coefficients and local time.