Kris's Research Notes

October 27, 2010

Model Parameters

Filed under: GaAs Simulations — Kris Reyes @ 2:16 pm

In this post, we examine the effect of the model parameters. Many parameters depend on the definition of an atom’s local neighborhood, and throughout this note, we define the local neighborhood of an atom to be the vector of neighbor counts (n_V, n_{Ga}, n_{As}), where e.g. n_{Ga} is the number of Gallium nearest neighbors. However, we often simplify things by considering only the neighbor count for the species the atom is not. That is, we often shall only be concerned with, for example, n_{Ga} for As atoms, and vice versa. Yet there are instance where we do use the other counts, as we will detail below. Regardless, the atom species along with this counts vector determine the atom’s extended species and we shall write As(n_V, n_{Ga}, n_{As}), Ga(n_V, n_{Ga}, n_{As}) to indicate the extended species. We will also write As(n_{Ga}), Ga(n_{As}) when we only care about those counts — meaning the other counts may be anything. (more…)


October 19, 2010

Droplet Crystallization – Part 3

Filed under: GaAs Simulations — Kris Reyes @ 11:32 pm

This is a follow up of this post. We had observed that, in the movies, once an As atom came close to the droplet/crystal interface, the crystal surface became unstable. By unstable, we mean that the As atoms on the surface of the crystal would start diffusing rapidly into the droplet — causing the droplet to etch into the crystal. We had concluded that this must have been  a bug. It was indeed is a bug, one that involves next nearest neighbor swapping.

We want next-nearest neighbor swapping in general. For example, consider one atom diffusing on the surface of the crystal:

Without next-nearest neighbor swapping, the atom would have to do three swaps with nearest neighbors in order to move one position to the right. This involves breaking several bonds. This probably doesn’t happen in real life.

Having next-nearest neighbor swaps contributes to the bug. Consider the following situation where a As atom (green) diffuses through a Ga droplet and approaches the droplet/crystal interface:

Consider the darker Ga^{(1)} atom (recall Ga^{(1)} means the Gallium atom has one Arsenic neighbor). The diffusion coefficient assigned to the indicated swap, \alpha( Ga^{(1)}, As^{(4)}) = \epsilon is small. This is how the surface becomes unstable.

To fix this, we assign a separate diffusion coefficient for nearest neighbor swaps and next-nearest neighbor swaps. Then we assign (for now) all next-nearest neighbor swaps  to \infty. Here is a movie if we set \epsilon = 0.11. The temperature was set at T=498 K, \mu_{As} = 0.35. We deposit Arsenic at the rate of r_{\downarrow As} = 1 monolayer/second for 15 seconds. (Here the Gallium droplet is blue, deposited Arsenic is purple, Gallium and Arsenic in the crystal are red and green, respectively):


Note, there is no crystal nucleation within the droplet, not even at the droplet/crystal interface (and hence, the “hills” of GaAs form away from the droplet). This is easy to explain. Consider the case where an As atom has diffused to the droplet/crystal interface:

To the atom, it is surrounded by Gallium atoms, and hence it is in an identical situation to diffusing in the droplet. In other word, the As atom does not know it is at the droplet/crystal interface, it only knows about its nearest neighbors. Therefore, it is likely to swap with either of the two Gallium neighbors above it.

We can attempt to address this by increasing the Ga^{(2)}-As^{(4)} and Ga^{(3)}-As^{(4)} bonds strengths. We can make them both 1.0 eV (the same as Ga^{(4)}-As^{(4)}). Here are movies where I vary \epsilon \in \left\{ 0.11, \hdots, 0.20\right\} (organized so that \epsilon is increasing from left to right, top to bottom). I simulated 10 seconds for this run:

I believe these trials show that GaAs is nucleating within the droplet. Compare the distribution of As within the droplet in the previous case (before we adjusted the bond strengths) to the As distribution in these movies. Before, the As seemed to be uniformly distributed within the droplet, whereas in these movies, the As is concentrated in the bottom.

Next, I simulate a scenario where we anneal for 10 seconds prior to any As depositions. This will allow the droplet to come to its equilibrium shape (with respect to our simulations). I then deposited As for 10 additional seconds. Here is the case where \epsilon = 0.15:

Here are is a run where the droplet size is increased (with and without annealing):

Note, we do not get the distribution of GaAs crystal inside the droplet that we expected in this post. That is, we had previously calculated that there should be more As atoms diffusing to the edges of the droplet than in the middle, and we don’t see this in the simulations.

To see why this happens, we slow down the movie. Here is one second of simulation time, sampled every 0.01 seconds.

We see that the As atoms on the droplet/crystal interface still diffuse readily on that interface. This would explain why we see a uniform distribution of GaAs nucleating at the inteface.  In fact, this may not be a bad thing (but I really don’t know).

October 15, 2010

Droplet Crystallization Part 2

Filed under: GaAs Simulations — Kris Reyes @ 6:48 pm

