Kris's Research Notes

November 15, 2010

Rereading Krug

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 8:54 am

I am currently rereading J. Krug’s paper “Atom Mobility for the Solid-on-Solid Model.” I think I understand more of it than last time, so I will now attempt to discuss the continuum limit for the evolution of a height profile. (more…)

October 12, 2010

The H-Theorem for a Time-Reversible Markov Chain.

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 6:16 pm

For a countable set of states \mathcal X and an evolving distribution P(t; x) on \mathcal X, we may define the entropy S at time t by

\displaystyle S(t) = -\sum_{x \in \mathcal X} P(t;x) \log P(t;x).

The H-Theorem asserts that (with some assumptions) \frac{dS}{dt} \geq 0; the entropy is always non-decreasing. Calculating this quantity, we see

\displaystyle \frac{dS}{dt} = -\sum_{x \in \mathcal X} \frac{dP(t;x)}{dt}\log P(t;x) + \frac{dP(t;x)}{dt}.

Since \sum_{x \in \mathcal X} \frac{dP(t;x)}{dt} = 0 we see we have

\displaystyle \frac{dS}{dt} = -\sum_{x \in \mathcal X} \frac{dP(t;x)}{dt} \log P(t;x).

We assume P(t) describes a Continuous Time Markov Chain. That is, we assume there exists transition rates R(x,y) for distinct x, y \in \mathcal X such that P(t) evolves according to:

\displaystyle \frac{dP(t;x)}{dt} = \sum_{\stackrel{y\in \mathcal X}{y\neq x}} P(t;y)R(y,x) - P(t;x)R(x,y).

If we write

\displaystyle R(x,x) = -\sum_{\stackrel{y\in \mathcal X}{y\neq x}} R(x,y),

we may write the above as

\displaystyle \frac{dP(t;x)}{dt} = \sum_{y\in \mathcal X} P(t;y)R(y,x).

If we consider the R as an infinite dimensional |\mathcal X| \times |\mathcal X| matrix and similarly consider P(t) as an infinite dimensional vector, we may write the above compactly in matrix notation:

\displaystyle \frac{dP}{dt} = P^TR.

We also assume the chain is time reversible. That is, there exists a distribution \pi(x) on \mathcal X such that the  detailed balance property holds:

\pi(y)R(y,x) = \pi(x)R(x,y).

For example, suppose \mathcal X was finite and the rates R were symmetric. Then the corresponding chain is time-reversible and the requisite \pi is the uniform distribution on \mathcal X. It is easy to see the H-theorem holds in this context, since we have

\displaystyle \frac{dP(t;x)}{dt} = \sum_{y\neq x \in \mathcal X} R(x,y)(P(t;y) - P(t;x)),

so that

\displaystyle \frac{dS}{dt} = -\sum_{\stackrel{x,y\in \mathcal X}{y\neq x}} \log P(t;x) R(x,y)(P(t;y) - P(t;x)).

Note that in this sum, both the terms

-R(x,y)\log P(t;x)(P(t;x) - P(t;y)),


R(x,y)\log P(t;y)(P(t;x) - P(t;y)).

occur. Grouping them together, and summing over x,y, we get:

\displaystyle \frac{dS}{dt}= \frac{1}{2}\sum_{x,y \in \mathcal X}R(x,y)\left(\log P(t;x)-\log P(t;y)\right)\left(P(t;x) - P(t;y)\right).

Each term in the sum is non-negative, so that \frac{dS}{dt} \geq 0.

We wish to consider the general time reversible case. Instead of being symmetric, we have

\displaystyle R(y,x) = \frac{\pi(x)}{\pi(y)} R(x,y).


\displaystyle \frac{dP(t;x)}{dt} = \sum_{\stackrel{y\in \mathcal X}{y\neq x}}\frac{R(x,y)}{\pi(y)}\left(\pi(x)P(t;y) - \pi(y)P(t;x)\right).

so that

\displaystyle \frac{dS}{dt} = -\sum_{\stackrel{x,y\in \mathcal X}{y\neq x}} \log P(t;x)\frac{R(x,y)}{\pi(y)}\left(\pi(x)P(t;y) - \pi(y)P(t;x)\right)

and as before, we observe the terms

\displaystyle -\log P(t;x)\frac{R(x,y)}{\pi(y)}\left(\pi(y)P(t;x) - \pi(x)P(t;y)\right),


