# The Azimuth Project Supersymmetry and the Fokker-Planck equation (Rev #1)

This page is a blog article in progress, written by Blake Stacey.

### Idea

We provide a pedagogical introduction to Fokker--Planck diffusion equations with spatiotemporally varying drift and diffusion coefficients. The treatment of Fokker–Planck equations with changes of variable is reviewed, followed by the transformation of diffusion equations into Schrödinger-like form, the application of supersymmetric quantum mechanics (SUSY QM) and the use of Floquet theory. The effective diffusion coefficient is defined in terms of first-passage times, and we argue that exchanging a potential for its SUSY partner should leave this coefficient unchanged, even in the case of periodic forcing.

### Introduction: Biological Applications

We investigate solutions of the Fokker–Planck diffusion equation with spatiotemporally varying drift and diffusion coefficients,

$\partial_t p(x,t) = -\partial_x \left(F(x,t) p(x,t)\right) + \partial_x^2 \left(D(x,t) p(x,t)\right).$

This equation is of biological significance in that it can describe, for example, the time evolution of a gene pool. Given a diploid population of $N$ individuals, we let $x$ stand for the frequency of one allelic form of a gene; the diffusion coefficient captures the effect of random mating,

$D(x,t) = D(x) = \frac{x (1 - x)}{4N},$

while the force field $F(x,t)$ incorporates both mutation and selection pressure:

$F(x,t) = \mu(1-2x) + f(t)x(1-x).$

Here, selection is implemented via the differences in Malthusian fitnesses $f(t) = f_2(t) - f_1(t)$, which we consider to vary on a timescale longer than a generation but shorter than the experiment.

The Fokker–Planck form has also been used to describe the stochastic dynamics of neural activity. Brunel uses the form

$\tau \partial_t p(V,t) = -\partial_V\left[(\mu(t) - V)p(V,t)\right] + \frac{\sigma^2(t)}{2} \partial_V^2 p(V,t),$

in which $V$ stands for the membrane depolarization of a randomly chosen neuron, $\sigma^2(t)$ is the fluctuation in excitatory and inhibitory inputs, and $\mu(t)$ is the average input signal.

### Changes of Variable

The Fokker–Planck Equation can be put into more convenient form through changes of variable.

Following section 5.1 of Risken,

$y = \int_0^x \sqrt{\frac{D}{D(\xi)}} d\xi.$

For the particular case of population genetics with random mating,

$y = \sqrt{4N D}\int_0^x\frac{d\xi}{\sqrt{\xi(1-\xi)}}.$

Using the standard formula

$\int\frac{dx}{\sqrt{a x^2 + b x + c}} = -\frac{1}{\sqrt{-a}} \arcsin \left(\frac{2 a x + b}{\sqrt{b^2 - 4 a c}}\right),$

this can be reduced to

$\frac{y}{\sqrt{4N D}} = \frac{\pi}{2} - \arcsin(1 - 2x) \equiv \theta.$

The force term is also transformed:

$G(y,t) = \sqrt{\frac{D}{D(x)}} \left[F(x,t) - \frac{1}{2}\partial_x D(x)\right].$

Substituting in our expressions for $D(x)$ and for the force field $F(x,t)$ gives

$G(y,t) = \sqrt{\frac{4N D}{x(1-x)}} \left[(1-2x)\left(\mu - \frac{1}{8N}\right) + f(t)x(1-x)\right].$

We eliminate $x$ in favor of $y$ or, more conveniently, the normalized variable $\theta$. Noting that

$1 - 2x = \cos\theta$

and

$x(1-x) = \frac{\sin^2\theta}{4},$

we can after some algebra show that

$G(\theta,t) = 4\sqrt{N D} \left[\left(\mu - \frac{1}{8N}\right)\tan\theta + \frac{f(t)}{4} \sin\theta\right].$

We can avail ourselves of the mathematical tools developed in quantum mechanics if we manipulate the Fokker–Planck equation into the form of a Schrödinger equation. This is a special case of a general trick employed in the study of differential equations, by which one studies non-self-adjoint differential operators by transmogrifying them into self-adjoint ones. See chapter 18 of Gbur’s textbook. To see how this is done, we start with a Fokker–Planck equation written operatorially as

