\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 6791--1 -- Spring 2013 \\
Homework 3}
\date{Due Friday, April 26 at 2:00 PM}
\newcommand{\oh}{\frac{1}{2}}
\newcommand{\grad}{\boldsymbol{\nabla}}
\newcommand{\dive}{\grad \cdot}
\newcommand{\divev}{\grad_{\bv} \cdot}
\newcommand{\divex}{\grad_{\bx} \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{\bv}{\mathbf{v}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\br}{\mathbf{r}}
\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{g}_{\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)}}
\newcommand{\lF}{\ell_F}
\newcommand{\Fmag}{\bar{F}}
\newcommand{\bF}{\mathbf{F}}
\newcommand{\pxv}{p_{\bX,\bV}}
\newcommand{\px}{p_{\bX}}
\newcommand{\kB}{k_B}
\newcommand{\lkk}{\ell_{KK}}
\newcommand{\tkk}{\tau_{KK}}
\newcommand{\lobs}{\ell_{o}}
\newcommand{\tobs}{\tau_{o}}
\newcommand{\bxua}{\bx^{(1)}}
\newcommand{\bxub}{\bx^{(2)}}
\newcommand{\bXua}{\bX^{(1)}}
\newcommand{\bXub}{\bX^{(2)}}
\newcommand{\bVua}{\bV^{(1)}}
\newcommand{\bVub}{\bV^{(2)}}
\newcommand{\ma}{m^{(1)}}
\newcommand{\mb}{m^{(2)}}
\newcommand{\Gmat}{\mathit{\mathsf{\Gamma}}}
\newcommand{\Gaa}{\Gmat^{(11)}}
\newcommand{\Gab}{\Gmat^{(12)}}
\newcommand{\Gba}{\Gmat^{(21)}}
\newcommand{\Gbb}{\Gmat^{(22)}}
\newcommand{\Gij}{\Gmat^{(ij)}}
\newcommand{\sigmat}{\mathit{\mathsf{S}}}
\newcommand{\sigaa}{\sigmat^{(11)}}
\newcommand{\sigab}{\sigmat^{(12)}}
\newcommand{\sigba}{\sigmat^{(21)}}
\newcommand{\sigbb}{\sigmat^{(22)}}
\newcommand{\sigij}{\sigmat^{(ij)}}
\newcommand{\bWua}{\bW^{(1)}}
\newcommand{\bWub}{\bW^{(2)}}
\newcommand{\rada}{a_1}
\newcommand{\radb}{a_2}
\newcommand{\Dlong}{\mathit{\mathsf{D}}}
\begin{document}
\maketitle
This homework has 165 points plus more than 155 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{Short-Time Analysis of Brownian Pinball Wizard (35 points plus 40 bonus points)}
The discrete model studied in the first homework can be easily reinterpreted as a continuous-time model, by simply keeping track of the particle position between collisions as well as immediately after collisions. Recall
that the particle moves always at a constant speed $ v $ but in a random direction. Moreover, the particle moves along a given direction for a random amount of time $ T $ according to a prescribed probability density function $ \pT (t) $. At the end of this time, the particle moves in a new random direction, independent of all previous motion, for another random amount of time given by the PDF $ \pT (t) $, and so on. Each time a random direction is chosen, it is to be uniformly distributed.
In the first homework you showed that the particle behaves diffusively at large time scales, meaning that the covariance matrix of the change in particle position over a long time interval $ \tau $ scales linearly with time:
\begin{equation}
\langle (\bX (t+ \tau) - \bX (t)) \otimes (\bX (t+\tau) - \bX (t) \rangle \sim 2 \Dlong \tau \text{ for } \tau \rightarrow \infty,
\label{eq:difflong}
\end{equation}
with (global) diffusivity matrix $ \Dlong $ computed in the first homework.
Here you are asked to consider its behavior at short to intermediate time scales.
\begin{enumerate}
\item (\textbf{15 points}) Compute a leading order approximation for the behavior of the covariance matrix for the change in particle position over a very short time interval $ \tau $. For what range of small $ \tau $ values should your result be expected to be fairly accurate? \label{part:short}
\item (\textbf{15 points}) Through taking appropriate statistics of simulations of continuous-time trajectories of this model, plot both diagonal and off-diagonal components of the covariance matrix for the change in particle position over a wide range of choices of time interval $ \tau $. Be sure to indicate an appropriate measure of statistical sampling error in your numerical results. \label{part:num}
\item (\textbf{5 points}) Compare your analytical result from part~\ref{part:short} with your numerical statistics in part~\ref{part:num}.
\item (\textbf{40 bonus points}) Develop an analytical expression for the covariance matrix for the change in particle position over an arbitrary time interval $ \tau $. To do this properly, you probably need to know or learn renewal theory for stochastic processes. For 20 bonus points, it's sufficient to make a decent hand-waving argument to obtain your analytical expression, so long as it is consistent with both Eq.~\eqref{eq:difflong} and your result from part~\ref{part:short}. Either way you derive the theoretical expression, compare it with your numerical results from part~\ref{part:num}.
\end{enumerate}
\section{Stochastic Mode Reduction, Free of Force (60 points plus 55 bonus points)}
In class we derived the Smoluchowski equation for Brownian motion (describing dynamics only of the position variable) from the Klein-Kramers equation (describing dynamics jointly of position and velocity) provided a natural time scale for position dynamics is long compared to a natural time scale for the velocity dynamics. The time scale for the position dynamics involved the length scale of variation of the external force. This raises the question of how this reduction procedure should be conducted when no external force is present. One might say that it should just be a special case of the reduction with external force, with some arbitrary large choice for the force length scale $ \lF $ and the amplitude of the external force set equal to zero: $ \Fmag =0 $. While this is technically correct, it is inelegant and suggests that this procedure is not really describing the relationship between the force-free Klein-Kramers equation
\begin{equation}
\frac{\partial \pxv (\bx,\bv;t)}{\partial t} = \divev \left(\frac{\gamma}{m} \bv \pxv (\bx,\bv;t)\right)
- \divex \left(\bv \pxv(\bx,\bv;t)\right) + \frac{\gamma \kB T}{m^2}\Delta_\bv \pxv (\bx,\bv;t)
\label{eq:kk}
\end{equation}
and the force-free Smoluchowski equation
\begin{equation*}
\frac{\partial \px (\bx;t)}{\partial t} = \frac{\kB T}{\gamma} \Delta_\bx \px (\bx;t)
\end{equation*}
in the most complete or natural way. Therefore, you are asked to begin again from scratch in this problem, developing a stochastic mode reduction procedure which is well suited to the force-free case.
Note that the governing parameters are the friction coefficient $ \gamma $, the particle mass $ m $, and the product of Boltzmann's constant and temperature $ \kB T $. As in class, we ignore parameters describing the initial data because those will become irrelevant after a sufficiently long time (or if we begin observations in a statistically stationary state.)
Note that the Klein-Kramers equation without force is (with some effort) exactly solvable, but you should not use this exact solution in what follows because the procedure described below can be generalized to non-solvable equations.
\begin{enumerate}
\item (\textbf{10 points}) Show that (up to multiplication by a constant) only one length scale $ \lkk $ and one time scale $\tkk $ can be formed by combinations of the governing parameters, and state what this length scale and time scale are.
\item (\textbf{10 points}) Nondimensionalize the Klein-Kramers equation with respect to the length scale $\lkk $ and time scale $\tkk $ you obtained in the previous part.
\item (\textbf{10 points}) You should have obtained a completely universal equation without any nondimensional parameters, meaning that for the force-free case, the Langevin model has identical qualitative behavior for any choice of governing parameters -- only the time scale and length scale differ. No real way presents itself to ask a question about how the equation behaves in the limit in which two time scales are well separated because the problem only has one time scale. To derive the connection with the force-free Smoluchowski equation, then, we instead introduce an observation time scale $ \tobs $ and observation length scale $ \lobs $ which we can think of as describing the resolution of our experimental apparatus.
The intuition is that the Smoluchowski equation should result when the Brownian particle is observed on length and time scales large compared with the fundamental length and time scales of the Brownian particle.
To make this mathematical, we define $ \mu = \tkk/\tobs $ and $ \epsilon = \lkk/\lobs $ as the ratios of the system's natural time and length scales to those of the observer, and will consider both of these ratios to be small. Handling two small parameters simultaneously is generally very difficult, so you should begin by considering a distinguished limit $ \epsilon \ll 1 $ and $ \mu = \epsilon^{2} $. The motivation for choosing $ \mu = \epsilon^{2} $ arises from an expected diffusive behavior where the observationally interesting length scale should scale as the square root of the observational time scale. Without this intuition, one would proceed somewhat by trial and error (see Part~\ref{part:arbpow}) below.
Express your Klein-Kramers equation now with respect to the observational length and time scales. (You can do this either by nondimensionalizing the original equation (\ref{eq:kk}) from scratch with respect to the observational length and time scales, or better, by simply rescaling your nondimensionalized equation by using the relationship between the length and time scales.) Note that you should not rescale the velocity variable along with the length and time variables; it should still be measured in terms of the natural velocity scale $ \lkk/\tkk $. Because the velocity is not being directly observed, a ``velocity observation scale" doesn't make sense.
\item (\textbf{30 points}) Develop now a perturbation expansion for the probability density with respect to the small parameter $ \epsilon $ to determine its leading order behavior, and show that it is consistently described by the (force-free) Smoluchowski equation. That is, show that in the absence of an external force, the Wiener process model is a good coarse-grained approximation of the more detailed Langevin model whenever the particle is observed on a coarse enough space and time scale. Note that this is a more general setting than that for which the coarse-graining works in the presence of the external force (a certain intrinsic time scale ratio had to be large; it's not enough to simply observe the system at long times and large scales).
\item (\textbf{15 bonus points}) \label{part:arbpow} Suppose you had instead chosen $ \mu = \epsilon^{\alpha} $ for some arbitrary choice of $ \alpha $. How would your calculation differ? Is there a reason why $ \alpha = 2 $ is somehow the ``best" choice?
\item (\textbf{20 bonus points}) Consider now the case in which a constant external force $ \bF $ (such as gravity or an external large-scale electric field) is applied, so that the Klein-Kramers equation becomes
\begin{align*}
\frac{\partial \pxv (\bx,\bv;t)}{\partial t} &= \divev \left[
\left(\frac{\gamma}{m} \bv - \frac{\bF}{m}\right)\pxv (\bx,\bv;t)\right] \nonumber \\
& \quad \qquad - \divex \left(\bv \pxv(\bx,\bv;t)\right) + \frac{\gamma \kB T}{m^2}\Delta_\bv \pxv (\bx,\bv;t)
\end{align*}
and the Smoluchowski equation becomes
\begin{equation*}
\frac{\partial \px (\bx;t)}{\partial t} = \divex \left( \frac{\bF}{\gamma} \px (\bx;t)\right) +
\frac{\kB T}{\gamma} \Delta_\bx \px (\bx;t).
\end{equation*}
Again if we try to use the approach from class, we'd have to introduce a fictitious length scale $ \lF $, which is rather unnatural. Develop a more appropriate relationship between the Klein-Kramers equation and Smoluchowski equation in the presence of a constant external force which only refers to the governing parameters actually in the problem, as well as possibly the observational length and time scale.
\item (\textbf{20 bonus points}) Finally consider the case in which the force is derived from an external potential $ \bF (\bx) = - \grad \phi (\bx) $ which has a pure power-law form $ \phi (\bx) = k |\bx|^\beta $ for some real constants $ k $ and $ \beta $ with appropriate physical dimensions. Reduce the Klein-Kramers equation to the Smoluchowski equation without introducing a fictitious force length scale $ \lF $.
\end{enumerate}
\section{Thermal Forcing of Interacting Particles (70 points plus more than 60 bonus points)}
Consider including the effects of Brownian motion to a simulation of a collection of particles (such as in colloidal suspensions and various molecular biology simulations with structures represented in terms of interacting particles). If they are close enough to interact, the correct thermal forcing scheme is not simply adding independent Brownian motion to each particle. To see this, we consider what thermal forcing of the governing equations is consistent with the laws of equilibrium statistical mechanics.
The main ideas can be seen from the case of two particles, with positions $ \bXua $, $ \bXub $ and velocities $ \bVua $, $ \bVub $ and masses $ \ma $, $ \mb $ respectively. The equations of motion
read:
\begin{align}
\difd \bXua &= \bVua, \nonumber \\
\ma \difd \bVua &= - \Gaa (\bXua-\bXub) \cdot \bVua \, \difd t - \Gab (\bXua - \bXub) \cdot \bVub \, \difd t
\nonumber \\
& \qquad \qquad - \grad_{\bXua} \phi (\bXua, \bXub) \, \difd t \nonumber \\
& \qquad \qquad + \sigaa (\bXua,\bXub) \, \difd \bWua + \sigab (\bXua,\bXub) \,\difd \bWub, \nonumber \\
\difd \bXub &= \bVub \, \difd t, \nonumber \\
\mb \difd \bVub &= - \Gba (\bXua-\bXub) \cdot \bVua \, \difd t - \Gbb (\bXua - \bXub) \cdot \bVub \, \difd t
\nonumber \\
& \qquad \qquad - \grad_{\bXub} \phi (\bXua, \bXub) \, \difd t \nonumber \\
& \qquad \qquad + \sigba (\bXua,\bXub) \, \difd \bWua + \sigbb (\bXua,\bXub) \,\difd \bWub,
\label{eq:couppart}
\end{align}
and the total energy of the system may be represented as
\begin{equation}
E (\bXua,\bXub) = \oh \ma |\bVua|^2 + \oh \mb |\bVub|^2 + \phi (\bXua,\bXub).
\end{equation}
Here $ \phi (\bxua,\bxub) $ denotes the specified potential energy as a function of the positions of the two particles, and $ \Gij (\br) $ are specified friction tensors (matrix functions) that generalize the simple friction coefficient for spherical particles. Note in particular that the friction tensors depend on the vector $ \bXua - \bXub $ describing the current separation of the particles. Also, the motion of one particle creates drag on the other because the first particle moves the fluid which in turn pushes the other particle. Because of the coupling of the motion of the particles, both through the conservative forces described by the potential energy and the frictional force, we allow the thermal forces acting on the particles to be coupled in principle. The way we do this is to generate two independent three-dimensional Wiener processes
$ \bWua (t) $ and $ \bWub (t) $, and allow these to act on both particles with suitable weightings and projections described by the tensors $ \sigij (\bXua,\bXub) $. Note that the thermal forcing is also allowed to depend on the current particle positions (so we must stipulate that the interpretation of the equations (\ref{eq:couppart}) is in the \Ito\ sense).
\begin{enumerate}
\item (\textbf{25 points}) Determine an appropriate choice of tensors $ \sigij (\bXua,\bXub) $ that render the thermal forcing scheme above consistent with the laws of statistical mechanics. How is it related to the friction tensor and potential energy and other parameters in the problem? For partial credit (15 points), work out this calculation for the case in which each particle is treated as moving in one dimension (so all the vectors in Eq.~(\ref{eq:couppart}) become scalars).
\item (\textbf{30 points plus 20 bonus points}) Simulate the system (\ref{eq:couppart}) in 3 dimensions using the following approximation to the friction tensors for the case of two spherical particles, with radii $ \rada $ and $ \radb$:
\begin{equation*}
\left[\begin{matrix}
\Gaa (\br) & \Gab (\br) \cr
\Gba (\br) & \Gbb (\br)
\end{matrix}\right]
= \left[\begin{matrix}
\frac{1}{6 \pi \eta \rada} \Id & \frac{1}{8 \pi \eta \max(r,\rada + \radb)}
\left(\Id + \frac{\br \otimes \br}{r^2}\right), \\
\frac{1}{8 \pi \eta \max(r,\rada + \radb)}
\left(\Id + \frac{\br \otimes \br}{r^2}\right) & \frac{1}{6 \pi \eta \radb} \Id
\end{matrix}\right]^{-1};
\end{equation*}
where $ r = |\br| $, $\eta $ is the dynamic viscosity of the fluid,
all blocks here are $ 3 \times 3 $ matrices, and $ \Id $ is the identity matrix. Sorry, this is the simplest reasonable approximation I know that involves any nontrivial frictional coupling of the particles. It's decent for $ r \gg \rada, \radb $ and not so good for $ r \lesssim \rada,\radb $. Choose your own values for the remaining parameters in the problem, including the potential energy $ \phi $. For normal credit, any reasonable nontrivial choice is fine. For bonus credit, choose some realistic potential energy and physically appropriate values for parameters. Plot some samples of pairs of trajectories.
\item (\textbf{15 points plus bonus}) Think of at least one interesting statistic involving the dynamics of the pair of particles that really involves both particles, and use your simulation code to estimate it. Your statistic could be a function of time and/or particle separation.
\item (\textbf{40 bonus points}) Extend the stochastic mode reduction procedure from class to obtain a reduced system from Eq.~(\ref{eq:couppart}) which only involves dynamics of the particle positions $ \bXua $ and $ \bXub $,and coarse-grains over their velocity values. Again, this will require the formulation of the limit involving a nondimensional parameter which, roughly speaking, says the time scale of the position dynamics is large compared to the time scale of the velocity (momentum) dynamics. Set up your nondimensionalization in a way that it is reasonable to assume all nondimensional parameters besides the key one can be assumed to be order unity. For somewhat reduced bonus credit (30 bonus points), you may assume both particles are identical. Also simply assume that the separation between the two particles always remains within a factor of 10 or so of the sum of their radii -- that is, don't worry about nonuniformities that emerge if the particles become too close or too far apart.
\end{enumerate}
\end{document}