\displaystyle \log P(t;y)\frac{R(y,x)}{\pi(x)}\left(\pi(y)P(t;x) - \pi(x)P(t;y)\right) = \log P(t;y)\frac{R(x,y)}{\pi(y)}\left(\pi(y)P(t;x) - \pi(x)P(t;y)\right)

occur in the above summation. Therefore

\displaystyle \frac{dS}{dt} = \sum_{\stackrel{x,y\in \mathcal X}{y\neq  x}}\frac{R(x,y)}{\pi(y)}\left(\log P(t;x) - \log P(t;y)\right)\left(\pi(y)P(t;x) -  \pi(x)P(t;y)\right)

October 8, 2010

Marginal Distributions – Part 6

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 4:14 pm

Now we use an initial profile

h_i = \begin{cases} 0 & i < 16 \\ 32\sin\left(\frac{2\pi}{32}(i-16)\right) & 16 \leq i \leq 48 \\ 0   & i > 48,\end{cases}

which differs from this post in that the amplitude of the sine function is smaller. (more…)

Marginal Distributions – Part 5

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 4:05 pm

Now I compare marginal distributions for a system with initial profile

h_i = \begin{cases} 0 & i < 16 \\  64\sin\left(\frac{2\pi}{32}(i-16)\right) & 16 \leq i \leq 48\\ 0  & i > 48.\end{cases}

That is, like in the previous post, but with an even more compressed sine function. The motivation behind this is to provide a buffer between the boundary and the mass in the middle, which we hope will lead to a continuous profile at the boundary. (more…)

Marginal Distributions – Part 4

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 3:55 pm

I reran the trials detailed in this post, but with an initial profile of the form

h_i = \begin{cases} 0 & i < 8 \\ 64\sin\left(\frac{2\pi}{48}(i-8)\right) & 8 \leq i \leq 56 \\ 0 & i > 56.\end{cases}

That is it is  flat near the boundaries and we have a compressed sine function in the middle (so that total mass is 0). Note we still get a discontinuity between i = 0 and i=1, but it is not as pronounced as before. (more…)

October 6, 2010

Comparing Marginal Distributions – Part 3

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 4:53 pm

This is a follow up of this post.

Recall we are evolving a height profile \bar h(t) that is initially sinusoidal \bar h_i(0) = 64 \sin\left(\frac{2\pi}{64}(i-1)\right) for i = 1, \hdots, 64 and we have fixed  h_0(t) = 0 for all t. We simulate the evolution with atom hopping rates

\displaystyle r(i, i\pm1) = \frac{\Omega^\prime}{2}e^{-\beta\gamma n_i}

where n_i is the number of lateral neighbors of the top atom at site i (always counting half bonds with the wall for n_M) and \Omega^\prime = 5\times 10^7. We also set the rates r(1, 0) and r(M, M+1) to zero. We noted two things in the previous post. First, we had observed a discontinuity between h_1(t) and h_0(t) = 0 for small t.  Second, comparing the expected and empirical marginal distributions  \mathbb P_{emp}(t; h_1) and

\displaystyle \mathbb P_{leq} (t;h_1) = \frac{\sinh \theta - \sinh \beta m_1(t)}{\cosh \theta} e^{\theta|h_1(t)| + m_1(t)h_1(t)},

we observed that the two distributions did not agree for small t.

We had discussed that this was probably due to the discontinuity between h_0 and h_1.  To avoid this, we decided to use an initial height profile that was very smooth at x=0, e.g. sin^5(x). So we repeat the previous experiment with initial profile

\bar h_i(0) = 64 \sin^5\left(\frac{2\pi}{64}(i-1)\right)

Here is a movie of the simulated evolution of \bar h(t) with rates specified above and total simulated time of 0.1 seconds. The predicted average height profile was solved using Euler method with a time step of \Delta t = 10^{-9} seconds.

As we can see, we still get a discontinuity between h_0 = 0 and h_1. Here are the predicted vs. observed height profiles along with marginals at several times in the interval [0, 0.1]. (more…)

October 4, 2010

Comparing Marginal Height Profiles – Part 2

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 3:58 pm

This is a follow up to this post.  Recall, in our KMC simulations, we had originally assigned atom hopping rates

