\documentclass[12pt]{article}
\usepackage{mysymbols}
\newcommand\En{\mathcal{E}}
\newcommand\Trace{\textrm{Tr}\,}
\newcommand\Diff{\tensor{\mathsf{D}}^{(2)}}
\newcommand\Unit{\tensor{\mathsf{I}}}
\newcommand\p{\mathbf{p}} % momentum vector
\newcommand\vvec{\mathbf{v}}
\newcommand\vrel{v_{\rm rel}}
\newcommand\vrelv{\vvec_{\rm rel}}
\newcommand\V{\mathbf{V}}
\newcommand\pu{\hat{p}}
\newcommand\vu{\hat{v}}
\newcommand\xu{\hat{x}}
\newcommand\yu{\hat{y}}
\newcommand\zu{\hat{z}}
\newcommand\unit{\hat{n}}
\newcommand\VCM{\mathbf{U}_\mathrm{CM}}
\newcommand\pCM{\p_\mathrm{CM}}
\newcommand\vCM{v_\mathrm{CM}}
\newcommand\vCMv{\mathbf{v}_\mathrm{CM}}
\newcommand\Pperp{\tensor{\mathsf{P}}_\perp}
\newcommand\Iperp{\tensor{\mathsf{I}}_{\perp}}
\newcommand\Gvec{\mathbf{G}}
\newcommand\Dvec{\mathbf{D}^{(1)}}
\title{Fokker-Planck Equation (FPE)}
\author{Nikolai G. Lehtinen}
\begin{document}
\maketitle
\section{Introduction}
Consider a distribution of particles $f(x_1,\dots,x_n,t)$ in $n$-dimensional
space with coordinates $x_k$ which also is a function of time $t$. Its
evolution  
is described by the Fokker-Planck equation (FPE), the derivation of
which is beautifully explained in \cite{nicholson}.
Generalizing this derivation to the case when distribution function depends on
several variables, we get
\begin{equation}
\partial_t f=\sum_{l=1}^\infty (-1)^l \sum_{k_1,\dots,k_l}
\partial_{x_{k_1}}\cdots\partial_{x_{k_l}} (f D^{(l)}_{k_1,\dots,k_l})
\label{eq:genfp}
\end{equation}
where $\partial$ denotes the partial derivative in respect to the
variable in the subscript. The tensors $D^{(l)}_{k_1,\dots,k_l}$ are
\[
D^{(l)}_{k_1,\dots,k_l}=\frac{1}{l!}\lim_{\Delta t \rightarrow \mbox{`0'}}
\frac{\langle\Delta x_{k_1}\cdots\Delta x_{k_l}\rangle}{\Delta t}
\]
where $\Delta x_{i}$ are the changes in variables $x_{i}$ in time
$\Delta t$ due to random collisions and other processes.
The limit  $\Delta t \rightarrow \mbox{`0'}$ means that we let $\Delta t$ 
become small  compared to any macroscopic time scale, while it still
remains large 
compared with time intervals between collisions.
The angular brackets $\langle\cdot\rangle$ denote the statistical (or ensemble)
average. The FPE is understood as (\ref{eq:genfp}) truncated after $l=2$.

\section{FPE in momentum space and interpretation of coefficients $D$}

We consider FPE in momentum space,
\begin{equation}
 \pderiv{f(\p)}{t}=-\nabla_\p[\Dvec f(\p)]+\nabla_\p[\nabla_\p(\Diff f(\p))]
\label{eq:fp}
\end{equation}
where
\[
  \Dvec= \lim_{\Delta t \rightarrow \mbox{`0'}} 
\frac{\langle\Delta \p\rangle}{\Delta t}
\]
and $\Diff$ is a second rank tensor
\[
  \Diff= \frac{1}{2}\lim_{\Delta t \rightarrow \mbox{`0'}} 
\frac{\langle\Delta \p\Delta\p^T\rangle}{\Delta t}
\]
Here $\Delta \p \Delta\p^T$ is a dyadic product of vector $\Delta\p$ 
with itself. In the second term of (\ref{eq:fp}), the inner
$\nabla_\p$ is the divergence of a tensor (the result of which is a
vector).

In general, $\Dvec$ is due to changes in particle momenta due to both
external force and random collisions. Without external force, it is
usually directed along $(-\p)$, and slows the particle down. However,
the decrease of average vector momentum can be due to spreading in
particle directions (which is described by the term with $\Diff$) and
thus not equal to the decrease of the absolute value of momentum
vector. Since
$\Delta(p^2)=2\p\cdot\Delta\p+(\Delta p_x)^2+(\Delta p_x)^2+(\Delta p_z)^2$,
we have an effective ``force'' (the dynamic friction coefficient)
\[
\mathbf{F}_1=\pu\frac{\langle\Delta\abs{\p}\rangle}{\Delta t}
=\frac{\pu}{2\abs{\p}}\frac{\langle\Delta (p^2)\rangle}{\Delta t}
=\Dvec+\frac{\pu}{\abs{\p}}\Trace{\Diff}
\]
where $\pu$ is the unit vector along $\p$ and $\Trace$ denotes the
trace operation.

Equation (\ref{eq:fp}) is usually rewritten as
\begin{equation}
 \pderiv{f(\p)}{t}=-\nabla_\p[(\Dvec-\nabla_\p\Diff) f(\p)]
  +\nabla_\p[\Diff\nabla_\p f(\p)]
\label{eq:fp2}
\end{equation}
Thus, $\Diff\nabla_\p f(\p)$ can be interpreted as a flux of particles
due to diffusion (in momentum space) with a diffusion coefficient $\Diff$, and
\[
\mathbf{F}_2=\Dvec-\nabla_\p\Diff
\]
as the force acting on particles.
If the diffusion does not
change the absolute value $p=\abs{\p}$, then both interpretations of
the force are the same, $\mathbf{F}_1=\mathbf{F}_2$, or
\begin{equation}
-\nabla_\p\Diff=\frac{\pu}{\abs{\p}}\Trace{\Diff}
\label{eq:diffconserve}
\end{equation}


\section{Calculation of $D$ due to small-angle collisions with
elastic fixed centers}
Let the momentum before collision be $\p=p\unit$ and after collision
$\p'=p\unit'$. The absolute value stays the same in elastic
collisions. The change of momentum is
\[
\Delta\p=p(\unit'-\unit)
\]
Then the needed coefficients are
\begin{eqnarray*}
\Dvec&=&p\frac{\langle(\unit'-\unit)\rangle}{\Delta t} \\
\Diff&=&\frac{p^2}{2}\frac{\langle(\unit'-\unit)(\unit'-\unit)^T\rangle}
{\Delta t}
\end{eqnarray*}
If we choose coordinate system so that $\p=p\zu$, then
\[
\unit'-\unit=\zu(1-\cos\theta)-(\xu\cos\phi+\yu\sin\phi)\sin\theta
\]
Here $\theta$, $\phi$ give the direction of the scattering
$\Omega$. To find the average over $\Omega$, we take into account
axial symmetry to average over $\phi$, so that
$\langle\sin^\phi\rangle=\langle\cos^2\phi\rangle=1/2$ and
$\langle\sin\phi\cos\phi\rangle=0$. We get
\begin{eqnarray}
\langle\unit'-\unit\rangle_\Omega&=&-\langle 1-\cos\theta\rangle \unit
\label{eq:angavunit}
\\
\langle(\unit'-\unit)(\unit'-\unit)^T\rangle_\Omega
&=&\frac{1}{2}\langle\sin^2\theta\rangle\Iperp +
\langle (1-\cos\theta)^2\rangle
\unit\unit^T
\label{eq:angavunit2}
\end{eqnarray}
where we introduced the perpendicular projection operator
\[
\Iperp=\Unit-\unit\unit^T
\]
($\Unit$ is the unit operator). Note that it has an interesting property
\begin{equation}
(\Iperp)_{kl}=p\partial_k\pu_l
\label{eq:Iperpproperty}
\end{equation}
Under the assumption of small-angle collisions equations
(\ref{eq:angavunit}--\ref{eq:angavunit2}) become
\begin{eqnarray}
\langle\unit'-\unit\rangle_\Omega & \approx &
-\frac{\langle\theta^2\rangle}{2}\unit
\label{eq:smangav}
\\
\langle(\unit'-\unit)(\unit'-\unit)^T\rangle_\Omega
& \approx & \frac{\langle\theta^2\rangle}{2}\Iperp
\label{eq:smangav2}
\end{eqnarray}
where we neglected 4-th power of $\theta$ and higher.

The small-angle collisions are best described in terms of {\em
momentum transfer frequency} $\nu_m$, since the total cross-section,
and therefore the collision frequency $\nu=1/\Delta t$ can be infinite
(e.g., for Coulomb potential). It is defined as
\[
\nu_m=\frac{\langle1-\cos\theta\rangle}{\Delta t}
\approx\frac{1}{2}\frac{\langle\theta^2\rangle}{\Delta t}
\]

From (\ref{eq:smangav}--\ref{eq:smangav2}) we get
\begin{eqnarray*}
\Dvec&=&-\nu_m\p
\\
\Diff&=&D(p)p^2\Iperp
\end{eqnarray*}
where $D(p)=\nu_m/2$ can be considered as the angular diffusion coefficient.
Note that any $\Diff\sim\Iperp$, such as this one,
conserves the absolute value of the
momentum. Therefore, (\ref{eq:diffconserve}) is valid:
\[
\nabla_\p\Diff=-2\p D(p)=-\frac{\p}{p^2}\Trace\Diff=-\nu\p
\]
and the dynamic friction is $\mathbf{F}=\Dvec-\nabla_\p\Diff=0$, as
expected for elastic collisions.


\section{Nonrelativistic elastic collisions with centers of finite mass}
Consider collisions with centers of mass $M$, assuming that $f$
describes particles of mass $m$. Example of this is electron
distribution in an ambient gas of molecules. We also will assume that
molecules move with speed $\V$, which has properties
$\langle\V\rangle=0$, and $\langle V^2\rangle=3T/M$, where $T$ is an
effective molecule temperature. We assume nonrelativistic motion here,
so that to get momenta in the a differente reference frame (RF), e.g.,
the center-of-mass reference frame (CM),
we use vector addition of velocities. We will write down the FP
equation for the velocities $\vvec$ and $v=\abs{\vvec}$.

The center of mass velocity is
\[\VCM=\frac{m\vvec+M\V}{M+m}\]

In the CM RF, the initial velocity is
\[
\vCMv=\vvec-\VCM=\vCM\unit
\]
where now $\unit$ is the unit vector along $\vCMv$.
We can calculate that
\[
\vCMv=\frac{M}{M+m}(\vvec-\V)=\frac{M}{M+m}\vrelv
\]
where $\vrelv$ is the relative velocity.

After a collision, the CM velocity becomes $\vvec'=\vCMv'+\VCM$, where
\[
\vCMv'=\vCM\unit'
\]
so that the elastic scattering in CM RF frame occurs from direction
$\unit$ into direction $\unit'$.

The velocity change is the same in all RF and is equal to
\[
\Delta\vvec=\Delta\vCMv=\vCM(\unit'-\unit)
=\frac{M}{M+m}\vrel(\unit'-\unit)
\]

The averages $\langle\Delta\vvec\rangle_\Omega$ and
$\langle\Delta\vvec\Delta\vvec^T\rangle_\Omega$ over the scattering angle
(in CM RF) are easily obtained using
(\ref{eq:angavunit}--\ref{eq:angavunit2}). Then we still have to
average over $\V$.
Now, since the scattering centers are moving, the
collision frequency $\nu$ depends on $\V$, because
$\nu=\nu(\vrel)$. Therefore, when averaging over $\V$, we must
average the products $\langle\nu\Delta\vvec\rangle$ and
$\langle\nu\Delta\vvec\Delta\vvec^T\rangle$:
\begin{eqnarray}
\langle\nu\Delta\vvec\rangle&=&-\langle\nu
\langle 1-\cos\theta\rangle_\Omega\vCMv\rangle_\V
=-\langle\nu_m\vCMv\rangle_\V
\label{eq:avexact}
\\
\langle\nu\Delta\vvec\Delta\vvec^T\rangle&=&
\left\langle\nu\vCM^2
\left[
\frac{1}{2}\langle\sin^2\theta\rangle_\Omega\Iperp
+\langle (1-\cos\theta)^2\rangle_\Omega\unit\unit^T
\right]\right\rangle_\V
\label{eq:avexact2}
\end{eqnarray}
where $\Iperp=\Unit-\unit\unit^T$ is the perpendicular projection
operator and $\nu_m=\nu\langle 1-\cos\theta\rangle_\Omega$ is the
momentum transfer rate (which was also discussed earlier in this
tutorial). Note that
\[
\langle\nu(\Delta \vvec)^2\rangle
=\Trace\langle\nu\Delta\vvec\Delta\vvec^T\rangle=
2\langle\nu_m(\vrel)\vCM^2\rangle_\V
\]

\subsection{Heavy scatterers}
Equations (\ref{eq:avexact}--\ref{eq:avexact2}) are rather useless in
general case because $\Delta\vvec$ is not small and FP equation cannot
be used. Let us consider a limiting case of heavy scattering centers,
$M\gg m$. Then, although $\Delta\vvec$ is not small in general case
small changes, at least for $v=\abs{\vvec}$ we can use the FP
approach.  To get $\langle\nu\Delta v\rangle$, we have to expand
$v'=\abs{\vvec'}$ in series over the small parameter
$\alpha=\sqrt{m/M}$. From the energy equipartition theorem, $mv^2\sim
MV^2$, and therefore $V/v\sim\alpha$ and $V/v'\sim\alpha$. We use the
fact that $\abs{\vCMv'}=\abs{\vCMv}$, or
$(\vvec'-\VCM)^2=(\vvec-\VCM)^2$. Multiplying both sides by
$(M+m)/(Mv)$, we get:
\begin{eqnarray*}
\left(\left(1+\frac{m}{M}\right)\frac{\vvec'}{v}
-\left(\frac{m}{M}\vu+\frac{\V}{v}\right)\right)^2&=&
\left(\left(1+\frac{m}{M}\right)\vu
-\left(\frac{m}{M}\vu+\frac{\V}{v}\right)\right)^2\\
\left(\frac{\vvec'}{v}
+\frac{m}{M}\left(\frac{\vvec'}{v}-\vu\right)-\frac{\V}{v}\right)^2&=&
\left(\vu-\frac{\V}{v}\right)^2\\
\end{eqnarray*}
Now, carefully expand in series over $\alpha$, keeping only terms up
to $O(\alpha^2)$:
\[
\frac{v'^2}{v^2}+2\frac{m}{M}\unit'\cdot(\unit'-\unit)
-2\frac{\vvec'}{v}\cdot\frac{\V}{v}=
1-2\vu\cdot\frac{\V}{v}
\]
\[
\frac{\Delta(v^2)}{v^2}=\frac{v'^2}{v^2}-1=-2\frac{m}{M}(1-\unit'\cdot\unit)
+2\frac{(\vvec'-\vvec)}{v}\cdot\frac{\V}{v}
\]
We use $\vvec'-\vvec=\vCMv'-\vCMv=\vCM(\unit'-\unit)$ and
\begin{equation}
\vCM=\vrel+O(\alpha^2v)=v-\unit\cdot\V+O(\alpha^2v)
\label{eq:vCMexpansion}
\end{equation}
and get
\begin{equation}
\frac{\Delta(v^2)}{v^2}=-2\frac{m}{M}(1-\unit'\cdot\unit)
+2\left(
\left(\unit\cdot\frac{\V}{v}\right)^2
-\left(\unit\cdot\frac{\V}{v}\right)\left(\unit'\cdot\frac{\V}{v}\right)
\right)
-2(\unit-\unit')\cdot\frac{\V}{v}
\label{eq:Dv2}
\end{equation}
where the first two terms on the RHS are $O(\alpha^2)$ and the last is
$O(\alpha)$. Let us average over the scattering direction. We use
$\unit'\cdot\unit=\cos\theta$ and
$\langle\unit'\rangle_\Omega=\unit\langle\cos\theta\rangle_\Omega$:
\[
\frac{\langle\Delta(v^2)\rangle_\Omega}{v^2}=
2\langle 1-\cos\theta\rangle_\Omega
\left(
-\frac{m}{M}+\left(\unit\cdot\frac{\V}{v}\right)^2-\unit\cdot\frac{\V}{v}
\right)
\]

Before we average over $\V$, let us switch from
$\unit$ to $\vu$, using
\[
\unit=\left(\vu-\frac{\V}{v}\right)\abs{\vu-\frac{\V}{v}}^{-1}
\]
To order $O(\alpha)$,
\[
\unit=\vu\left(1+\vu\cdot\frac{\V}{v}\right)-\frac{\V}{v}
=\vu-\frac{\V_\perp}{v}
\]
where $\V_\perp=\V-\vu(\vu\cdot\V)$. Thus, to the same order of
$O(\alpha^2)$,
\[
\frac{\langle\Delta(v^2)\rangle_\Omega}{v^2}=
2\langle 1-\cos\theta\rangle_\Omega
\left(
-\frac{m}{M}+\left(\vu\cdot\frac{\V}{v}\right)^2-\vu\cdot\frac{\V}{v}
+\frac{\V_\perp^2}{v^2}
\right)
\]
Now it is time to average over $\V$: 
\[
\frac{\langle\nu\Delta(v^2)\rangle}{v^2}=
2\left\langle\nu_m(\vrel)
\left(
-\frac{m}{M}+\left(\vu\cdot\frac{\V}{v}\right)^2-\vu\cdot\frac{\V}{v}
+\frac{\V_\perp^2}{v^2}
\right)\right\rangle_\V
\]
We use (\ref{eq:vCMexpansion}) to expand
$\nu_m(\vrel)=\nu_m(v)-(\vu\cdot\V)\nu'_m(v)+O(\alpha^2\nu)$.
Substituting and keeping terms up to $O(\alpha^2\nu)$, we get
\begin{eqnarray*}
\frac{\langle\nu\Delta(v^2)\rangle}{v^2}&=&
-2\frac{m}{M}\nu_m
+2\nu_m\frac{\langle\V_\perp^2\rangle}{v^2}
+2\left\langle\left(\vu\cdot\frac{\V}{v}\right)^2\right\rangle_\V
(\nu_m+v\nu'_m)\\
&=&-2\frac{m}{M}\nu_m+
\frac{2}{3}\frac{\langle\V^2\rangle}{v^2}(3\nu_m+v\nu'_m)
\end{eqnarray*}
where we used $\langle\V_\perp^2\rangle=\frac{2}{3}\langle\V^2\rangle$
and $\langle(\vu\cdot\V)^2\rangle=\frac{1}{3}\langle\V^2\rangle$ since
$\V$ is isotropic. Note that the term of order $O(\alpha)$ after
averaging became $O(\alpha^2)$.

To write FP equation for the energy $\En=\frac{mv^2}{2}$, beside
$\langle\nu\Delta\En\rangle=\En\langle\nu\Delta(v^2)\rangle$ we also
need $\langle\nu(\Delta\En)^2\rangle$. Using only the last (biggest)
term in (\ref{eq:Dv2}), we get:
\[
\frac{(\Delta\En)^2}{\En^2}=\left(
(\unit-\unit')\cdot\frac{\V}{v}
\right)^2+O(\alpha^3)
\]
It is more convenient to average first over $\V$, keeping $\unit$ and
$\unit'$ constant. We can do it here, unlike the calculation for
$\Delta(v^2)$ because our error in doing this (due to the fact that
$V$ is not completely isotropic relative to $\unit$ and $\unit'$ in
the CM RF) will be of $O(\alpha^3)$. So
\[
\frac{\langle(\Delta\En)^2\rangle_\V}{\En^2}=
\frac{4}{3}
\abs{\unit-\unit'}^2
\frac{\langle\V^2\rangle}{v^2}
\]
Finally, averaging over the scattering direction,
\[
\frac{\langle\nu(\Delta\En)^2\rangle}{\En^2}=
\frac{8}{3}\nu_m\frac{\langle\V^2\rangle}{v^2}
\]
Summarizing for a thermal gas of heavy particles,
$\langle\V^2\rangle=3T/M$:
\begin{eqnarray*}
\langle\nu\Delta\En\rangle&=&
-2\frac{m}{M}\nu_m\En+
\frac{m}{M}\left(3\nu_m+2\En\frac{d\nu_m}{d\En}\right)T
\\
\langle\nu(\Delta\En)^2\rangle&=&
4\frac{m}{M}\nu_mT\En
\end{eqnarray*}
where we switched from $\nu'_m=d\nu_m/dv$ to derivative over $\En$.
The FP equation in the form (\ref{eq:fp2}) is
\[
\pderiv{n}{t}=-\pderiv{}{\En}\left(Fn-D\pderiv{n}{\En}\right)
\]
where
\begin{eqnarray*}
F(\En)&=&\langle\nu\Delta\En\rangle-
\frac{1}{2}\pderiv{\langle\nu(\Delta\En)^2\rangle}{\En}
=-2\frac{m}{M}\nu_m\left(\En-\frac{T}{2}\right)\\
D(\En)&=&\frac{1}{2}\langle\nu(\Delta\En)^2\rangle=2\frac{m}{M}\nu_mT\En
\end{eqnarray*}
Note that in an $n$-dimensional gas there is a requirement that 
\[
\frac{F(\En)}{D(\En)}=-\frac{1}{T}+\frac{n-2}{2\En}
\]
in order to get a stationary Maxwellian distribution
$n(\En)\propto\En^{(n-2)/2}e^{-\En/T}$.

\begin{thebibliography}{9}
\bibitem{nicholson}
D.~R.~Nicholson, ``Introduction to Plasma theory'', Appendix~B,
J.~Wiley and Sons, New York (1983).
\end{thebibliography}
\end{document}

