Exact power spectra of Brownian motion with solid friction

We continue in this paper our study of Brownian motion involving solid (dry or Coulomb) friction in addition to viscous friction; see [1–3]. As in these works, we consider the piecewise linear Langevin equation

$$\begin{equation} \dot{v}=-\gamma v - \Delta \sigma (v) +F + \sqrt{\Gamma } \xi (t), \end{equation} \tag{ 1 }$$

where γ > 0 denotes the viscous coefficient, Δ > 0 the dry friction coefficient, F an external constant forcing, ξ(t) a Gaussian white noise and Γ the related diffusion constant. The term −Δσ(v), where σ(v) denotes the sign of v with the convention σ(0) = 0, represents the dry friction force. Its physical interpretation follows by considering equation (1) in the deterministic limit Γ = 0: for |F| < Δ, the stationary state of this equation is the 'sticking' state v = 0, whereas for |F| > Δ, the stationary state is a 'sliding' state with v ≠ 0. Thus, to induce motion (v ≠ 0) from rest (v = 0), F has to be larger than Δ, the dry friction contact force.¹

The effect of dry friction on the properties of Brownian motion was studied by de Gennes [4], who showed, for the special case γ = F = 0, that the velocity–velocity correlation function 〈v(t)v(0)〉 acquires a dependence on the noise power when Δ > 0. In [3], we extended his study by obtaining eigenfunction expansions of the propagator p(v, t|v₀, 0) for the general case γ > 0 and F > 0. From these expansions, which involve a special function known as the parabolic cylinder function, we were also able to compute 〈v(t)v(0)〉. Our main finding was that the stick and slip states of the deterministic system (Γ = 0) translate into stick and slip regimes in the noisy system (Γ ≠ 0), which are characterized by a strong and weak dependence, respectively, of the correlation time of the velocity–velocity correlation function with the external force F.

Our goal in this paper is to show that our exact results of [3] can be expressed in a more convenient way by solving the Fokker–Planck equation associated with equation (1) in Laplace space rather than in direct space. The resulting expressions for the propagator and velocity–velocity correlation function are indeed somewhat more compact and more elegant than the eigenfunction expansions presented earlier. Our solution in Laplace space also enables us to complete the study of equation (1) by deriving exact, closed-form expressions for the power spectrum of this equation, i.e. the Fourier transform of the velocity–velocity correlation function. These expressions extend early articles by Caughey and Dienes [5] and by Atkinson and Caughey [6, 7], recently brought to our attention, which considered the power spectrum of the Langevin equation with pure Coulomb friction, i.e. the same case (γ = F = 0) considered by de Gennes.