r_{i, i\pm1} =\begin{cases} \frac{1}{2}R_i & i = 2, \hdots , N-1 \\ R_i & i = 1, N. \end{cases}

where R_i = \Omega^\prime e^{-\beta \gamma n_i}, and n_i is the number of lateral bonds of the top atom at site i (and we count the atom-wall bond at site N as half a bond). Here \Omega^\prime = 5\times 10^7 and \gamma = 0.25 eV. The different cases here in r_{i, i\pm1} arise due to how the simulations sample a hop event — the code first samples a site according to rate R_i, then hops the atom to the left or right with equal probability (hence the 1/2) if i \neq 1, N and directly to the right or left if i = 1 or N, respectively (and hence no 1/2 factor).  This is incongruous to our analysis of the average height profile and its evolution.

We can rectify this in two ways. First, we can change the simulation so that r_{N, N-1} = \frac{1}{2} R_N and r_{1, 2} = \frac{1}{2} R_1.  Second, we can change the evolution of \bar h(t) by replacing expressions \frac{1}{2} R_1 with R_1 and similarly \frac{1}{2} R_N with R_N. We opt to make the first change. Then we have the usual expressions for the evolution of \bar h_k:

\displaystyle \frac{d \bar h_k}{dt} = \begin{cases}\frac{\Omega}{2}\left(-e^{\beta\mu_1} + e^{\beta\mu_2}\right) & k = 1 \\ \frac{\Omega}{2}\left(e^{\beta\mu_{i-1}} - 2e^{\beta\mu_i} e^{\beta\mu_{i+1}}\right) & k = 2, \hdots, N-1 \\ \frac{\Omega}{2}\left(e^{\beta\mu_{N-1}} - e^{\beta\mu_N}\right) & k =N \\ \end{cases},

where \Omega = \Omega^\prime e^{-\beta\gamma}.

This allows us to approximate the evolution of the average height profile. In this run, we have 4000 KMC simulations of an initially sinusoidal height profile with period 64 and amplitude 64. We measure the height profile every 0.01 second. This allows to calculate the empirical average height profile h_{emp}(t) and the empirical distributions for the height differences. We compare this to the expected height profile \bar h, which we compute using the evolution equation above and initial conditions \bar h(0.01) = \bar h_{emp}(0.01). We also compare empirical vs. predicted marginal distributions.

We plot the empirical data along with the predicted data below. Each row corresponds to t \in \left\{0.01, 0.04, 0.09, 0.14, 0.19, 0.24\right\}. The first plot in a row is the plot of \bar h and \bar h_{emp}. Subsequent plots in a row are expected and observed plots of the marginal distribution for h_i(t) - h_{i-1}(t) for i \in \left\{1, 16, 32, 48, 64\right\}. In these graphs, blue plots are data observed from the simulations while red plots are the expected average height profile and marginal distributions. We note two things. First, the average height profile predicts the simulation data very well. Second, except for the case where i=1, we get a good match in the marginal distributions as well. When t = 1, 4 or 9, the marginal distribution for h_1 is not good. Note these correspond to the case when h_1 is much larger than h_0 =0. (more…)

September 29, 2010

How the average height profile evolves over time.

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 6:47 am

Recall we defined the average height profile

\bar h_i(t) = \mathbb E_{P(t)} [h_i] = \displaystyle \sum_{h \in \mathcal X} h_i P(t;h).

If we approximated P(t) with the local equilibrium distribution