$\partial_t p(x,t) = \mathcal{L}_0 p(x,t),$

where $\mathcal{L}_0$ has a drift term linear in the derivative $\partial_x$ and a diffusion term quadratic in $\partial_x$:

$\mathcal{L}_0 = D\partial_x^2 -\partial_x h(x).$

The operator $\mathcal{L}_0$ is not Hermitian, since taking the conjugate changes the sign of the derivative $\partial_x$. (Recall from basic quantum mechanics that the momentum operator in the position basis is $-i\partial_x$, and that momentum, being an observable quantity, needs a Hermitian operator.) To make our Fokker–Planck time evolution look more like a Schrödinger equation, we shall try to invent a Hermitian analogue of $\mathcal{L}_0$. The trick we employ is to define

$\phi_0(x) \equiv \exp\left(\frac{1}{2D} \int_0^x h(y) d y\right),$

and to rewrite the probability distribution $p(x,t)$ as

$p(x,t) = \phi_0(x) \psi(x,t).$

As we shall soon see, the function $\psi(x,t)$ plays the role of the wavefunction in the Schrödinger-like version of our time-evolution equation. Because the Fokker–Planck operator involves first and second derivatives with respect to $x$, we note that

$\partial_x \phi_0(x) = \frac{h}{2D}\phi_0(x)$

and thus