As in [3], we study the propagator P(x, t|x', 0) of equation (1), expressed in terms of the non-dimensional variables $\sqrt{2\gamma /\Gamma } v \rightarrow x$, γt → t. The Fokker–Planck equation governing the evolution of this propagator has the form

$$\begin{equation} \frac{\partial P}{\partial t} = \frac{\partial }{\partial x} (x+\delta \sigma (x)-f) P+ \frac{\partial ^2 P}{\partial x^2}, \end{equation} \tag{ 2 }$$

where $\delta = \Delta \sqrt{2/(\gamma \Gamma )}$ measures the magnitude of the dry friction relative to the viscous damping while $f=F/\sqrt{2/(\gamma \Gamma )}$ stands for the external constant force, measured also against the viscous damping.

The time-independent or stationary solution of the Fokker–Planck equation has the standard form

$$\begin{equation} \rho _f(x)=\frac{{\rm e}^{-\Phi (x)}}{Z} \end{equation} \tag{ 3 }$$

in terms of the potential

$$\begin{equation} \Phi (x)=\frac{(|x|+\delta )^2}{2}-f x. \end{equation} \tag{ 4 }$$

Depending on the sign of the parameter δ, the potential refers to a dry friction problem (δ > 0) or a Kramer-type tunneling problem (δ < 0).²

For piecewise-linear Fokker–Planck equations, the propagator can be obtained exactly, as was observed in [6, 7], by considering its Laplace transform:

$$\begin{equation} \tilde{P}(x,s|x^{\prime },0)=\int _0^\infty {\rm e}^{-s t} P(x,t|x^{\prime },0) \,{\rm d}t. \end{equation} \tag{ 5 }$$

With this transform, equation (2) becomes a second-order ordinary differential equation

$$\begin{equation} \fl s \tilde{P}(x,s|x^{\prime },0) -\delta (x-x^{\prime }) = \frac{{\rm d}}{{\rm d} x} (x+\delta \sigma (x)-f) \tilde{P}(x,s|x^{\prime },0)+ \frac{{\rm d}^2 \tilde{P}}{{\rm d} x^2}, \end{equation} \tag{ 6 }$$

which can be solved in terms of parabolic cylinder functions. The detail of this solution is given in appendix A. The final result obtained for positive values of the initial condition, i.e. x' > 0, has for expression

$$\begin{eqnarray} &&\fl \tilde{P}(x,s|x^{\prime },0) = g_<(s,\delta ,f)\, \frac{\Gamma (s)}{\sqrt{2\pi }}\, {\rm e}^{(x^{\prime }+\delta -f)^2/4}\, D_{-s}(x^{\prime }+\delta -f)\nonumber \\ &&\times {\rm e}^{-(x-\delta -f)^2/4}\, D_{-s}(-x+\delta +f) \end{eqnarray} \tag{ 7 }$$

for x < 0,

$$\begin{eqnarray} &&\fl \tilde{P}(x,s|x^{\prime },0) = g_>(s,\delta ,f) \, \frac{\Gamma (s)}{\sqrt{2\pi }} \, {\rm e}^{(x^{\prime }+\delta -f)^2/4}\, D_{-s}(x^{\prime }+\delta -f)\nonumber \\ &&\times \,{\rm e}^{-(x+\delta -f)^2/4}\, D_{-s}(x+\delta -f)\nonumber \\ &&+\,\frac{\Gamma (s)}{\sqrt{2\pi }}\, {\rm e}^{(x^{\prime }+\delta -f)^2/4}\, D_{-s}(x^{\prime }+\delta -f)\nonumber \\ &&\times\, {\rm e}^{-(x+\delta -f)^2/4}\, D_{-s}(-x-\delta +f) \end{eqnarray} \tag{ 8 }$$

for 0 < x < x', and

$$\begin{eqnarray} &&\fl \tilde{P}(x,s|x^{\prime },0) = g_>(s,\delta ,f) \, \frac{\Gamma (s)}{\sqrt{2\pi }}\, {\rm e}^{(x^{\prime }+\delta -f)^2/4}\, D_{-s}(x^{\prime }+\delta -f) \nonumber \\ &&\times\, {\rm e}^{-(x+\delta -f)^2/4} D_{-s}(x+\delta -f)\nonumber \\ &&+\,\frac{\Gamma (s)}{\sqrt{2\pi }}\, {\rm e}^{(x^{\prime }+\delta -f)^2/4}\, D_{-s}(-x^{\prime }-\delta +f)\nonumber \\ &&\times \,{\rm e}^{-(x+\delta -f)^2/4}\, D_{-s}(x+\delta -f) \end{eqnarray} \tag{ 9 }$$

for x > x'. The coefficients g_< and g_> are defined in appendix A by equations (A.21) and (A.23), respectively. The propagator for negative values of the initial condition, i.e. x' < 0, is obtained from the equations above simply by replacing x, x', and f by −x, −x', and −f, respectively.

The set of equations (7)–(9) is the main result of this paper. To gain some insight into the form of this solution, consider the case without dry friction and external forcing, i.e. δ = 0 and f = 0, which corresponds to the Ornstein–Uhlenbeck process. Then g_<(s, 0, 0) = 1 and g_>(s, 0, 0) = 0, according to equations (A.21)–(A.23) of appendix A, and we are led to³

$$\begin{equation} \fl \tilde{P}(x,s|x^{\prime },0) = \frac{\Gamma (s)}{\sqrt{2 \pi }} \left\lbrace \begin{array}{@{}lcl@{}}{\rm e}^{(x^{\prime })^2/4} D_{-s}(x^{\prime }) \,{\rm e}^{-x^2/4} D_{-s}(-x), & & x<x^{\prime } \\ {\rm e}^{(x^{\prime })^2/4} D_{-s}(-x^{\prime }) {\rm e}^{-x^2/4} D_{-s}(x), & & x>x^{\prime }, \end{array} \right. \end{equation} \tag{ 10 }$$

which is the Laplace transform of the Ornstein–Uhlenbeck propagator. Comparing this particular form of $\tilde{P}$ with the general solution above, we see that the products of parabolic cylinder functions occurring in equations (7)–(9) are essentially the Laplace transform of the Ornstein–Uhlenbeck process, modified to include the dry friction and external forces. Following the method of images, the additional term with the coefficients g_< and g_> can then be considered, at a superficial level, as convolution integrals in the time domain, with g_< and g_> playing the role of the source terms.

From the relatively compact solution of the propagator in Laplace space, we are not able to obtain the propagator itself in closed form, except for special cases, such as the Ornstein–Uhlenbeck process (δ = f = 0) and the pure dry friction case corresponding to γ = F = 0. However, for all cases it is possible to obtain the propagator numerically by inverting the Laplace transform. Figure 1 shows the result of this procedure using the so-called Talbot method [8–10].⁴ This figure reproduces exactly our previous results for the propagator based on the spectral decomposition of the Fokker–Planck solution (see figure 10 of [3]). In general, we have found that the numerical computation of P(x, τ|x', 0) from the inverse Laplace transform of $\tilde{P}(x,s|x^{\prime },0)$ is stable and can be carried out to arbitrary level of accuracy for a large range of physically-relevant parameter values. Our Laplace solution can therefore be considered a useful complement to the spectral solution presented in [3].

Zoom In Zoom Out Reset image size

Using the closed-form solution for $\tilde{P}(x,s|x^{\prime },0)$, we now derive an expression for the Laplace transform of the auto-correlation function:

$$\begin{equation} \tilde{C}(s) =\int _0^\infty {\rm e}^{-st} \langle x(t)x(0)\rangle \,{\rm d}t. \end{equation} \tag{ 11 }$$

In terms of $\tilde{P}(x,s|x^{\prime },0)$, we thus have

$$\begin{equation} \tilde{C}(s)=\int _{-\infty }^\infty {\rm d}x^{\prime } \int _{-\infty }^\infty {\rm d}x\, x\, x^{\prime } \tilde{P}(x,s|x^{\prime },0) \rho _f(x^{\prime }), \end{equation} \tag{ 12 }$$

where ρ_f(x) is the stationary density of the Fokker–Planck equation, given by

$$\begin{equation} \rho _f(x) = \frac{1}{Z} \left\lbrace \begin{array}{@{}lcl@{}}{\rm e}^{-(x+\delta -f)^2/2}\, {\rm e}^{(\delta -f)^2/2} & \mbox{ if } & x>0 \\ {\rm e}^{-(x-\delta -f)^2/2)}\, {\rm e}^{(\delta +f)^2/2} & \mbox{ if } & x<0 \end{array} \right. \end{equation} \tag{ 13 }$$

with normalization

$$\begin{equation} Z=\frac{D_{-1} (\delta +f)}{D_0(\delta +f)}+\frac{D_{-1} (\delta -f)}{D_0(\delta -f)}. \end{equation} \tag{ 14 }$$

This expression of ρ_f(x) is invariant with respect to the inversions x → −x and f → −f, as is the propagator. As a result, we can rewrite $\tilde{C}(s)$ as

$$\begin{equation} \tilde{C}(s)= \int _0^\infty {\rm d}x^{\prime } \int _{-\infty }^\infty {\rm d}x ( x\, x^{\prime } \tilde{P}(x,s|x^{\prime },0) \rho _f(x^{\prime }) + `f \rightarrow -f{\hbox{'}}), \end{equation} \tag{ 15 }$$

where the symbol 'f → −f' indicates the contribution to the kernel obtained by replacing f by −f.

We proceed to evaluate the integrals in equation (15). The first one in x can be performed using a known relation for the derivative of parabolic cylinder functions, which can be obtained from the relation shown in equation (A.17) of appendix A. The result is

$$\begin{eqnarray} &&\fl \int _0^\infty {\rm d}x^{\prime } \int _{-\infty }^\infty {\rm d}x \, x\, x^{\prime }\, \tilde{P}(x,s|x^{\prime },0) \rho _f(x^{\prime }) = \int _0^\infty {\rm d}x^{\prime } \,\left(\frac{(x^{\prime })^2}{s} - \frac{x^{\prime }(x^{\prime }+\delta -f)}{s(s+1)}\right) \rho _f(x^{\prime }) \nonumber \\ &&+\, \frac{D_{-s-2}(\delta -f) \Gamma (s)}{Z \sqrt{2 \pi }} \bigg ( -\frac{{\rm e}^{-(\delta +f)^2/4}}{{{\rm e}^{-(\delta -f)^2/4}}} g_<(s,\delta ,f) D_{-s-2}(\delta +f) \nonumber \\ &&+\, g_>(s,\delta ,f) D_{-s-2}(\delta -f) + D_{-s-2}(-\delta +f)\bigg ). \end{eqnarray} \tag{ 16 }$$

The second term on the right-hand side of this expression can be simplified using the definitions (A.21) and (A.23) as well as the identities (A.17) and (A.19):

$$\begin{eqnarray} &&\fl D_{-s-2}(\delta -f)\Bigg ( -\frac{{\rm e}^{-(\delta +f)^2/4}}{{{\rm e}^{-(\delta -f)^2/4}}} g_<(s,\delta ,f) D_{-s-2}(\delta +f) \nonumber \\ &&\qquad + g_>(s,\delta ,f) D_{-s-2}(\delta -f)+ D_{-s-2}(-\delta +f)\Bigg )\nonumber \\ &&\fl \qquad = \frac{2 \delta \sqrt{2 \pi }}{s(s+1) \Gamma (s)} \frac{D_{-s-1}(\delta +f)}{D_{-s}(\delta -f) D_{-s-1}(\delta +f) + D_{-s}(\delta +f) D_{-s-1}(\delta -f)}. \end{eqnarray} \tag{ 17 }$$

By combining this result in equation (15), and by performing the integral over x', we then obtain

$$\begin{eqnarray} &&\fl \tilde{C}(s) = \frac{1}{s} \langle x^2\rangle _f - \frac{1}{s(s+1)} ( \langle x^2\rangle _f + \delta \langle | x|\rangle _f - f \langle x \rangle _f ) \nonumber \\ &&+ \frac{2\delta }{s(s+1)} \frac{1}{Z} \frac{D_{-s-1}(\delta +f) D_{-s-2}(\delta -f) + D_{-s-1}(\delta -f) D_{-s-2}(\delta +f)}{D_{-s}(\delta +f) D_{-s-1}(\delta -f) + D_{-s}(\delta -f) D_{-s-1}(\delta +f)}, \nonumber\\ \end{eqnarray} \tag{ 18 }$$

where 〈 · 〉_f denotes the expected value with respect to the stationary distribution ρ_f. Surprisingly, all of the stationary expected values appearing above can be conveniently written in terms of parabolic cylinder functions:

$$\begin{eqnarray} \langle x^2 \rangle _f &=& 2 \frac{D_0(\delta +f) D_{-3}(\delta -f) + D_0(\delta -f) D_{-3}(\delta +f)}{D_0(\delta +f) D_{-1}(\delta -f) + D_0(\delta -f) D_{-1}(\delta +f)} \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \langle |x| \rangle _f &=& \frac{D_0(\delta +f) D_{-2}(\delta -f) + D_0(\delta -f) D_{-2}(\delta +f)}{D_0(\delta +f) D_{-1}(\delta -f) + D_0(\delta -f) D_{-1}(\delta +f)} \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} \langle x \rangle _f &=& \frac{D_0(\delta +f) D_{-2}(\delta -f) - D_0(\delta -f) D_{-2}(\delta +f)}{D_0(\delta +f) D_{-1}(\delta -f) + D_0(\delta -f) D_{-1}(\delta +f)} . \end{eqnarray} \tag{ 21 }$$

The expression shown in equation (18) can be re-written in other ways to make some of its properties more explicit. Using the identity (A.17), it is easy to check, in particular, that the stationary expectation values obey the relation

$$\begin{equation} \langle x^2\rangle _f +\delta \langle |x|\rangle _f-f \langle x \rangle _f=1, \end{equation} \tag{ 22 }$$

so that

$$\begin{eqnarray} &&\fl \tilde{C}(s) = \frac{1}{s}\left(\langle x^2\rangle _f -1 +\frac{2 \delta }{Z} \frac{D_{-s-1}(\delta +f) D_{-s-2}(\delta -f) + D_{-s-1}(\delta -f) D_{-s-2}(\delta +f)}{D_{-s}(\delta +f) D_{-s-1}(\delta -f) + D_{-s}(\delta -f) D_{-s-1}(\delta +f)} \right) \nonumber \\ &&\fl \qquad + \frac{1}{s+1}\left( 1-\frac{2 \delta }{Z} \frac{D_{-s-1}(\delta +f) D_{-s-2}(\delta -f) + D_{-s-1}(\delta -f) D_{-s-2}(\delta +f)}{D_{-s}(\delta +f) D_{-s-1}(\delta -f) + D_{-s}(\delta -f) D_{-s-1}(\delta +f)} \right). \end{eqnarray} \tag{ 23 }$$

This form of the Laplace transform of the correlation function, or resolvent, is quite useful to uncover its analytical structure. Because of equation (14) and the identity (A.17) (for ν = 0), the coefficient of 1/(s + 1) vanishes at s = −1 for δ ≠ 0. Thus, the apparent singularity at s = −1 is removable. For a similar reason, the coefficient of 1/s at s = 0, which corresponds to the residue at s = 0, is given by 〈x〉²_f. Thus, the long-time limit of the correlation function is given, as expected, by the square of the mean velocity. All the other poles of the resolvent $\tilde{C}(s)$ are determined by the denominator of the expression above, i.e. by the zeros of

$$\begin{equation} D_{-s}(\delta +f) D_{-s-1}(\delta -f) + D_{-s}(\delta -f) D_{-s-1}(\delta +f), \end{equation} \tag{ 24 }$$

which is precisely the characteristic equation derived from the Fokker–Planck operator [3].

Another form of $\tilde{C}(s)$ can be obtained by rearranging equation (23) using the aforementioned identities to obtain

$$\begin{eqnarray} &&\fl \tilde{C}(s) = \frac{\langle x^2 \rangle _f-1}{s} + \frac{1}{s+1} + \frac{2 \delta /Z}{s(s+1)^2} -\frac{4 \delta ^2/Z}{s(s+1)^2} \left( \frac{D_{-s}(\delta +f)}{D_{-s-1}(\delta +f)} + \frac{D_{-s}(\delta -f)}{D_{-s-1}(\delta -f)} \right)^{-1}. \nonumber \\ \end{eqnarray} \tag{ 25 }$$

This expression is better suited for numerical calculations than either equation (21) or (23). The limit δ = f = 0 is also clearer at the level of this expression. Noting from equation (22) that 〈x²〉_f = 1 in this limit, we recover $\tilde{C}(s)=1/(s+1)$, which characterizes the simple exponential decay of the correlation function of the Ornstein–Uhlenbeck process. Finally, we obtain from equation (25) a simple analytic expression for the power spectrum $p(\omega )=\mbox{Re}\, \tilde{C}(s={\rm i}\omega )$, i.e. for the Fourier transform of the auto-correlation function, namely,

$$\begin{eqnarray} \fl p(\omega ) = \frac{1}{1+\omega ^2} &-& \frac{4 \delta /Z}{(1+\omega ^2)^2} \nonumber \\ &-&\frac{4 \delta ^2/Z}{\omega } \mbox{Im} \left(\frac{1}{(1+{\rm i}\omega )^2} \left( \frac{D_{-{\rm i} \omega }(\delta +f)}{D_{-{\rm i}\omega -1}(\delta +f)} + \frac{D_{-{\rm i} \omega }(\delta -f)}{D_{-{\rm i}\omega -1}(\delta -f)} \right)^{-1} \right). \end{eqnarray} \tag{ 26 }$$

This result is used in the next section to discuss the stick–slip transition occurring at f = δ.

We have discussed in detail the behavior of the correlation function 〈x(t)x(0)〉 as a function of δ and f in [3] and, more precisely, how the stick–slip transition that appears in the deterministic (Γ = 0) equation when f = δ is modified in the presence of noise (Γ > 0) to a smooth crossover between a stick and a slip regimes, characterized by different exponential decay of 〈x(t)x(0)〉. Figure 2 shows how this crossover shows up at the level of the power spectrum. We see that p(ω) is rather flat in the stick regime (f < δ), and that it starts to develop a sharp zero-frequency peak at the stick–slip transition f = δ. This peak is the translation in frequency of the increase of the correlation time associated with 〈x(t)x(0)〉 as we go from the slip to the stick regimes [3].

Zoom In Zoom Out Reset image size

To gain further insight into the behavior of p(ω), we plot this function on a log–linear scale in figure 3 together with an asymptotic expansion of this function obtained in the limit of large dry friction δ and large forcing f. Mathematically, this expansion is equivalent to the small noise limit of the Langevin equation, and is obtained by using the representation (23) for $\tilde{C}(s)$ to express the power spectrum as

$$\begin{eqnarray} &&\fl p(\omega ) = \frac{1}{\omega } \mbox{Im}\Bigg (\frac{1}{1+{\rm i} \omega }\nonumber \\ &&\fl \quad \times \left ( \frac{2 \delta }{F_0(\delta +f)+F_0(\delta -f)}\frac{F_{1+{\rm i}\omega }(\delta +f) +F_{1+{\rm i}\omega }(\delta -f)}{ 2\delta +(1+{\rm i}\omega )(F_{1+{\rm i}\omega }(\delta +f)+F_{1+{\rm i}\omega }(\delta -f))}-1 \right )\Bigg ), \end{eqnarray} \tag{ 27 }$$

where

$$\begin{equation} F_s(\delta )=\frac{D_{-s-1}(\delta )}{D_{-s}(\delta )}. \end{equation} \tag{ 28 }$$

With some results of asymptotic analysis, presented in appendix B, we can then obtain the following approximation, which is valid for large values of the parameters δ and f:

$$\begin{eqnarray} \fl\displaystyle p(\omega ) \approx \left\lbrace \begin{array}{@{}ll@{}}\dsty\frac{1}{\omega } \mbox{Im} \left(\dsty\frac{\delta ^2-f^2}{(1+{\rm i} \omega )^2} \frac{\sqrt{\frac{(\delta +f)^2}{4}+1 +{\rm i}\omega }+\sqrt{\frac{(\delta -f)^2}{4}+1 +{\rm i}\omega }-\delta }{\sqrt{\frac{(\delta +f)^2}{4}+1 +{\rm i}\omega } +\sqrt{\frac{(\delta -f)^2}{4}+1 +{\rm i}\omega }+\delta } -\frac{1}{1+{\rm i} \omega }\! \right) & 0\le f <\delta \\ 1/(1+\omega ^2) \qquad & 0 \le \delta < f; \end{array} \right. \end{eqnarray} \tag{ 29 }$$

(see [11] for related results derived via continued fraction expansions). Figure 3 shows that this asymptotic formula is relatively accurate, even for rather small parameter values. Larger deviations are visible close to the transition point δ = f. The inset of figure 3 also shows that the tail of p(ω) decays at large frequencies as ω⁻², which is the sign that temporal correlations decay exponentially, as in the Ornstein–Uhlenbeck process.

Zoom In Zoom Out Reset image size

To give an idea of how these results might compare in practice with experimental results, we show in figure 4 simulated paths of the Brownian motion equation with dry friction for the rescaled variable x(t) together with their corresponding power spectrum. Two paths are shown on the left-hand side of figure 4: one in the stick regime (δ = 12 and f = 10) and one in the slip regime (δ = 12 and f = 14). The power spectrum characterizing each of these regimes, shown on the right-hand side of figure 4, is computed numerically by averaging the spectra of many random paths (here 50) over a relatively long time (here T = 100). The numerical results compare well with the theory, as can be seen. Instead of averaging different spectra, one can also use, as is well known, a frequency window larger than the frequency spacing Δω = 2π/T to obtain a relatively smooth spectrum.

Zoom In Zoom Out Reset image size

To close this section, let us now briefly discuss the behavior of two quantities derived from the power spectrum. The first is the stationary variance Δx², which is proportional to the total spectral weight:

$$\begin{equation} \Delta x^2=\langle x^2\rangle _f-\langle x \rangle _f^2=\frac{1}{\pi }\int _{-\infty }^\infty p(\omega ) \,{\rm d} \omega . \end{equation} \tag{ 30 }$$

As shown in figure 5, there is a sharp increase of Δx² at the stick–slip transition f = δ, separating a low variance (stick) regime from a high variance (slip) regime. The second quantity is the (non-dimensional) diffusion constant:

$$\begin{equation} D=\int _0^\infty \langle x(t)x(0)\rangle \,{\rm d}t = \lim _{\omega \rightarrow 0} p(\omega ). \end{equation} \tag{ 31 }$$

This quantity is of particular interest, since it can be measured in experiments. Its expression is obtained from equations (27) or (29) and is plotted in figure 5. We see that in the slip regime, the diffusion constant does not depend much on the external force, as is the case for the Ornstein–Uhlenbeck process. This is consistent with the observation that the slip regime is essentially a regime of normal Brownian motion in which dry friction force plays little role; see [3] for more details. In the stick regime, on the other hand, the diffusion constant is relatively small, and sharply increases when the stick–slip transition is approached. This behavior is also seen if we derive D from the asymptotic expression (29). This is illustrated with the dashed lines in figure 5.

Zoom In Zoom Out Reset image size

We have presented the exact solution of the propagator of a Langevin equation modeling Brownian motion in the presence of solid friction, viscous damping, and an external constant force. The main feature of this equation is that it shows a stick–slip transition often encountered in real systems involving solid friction and external forcing. The solution of this equation follows by considering the Laplace transform of its associated Fokker–Planck equation, and serves as a complement to a previous exact solution, derived in [3] using the spectral decomposition of the Fokker–Planck operator. It also extends a similar Laplace solution, previously derived in [5–7] for the special case where only solid friction is present.

A clear advantage of the Laplace solution over the spectral solution is that the former is given as an explicit and compact formula, which can be used to obtain the propagator by inverse Laplace transform. The Laplace solution is also useful as it enables us to obtain the power spectrum of the system, in addition to the diffusion constant of the Brownian motion affected by solid friction. Both of these characteristics are easily accessible experimentally, and so might be useful to compare the model with experiments involving noise and solid friction, such as those recently reported for example in [12, 13].

HT is grateful for the support and hospitality of the National Institute of Theoretical Physics at the University of Stellenbosch, South Africa, where part of this work was written. The work of WJ is partly supported by EPSRC through grant no. EP/H04812X/1.

We provide here some details of the derivation of the expressions (7)–(9) for the Laplace transform $\tilde{P}(x,s|x^{\prime },0)$ of the propagator P(x, t|x', 0).

The derivation starts with the Fokker–Planck equation (5), which for positive values for the initial condition x' > 0 reads

$$\begin{equation} \fl s \tilde{P}(x,s|x^{\prime },0) = \frac{{\rm d}}{{\rm d} x} (x-\delta -f) \tilde{P}(x,s|x^{\prime },0)+ \frac{{\rm d}^2 \tilde{P}}{{\rm d} x^2}, \quad x<0 \end{equation} \tag{ A.1 }$$

$$\begin{equation} \fl s \tilde{P}(x,s|x^{\prime },0) = \frac{{\rm d}}{{\rm d} x} (x+\delta -f) \tilde{P}(x,s|x^{\prime },0)+ \frac{{\rm d}^2 \tilde{P}}{{\rm d} x^2},\quad 0<x<x^{\prime } \end{equation} \tag{ A.2 }$$

$$\begin{equation} \fl s \tilde{P}(x,s|x^{\prime },0) = \frac{{\rm d}}{{\rm d} x} (x+\delta -f) \tilde{P}(x,s|x^{\prime },0)+ \frac{{\rm d}^2 \tilde{P}}{{\rm d} x^2}, \quad x^{\prime }<x. \end{equation} \tag{ A.3 }$$

These equations must be solved with the following matching conditions at x = 0:

$$\begin{equation} \fl \tilde{P}(x=0^-,s|x^{\prime },0) = \tilde{P}(x=0^+,s|x^{\prime },0) \end{equation} \tag{ A.4 }$$

$$\begin{equation} \fl (-\delta -f) \tilde{P}(x=0^-,s|x^{\prime },0) + \left.\frac{{\rm d} \tilde{P}}{{\rm d}x}\right|_{x=0^-} = (\delta -f) \tilde{P}(x=0^+,s|x^{\prime },0) + \left.\frac{{\rm d} \tilde{P}}{{\rm d}x}\right|_{x=0^+} \end{equation} \tag{ A.5 }$$

$$\begin{equation} \fl \tilde{P}(x=x^{\prime }-0,s|x^{\prime },0) = \tilde{P}(x=x^{\prime }+0,s|x^{\prime },0) \end{equation} \tag{ A.6 }$$

$$\begin{equation} \fl -1 = \left.\frac{{\rm d} \tilde{P}}{{\rm d}x}\right|_{x=x^{\prime }+0}-\left.\frac{{\rm d} \tilde{P}}{{\rm d}x}\right|_{x=x^{\prime }-0} , \end{equation} \tag{ A.7 }$$

which results from the continuity of the probability current, as well as the usual decaying boundary conditions at x = ±∞.

The equations (A.1)–(A.3) have the form of the Hermite differential equation

$$\begin{equation} u^{\prime \prime }(z)+(zu(z))^{\prime }+\nu u(z)=0. \end{equation} \tag{ A.8 }$$

The solution can be written in terms of parabolic cylinder functions either as ${\rm e}^{-z^2/4} D_\nu (z)$ or ${\rm e}^{-z^2/4} D_\nu (-z)$.⁵ In view of the asymptotic property of the parabolic cylinder function,

$$\begin{equation} D_\nu (z) \sim z^\nu {\rm e}^{-z^2/4},\quad z\rightarrow \infty , \end{equation} \tag{ A.9 }$$

the solution of equations (A.1) and (A.3) having vanishing currents at infinity can be written as

$$\begin{equation} \tilde{P}(x,s|x^{\prime },0) = C_-(x^{\prime },s) \,{\rm e}^{-(x-\delta -f)^2/4} D_{-s}(-x+\delta +f) \end{equation} \tag{ A.10 }$$

for x < 0 and

$$\begin{equation} \tilde{P}(x,s|x^{\prime },0) =C_+(x^{\prime },s) \,{\rm e}^{-(x+\delta -f)^2/4} D_{-s}(x+\delta -f) \end{equation} \tag{ A.11 }$$

for x > x'. On the other hand, the solution of equation (A.2) is given by a linear combination of the two fundamental solutions:

$$\begin{eqnarray} &&\fl \tilde{P}(x,s|x^{\prime },0) = B_+(x^{\prime },s)\, {\rm e}^{-(x+\delta -f)^2/4}\, D_{-s}(x+\delta -f) \nonumber \\ &&+ B_{-}(s,x^{\prime })\, {\rm e}^{-(x+\delta -f)^2/4}\, D_{-s}(-x-\delta +f), \qquad 0<x<x^{\prime }. \end{eqnarray} \tag{ A.12 }$$

The amplitudes B_± and C_± are determined by the matching conditions (A.4)–(A.7), which result in a set of inhomogeneous linear equations:

$$\begin{eqnarray} &&\fl C_-(x^{\prime },s)\, {\rm e}^{-(\delta +f)^2/4}\, D_{-s}(\delta +f) = B_+(x^{\prime },s)\, {\rm e}^{-(\delta -f)^2/4}\, D_{-s}(\delta -f) \nonumber \\ &&+ B_{-}(s,x^{\prime })\, {\rm e}^{-(\delta -f)^2/4}\, D_{-s}(-\delta +f) \end{eqnarray} \tag{ A.13 }$$

$$\begin{eqnarray} &&\fl -C_-(x^{\prime },s)\, {\rm e}^{-(\delta +f)^2/4}\, D_{-s-1}(\delta +f) = B_+(x^{\prime },s)\, {\rm e}^{-(\delta -f)^2/4}\, D_{-s-1}(\delta -f) \nonumber \\ &&- B_-(x^{\prime },s)\, {\rm e}^{-(\delta -f)^2/4}\, D_{-s-1}(-\delta +f) \end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray} &&\fl C_+(x^{\prime },s) D_{-s}(x^{\prime }+\delta -f) = B_+(x^{\prime },s) D_{-s}(x^{\prime }+\delta -f) \nonumber \\ &&+ B_-(x^{\prime },s) D_{-s}(-x^{\prime }-\delta +f) \end{eqnarray} \tag{ A.15 }$$

$$\begin{eqnarray} &&\fl {\rm e}^{(x^{\prime }+\delta -f)^2/4} = C_+(x^{\prime },s) D_{-s+1}(x^{\prime }+\delta -f)\nonumber \\ &&- B_+(x^{\prime },s) D_{-s+1}(x^{\prime }+\delta -f) \nonumber \\ &&+ B_-(x^{\prime },s) D_{-s+1}(-x^{\prime }-\delta +f). \end{eqnarray} \tag{ A.16 }$$

In writing these equations, we have used the following identities for the parabolic cylinder functions:

$$\begin{eqnarray} &&\nu \, {\rm e}^{-z^2/4}\, D_{\nu -1}(z) = {\rm e}^{-z^2/4} (z D_\nu (z)-D_{\nu +1}(z)) \nonumber \\ &&{\hphantom{\fl \nu \, {\rm e}^{-z^2/4}\, D_{\nu -1}(z)}}= z\, {\rm e}^{-z^2/4}\, D_\nu (z) +({\rm e}^{-z^2/4}\, D_\nu (z))^{\prime } \end{eqnarray} \tag{ A.17 }$$

to evaluate and simplify the derivatives.

Equations (A.13)–(A.16) can be solved directly to find B_± and C_±. The difference of equations (A.15) and (A.16) yields

$$\begin{equation} B_-(x^{\prime },s) = \frac{\Gamma (s)}{2 \sqrt{\pi }} {\rm e}^{(x^{\prime }+\delta -f)^2/4} D_{-s}(x^{\prime }+\delta -f), \end{equation} \tag{ A.18 }$$

if we take into account the product identity

$$\begin{equation} D_\nu (z) D_{\nu -1}(-z)+D_\nu (-z) D_{\nu -1}(z)=\frac{\sqrt{2\pi }}{\Gamma (-\nu +1)}, \end{equation} \tag{ A.19 }$$

which follows from the Wronskian of the fundamental system.⁶ Then the difference of equations (A.13) and (A.14) yields

$$\begin{equation} C_-(x^{\prime },s)=g_<(s,\delta ,f) B_-(x^{\prime },s), \end{equation} \tag{ A.20 }$$

where

$$\begin{equation} g_<(s,\delta ,f) = \frac{ {\rm e}^{-(\delta -f)^2/4} \left(\frac{D_{-s}(-\delta +f)}{D_{-s}(\delta -f)}+ \frac{D_{-s-1}(-\delta +f)}{D_{-s-1}(\delta -f)}\right)} { {\rm e}^{-(\delta +f)^2/4} \left(\frac{D_{-s}(\delta +f)}{D_{-s}(\delta -f)}+ \frac{D_{-s-1}(\delta +f)}{D_{-s-1}(\delta -f)} \right)}. \end{equation} \tag{ A.21 }$$

Summing equations (A.13) and (A.14), we also find

$$\begin{equation} B_+(x^{\prime },s)=g_>(s,\delta ,f) B_-(x^{\prime },s) \end{equation} \tag{ A.22 }$$

with

$$\begin{equation} g_>(s,\delta ,f)=\frac{ \frac{D_{-s-1}(-\delta +f)}{D_{-s-1}(\delta +f)} -\frac{D_{-s}(-\delta +f)}{D_{-s}(\delta +f)}}{ \frac{D_{-s}(\delta -f)}{D_{-s}(\delta +f)}+ \frac{D_{-s-1}(\delta -f)}{D_{-s-1}(\delta +f)}}. \end{equation} \tag{ A.23 }$$

Finally, summing equations (A.15) and (A.16), we find

$$\begin{equation} C_+(x^{\prime },s)=B_+(x^{\prime },s)+\frac{\Gamma (s)}{\sqrt{2\pi }} {\rm e}^{(x^{\prime }+\delta -f)^2/4} D_{-s}(-x^{\prime }-\delta +f). \end{equation} \tag{ A.24 }$$

The substitution of all these expressions back into equations (A.10)–(A.12) leads us to equations (7)–(9).

The power spectrum, defined in equation (26), is entirely determined by the ratio F_s(δ) of parabolic cylinder functions defined in equation (28). Using the linear recurrence relation (A.17), this ratio is found to obey the relation

$$\begin{equation} F_s(\delta )=\frac{1}{\delta + (s+1) F_{s+1}(\delta )}, \end{equation} \tag{ B.1 }$$

iteration of which leads to continued fraction expansions [11].

Asymptotic expressions for F_s(δ) for large values of δ and s can be obtained from the contour integral representation of parabolic cylinder functions,

$$\begin{equation} D_\nu (z)=\frac{{\rm e}^{-z^2/4}}{{\rm i}\sqrt{2\pi }}\int _{c-{\rm i}\infty }^{c+{\rm i}\infty }{\rm e}^{-zt+t^2/2}t^\nu {\rm d}t, \end{equation} \tag{ B.2 }$$

which is valid for all $z,\nu \in \mathbb C$, and any integral contour satisfying $|\arg t|<\pi /2$ and c > 0. To this end, we introduce a large parameter N and change the variable of integration to $u=t/\sqrt{N}$. This leads to

$$\begin{equation} F_{Ns}(\sqrt{N}\delta )=\frac{1}{\sqrt{N}} \frac{ \int _{c-{\rm i}\infty }^{c+{\rm i}\infty }{\rm e}^{N(-\delta u+u^2/2-s\log u)}\frac{{\rm d}u}{u}}{ \int _{c-{\rm i}\infty }^{c+{\rm i}\infty }{\rm e}^{N(-\delta u+u^2/2-s\log u)}{\rm d}u}. \end{equation} \tag{ B.3 }$$

Each of the integrals above has the Laplace form

$$\begin{equation} \int _{\mathcal{C}}{\rm e}^{Nf(z)}g(z)\,{\rm d}z, \end{equation} \tag{ B.4 }$$

which can be approximated following the saddlepoint method as

$$\begin{equation} \int _{\mathcal{C}}{\rm e}^{Nf(z)}g(z){\rm d}z\sim \pm {\rm i} {\rm e}^{Nf(z_s)}g(z_s)\sqrt{\frac{2\pi }{N f^{\prime \prime }(z_s)}}, \end{equation} \tag{ B.5 }$$

where z_s is the saddlepoint of the integral given by f'(z_s) = 0, and the sign depends on the orientation of the contour C, which is deformed to a steepest-descent contour passing through the saddlepoint. This approximation is valid to leading order in N provided g does not vanish at z_s.

In our case, we obtain the saddlepoints from

$$\begin{equation} 0=\frac{\partial }{\partial u}\left(-\delta u+\frac{1}{2}u^2-s\log u\right)=-\delta +u-\frac{s}{u}, \end{equation} \tag{ B.6 }$$

which has two solutions:

$$\begin{equation} u_1=\case{1}{2}(\delta +\sqrt{\delta ^2+4s}),\quad u_2=\case{1}{2}(\delta -\sqrt{\delta ^2+4s}). \end{equation} \tag{ B.7 }$$

On the one hand, if δ² + 4s > 0, then u₁ > 0 > u₂ and the steepest-descent contour passes through u₁ only. In this case, much of the saddlepoint approximation cancels, leaving us with

$$\begin{equation} F_{Ns}(\sqrt{N}\delta )\sim \frac{1}{\sqrt{N}}\frac{1}{u_1}=\frac{2}{\sqrt{N}\delta +\sqrt{N\delta ^2+4Ns}}. \end{equation} \tag{ B.8 }$$

On the other hand, if δ² + 4s < 0, then u₁ and u₂ are complex conjugate to each other, and the steepest-descent contour passes through both saddlepoints. The crossover between these cases is given by coalescing saddlepoints; a uniform asymptotic expansion for this case can be obtained in terms of Airy functions [15].

A second case of interest is given when the argument of F is kept constant, i.e.

$$\begin{equation} F_s(\sqrt{N}\delta )=\frac{1}{\sqrt{N}}\frac{ \int _{c-{\rm i}\infty }^{c+{\rm i}\infty }{\rm e}^{N(-\delta u+u^2/2)}\frac{{\rm d}u}{u^{s+1}}}{ \int _{c-{\rm i}\infty }^{c+{\rm i}\infty }{\rm e}^{N(-\delta u+u^2/2)}\frac{{\rm d}u}{u^s}}. \end{equation} \tag{ B.9 }$$

In this case, there is only one saddlepoint at u₀ = δ, which collides with a singularity of the integrand when δ approaches zero. For δ > 0, we then find

$$\begin{equation} F_s(\sqrt{N}\delta )\sim \frac{1}{\sqrt{N} u_0}=\frac{1}{\sqrt{N}\delta }, \end{equation} \tag{ B.10 }$$

which matches the previous asymptotic evaluation. For δ < 0, i.e. for Kramer's problem, the integral is dominated by the singularity at the origin, and its contribution gives

$$\begin{equation} F_s(\sqrt{N}\delta )\sim \frac{-\sqrt{N}\delta }{s} \end{equation} \tag{ B.11 }$$

for s ≠ 0, −1, −2, .... If s = −1, −2, ..., we find instead $F_s(\sqrt{N}\delta )\sim 1/(\sqrt{N}\delta )$.