Recall in the previous post, we simulated Ga droplet crystallization. In the simulation, we formed Ga droplets, and then turned on a flux of As. The droplets immediately formed shells of GaAs, which we do not observe in experiments. We had concluded in the previous meeting that it would be interesting to see what happens if we allow As atoms to diffuse through the Ga droplet at a fast rate.

First, I refined the notions of “extended species”. We define the extended species Ga^{(i)} to be a Ga has exactly i neighbors of type As, and similarly for As^{(i)}. In fact, the code has been modified to handle species of the form Ga^n, where n = (n_V, n_{Ga}, n_{As}) is a vector of counts, e.g. n_{Ga} is the number of Gallium neighbors. But as a first step, we forget about the other counts besides n_{As} so that, for example Ga^{(2,1,1)} is more or less identical to  Ga^{(1,2,1)}. That is, we obtain the extended type definition above from this more general setting. For now, I used the same bonding energies after we map to the coarser definitions for extended type.

Next, we define the swapping rate of two atoms A, B as

r_{swap}(A, B) = \Omega \exp\left[ \alpha(A, B) ( E(A) + E(B) )\right].

where as before \Omega is a constant prefactor, E(A), E(B) are the local energies of A and B, respectively. The diffusion coefficient \alpha(A, B), is given by:

  • \alpha(A,B) = 1 if one, but not both of the atoms is a vacuum (surface diffusion);
  • \alpha(A,B) = \epsilon if one of the atoms is Ga^{(1)} and the other is As (fast diffusion through droplet);
  • \alpha(A,B) = \infty otherwise.

I varied \epsilon \in \left\{0.5, 0.8, \hdots, 0.15, 1 \right\}. I started with a hemispherical droplet of radius 32 lattice sites. I fixed T = 498, \mu_{As} = 0.35. I used deposition rates r_{\downarrow As} = 1 monolayer/second. and r_{\downarrow Ga} = 0. For each \epsilon I ran two trials. First I simulated 0.1 seconds, measuring the configuration every 0.001 seconds.  Second I simulated 5 seconds measuring every 0.05 seconds.


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)

Diffusion through a Droplet

Filed under: GaAs Simulations — Kris Reyes @ 6:14 pm

Recall we wish to simulate the crystallization of Ga droplets with the presence of As flux. Let us simplify a bit, and suppose the Ga droplet is hemispherical of radius R. Let’s describe the concentration of  As in the droplet by a function u(r, \theta), where we have used spherical coordinates to describe the location inside the hemisphere (so that our domain is r \in [0, R], \theta \in [0, \pi]).

We wish to solve Laplace’s equation (WHY??):

\displaystyle \nabla^2 u =\frac{1}{r}\frac{\partial}{\partial r}\left( r \frac{\partial u}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u}{\partial \theta ^2} = 0

subject to boundary conditions

  • \displaystyle \frac{\partial u}{\partial r}\big|_{r = R} = c, i.e. constant flux on the surface of the droplet;
  • u(r, 0) = u(r,\pi)= 0, i.e. no concentration on the bottom of the droplet.

To solve this, we assume a separable solution of the form:

u(r, \theta) = F(r)G(\theta).

Applying the Laplacian we get

\displaystyle \nabla^2\left[F(r)G(\theta)\right] = \frac{1}{r} \frac{dF}{dr} G + \frac{d^2F}{dr^2}G + \frac{1}{r^2}F\frac{d^2G}{d\theta^2} = 0.

Separating variables, we get:

\displaystyle \frac{r^2}{F}\left(\frac{1}{r}\frac{dF}{dR} + \frac{d^2F}{dr^2}\right) = \lambda = -\frac{1}{G} \frac{d^2G}{d\theta^2},

for some constant \lambda. This gives us our ODEs:

\displaystyle \lambda F - r\frac{dF}{dr} - r^2\frac{d^2F}{dr^2} = 0.

\displaystyle \lambda G + \frac{d^2G}{d\theta^2} = 0.

We consider the second equation first, which is an eigenvalue problem.

If \lambda > 0, we have

G(\theta) = c_1 \cos(\sqrt{\lambda} \theta) + c_2 \sin(\sqrt{\lambda}\theta)

and given the boundary conditions G(0) = G(\pi) = 0, we see that \lambda = n^2 for n = \mathbb Z^+ and the corresponding eigenvector is of the form:

G_n(\theta) = c_n \sin(n \theta).

Now if \lambda < 0, we have a general solution of the form

G(\theta) = b_1\cosh(\sqrt{|\lambda|}\theta) + b_2\sinh(\sqrt{|\lambda|}\theta)

However, only the trivial solution satisfies the boundary conditions.

Lastly if \lambda = 0, we have G(\theta) = a_1 + a_2\theta, but once again only the trivial solution satisfies boundary conditions.

Therefore the only eigenvalues of the above operator are of the form \lambda = n^2 for n \in \mathbb Z^+ with corresponding eigenvectors G_n(\theta) = c_n \sin(n\theta).

