\documentclass[12pt]{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{equations}
\usepackage{own}
%\setlength{\textwidth}{6.5in}
%\renewcommand{\baselinestretch}{1.1}
%\renewcommand{\theequation}{\thesection.\arabic{equation}}
% see page 19 of Advanced LaTeX course by Michel Goossens for nifty picture
%\voffset -1in
%\hoffset -1in
\oddsidemargin 0.55625in % Left margin on odd-numbered pages.
\evensidemargin 0in % Left margin on even-numbered pages.
\reversemarginpar % put notes in opposite margin than defaults
\textheight 8.70in
\addtolength{\textheight}{-\baselineskip}
\textwidth 5.98in
\topmargin -.28in % distance from top of page to header - 1in
\headheight 0.5in % height of box containing the header
\headsep \baselineskip % distance between head and body
%\footheight 0.5in % height of foot
\footskip 0.5in % bottom of last line to bottom of foot
\topsep=0in % get rid of extra vertical space in list modes
\partopsep=0in %
\parsep=0.25in %
\parskip 0.1 in
\title{Stochastic Modeling in Physics and Microbiology\\
MATH 6960--1 -- Spring 2013 \\
Homework 2}
\date{Due Tuesday, March 19 at 2:00 PM}
\newcommand{\grad}{\boldsymbol{\nabla}}
\newcommand{\dive}{\grad \cdot}
\renewcommand{\theenumi}{\alph{enumi}}
\newcommand{\mathi}{\mathrm{i}}
\newcommand{\expe}{\mathrm{e}}
\newcommand{\difd}{\mathrm{d}}
\newcommand{\Prob}{\,\mathrm{Prob}\,}
\newcommand{\bbE}{\mathbb{E}}
\newcommand{\Ito}{It\^{o}}
\newcommand{\pT}{p_T}
\newcommand{\pTeps}{p_T^{(\epsilon)}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\bW}{\mathbf{W}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\bxp}{\mathbf{x}^{\prime}}
\newcommand{\bxprobe}{\mathbf{x}^{(p)}}
\newcommand{\bk}{\mathbf{k}}
\newcommand{\hbuk}{\hat{\mathbf{u}}_{\bk}}
\newcommand{\hbumk}{\hat{\mathbf{u}}_{\bk}}
\newcommand{\hpk}{\hat{p}_{\bk}}
\newcommand{\hqk}{\hat{q}_{\bk}}
\newcommand{\hgk}{\hat{h}_{\bk}}
\newcommand{\bFtherm}{\mathbf{F}_T}
\newcommand{\bfthermk}{\hat{\mathbf{f}}_{T,\bk}}
\newcommand{\bxinit}{\mathbf{x}_{\mathrm{in}}}
\newcommand{\bftherm}{\mathbf{f}_T}
\newcommand{\tp}{t^{\prime}}
\newcommand{\Id}{\mathit{\mathsf{I}}}
\newcommand{\bumeas}{\mathbf{U}_{m}}
\newcommand{\lprobe}{\ell_p}
\newcommand{\lpart}{a}
\newcommand{\bVpart}{\mathbf{V}}
\newcommand{\XXtens}{\mathit{\mathsf{C}}^{(\bX \bX)}}
\newcommand{\XVtens}{\mathit{\mathsf{C}}^{(\bX \bV)}}
\begin{document}
\maketitle
This homework has 160 points plus more than 30 bonus points available but, as always, homeworks are graded out of 100 points. Full credit will generally be awarded for a solution only if it is both correctly and efficiently presented using the techniques covered in the lecture and readings, and if the reasoning is properly explained. If you used software or simulations in solving a problem, be sure to include your code, simulation results, and/or worksheets documenting your work. If you score more
than 100 points, the extra points do count toward your homework total.
\section{Brownian Tails (\textbf{105 points plus 30 bonus points})}
One shortcoming of the Langevin description of Brownian motion is that the prediction of an exponential correlation function
\begin{equation*}
\langle \bV (t) \otimes \bV(\tp) \rangle = \frac{k_B T}{m} \Id\, \expe^{- (\gamma/m) |t-\tp|}
\end{equation*}
does not agree with numerical simulations (since the 1960s) and high resolution experimental observations (since the 1980s). The theoretical explanation for the discrepancy is that while the individual particles in the surrounding fluid have a very fast time scale, their collective motion as a hydrodynamic fluid has a somewhat slower response to the motion of the particle which feeds back to affect the particle motion. We will consider a simplified cartoon model for this effect which actually does incorporate the basic concept correctly, though several physical details are too crudely considered.
Suppose we model the motion of the particle of interest as follows:
\begin{align}
\frac{\difd \bX (t)}{\difd t} &= \bu (\bX (t),t), \label{eq:Xlag} \\
\bX (t=0) & = \bxinit. \nonumber
\end{align}
where $ \bu (\bx,t) $ denotes the velocity of the fluid in which the particle moves, as a function of space and time. (This is an approximation which neglects particle inertia and shape.) We further assume the fluid is in a thermal equilibrium state governed by the time-dependent incompressible Stokes equations with random thermal force density $ \bftherm (\bx,t) $ (representing the internal transfer of momentum due to collisions of the fluid molecules):
\begin{align}
\rho \frac{\partial \bu (\bx,t)}{\partial t} &= - \mu \Delta \bu (\bx,t) - \grad p (\bx,t) + \bftherm (\bx,t),
\nonumber \\
\dive \bu (\bx,t) &= 0. \label{eq:upde}
\end{align}
(This is a fairly good fluid model.)
Here $ \mu $ is the dynamic viscosity, $ p $ is the pressure, and $ \rho $ is the density of the fluid.
Note that $ \bu $ and $p $ are random functions but are indicated using lower case letters as is common in fluid mechanics. Note that thermal equilibrium implies that the \emph{statistical average} of the pressure gradient $ \grad p $ and the velocity $ \bu (\bx,t) $ must be zero; they both have random fluctuations (which turn out to be significant only on the microscale).
For simplicity, we assume the fluid fills a cubic domain $ D $ with length $2 \pi L $ along each side, and satisfies periodic boundary conditions. So long as $ L $ is large compared with the particle size, this gives reasonable results. We then can write the various fluid fields as Fourier series:
\begin{align}
\bu (\bx,t) &= \sum_{\bk \in \BBZ^3} \hbuk (t) \expe^{\mathi \bk \dotprod \bx/L}, \label{eq:ufour} \\
p (\bx,t) &= \sum_{\bk \in \BBZ^3} \hpk (t) \expe^{\mathi \bk \dotprod \bx/L}, \\
\bftherm (t) &= \sum_{\bk \in \BBZ^3} \bfthermk (t) \expe^{\mathi \bk \dotprod \bx/L},
\end{align}
where $ \BBZ^{3} $ is the three-dimensional lattice of integers.
This reduces the partial differential equation (\ref{eq:upde}) to the constrained system of decoupled ordinary differential equations, indexed by wavenumber $ \bk \in \BBZ^3$
\begin{align}
\rho \frac{\difd \hbuk}{\difd t} & = - \mu L^{-2} |\bk|^2 \hbuk - \mathi L^{-1} \bk \hpk + \bfthermk (t), \nonumber \\
\bk \dotprod \hbuk &= 0. \label{eq:incompress}
\end{align}
Actually another constraint is that since the velocity field is real-valued, we must have
$ \hbumk (t) = \hbuk^{\ast} (t) $, where $ \ast $ denotes complex conjugate. This creates some annoying
complications in the calculation, which you are welcome to just ignore for normal credit. (If you pursue the consequences, you can earn bonus credit, but note that you'll have to fix up some of the following equations where I don't take into account the complex conjugacy relations.)
For the purposes of this problem, we take the energy of the system as
\begin{equation*}
E = \rho \frac{1}{2} \int_D |\bu |^2 \, \difd \bx = \frac{\rho L^3}{2} \sum_{\bk \in \BBZ^3} |\hbuk|^2;
\end{equation*}
the particle does not contribute to this energy because we have neglected its mass. (The fact that the energy has an infinite number of modes creates some mathematical headaches if you think about it too carefully, but in practice, there is no problem in treating the sum as if it were over a finite number of terms. That's because for any prescribed length scale of resolution $ \lambda $, the contributions of modes with $ |\bk| \gg L/\lambda $ are negligible.)
\begin{enumerate}
\item (\textbf{20 points}) Use similar arguments from class to determine an appropriate structure for the thermal force density $ \bfthermk (t) $ on each Fourier mode of the velocity field. You may assume the thermal force density acts independently on each Fourier mode; this can be shown rigorously to be a consequence of spatial homogeneity and other modeling assumptions you should be making. Specify all constants in terms of given physical parameters and temperature and Boltzmann's constant.
You will find that the thermal force density is not quite uniquely defined because part of it will simply get absorbed by the pressure term, whose job is to enforce the incompressibility constraint (\ref{eq:incompress}). Don't worry about this, you do not need to write down the most general thermal force density that works. Simply specify one which satisfies the standard principles discussed in class, and it turns out that all possible choices actually give rise to the same velocity but different pressures (which we don't care much about here).
If you are completely confused about how to solve the equations, just drop the pressure and incompressibility constraint (\ref{eq:incompress}) to still get substantial partial credit.
\item (\textbf{20 points}) Derive the mean and autocorrelation function for the Fourier modes $ \hbuk (t) $ of the velocity field in thermal equilibrium.
\item (\textbf{10 points}) Derive the mean and autocorrelation function for the fluid velocity $ \bu (\bxprobe,t) $ at a fixed probe location $ \bxprobe $, in thermal equilibrium.
\item (\textbf{5 points}) Show that the variance of $ \bu (\bxprobe,t) $ is infinite. You might worry about this unless you have been assuaged by a good deal of experience with theoretical physics.
\item (\textbf{10 points}) To make practical use of your results, the observations of the velocity field must be regularized. This happens naturally from the finite size of the probe, which implies that the measured fluid velocity is actually
\begin{equation*}
\bumeas (t) \equiv \int_{D} \bu (\bxp,t) h (\bxprobe-\bxp) \, \difd \bxp
\end{equation*}
where $ h (\bx) $ describes the measurement profile of the probe. (If it's a good probe, it looks something like a Dirac delta function smoothed out over the probe size $ \lprobe $. For simplicity we assume the probe doesn't bleed different velocity components onto each other.)
Substituting the expression (\ref{eq:ufour}) and Fourier expansion of the probe's measurement profile:
\begin{equation*}
h (\bx) = \sum_{\bk \in \BBZ^3} \hgk \expe^{i \bk \dotprod \bx/L}
\end{equation*}
and computing the integral over $ \bxp $, we obtain the following representation for the fluid velocity measured by a probe at $ \bxp $:
\begin{equation*}
\bumeas (t) \equiv L^3 \sum_{\bk \in \BBZ^3} \hgk \hbuk (t) \expe^{\mathi \bk \dotprod \bxprobe/L}.
\end{equation*}
This sum will converge because the Fourier coefficients $ \hgk $ decay rapidly for $ |\bk| \gg L/\lprobe $.
Obtain then an expression for the mean and autocorrelation function of $ \bumeas (t) $ in thermal equilibrium.
\item (\textbf{10 points plus 10 bonus points}) Plot your autocorrelation function for $ \bumeas (t) $ using some reasonable choice of probe filter. One easy one is to take
\begin{equation*}
\hgk \equiv \begin{cases}
1 & $|\bk| \leq L/\lprobe$, \\
0 & $|\bk| > L/\lprobe$.
\end{cases}
\end{equation*}
For bonus credit, you may compare the results obtained using more realistic probe filter profiles.
\item (\textbf{20 points})
According to (\ref{eq:Xlag}), the velocity of the particle is given by the fluid velocity evaluated at the current particle position $ \bu (\bX (t),t) $. Because of the compounded randomness, this function is difficult to analyze. But under certain parameter regimes (``small Kubo number") of relevance in microscale physics, the statistics of this particle velocity are well approximated by $ \bu (\bxinit,t) $
obtained by simply freezing the particle location at its initial location. Again, one would get a divergent result unless one regularizes the expression by realizing the particle, like the probe, has a finite size, so that the particle velocity $ \bVpart (t) $ can be expressed in the form:
\begin{equation*}
\bVpart (t) \equiv L^3 \sum_{\bk \in \BBZ^3} \hqk \hbuk (t) \expe^{\mathi \bk \dotprod \bxinit/L}.
\end{equation*}
where $ \hqk $ are Fourier coefficients of a shape function $ q $ describing the interaction of the particle with the fluid. ($ \hqk $ decays rapidly for $ |\hqk| \gg L/\lpart $ where $ \lpart $ is the size of the particle).
Compute then the mean and autocorrelation function of the particle velocity $ \bVpart (t) $ and from this compute the mean and covariance (not the entire autocorrelation function) for the particle position $ \bX (t) $, which under the above approximations, satisfies the differential equation
\begin{align*}
\frac{\difd \bX (t)}{\difd t} &= \bVpart (t), \\
\bX (0) &= \bxinit.
\end{align*}
Again use some reasonable choice for $ \hqk $ as for the case of the velocity measured by the probe.
\item (\textbf{10 points}) Plot your formula for the covariance of $ \bX (t) $ and from it assess whether the change in the model affects the conclusion that at long times the particle obeys a simply diffusive behavior:
\begin{equation*}
\difd \bX (t) \sim \sqrt{2 \kappa} \difd \bW (t)
\end{equation*}
for some constant $ \kappa $.
\item (\textbf{20 bonus points}) Investigate the long-time behavior of the covariance of $ \bX (t) $ analytically by considering the large domain limit $ L \rightarrow \infty$. This should allow the sum over wavenumber $ \bk $ in your expression to be viewed as a Riemann sum approximation of an integral, for which you can more easily evaluate the large $ t $ asymptotics.
\end{enumerate}
\section{Stochastically Sane Shortcut (55 points)}
Having worked out the statistics of the behavior of a particle undergoing Brownian motion according to the Langevin equation, you may wish to reread Section 1.2.2 of Gardiner's \emph{Handbook of Stochastic Methods}, which presents a quick physicist-style derivation of the diffusivity of Brownian motion. It begins by writing the Langevin equation as a second order equation for position $ \bX (t) $ as a function of time $ t $ (Gardiner 1.2.14), which in our notation and in $d$ dimensions reads:
\begin{equation}
m \frac{\difd^2 \bX}{\difd t^2} = - \gamma \frac{\difd \bX}{\difd t} + \bFtherm \label{eq:secord}
\end{equation}
where $ m $ is mass, $ \gamma $ is the friction coefficient, and $ \bFtherm (t) $ is the random thermal force. Taking half the tensor product of Eq.~(\ref{eq:secord}) with $ \bX $ on the left and on the right, Gardiner obtains in (1.2.15):
\begin{equation*}
\frac{m}{2} \frac{\difd^2 (\bX \otimes \bX)}{\difd t^2} - m \frac{\difd \bX}{\difd t} \otimes \frac{\difd \bX}{\difd t}
= - \frac{\gamma}{2} \frac{\difd ( \bX \otimes \bX)}{\difd t} + \frac{1}{2}\left(\bFtherm \otimes \bX + \bX \otimes \bFtherm\right).
\end{equation*}
This equation is now averaged. The first and third term become derivatives of the function
$ \langle \bX \otimes \bX \rangle $. The second term involves $ \langle \bV \otimes \bV \rangle
$, which can be expressed in terms of the statistics of the particle velocity $ \bV (t) = \difd \bX (t)/\difd t $ and can be separately computed (as we did in class). The term involving the product of the force and the position
\begin{equation}
\langle \bFtherm \otimes \bX + \bX \otimes \bFtherm \rangle \label{eq:fxvanish}
\end{equation}
is argued to vanish by some intuitive argument. This leaves a linear differential equation for $ \langle \bX \otimes \bX \rangle$ (Gardiner 1.2.16), which can be solved through elementary methods.
This winds up giving the correct answer, but the motivation for assuming (\ref{eq:fxvanish}) vanishes is rather subtle. As we pointed out in class, assuming the velocity and thermal force are uncorrelated gives a patently wrong answer, whereas here the assumption that the position and thermal force are uncorrelated gives a correct answer. Knowing when and when not to make this assumption is a bit murky to me -- which is why I recommend going back to the careful formulation in terms of stochastic differential equations where the basic idea of this calculation can be pursued without ambiguity. You are asked here to put this physically based derivation on such a sound footing in a manner similar to how we derived an equation for the correlation tensor of the velocity of the Brownian particle in class. You may make use of the results of that calculation but you should not insert results from our exact solution for $ \bX (t) $ from class. (That would obviate the whole point of the shortcut calculation.)
In the following procedure, if the tensor calculations confuse you, you may work in one dimension for substantial partial credit.
\begin{enumerate}
\item (\textbf{20 points}) I don't know how to write down a second order stochastic differential equation properly so work with the second order system of stochastic differential equations which characterize the Langevin equation model for Brownian motion:
\begin{align*}
\difd \bX (t) &= \bV (t) \, \difd t, \\
m \difd \bV (t) &= - \gamma \bV (t) \, \difd t + \sqrt{2 \gamma k_B T} \, \difd \bW (t),
\end{align*}
where $ k_B $ is Boltzmann's constant, $ T $ is absolute temperature, and $ \bW (t) $ is a standard vector-valued Brownian motion.
From this system, derive stochastic differential expressions for the tensors
$ \XXtens (t) \equiv \bX (t) \otimes \bX (t)$ and $ \XVtens (t) \equiv
\frac{1}{2} \left( \bX (t) \otimes \bV (t) + \bV (t) \otimes \bX (t) \right)$.
I strongly recommend using the \Ito\ stochastic calculus rules. Be very clear in justifying how you are performing the calculations.
\item (\textbf{5 points}) Average your stochastic differential expressions over the random thermal forcing to
obtain ordinary differential equations for the averages $ \langle \XXtens (t) \rangle $ and
$ \langle \XVtens (t) \rangle $ \label{eq:momeqs}
\item (\textbf{10 points}) Explain carefully how the assumed vanishing of the expression (\ref{eq:fxvanish}) is manifested in your stochastic calculations. Why is it clearly justified from the stochastic approach?
\item (\textbf{10 points}) Your equations from Eq.~(\ref{eq:momeqs}) should become a closed system of two linear ordinary differential equations for $ \langle \XXtens (t) \rangle $ and $ \langle \XVtens (t) \rangle $ once you substitute the results from class for the covariance of the particle velocity. Solve these differential equations.
\item (\textbf{10 points}) Compute the (long-time) diffusion matrix for the Brownian particle from your results, and show it agrees with the results obtained by solving directly for $ \bX (t) $ and then averaging. Your current calculation is based on averaging first, then integrating.
\end{enumerate}
\end{document}