\displaystyle \pi_{leq}(t;h) = \frac{1}{\Xi} \exp\left[\beta (H(h) - \sum_{i=1}^N \mu_i h_i\right],

we showed that \bar h evolves approximately like

\displaystyle \frac{d\bar h_k}{dt} =\frac{\Omega}{2} \left(e^{\beta\mu_{k-1}} - 2 e^{\beta\mu_k} + e^{\beta \mu_{k+1}}\right),

for k=2, \hdots, N-1 and

(1) \displaystyle \frac{d\bar h_1}{dt} =\frac{\Omega}{2}  \left( - e^{\beta\mu_1} + e^{\beta  \mu_2}\right),


(2) \displaystyle \frac{d\bar h_N}{dt} = \frac{\Omega}{2} \left( e^{\beta\mu_{N-1}} - e^{\beta\mu_N}\right).

Here, \Omega is the prefactor used in the hopping rate

(3) \displaystyle r (h, h-e_i+e_{i\pm 1}) = \Omega e^{\beta[H(h) - H(h(i))]} = \Omega e^{-\beta\gamma(n_i(h)-2)},

where n_i(h) is the number of bonds the top atom at site i has. If instead we wished to use the number of lateral neighbors, \ell_i(h), then n_i(h) = \ell_i(h)+1 and the rates are given by

\displaystyle r (h, h-e_i+e_{i\pm 1}) = \Omega e^{-\beta\gamma(\ell_i(h)-1)} = \Omega e^{\beta\gamma}e^{-\beta\gamma \ell_i(h)}.

We make this point because in KMC simulations, we often define rates in terms of \ell_i(h). For example, in our current simulation we set our rates such that \Omega e^{\beta\gamma} = 5\times 10^7, hence in our analysis we must fix \Omega = e^{-\beta\gamma}5\times 10^7. With \gamma = 0.25, this means \Omega = 2.748\times 10^6. Note this is a large number.

Consider the average profile near equilibrium. Here the average profile does not change much in space (\bar h_i - \bar h_{i-1} \approx 0) or time (\frac{d\bar h}{dt} \approx 0). Then

m_i = kT \mathcal M(\bar h_i - \bar h_{i-1}) = kt \mathcal M(\epsilon).

Consider kT \mathcal M(x) near x = 0:

Observe the m_i are small even for (relatively) large difference in height profile. Then the mu_i = -(m_{i+1} - m_i)  are also small. Using the above graph as an example,  we see that near equilibrium \mu_i \in [-0.3,0.3] with high probability.  In equilibrium all the $\mu_i$ are equal and so, by examining equations (1), (2), and (3) we see that \frac{d\bar{h_i}}{dt} = 0 in this case.

Now consider the system not near equilibrium. In particular, suppose we had the following height profile:

We wish to consider how our model predicts it will evolve with respect to the equations (1), (2), and (3) above. To that end, consider the plot of e^{\beta \mu_i} where the $\mu_i$ are calculated from this average profileAs we see, the values are somewhat close together — but not close enough! That is, consider \displaystyle \frac{1}{\Omega/2} \frac{d\bar  h}{dt}, the unnormalized rates:

When we scale by \Omega/2 \approx 10^6 we see that \frac{d\bar h}{dt} is very large! This could leads to some unstable behavior if we try to evolve \frac{d\bar h}{dt} using e.g Euler’s Method.

September 24, 2010

Comparing Marginal Distributions

Filed under: Stat. Mech for SOS Model — Kris Reyes @ 2:49 pm

Recall: We wish to compare the distributions \pi_{leq}(t) and P(t), where

\pi_{leq}(t;h) = \frac{1}{\Xi} e^{-\beta H_{leq}(t;h)},


P(t;h) = \frac{1}{Z} e^{\beta H(t;h)},

for energy functions H, H_{leq} and partition functions Z, \Xi. Because of the form of H_{leq} and H, we may instead (and perhaps more naturally) write these probabilities in terms of the height differences latex g_i = h_i - h_{i-1}. We do so below and write \pi^\Delta_{leq}(t; g) for this probability, where g = (g_1, \hdots, g_n) . For example

\displaystyle \pi^\Delta_{leq}(t;g) = \frac{1}{\Xi}e^{-\beta H_{leq}(t;g)} = \frac{1}{\Xi} \exp\left[-\beta(\sum_{i=1}^N |g_i| + m_ig_i)\right]

One way to visualize this is to consider the marginals

\displaystyle \pi^{\Delta_i}_{leq} (t; g_i) = \sum_{g_1, \hdots, \hat g_i, \hdots g_N} \pi^\Delta_{leq}(t; (g_1,\hdots,g_N)),

and similarly for P^{\Delta_i}(t;g_i). The marginal distribution \pi{^\Delta_i}_{leq}(t;g_i)  may be simplified as

\displaystyle \pi^{\Delta_i}_{leq}(t; g_i) = e^{\theta|g_i| + \beta m_i(t)g_i}\frac{\cosh\theta - \cosh \beta m_i(t)}{\sinh \theta}

where \theta = \beta\gamma/2 and m_i is a function of the average height profile at time t. We compare this distribution with the empirical estimation for P^{\Delta_i}(t, g_i), which we obtained by KMC simulation.