$\partial_x^2 \phi_0(x) = \frac{h'}{2D}\phi_0(x) + \frac{h^2}{4D^2} \phi_0(x).$

The left-hand side of our Fokker–Planck equation is, using our new “wavefunction”,

$\partial_t p(x,t) = \phi_0(x) \partial_t \psi(x,t).$

If we move $\phi_0(x)$ to the other side of the equation, we get that

$\partial_t \psi(x,t) = \phi_0^{-1}\mathcal{L}_0\phi_0\psi(x,t).$

Working out $\mathcal{L}_0\phi_0\psi$ yields several terms, all of which contain a factor $\phi_0$, since every derivative of $\phi_0$ yields a $\phi_0$ again. Therefore, when we divide by $\phi_0$ we find

$\phi_0^{-1}\mathcal{L}_0\phi_0\psi = D\partial_x^2\psi + \left(-\frac{h'}{2} - \frac{h^2}{4D}\right)\psi.$

Both operators acting on $\psi$ are Hermitian, the first because it is quadratic in the derivative and the second because it contains only real functions of $x$.

Consequently, we define the new operator

$\mathcal{H} \equiv \phi_0^{-1}\mathcal{L}_0\phi_0,$

and we rewrite our equation in the form

$\partial_t \psi(x,t) = \mathcal{H}\psi(x,t).$

This is the Schrödinger-like formulation of the Fokker–Planck equation.

It is interesting to note what happens when the diffusion coefficient is not constant, i.e., when $D(x)$ is a non-constant function of $x$. If we define

$\phi_0^*(x) \equiv \exp\left(\int_0^x\frac{h(y)}{2D(y)} dy\right),$

then a series of manipulations following the steps above leads to the result

$\mathcal{H}^* = D\partial_x^2 + 2D'\partial_x + \left(-\frac{h'}{2} - \frac{h^2}{4D}\right) + \frac{D'h}{2D},$

which clearly reduces to the $\mathcal{H}$ given above when $D(x) = D$. Note that $\mathcal{H}^*$ contains a term linear in $\partial_x$, and thus is not Hermitian.

### Supersymmetry

Cooper, Klein and Sukhatme (1995) give a good general introduction to supersymmetry methods in quantum mechanics. Basically, it’s the solution to the simple harmonic oscillator using creation and annihilation operators, but all grown up. We start with a Hamiltonian $\mathcal{H}_1$, which we say we can factor into an operator and its adjoint:

$\mathcal{H}_1 = A^\dagger A,$

By exchanging the order of these factors, we obtain a new Hamiltonian, which we call the partner of the first:

$\mathcal{H}_2 = A A^\dagger.$

A typical Hamiltonian in one-dimensional quantum mechanics has a momentum term which in the position basis becomes a second derivative and a position-dependent term indicating the potential energy of the system. Scaling to discard unnecessary units, we can choose the operator $A$ to take the form

$A = -\partial_x + W(x),$

where $W(x)$ is known as the superpotential. As indicated above, taking the adjoint changes the sign on the spatial derivative:

$A^\dagger = \partial_x + W(x).$

In terms of the superpotential,

$\mathcal{H}_1 = -\partial_x^2 - W'(x) + W^2(x).$

Exchanging the order of the operators yields

$\mathcal{H}_2 = -\partial_x^2 + W'(x) + W^2(x).$

The Hamiltonians $\mathcal{H}_1$ and $\mathcal{H}_2$ are isospectral, as one can easily demonstrate. SUSY plus the extra condition known as shape invariance allows one to solve for eigenfunctions and eigenvalues using operator methods, like one does for the harmonic oscillator.

The Schr"odinger form of the Fokker–Planck equation, exhibits a supersymmetry of this form, between the Hamiltonian with force coefficient $F(x)$ and its partner with force coefficient $-F(x)$. We can factor the Fokker–Planck Hamiltonian given earlier using the two operators

$a = \sqrt{D}\partial_x - \frac{1}{2\sqrt{D}} F,$

$a^\dagger = -\sqrt{D}\partial_x - \frac{1}{2\sqrt{D}} F.$

With this choice of constant prefactors,

$\mathcal{H} = -a^\dagger a.$

The constants here were chosen for compatibility with Jung and Hänggi.

If we use the same operators $a$ and $a^\dagger$ but with $D$ a function of $x$, then the non-Hermitian operator $\mathcal{H}^*$ defined above can be written

$\mathcal{H}^* = -a^\dagger a + \frac{3D}{4}'\partial_x + \frac{D' F}{4D}.$

The derivative can be expressed in terms of the SUSY operators as

$\mathcal{H}^* = -a^\dagger a - \left(\frac{3D'}{4\sqrt{D}}\right)a^\dagger + \frac{5}{8} \frac{D'F}{D}.$

### Floquet Theory and Time Dependence

We wish to understand what happens when the influences on a diffusion process vary periodically over time. A particle trapped in a potential well may be subjected to an extra oscillating field; a population of microbes in a dish might be dividing once every half-hour in a laboratory which is dark for twelve hours out of every twenty-four. Differential equations with periodic driving forces are treated with Floquet theory. One way to apply Floquet theory to our diffusion equations is to rewrite a problem involving diffusion along a one-dimensional line as one of diffusion in a 2D space, with a new coordinate $\theta$ taking over the role of the time $t$.

Following Jung and Hänggi, we write

$\partial_t p = -\partial_x [h(x) + A f(\Omega t)]p + D\partial_x^2 p,$

and transform it into the 2D process

$\partial_t W(x,\theta;t) = (\mathcal{L}_0 - A f(\theta) \partial_x - \Omega \partial_\theta) W(x,\theta;t),$

in which the operator $\mathcal{L}_0$ is

$\mathcal{L}_0 = -\partial_x h + D\partial_x^2,$

capturing the tao of the Fokker–Planck. If we set the right-hand side of the Fokker–Planck equation to zero, we find a differential equation for the stationary probability density $W_{st}(x,\theta)$:

$(\mathcal{L}_0 - f(\theta)\partial_x - \Omega\partial_\theta) W_{st}(x,\theta) = 0.$

The probability distribution in the 1D view of the problem, $p(x,t)$, does not tend to a steady-state solution, but its asymptotic behavior is just given by $W_{st}(x,\theta)$.

$p_{as}(x,t) = 2\pi W_{st}(x,\theta),$

with

$t = \frac{\theta}{\Omega}.$

The factor of $2\pi$ arises because we have compactified the $\theta$ coordinate in order to represent the temporally periodic influence $f(t)$.

If we soup up the 1D process just a bit,

$\partial_t p = -\partial_x [h(x) + g(x) f(\Omega t)] p + D\partial_x^2 p,$

then the derivative with respect to $x$ gives us one more term:

$\partial_t W = [\mathcal{L}_0 - g' f(\theta) - g f(\theta)\partial_x - \Omega\partial_\theta]W,$

where we might as well define the expression in brackets to be a new operator $\mathcal{L}$.

Echoing the formalism developed in the section above, we define the operator

$a = \sqrt{D}\partial_x - \frac{1}{2\sqrt{D}} h,$

whose conjugate is

$a^\dagger = -\sqrt{D}\partial_x - \frac{1}{2\sqrt{D}} h.$

The product of these operators acting on a test function $\psi$ is

$a^\dagger a \psi = -D \partial_x^2 \psi + \left(\frac{h'}{2} + \frac{h^2}{4D}\right)\psi.$

Exchanging the order of $a^\dagger$ and $a$ merely changes the sign on the $h'$ term, allowing us to deduce that

$[a, a^\dagger] = -h'.$

If we attempt to construct a Hermitian version of $\mathcal{L}$,

$\mathcal{H} \equiv \phi_0^{-1}\mathcal{L}\phi_0,$

wherein

$\phi_0(x) = \exp\left(\frac{1}{2D} \int_0^x h(y) d y \right),$

we find that

$\mathcal{H} = -a^\dagger a - f g\partial_x -\Omega\partial_\theta - f g' - \frac{f g h}{2D}.$

The terms involving $f$ and $g$ can be combined by noting that $a^\dagger$ contains a derivative with respect to $x$.

$\mathcal{H} = -a^\dagger a + \frac{f a^\dagger}{\sqrt{D}}g - \Omega\partial_\theta.$

If the force $h(x)$ is exchanged with $\tilde{h}(x) \equiv -h(x)$, and these steps are repeated, one deduces that

$\tilde{\mathcal{H}} = -a a^\dagger - \frac{f a}{\sqrt{D}}g - \Omega\partial_\theta.$

In SUSY QM, the operators $a$ and $a^\dagger$ map eigenstates of one Hamiltonian into eigenstates of its partner. We consider, therefore, products of these operators with our Hamiltonians. By straightforward algebra,

$a\mathcal{H}^\dagger = -a a^\dagger a + \frac{a f}{\sqrt{D}} a g + a\Omega\partial_\theta,$

and

$\tilde{\mathcal{H}} a = -a a^\dagger a + \frac{f a}{\sqrt{D}}g a - a\Omega\partial_\theta.$

Using $\mathfrak{T}$ to denote the operation of $\theta$ reversal,

$a\mathfrak{T}\mathcal{H}^\dagger = -a a^\dagger a + \frac{(\mathfrak{T}f)}{\sqrt{D}} a g - a\Omega\partial_\theta.$

Subtracting the previous two equations gives

$a(\mathfrak{T}\mathcal{H}^\dagger) - \tilde{\mathcal{H}}a = \frac{a}{\sqrt{D}}\left[(\mathfrak{T}f)a g - f g a\right].$

Because $f(\theta)$ is periodic, we might as well call it even.

$a(\mathfrak{T}\mathcal{H}^\dagger) - \tilde{\mathcal{H}}a = \frac{a f}{\sqrt{D}} [a, g].$

If $[a, g] = 0$, then we recover the case studied by Jung and Hänggi, in which a SUSY exists between $\mathcal{H}^\dagger$ and $\tilde{\mathcal{H}}$. To illustrate, consider the eigenfunctions of $\mathcal{H}$ and $\mathcal{H}^\dagger$, defined by

$\mathcal{H}\Sigma(x,\theta) = \lambda\Sigma(x,\theta)$

and

$\mathcal{H}^\dagger\Lambda(x,\theta) = \lambda\Lambda(x,\theta)$

Left-multiplying the previous equation by $a$ and acting with $\mathfrak{T}$ shows that

$a(\mathfrak{T}\mathcal{H}^\dagger)\Lambda(x,-\theta) = \lambda a\Lambda(x,-\theta),$

and making the substitution $a(\mathfrak{T}\mathcal{H}^\dagger) = \tilde{\mathcal{H}}a$, we can write,

$\tilde{\mathcal{H}}a \Lambda(x,-\theta) = \lambda a\Lambda(x,-\theta).$

The operator $a$ maps eigenstates of $\mathfrak{T}\mathcal{H}^\dagger$ into those of $\tilde{\mathcal{H}}$.

To understand the effect of the function $g(x)$ not being constant, we need to understand the commutator $[a, g]$. Because $[g, h] = 0$, we only need to consider the term in $a$ containing the derivative $\partial_x$. That is,

$[a, g] = \sqrt{D}[\partial_x, g].$

When we let this operator act on a wavefunction $\psi$, the surviving term is proportional to $g'\psi$. We get that

$a(\mathfrak{T}\mathcal{H}^\dagger - \sqrt{D}f g') - \tilde{\mathcal{H}}a = 0.$

Thus, it appears that we can employ SUSY for the Brunel description of neural dynamics, at least for the special case of constant variation $\sigma^2$, but not for population genetics with periodically varying fitness, at least not easily.

### Comparison to Simulation

To understand the long-term behavior of solutions of the Fokker–Planck equation, we employ measures which summarize salient features of the probability distribution $p(x,t)$. For example, we are curious to know how much $p(x,t)$ will have spread out when $t$ has grown rather large. A natural choice is to examine the variance, $\langle x^2(t) \rangle_c = \langle x^2(t) \rangle - \langle x(t)\rangle^2$, scaled by the elapsed time:

$D_{eff} \equiv \lim_{t \rightarrow\infty} \frac{\langle x^2(t)\rangle - \langle x(t)\rangle^2} {2t}.$

This is the effective diffusion coefficient defined by Riemann et al. and used by Sahoo et al. in the case of periodic partner potentials. It can be shown analytically that the effective diffusion coefficients for such partner potentials $\pm V(x)$ are in fact equal. Furthermore, Sahoo et al. find through numerical simulation that the effective diffusion coefficients are also equal in the case of periodic forcing.

Let $t(x_0\rightarrow b)$ denote the tiime required to reach $b$ from $x_0$. Then

$T_n(x_0\rightarrow b) \equiv \langle t^n(x_0\rightarrow b) \rangle.$

For the case of a particle in a tilted periodic potential with

$V(x) \equiv V_0(x) - x F$

and

$V_0(x) = V_0(x + L),$

the moments $T_n$ satisfy the recursion relation

$T_n(x_0\rightarrow b) = \frac{n}{D_0} \int_{x_0}^b dx\, \exp\left(\frac{V(x)}{k_B T}\right) \int_{-\infty}^x dy\, \exp\left(-\frac{V(y)}{k_B T}\right) T_{n-1}(x_0\rightarrow b),$

where $T_0(x_0\rightarrow b)$ is defined to equal unity. The effective diffusion $D_{eff}$ can be derived in terms of the “first passage time dispersion”

$\Delta T_2(x_0\rightarrow b) \equiv \langle t^2(x_0\rightarrow b) \rangle - \langle t(x_0\rightarrow b)\rangle^2,$

which is just the second cumulant, of course.

The units of $D_{eff}$ are length squared over time. We have one obvious quantity having units of length, namely the potential periodicity $L$. The dispersion $\Delta T_2$ carries units of time squared, so because we expect $D_{eff}$ to grow with increasing dispersion, we have to divide by a quantity with units of time cubed. This dimensional analysis is our heuristic justification for the result of Riemann et al.,

$D_{eff} = \frac{L^2}{2} \frac{\Delta T_2(x_0 \rightarrow x_0 + L)} {T_1(x_0 \rightarrow x_0 + L)},$

which they derive using a coarse-graining argument.

In physical units, the Fokker–Planck equation for this process is

$\partial_t p(x,t) = \partial_x V'(x) p(x,t) + k_B T\partial_x^2 p(x,t).$

Writing $D_0$ for the bare diffusion constant $k_B T/\gamma$, we define

$I_+(x) \equiv \frac{1}{D_0} e^{\frac{V(x) - Fx}{k_BT}} \int_{x-L}^x dy\, \exp\left(\frac{-V(y) + Fy}{k_BT}\right),$

and the complementary quantity

$I_-(x) \equiv \frac{1}{D_0} e^{\frac{-V(x) - Fx}{k_BT}} \int_x^{x+L} dy\, \exp\left(\frac{V(y) - Fy}{k_BT}\right).$

The effective diffusion is

$D_{eff} = D_0 \frac{\int_{x_0}^{x_0 + L} I_\pm(x) I_+(x) I_-(x)\, dx} {\left[\int_{x_0}^{x_0 + L} I_\pm(x)\, dx\right]^3}.$

Here, the subscript $\pm$ indicates that the index on that function may be chosen to be either $+$ or $-$, and $x_0$ is an arbitrary reference point of convenience.

The symmetry of the factors in this expression indicates that under the exchange $V(x) \rightarrow -V(x)$, $F \rightarrow -F$ the effective diffusion is unchanged.

Reichl and friends find that the mean first-passage time for a Brownian rotor is governed by the lowest Floquet quasi-eigenvalues of the rotor’s Fokker–Planck equation, and varies as the inverse of the lowest Floquet decay rate. That is to say, the Floquet procedure naturally provides quantities with the dimensions of time, and those quantities are the natural timescales for phenomena of interest. By analogy, we expect that the first-passage time in a periodic potential will behave in a similar manner; because the Floquet decay rates are the eigenvalues of the corresponding 2D Fokker–Planck problem, the SUSY partner of that problem will map to a 1D diffusion problem with the same Floquet decay rates. We expect, therefore, that Floquet timescales will be the same for SUSY partners. Furthermore, quantities derived from these timescales, such as the effective diffusion, should also be invariant under SUSY exchange.

### Open Questions

Two questions call themselves to mind, both concerning generalizations of the diffusion problems treated above. First, there is the situation of population genetics with time-varying fitness, and second, we have the case of stochastic neural dynamics with time-varying inputs. In the former problem, we aim to understand the eigenvalues of the operator

$\mathcal{L}_{pop} = \mathcal{L}_0 - g'(x)f(\theta) - g(x)f(\theta)\partial_x - \Omega\partial_\theta.$

As deduced above, it’s the cross term $g(x)f(\theta)\partial_x$ which makes this generalization tricky.

The second situation also involves a generalized time-evolution operator. Suppose that the external drive $\mu(t)$ is periodically varying, and that the fluctuation $\sigma^2(t)$ which provides the diffusion coefficient has a periodic variation on top of a constant:

$\frac{\sigma^2(t)}{2\tau} = D + f_2(t).$

Then, we can rewrite the Fokker–Planck equation in the form $\partial_t p(x,t) = \mathcal{L}p(x,t)$, and changing over to the 2D perspective, we can cook up the new time-evolution operator,

$\mathcal{L}_{neuro} = \mathcal{L}_0 - f_1(\theta)\partial_x + f_2(\theta) \partial_x^2 - \Omega\partial_\theta.$

Again, we are concerned with the eigenvalues of this operator, which provide the natural timescales of interest for our problem, and in particular, we’d like to know if this operator is isospectral with a SUSY partner.

Finally, there is the question of shape invariance. Supposing that a SUSY exists between partner time-evolution operators, does this extra condition also hold, implying analytic solubility by algebraic methods?

• P. Alpatov and L. E. Reichl (1994), “Spectral properties of a time-periodic Fokker-Planck equation” Physical Review E 49: 2630–38.
• F. Cooper, A. Khare and U. Sukhatme (1995), “Supersymmetry and quantum mechanics” Physics Reports 251: 267–385, arXiv:hep-th/9405029.
• N. Brunel (2000), “Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons” Journal of Computational Neuroscience 8: 183–208.
• Gregory J. Gbur (2011), Mathematical Methods for Optical Physics and Engineering. Cambridge University Press.
• P. Jung and P. Hänggi (1991), “Amplification of small signals via stochastic resonance” Physical Review A 44, 12: 8032.
• M. Kardar (2005), lecture notes for 8.592J: Statistical physics in biology.
• P. Reimann et al. (2001), “Giant acceleration of free diffusion by use of tilted periodic potentials” Physical Review Letters 87, 1: 010602.
• P. Reimann et al. (2002), “Diffusion in tilted periodic potentials: Enhancement, universality, and scaling” Physical Review E 65, 3: 031104.
• H. Risken (1984), The Fokker–Planck Equation. Springer.
• M. Sahoo et al. (2006), “Supersymmetry and Fokker–Planck dynamics in periodic potentials”, arXiv:cond-mat/0610294.

category: blog