Then we have corresponding solutions to the first ODE: F_n(r) = a_n r^p_n, where p satisfies

n^2 r^{p_n} - p_ nr^{p_n} - p_n(p_n-1)r^{p_n} = 0,

which is what is obtained if we plug in this form of F into the above ODE. That is, p_n = \pm n and F(r) = a_nr^{\pm n}. We require F(r) is bounded near r = 0, and so we ignore solutions F(r) = a_nr^{-n}.

Therefore we have a general solution of the form (incorporating all undetermined coefficients)

\displaystyle u(r,\theta) = \sum_{n=1}^\infty A_n r^{n}\sin(n\theta).

We now apply the final boundary condition \frac{\partial u}{\partial r}\big|_{r = R} = c:

\displaystyle \frac{\partial u}{\partial r}\big|_{r = R} = \sum_{n=1}^\infty nA_n R^{n-1} \sin(n\theta) = c.

By Fourier Inversion, we can solve

\displaystyle n A_n R^{n-1} = \frac{c}{\pi} \int_0^\pi \sin(n\theta) d\theta = \frac{c}{n\pi} \left(1-\cos(n\pi)\right)

and so

\displaystyle A_n = \frac{c}{\pi}\frac{(1-\cos(n\pi))}{n^2R^{n-1}}.

Plugging this in to the expression for u(r, \theta) we get:

\displaystyle u(r, \theta) = \frac{c}{\pi} \sum_{n=1}^\infty \frac{1-\cos(n\pi)}{n^2R^{n-1}} r^{n}\sin(n\theta).

We may compute the flux through the bottom of the crystal by calculating:

\displaystyle \frac{\partial u}{\partial y} = \frac{\partial u}{\partial r}\frac{\partial r}{\partial y} + \frac{\partial u}{\partial \theta}\frac{\partial \theta}{\partial y}.

That is:

\displaystyle \frac{\partial u}{\partial y} = \frac{c}{\pi} \sum_{n=1}^\infty \frac{1-\cos(n\pi)}{n R^{n-1}} \left[\sin\theta \sin(n\theta) r^{n-1} + \frac{\cos(n\theta)}{\cos\theta} r^{n-1} - \frac{\cos(n\theta)\sin^2\theta}{\cos\theta}r^{n-1}\right]

and when so \theta = 0, i.e. when y=0, x \geq 0, we get:

\displaystyle \frac{\partial u}{\partial y}\bigg|_{\theta = 0} = \frac{c}{\pi} \sum_{n=1}^\infty \frac{1-\cos(n\pi)}{n R^{n-1}} r^{n-1} = \frac{2c}{\pi}\sum_{k=1}^\infty \frac{1}{2k+1}\left(\frac{r}{R}\right)^{2k}.

I don’t know of a closed form for this, but we may plot truncated sums (for e.g. n \leq 10000) and examine the flux vs r. We may also vary the flux c. Here are the plots of \frac{\partial u}{\partial y}\big|_{\theta = 0} for R = 1, 0 \leq r \leq R and c varying from 0.1 to 1.0.



Note this is what we want: the amount of material being deposited is greater near the edge of the droplet than the middle.

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…)

Ga Droplet crystallization

Filed under: GaAs Simulations — Kris Reyes @ 2:41 pm

Just for fun and to see how our current code would handle it, I tried running the droplet crystallization experiment. Denis provided me details of the experiment, and I wanted to simulate a portion of it. Specifically:

  1. I start out with 10 layers of GaAs Crystal and fix temperature to 498 K.
  2. I turn on a Ga Flux of 0.1 monolayers Ga/second for 37 seconds.
  3. I turn off the Ga Flux and let the crystal anneal for 60 seconds.
  4. I turn on an As Flux of 3 monolayers As/second for 90 seconds.

The extended species definitions are the same as the ones detailed in this post. The pairwise energies are also the same, except I changed the Ga^{(0)}-Ga^{(0)} bond strength to 0.28 eV (I thought this would give us shorter, wider droplets — it does not). I set the desorption potentials \mu_{Ga} = 3 (effectively \infty) and \mu_{As} = 0.35. I also had to implement bulk diffusion. I used a very simple (and so incorrect) model for the rates of bulk diffusion. The rate of swapping atoms A, B when both A, B are not vacuum and are of different species was given by

r_{swap}(A,B) = \Omega \exp[\alpha(E(A) + E(B))],

where \alpha is a constant and E(A), E(B) are the local bond energies at atom A, B, respectively. In this run, I set \alpha = 0.25.

Here are the movies sped up twice and four times as fast as real time. Both movies are samples of the system every 0.4 seconds:

  • (2\times real time): Movie (25 MB)
  • (4\times real time): Movie (25 MB)

I actually sampled the system every 0.1 seconds, but movies at this full resolution are huge — at least for animated gifs. I am looking in to making actual movies instead.

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…)

Older Posts »