Diffusion, Convolution
The following notes are
straight out of Howard Berg's ``Random Walks in Biology'' with a few ``fun''
additions.
## Walker on a Line
Let's look at a bunch of particles moving on a line in 1-D. Assume that after
every time interval of $\tau$ seconds, our particles move with equal
probability to the left or to the right, with some characteristic step size
$\delta=v_x \tau$. It's important to note that
(*) $v_x,\tau$ will depend on the size of our particle, and, since we're
talking about Statistical Mechanics, the medium in which it is immersed, as
well as the temperature. You should know this already since $kT \sim \langle m
v_x^2 \rangle$.
(*) The particle's steps at each moment should be thought of as independent,
or, uncorrelated random variables; otherwise things will get very, very
difficult. The snootier description of this point is the phrase i.i.d. :
independently identically distributed, which means, if we write $s_i(n)$ as
the $n^\mathrm{th}$ step of the $i^\mathrm{th}$ particle:
\begin{eqnarray}
\langle s_i(n) s_I(n+1) \rangle = 0
\end{eqnarray}
(*) We should also assume that each particle moves independently from all
others. This is a simple boiler plate ``dilute'' or ideal assumption.
If the probability of going right and left is equal -- let's say $p$ and $q$
respectively -- then we expect our particles to go nowhere on average. But,
interestingly enough, the average of the \textit{square} of our displacement
for each particle will be nonzero. If we denote $x_i(n)$ as the position of
the $i^\mathrm{th}$ particle after $n$ steps, then
\begin{eqnarray}
\langle x_i(n)^2 \rangle \sim n
\end{eqnarray}
Which is typical of a diffusion process.
## Linear Recurrence Relation
The following \textbf{ linear recurrence relation} is a statement of our
particle movement:
\begin{eqnarray}
x_i(n) &=& x_i(n-1) \pm \delta
\end{eqnarray}
After every time step $n-1 \to n$, we have some change to our position, which
is plus or minus $\delta$. If we average over all our particles, we get
something called the \textbf{Martingale} property:
\begin{eqnarray}
\langle x(n) \rangle &=& \frac{1}{N}\sum_{i=1}^N x_i(n)\\
&=& \frac{1}{N}\sum_{i=1}^N x_i(n-1) \pm \delta
\end{eqnarray}
Since half our particles in the sum will have the plus sign, and half the
negative our second term on the RHS sums to zero and we get:
\begin{eqnarray}
\langle x(n) \rangle &=& \langle x(n-1) \rangle
\end{eqnarray}
This is a fancier way of saying: the particles don't move, on average. But,
what about the average of the \textbf{square} of $x(n)$?
\begin{eqnarray}
\langle x(n)^2 \rangle &=& \frac{1}{N}\sum_{i=1}^N x_i(n)^2\\
&=& \frac{1}{N}\sum_{i=1}^N \left( x_i(n-1)^2 \pm 2 \delta x_i(n-1) +
\delta^2 \right)
\end{eqnarray}
Once again, our second term on the RHS will vanish as we sum, and we are left
with:
\begin{eqnarray}
\langle x(n)^2 \rangle &=& \langle x(n-1)^2 \rangle+\delta^2 \\
&=& n \delta^2
\end{eqnarray}
Where in the second line I've worked back through all the recursive steps,
assuming $\langle x^(0)\rangle=0$. So the variance of our particle position is
linearly increasing with time! If we define $\frac{\delta^2}{2 \tau}=D$, we
can write this as:
\begin{eqnarray}
\langle x^2(t) \rangle &=& \frac{t}{\tau} \delta^2 \\
\langle x^2(t) \rangle &=& 2 D t
\end{eqnarray}
## Convolution Approach
But, you guys already knew this! Because we can view $\delta$ as a simple
random variable drawn from some probability distribution $\delta \sim
P(\delta)$. Let's examine the binomial distribution
\begin{eqnarray}
P(k,n) &=& p^k q^{n-k} \left(\begin{array}{c}n \\k \end{array}\right)
\\
q &=& 1-p
\end{eqnarray}
This is the probability of $k$ ``successes'' after $n$ trials. It has the
first two cumulants:
\begin{eqnarray}
\langle k \rangle &=& np \\
\langle k^2 \rangle -\langle k \rangle^2 &=& npq
\end{eqnarray}
For the time being, let's say that after every time interval $\tau$ that
particle has only had one chance, one trial, to move left or right: this means
we set $n=1$. (We will take the $n \to \infty$ limit later, which corresponds
to Gaussian statistics and temporal coarse graining.) If we define a random
variable $b$:
\begin{eqnarray}
b &\sim& P(k,1) \\
b &=& 0,1
\end{eqnarray}
we can scale this variable such that we get our proper $\pm \delta$ with
probability $p$ and $q$ -- which are now no longer equal:
\begin{eqnarray}
s &=& 2b-1 \\
s &=& \pm 1
\end{eqnarray}
and so we write:
\begin{eqnarray}
x_i(n) &=& x_i(n-1)+s \delta
\end{eqnarray}
$b$ is now a random variable, expectation values:
\begin{eqnarray}
\langle b \rangle &=& p \\
\langle b^2 \rangle -\langle b \rangle^2 &=& pq \\
\langle b^2 \rangle &=& pq + p^2
\end{eqnarray}
which, means our statistics for $s$ are:
\begin{eqnarray}
\langle s \rangle &=& 2p-1=p-q \\
\langle s^2 \rangle &=& \langle 4b^2-4b+1 \rangle \\
&=& 4pq+4p^2-4p +1 \\
&=& 4pq+4p(p-1) + 1 \\
&=& 4pq - 4pq + 1\\
\langle s^2 \rangle &=&1\\
\langle s^2 \rangle -\langle s \rangle &=& 1-(p-q)^2
\end{eqnarray}
So to summarize:
\begin{eqnarray}
\langle s \rangle &=& 2p-1=p-q \\
\langle s^2 \rangle &=&1\\
\langle s^2 \rangle -\langle s \rangle &=& 1-(p-q)^2
\end{eqnarray}
We now have everything we need for our linear recurrence relation, and not
that since $p\neq q$, we're going to have drifting particle. Treating all
quantities now as random variables rather than ensemble averages over our 1-D
``gas'' of walkers:
\begin{eqnarray}
\langle x_i(n) \rangle &=& \langle x_i(n-1) + s\delta \rangle \\
\langle x_i(n) \rangle &=& \langle x_i(n-1) \rangle + (p-q)\delta \\
\langle x_i(n) \rangle &=& n (p-q) \delta
\end{eqnarray}
We see that for $p=q$ we have no drift, or no net movement (martingale). For
$p>q$ we have net movement to the right (super-martingale) and $p<q$ we
have net movement to the left. Now for variance:
\begin{eqnarray}
\langle x_i(n)^2 \rangle &=& \langle (x_i(n-1) + s\delta)^2 \rangle \\
\langle x_i(n)^2 \rangle &=& \langle x_i(n-1)^2 + s^2\delta^2 +
2sx_i(n-1) \rangle \\
\langle x_i(n)^2 \rangle &=& \langle x_i(n-1)^2 \rangle + \langle s^2
\rangle \delta^2 + \langle 2sx_i(n-1) \rangle \\
\langle x_i(n)^2 \rangle &=& \langle x_i(n-1)^2 \rangle + \delta^2 +
\langle 2s \rangle \langle x_i(n-1) \rangle
\end{eqnarray}
Using our former expression, we get:
\begin{eqnarray}
\langle x_i(n)^2 \rangle &=& \langle x_i(n-1)^2 \rangle + \delta^2 +
\langle 2s \rangle n (p-q) \delta \\
\langle x_i(n)^2 \rangle &=& \langle x_i(n-1)^2 \rangle + \delta^2 +
2(p-q)^2 n \delta
\end{eqnarray}
This last formula looks like a mess, and turns out to be much easier to write
down the variance of $x_i(n)$:
\begin{eqnarray}
\langle x_i(n)^2 \rangle-\langle x_i(n) \rangle^2 &=& n \delta^2
\left(1-(p-q)^2 \right)
\end{eqnarray}
Now, one might ask, how did I do this? You can do the grunt work, via ensemble
averaging $\sum_{i=1}^N$ or using our equation above if you like, or you can
\textbf{remember that} $x_i(n)$ \textbf{is drawn from a sum of random
variables, which means convolution, which means sum of cumulants}! As we
learned earlier in class,
\begin{eqnarray}
X_i & \sim & P(x) \\
Z &=& X_1+X_2 \\
P(Z) &\sim & (P \star P)(z) \\
\langle Z \rangle &=& \langle X_1 \rangle + \langle X_2 \rangle \\
\mathrm{Var}(Z) &=& \mathrm{Var}(X_1) + \mathrm{Var}(X_2)
\end{eqnarray}
This makes it easy to write, down, immediately:
\begin{eqnarray}
E\left(x_i(n)\right) &=& n (p-q) \delta \\
\mathrm{Var}\left(x_i(n)\right) &=& n \delta^2 \left(1-(p-q)^2 \right)
\end{eqnarray}
For each and every particle, $x_i$.
## The $n \to \infty$ limit
Now, let's assume that we are quite lazy in observing our particle, and that,
over each time interval $\tau$, it has in fact \textbf{many, many}
opportunities to move right or left. This means that our distribution on step
successes is a binomial:
\begin{eqnarray}
P(k,n) &=& p^k q^{n-k} \left(\begin{array}{c}n \\k \end{array}\right)
\end{eqnarray}
but, with $n$ very very large, and $np$ equal to some finite constant, $\mu$.
It is a somewhat easy proof to show that in this limit, as $n \to \infty$ and
$np=\mu$, $npq=\sigma^2$ our Binomial distribution becomes a Gaussian:
\begin{eqnarray}
P(k, `` \infty '') &=& \frac{1}{\sqrt{2\pi \sigma^2}}
e^{\frac{(k-\mu)^2}{2\sigma^2}}
\end{eqnarray}
So now, our step variable, $\delta$ is drawn from some Gaussian distribution,
with mean $\mu$ and variance $\sigma^2$. We won't do the tedious change of
variables again here, because you see how this is done. The point is, we have
gone from a binomial-choice process to a Gaussian one. And it's important to
note that our walker, after sequence of $n$ steps can be viewed as a sum of
independent $\delta$'s:
\begin{eqnarray}
x_i(n) &=& x_i(n-1)+\delta_n \\
x_i(n) &=& \sum_{j=1}^n \delta_j \\
\delta_j &\sim & \mathcal{N}(\mu,\sigma^2)= \frac{1}{\sqrt{2\pi
\sigma^2}} e^{\frac{(\delta-\mu)^2}{2\sigma^2}}
\end{eqnarray}
This means, by addition of cumulants alone, that:
\begin{eqnarray}
\mathrm{E}\left(x_i(n) \right) &=& n \mu \\
\mathrm{Var}\left(x_i(n) \right) &=& n \sigma^2
\end{eqnarray}
Once again, we find a non-zero drift -- possibly -- and a linear scaling of
variance with time. It is important thing to note that
(****)The convolution of a Gaussian is another Gaussian, with mean
$\mu_{12}=\mu_1+\mu_2$ and variance $\sigma_{12}=\sigma_1+\sigma_2$.
So, for example, if our particle is initially placed at position $x=0$, it's
probability density is a dirac delta function:
\begin{eqnarray}
x_i(0) &\sim & \delta(x)
\end{eqnarray}
After a single step $\delta$, the new probability density is the convolution:
\begin{eqnarray}
\delta &\sim & P(\delta)= \frac{1}{\sqrt{2\pi \sigma^2}}
e^{\frac{(\delta-\mu)^2}{2\sigma^2}}\\
x_i(1) &=& x_i(0)+\delta \\
x_i(1) &\sim & (\delta \star P)(x) \\
(\delta \star P)(x) &=& \int dx^\prime \delta(x-x^\prime) P(x^\prime)
\\
&=& P(x)=\mathcal{N}(\mu,\sigma^2)
\end{eqnarray}
Or, just a Gaussian. So that was pretty easy, since the convolution of a
function with Dirac delta is just the function itself. But what about after
the second step? Then we have:
\begin{eqnarray}
x_i(2) &\sim & (P \star P)(x) \\
(P \star P)(x) &=& \mathcal{N}(2\mu,2\sigma^2)= \frac{1}{\sqrt{2\pi
2\sigma^2}} e^{\frac{(x-2\mu)^2}{4\sigma^2}}
\end{eqnarray}
You can see how this goes, we are going to have:
\begin{eqnarray}
x_i(n) & \sim & (P \star P \star \cdots \star P)(x) \\
(P \star P \star \cdots \star P)(x) &=& \mathcal{N}(n\mu,n\sigma^2)=
\frac{1}{\sqrt{2\pi n\sigma^2}} e^{\frac{(x-n\mu)^2}{2n\sigma^2}}
\end{eqnarray}
This probability density is a Gaussian that is expanding its variance linearly
with time and drifting with some $\mu$! If we denote $\sigma^2n=\sigma^2
t/\tau=2Dt$, and $\frac{\mu}{\tau}=v$, we can write in more conventional form:
\begin{eqnarray}
P(x(t))&=& \mathcal{N}(vt,2Dt)= \frac{1}{\sqrt{4\pi Dt}}
e^{\frac{(x(t)-vt)^2}{4Dt}}
\end{eqnarray}
## What do you mean, Diffusion?
Now, one way for us to think about flux, or the dynamics of many particles is
to count how many are crossing a surface at each time step. If set a ``net''
in the $y-z$ plane at $x=\delta/2$, we can write the number of particles going
through this net as:
\begin{eqnarray}
\frac{1}{2}\left(-N(x+\delta)+N(x) \right)
\end{eqnarray}
Since, half the particles on the left of the net -- $N(x)$ -- will move
through it, to the right, and half the particles on the right of the net --
$N(x+\delta)$ -- will move through the opposite way, to the left. Dividing by
the area of net and the amount of time $\tau$, we get flux:
\begin{eqnarray}
J_x &=& \frac{1}{2 A\tau}\left(-N(x+\delta)+N(x) \right)
\end{eqnarray}
and, multiplying by $\delta^2/\delta^2$, we can write
\begin{eqnarray}
J_x &=& \frac{\delta^2}{2\tau
\delta}\left(\frac{-N(x+\delta)+N(x)}{\delta A} \right)\\
&=& \frac{\delta^2}{2\tau
\delta}\left(-\rho(x+\delta)+\rho(x)\right)\\
J_x &=& -\frac{\delta^2}{2 \tau \delta}\frac{\partial \rho}{\partial
x}
\end{eqnarray}
and know we notice our former constant $\delta^2/\tau$ as $D$,
\begin{eqnarray}
J_x &=& - D \frac{\partial \rho}{\partial x}
\end{eqnarray}
or in more general form, since we could placed our surface along any axis:
\begin{eqnarray}
\vec{J} &=& - D \vec{\nabla}\rho(\mathbf{x})
\end{eqnarray}
You see that flow will move down concentration gradients with some
``reaction'' coefficient $D$. This Fick's first law, and if we combine it with
the continuity equation -- which is really just a sanity check for
conservation of particles -- we get Fick's second law:
\begin{eqnarray}
\vec{\nabla} \cdot \vec{J} + \frac{\partial \rho}{\partial t} &=& 0 \\
- D \nabla^2 \rho(\mathbf{x})+ \frac{\partial \rho}{\partial t} &=&
0\\
\left(\frac{\partial}{\partial t}- D \nabla^2\right) \rho(\mathbf{x},t)
&=& 0
\end{eqnarray}
This is the \textbf{diffusion} equation, and it shows that our concentration
will only reach steady state when $\nabla^2 \rho = 0$.
## Green's Function for the Diffusion equation
Now, we would like to find something called a Green's function for the
diffusion equation, which basically tells us, how does the system react if we
put a tiny bit of particles at some point in space, at some point in time,
$x,t$:
\begin{eqnarray}
\left(\frac{\partial}{\partial t}- D \nabla^2\right) G(\mathbf{x},t)
&=& \delta^D(x,t)
\end{eqnarray}
The term on the write is our source function, and basically says that at time
$t=0$ we put a small amount of concentration at $x=0$. This equation can be
solved analytically, and in fact, it is no different the Green's function for
the Schrodinger equation of the free particle:
\begin{eqnarray}
\left(i\hbar \frac{\partial}{\partial t}+\frac{\hbar^2}{2m}\nabla^2
\right)\psi(x,t) &=& V(x) \psi(x)
\end{eqnarray}
If we set $V(x)=0$ and $\frac{i\hbar}{2m}=D$ we just have the diffusion
equation again.
The best way to find a green's function, or a delta function response is by
fourier transform. Let's write:
\begin{eqnarray}
G(x,t) &=& \int \frac{d^3 k d\omega}{(2\pi)^4} e^{i(kx-\omega
t)}\tilde{G}(k,\omega)
\end{eqnarray}
This means:
\begin{eqnarray}
\tilde{G}(k,w) &=& \frac{i}{w-ik^2D}\\
G(x,t) &=& \int \frac{d^3k dw}{(2\pi)^4} \frac{i
e^{i(kx-wt)}}{w-ik^2D}
\end{eqnarray}
The best way to approach this integral is via contour integration in the
complex plane. We wick rotate first, and send $w \to -iw$ to get:
\begin{eqnarray}
G(x,t) &=& \int \int_{i \infty}^{-i \infty}\frac{d^3 k
dw^\prime}{(2\pi)^4} \frac{e^{-w^\prime t}}{i(w^\prime-k^2 D)}
\end{eqnarray}
This integrand has a simple pole at $w^\prime=k^2 D$, and we now use the
residue theorem to write:
\begin{eqnarray}
G(x,t) &=& \int \frac{d^3 k}{(2\pi)^2} \frac{2\pi i}{2\pi}
\mathrm{Res}\left(\frac{e^{-w^\prime t}}{i(w^\prime-k^2 D)}\right)
\end{eqnarray}
We see that the $i$'s and $2\pi$'s cancle, and we are left with the fourier
transform of a gaussian:
\begin{eqnarray}
G(x,t)&=& \int \frac{d^3 k }{(2\pi)^3} e^{i \mathbf{k}\cdot
\mathbf{x}- k^2 D t}\\
G(x,t) &=& \frac{1}{\sqrt{4\pi D t}}e^{-\frac{x^2}{4Dt}}
\end{eqnarray}
This Green's function, which is the delta function response of a diffusive
system, is exactly the probability distribution we derived earlier, for a
particle that initially had some probability density:
\begin{eqnarray}
P(x,t=0) &=& \delta^D(x)
\end{eqnarray}
and at later times:
\begin{eqnarray}
P(x,t) &=& \frac{1}{\sqrt{4\pi D t}}e^{-\frac{x^2}{4Dt}}
\end{eqnarray}
So, the Green's function is in fact the fundamental solution, or the
probabilistic kernel for all diffusion problems. If we have some source of
particles in a system, or, some source of probability flux -- since we are now
treating both on the same grounds -- we have:
\begin{eqnarray}
\left( \frac{\partial}{\partial t}-D \nabla^2 \right)P(x,t) &=& j(x,t)
\end{eqnarray}
The solution to this in-homogeneous differential equation is a convolution of
our source with our Green's function:
\begin{eqnarray}
P(x,t) &=& (G \star j)(x,t) \\
&=& \int d^3 x^\prime dt^\prime G(x-x^\prime,
t-t^\prime)j(x^\prime,t^\prime)
\end{eqnarray}
We can now solve for the probability distribution of any system, for any
``forcing'' function, or source $j(x,t)$, subject to diffusion.