- #1
ergospherical
- 967
- 1,276
I wrote these "notes" up a couple of weeks ago when I was trying to see what I remembered and put a few different concepts in one place. There's a lot of material available for stochastic physics but there are very varied approaches between different sources (random, even).
Maybe somebody will find the content useful. It's not necessarily a topic which is covered in much depth in many courses. It's not rigorous, either. (Hence "physics-focussed"...)
Fokker-Planck
We'll look at a 1-dimensional case of a particle moving in a force-field ##f(x)## and subject to friction ##-\gamma \dot{x}##.$$\ddot{x} = - \gamma \dot{x} +f(x) + N(t)$$The term ##N(t)## is Gaussian white-noise, arising due to microscopic (thermal) effects. In the over-damped limit where ##|\ddot{x}| \ll \gamma |\dot{x}|##, this is simplified to a first order differential equation.$$\gamma \dot{x} = f(x) + N(t)$$You can re-write this by identifying the stochastic part ##N(t) dt## as proportional to increments ##dW## of a random Wiener process (in short, ##dW## follows a normal distribution with mean ##0## and variance ##dt##).$$dx = \frac{1}{\gamma} f(x) dt + \sqrt{\frac{2kT}{\gamma}} dW$$
The coefficient ##\sqrt{2kT/\gamma} \equiv \sigma## of ##dW## is the standard deviation of the stochastic noise. The conversion of the stochastic differential equation above into the Fokker-Planck equation for the full probability density ##P(x,t)## is a textbook derivation (see `Kramers-Moyal expansion').$$\frac{\partial P}{\partial t} = - \frac{\partial}{\partial x} \left( \frac{f(x)}{\gamma} P\right) + \frac{\sigma^2}{2} \frac{\partial^2 P}{\partial x^2}$$If ##f(x) = f_0## is a constant, say, then the solution is a Gaussien with mean position ##\bar{x}(t) = (f_0/\gamma) t## and constant variance ##\sigma^2##,$$P(x,t) = \frac{1}{\sqrt{2\pi \sigma^2 t}} \exp \left(- \frac{(x-\bar{x}(t))^2}{2\sigma^2 t}\right)$$
Harmonic well & Ornstein-Uhlenbeck
It is more complicated in something like a harmonic force ##f(x) = -k(x-x_0)##. In fact the form of the SDE is then in the form of an Ornstein-Uhlenbeck process$$dx = -\frac{k}{\gamma} (x-x_0) dt + \sigma dW$$We can get ##P(x,t)## with a trick. Consider a new stochastic variable ##y \equiv e^{kt} x##. With the help of the Itô Lemma, $$dy = (\partial_t y)dt + (\partial_x y) dx + \frac{1}{2} (\partial_x^2 y) dx^2$$ The ##dx^2## is included in the sense that ##dW^2 = dt## (definition of Wiener process), and therefore ##dx^2 = \sigma^2 dt## (ignoring ##dt## and ##dt dW## which are smaller order than ##dt##). In any case it makes no difference because in fact ##\partial_x^2 y = 0##.$$dy = kx_0 e^{kt} dt + e^{kt} \sigma dW$$So then by integration with the initial condition ##x(0) = 0##, and converting back in terms of ##x##,$$x = x_0(1-e^{-kt}) + \sigma \int_0^t e^{k(t'-t)} dW_{t'}$$The mean is ##\bar{x}(t) = x_0(1-e^{-kt})##, because ##\overline{dW} = 0##. For the variance check that you agree with$$\mathrm{Var}(x) = \overline{x^2} - \bar{x}^2 = \overline{\left( \sigma \int_0^t e^{k(t'-t)} dW_{t'} \right)^2}$$To evaluate the RHS, you can write it as a double integral$$\mathrm{Var}(x) = \sigma^2 \int_0^t \int_0^t e^{k(t'-t)} e^{k(t''-t)} \overline{dW_{t'} dW_{t''}}$$with ##\overline{dW_{t'} dW_{t''}} = \delta(t'-t'')## to get$$\mathrm{Var}(x) = \frac{\sigma^2}{2k}(1-e^{-2kt})$$The fundamental Gaussian solution for ##P(x,t)## is $$P(x,t) = \frac{1}{\sqrt{\pi \frac{\sigma^2}{k}(1-e^{-2kt}) t}} \exp \left(- \frac{(x-\bar{x}(t))^2}{\frac{\sigma^2}{k}(1-e^{-2kt}) t}\right)$$with ##\bar{x}(t) = x_0(1-e^{-kt})## as before.
First passage times
So you may have the probability ##P(x,t)##. It's of interest in a lot of applications to figure out when a particle is likely to escape from an enclosure ##\mathcal{D}## given that there are absorbing sinks/traps at the boundaries ##\partial \mathcal{D}##. Specifically, imagine an enclosure with a wall on the left at ##x_L## and a sink/trap on the right at ##x_{\mathrm{abs}}##. The particle is going to begin at ##x_0 \in \mathcal{D} - \partial \mathcal{D}##, somewhere between.
There are two ways to proceed. The probability that the particle remains in the enclosure at any given time is obtained by just integrating the probability over the region,$$R(t) = \int_{\mathcal{D}} P(x,t) dx$$Also, if ##\mathfrak{p}## is the probability density for the first passage time then$$R(t) = P(\text{first passage $> t$}) = \int_t^{\infty} \mathfrak{p}(t') dt'$$
You can see that ##\mathfrak{p} = - \partial_t R##. The expected time taken for the first passage is
\begin{align*}
\bar{\mathfrak{t}} &= \int_0^{\infty} t \mathfrak{p}(t) dt = - \int_0^{\infty} t (\partial_{t} R) dt = \int_0^{\infty} R(t) dt \\
&= \int_0^{\infty} \int_{\mathcal{D}} P(x,t) dx dt
\end{align*}
But this is not usually the easiest way, because you need to solve Fokker-Planck for the full ##P(x,t)##. Proceed another way. For an SDE of the form ##dx = \mu(x) dt + \sigma dW##, look at the operator ##L## defined b$$L\phi = \partial_x( \mu(x) \phi) - \frac{1}{2} \sigma^2 \partial_x^2 \phi$$In terms of ##L## the Fokker-Planck equation is ##\partial_t P = -LP##. Integration by parts gives it's adjoint (##\int \phi L^{\dagger} \varphi = \int \varphi L \phi##) as $$L^{\dagger}\phi = -\partial_x( \mu(x) \phi) - \frac{1}{2} \sigma^2 \partial_x^2 \phi$$You can start from Chapman-Kolmogorov. For brevity of the formulas we can just say `0' instead of writing ##(x_0,t_0)##.$$P(x_1,t_1|0) = \int P(x_1, t_1|x,t) P(x,t|0)$$
Differentiating with respect to ##t##
\begin{align*}
0 &= \int \left[ \partial_t P(x_1, t_1|x,t) P(x,t|0) + P(x_1,t_1|x,t) \partial_t P(x,t|0) \right] \ dx \\
&= \int \left[ \partial_t P(x_1, t_1|x,t) P(x,t|0) - P(x_1,t_1|x,t) \underbrace{L}_{\text{FP op.}} P(x,t|0) \right] \ dx \\
&= \int \left[ \partial_t P(x_1, t_1|x,t) - \underbrace{L^{\dagger}}_{\text{FP adj. op.}} P(x_1,t_1|x,t) \right] P(x,t|0) \ dx
\end{align*}
Therefore the time derivative with respect to the prior argument is$$\partial_t P(x_1,t_1|x,t) = L^{\dagger} P(x_1,t_1|x,t)$$
Operating with this on the mean passage time ##\bar{\mathfrak{t}}##, you can show that
\begin{align*}
L^{\dagger} \bar{\mathfrak{t}} = \int_0^{\infty} \int_{\mathcal{D}} \frac{\partial}{\partial t_0} P(x,t| x_0, t_0) \ dx dt &= \int_0^{\infty} \int_{\mathcal{D}} \frac{\partial}{\partial t_0} P(x,t-t_0| x_0, 0) \ dx dt \\
&= - \int_0^{\infty} \int_{\mathcal{D}} \frac{\partial}{\partial t} P(x,t-t_0| x_0, 0) \ dx dt \\
&= -\int_0^{\infty} \partial_t R \ dt \\
&= \int_0^{\infty} \mathfrak{p}(t) \ dt \\
&= 1
\end{align*}
In the top-line, time-shift invariance is needed to replace ##t \mapsto t - t_0##. The key point is that you get $$L^{\dagger} \bar{\mathfrak{t}} = -\partial_x(\mu(x)\bar{\mathfrak{t}}) - \frac{1}{2} \sigma^2 \partial_x^2 \bar{\mathfrak{t}} = 1$$This is a pretty helpful form, e.g. if ##\mu(x) = f(x)/\gamma## where ##f(x) = -U'(x)## is a potential force. Because in that case, see if you can solve this for$$\bar{\mathfrak{t}} = \frac{2}{\sigma^2} \int_{x_0}^{x_{\mathrm{abs}}} e^{2 U(x)/\sigma^2} \left[ \int_{x_L}^{x} e^{-2 U(x')/\sigma^2} dx' \right] dx$$
In the limits some care is needed. ##x_0## is the initial position of the particle and ##x_{\mathrm{abs}}## is the position of the absorbing trap or sink.
You can play around with this for some examples. For free diffusion with ##U=0##, the exponentials go over to ##1## and $$\bar{\mathfrak{t}} = \frac{1}{\sigma^2}\left[ (x_{\mathrm{abs}}-x_L)^2 - (x_0 - x_L)^2\right]$$
Kramers escape
A short final application. You have a potential with a local minimum at ##x_a = 0##, with value ##U(a) = 0## and a local maximum at ##x_b > x_a##, with value ##U_b##. (Assume also that ##U \rightarrow \infty## below ##x < 0##, so the particle is prevented from escaping to ##-\infty##.)
At ##x_a##, the potential goes locally like ##U(x) \sim \frac{1}{2} \omega_a^2 x^2##. At ##x_b##, the potential goes locally like ##U(x) \sim U_b -\frac{1}{2} \omega_b^2 x^2##.
If the particle starts at ##a##, i.e. ##x_a = 0##, how long will it take to escape over the barrier at ##b##? With the expression for the mean passage time, and taking care with the limits,$$\bar{\mathfrak{t}} = \frac{2}{\sigma^2} \int_0^{x_b} e^{2U(x)/\sigma^2} \left[\int_{-\infty}^x e^{-2U(x')/\sigma^2} dx' \right]dx$$The quick-and-dirty approach is to take ##e^{-2U(x')/\sigma^2}## to be dominated by the portion around the local minimum ##x \approx 0##(and vice versa for ##e^{2U(x)/\sigma^2}##, around ##x \approx x_b##). For ##e^{-2U(x')/\sigma^2}##, considering only the part of the integral in the region ##[-\epsilon, \epsilon]##,$$\int_{-\infty}^x e^{-2 U(x)/\sigma^2} \approx e^{-2 U_a/\sigma^2} \sqrt{\frac{\pi \sigma^2}{\omega_a^2}}$$For ##e^{2U(x)/\sigma^2}##, the upper limit ##L## means we only get to integrate over \textit{half} of the region around the maximum (i.e. only over ##[L-\epsilon, L]##. So we end up needing a factor of ##1/2##,$$\int_0^L e^{2 U(x)/\sigma^2} \approx \frac{1}{2}e^{2 U_b/\sigma^2} \sqrt{\frac{\pi \sigma^2}{\omega_b^2}}$$The overall passage time grows exponentially with the difference ##U_b - U_a##,$$\bar{\mathfrak{t}} = \frac{\pi \sigma^2}{2\omega_a \omega_b}e^{2(U_b-U_a)/\sigma^2}$$
Maybe somebody will find the content useful. It's not necessarily a topic which is covered in much depth in many courses. It's not rigorous, either. (Hence "physics-focussed"...)
Fokker-Planck
We'll look at a 1-dimensional case of a particle moving in a force-field ##f(x)## and subject to friction ##-\gamma \dot{x}##.$$\ddot{x} = - \gamma \dot{x} +f(x) + N(t)$$The term ##N(t)## is Gaussian white-noise, arising due to microscopic (thermal) effects. In the over-damped limit where ##|\ddot{x}| \ll \gamma |\dot{x}|##, this is simplified to a first order differential equation.$$\gamma \dot{x} = f(x) + N(t)$$You can re-write this by identifying the stochastic part ##N(t) dt## as proportional to increments ##dW## of a random Wiener process (in short, ##dW## follows a normal distribution with mean ##0## and variance ##dt##).$$dx = \frac{1}{\gamma} f(x) dt + \sqrt{\frac{2kT}{\gamma}} dW$$
The coefficient ##\sqrt{2kT/\gamma} \equiv \sigma## of ##dW## is the standard deviation of the stochastic noise. The conversion of the stochastic differential equation above into the Fokker-Planck equation for the full probability density ##P(x,t)## is a textbook derivation (see `Kramers-Moyal expansion').$$\frac{\partial P}{\partial t} = - \frac{\partial}{\partial x} \left( \frac{f(x)}{\gamma} P\right) + \frac{\sigma^2}{2} \frac{\partial^2 P}{\partial x^2}$$If ##f(x) = f_0## is a constant, say, then the solution is a Gaussien with mean position ##\bar{x}(t) = (f_0/\gamma) t## and constant variance ##\sigma^2##,$$P(x,t) = \frac{1}{\sqrt{2\pi \sigma^2 t}} \exp \left(- \frac{(x-\bar{x}(t))^2}{2\sigma^2 t}\right)$$
Harmonic well & Ornstein-Uhlenbeck
It is more complicated in something like a harmonic force ##f(x) = -k(x-x_0)##. In fact the form of the SDE is then in the form of an Ornstein-Uhlenbeck process$$dx = -\frac{k}{\gamma} (x-x_0) dt + \sigma dW$$We can get ##P(x,t)## with a trick. Consider a new stochastic variable ##y \equiv e^{kt} x##. With the help of the Itô Lemma, $$dy = (\partial_t y)dt + (\partial_x y) dx + \frac{1}{2} (\partial_x^2 y) dx^2$$ The ##dx^2## is included in the sense that ##dW^2 = dt## (definition of Wiener process), and therefore ##dx^2 = \sigma^2 dt## (ignoring ##dt## and ##dt dW## which are smaller order than ##dt##). In any case it makes no difference because in fact ##\partial_x^2 y = 0##.$$dy = kx_0 e^{kt} dt + e^{kt} \sigma dW$$So then by integration with the initial condition ##x(0) = 0##, and converting back in terms of ##x##,$$x = x_0(1-e^{-kt}) + \sigma \int_0^t e^{k(t'-t)} dW_{t'}$$The mean is ##\bar{x}(t) = x_0(1-e^{-kt})##, because ##\overline{dW} = 0##. For the variance check that you agree with$$\mathrm{Var}(x) = \overline{x^2} - \bar{x}^2 = \overline{\left( \sigma \int_0^t e^{k(t'-t)} dW_{t'} \right)^2}$$To evaluate the RHS, you can write it as a double integral$$\mathrm{Var}(x) = \sigma^2 \int_0^t \int_0^t e^{k(t'-t)} e^{k(t''-t)} \overline{dW_{t'} dW_{t''}}$$with ##\overline{dW_{t'} dW_{t''}} = \delta(t'-t'')## to get$$\mathrm{Var}(x) = \frac{\sigma^2}{2k}(1-e^{-2kt})$$The fundamental Gaussian solution for ##P(x,t)## is $$P(x,t) = \frac{1}{\sqrt{\pi \frac{\sigma^2}{k}(1-e^{-2kt}) t}} \exp \left(- \frac{(x-\bar{x}(t))^2}{\frac{\sigma^2}{k}(1-e^{-2kt}) t}\right)$$with ##\bar{x}(t) = x_0(1-e^{-kt})## as before.
First passage times
So you may have the probability ##P(x,t)##. It's of interest in a lot of applications to figure out when a particle is likely to escape from an enclosure ##\mathcal{D}## given that there are absorbing sinks/traps at the boundaries ##\partial \mathcal{D}##. Specifically, imagine an enclosure with a wall on the left at ##x_L## and a sink/trap on the right at ##x_{\mathrm{abs}}##. The particle is going to begin at ##x_0 \in \mathcal{D} - \partial \mathcal{D}##, somewhere between.
There are two ways to proceed. The probability that the particle remains in the enclosure at any given time is obtained by just integrating the probability over the region,$$R(t) = \int_{\mathcal{D}} P(x,t) dx$$Also, if ##\mathfrak{p}## is the probability density for the first passage time then$$R(t) = P(\text{first passage $> t$}) = \int_t^{\infty} \mathfrak{p}(t') dt'$$
You can see that ##\mathfrak{p} = - \partial_t R##. The expected time taken for the first passage is
\begin{align*}
\bar{\mathfrak{t}} &= \int_0^{\infty} t \mathfrak{p}(t) dt = - \int_0^{\infty} t (\partial_{t} R) dt = \int_0^{\infty} R(t) dt \\
&= \int_0^{\infty} \int_{\mathcal{D}} P(x,t) dx dt
\end{align*}
But this is not usually the easiest way, because you need to solve Fokker-Planck for the full ##P(x,t)##. Proceed another way. For an SDE of the form ##dx = \mu(x) dt + \sigma dW##, look at the operator ##L## defined b$$L\phi = \partial_x( \mu(x) \phi) - \frac{1}{2} \sigma^2 \partial_x^2 \phi$$In terms of ##L## the Fokker-Planck equation is ##\partial_t P = -LP##. Integration by parts gives it's adjoint (##\int \phi L^{\dagger} \varphi = \int \varphi L \phi##) as $$L^{\dagger}\phi = -\partial_x( \mu(x) \phi) - \frac{1}{2} \sigma^2 \partial_x^2 \phi$$You can start from Chapman-Kolmogorov. For brevity of the formulas we can just say `0' instead of writing ##(x_0,t_0)##.$$P(x_1,t_1|0) = \int P(x_1, t_1|x,t) P(x,t|0)$$
Differentiating with respect to ##t##
\begin{align*}
0 &= \int \left[ \partial_t P(x_1, t_1|x,t) P(x,t|0) + P(x_1,t_1|x,t) \partial_t P(x,t|0) \right] \ dx \\
&= \int \left[ \partial_t P(x_1, t_1|x,t) P(x,t|0) - P(x_1,t_1|x,t) \underbrace{L}_{\text{FP op.}} P(x,t|0) \right] \ dx \\
&= \int \left[ \partial_t P(x_1, t_1|x,t) - \underbrace{L^{\dagger}}_{\text{FP adj. op.}} P(x_1,t_1|x,t) \right] P(x,t|0) \ dx
\end{align*}
Therefore the time derivative with respect to the prior argument is$$\partial_t P(x_1,t_1|x,t) = L^{\dagger} P(x_1,t_1|x,t)$$
Operating with this on the mean passage time ##\bar{\mathfrak{t}}##, you can show that
\begin{align*}
L^{\dagger} \bar{\mathfrak{t}} = \int_0^{\infty} \int_{\mathcal{D}} \frac{\partial}{\partial t_0} P(x,t| x_0, t_0) \ dx dt &= \int_0^{\infty} \int_{\mathcal{D}} \frac{\partial}{\partial t_0} P(x,t-t_0| x_0, 0) \ dx dt \\
&= - \int_0^{\infty} \int_{\mathcal{D}} \frac{\partial}{\partial t} P(x,t-t_0| x_0, 0) \ dx dt \\
&= -\int_0^{\infty} \partial_t R \ dt \\
&= \int_0^{\infty} \mathfrak{p}(t) \ dt \\
&= 1
\end{align*}
In the top-line, time-shift invariance is needed to replace ##t \mapsto t - t_0##. The key point is that you get $$L^{\dagger} \bar{\mathfrak{t}} = -\partial_x(\mu(x)\bar{\mathfrak{t}}) - \frac{1}{2} \sigma^2 \partial_x^2 \bar{\mathfrak{t}} = 1$$This is a pretty helpful form, e.g. if ##\mu(x) = f(x)/\gamma## where ##f(x) = -U'(x)## is a potential force. Because in that case, see if you can solve this for$$\bar{\mathfrak{t}} = \frac{2}{\sigma^2} \int_{x_0}^{x_{\mathrm{abs}}} e^{2 U(x)/\sigma^2} \left[ \int_{x_L}^{x} e^{-2 U(x')/\sigma^2} dx' \right] dx$$
In the limits some care is needed. ##x_0## is the initial position of the particle and ##x_{\mathrm{abs}}## is the position of the absorbing trap or sink.
You can play around with this for some examples. For free diffusion with ##U=0##, the exponentials go over to ##1## and $$\bar{\mathfrak{t}} = \frac{1}{\sigma^2}\left[ (x_{\mathrm{abs}}-x_L)^2 - (x_0 - x_L)^2\right]$$
Kramers escape
A short final application. You have a potential with a local minimum at ##x_a = 0##, with value ##U(a) = 0## and a local maximum at ##x_b > x_a##, with value ##U_b##. (Assume also that ##U \rightarrow \infty## below ##x < 0##, so the particle is prevented from escaping to ##-\infty##.)
At ##x_a##, the potential goes locally like ##U(x) \sim \frac{1}{2} \omega_a^2 x^2##. At ##x_b##, the potential goes locally like ##U(x) \sim U_b -\frac{1}{2} \omega_b^2 x^2##.
If the particle starts at ##a##, i.e. ##x_a = 0##, how long will it take to escape over the barrier at ##b##? With the expression for the mean passage time, and taking care with the limits,$$\bar{\mathfrak{t}} = \frac{2}{\sigma^2} \int_0^{x_b} e^{2U(x)/\sigma^2} \left[\int_{-\infty}^x e^{-2U(x')/\sigma^2} dx' \right]dx$$The quick-and-dirty approach is to take ##e^{-2U(x')/\sigma^2}## to be dominated by the portion around the local minimum ##x \approx 0##(and vice versa for ##e^{2U(x)/\sigma^2}##, around ##x \approx x_b##). For ##e^{-2U(x')/\sigma^2}##, considering only the part of the integral in the region ##[-\epsilon, \epsilon]##,$$\int_{-\infty}^x e^{-2 U(x)/\sigma^2} \approx e^{-2 U_a/\sigma^2} \sqrt{\frac{\pi \sigma^2}{\omega_a^2}}$$For ##e^{2U(x)/\sigma^2}##, the upper limit ##L## means we only get to integrate over \textit{half} of the region around the maximum (i.e. only over ##[L-\epsilon, L]##. So we end up needing a factor of ##1/2##,$$\int_0^L e^{2 U(x)/\sigma^2} \approx \frac{1}{2}e^{2 U_b/\sigma^2} \sqrt{\frac{\pi \sigma^2}{\omega_b^2}}$$The overall passage time grows exponentially with the difference ##U_b - U_a##,$$\bar{\mathfrak{t}} = \frac{\pi \sigma^2}{2\omega_a \omega_b}e^{2(U_b-U_a)/\sigma^2}$$
Last edited: