Kris's Research Notes

April 12, 2011

Detailed Balance Violations

Filed under: GaAs Simulations — Kris Reyes @ 3:26 pm

Detail balance is violated in our simulations. In this note, I describe how this occurs and examine how often these violations take place.

Recall for a continuous time Markov chain on a state space \mathcal X, we specify rates r(x, y) for the transition x \rightarrow y for states x, y \in \mathcal X. We say that \mathcal X is a semi-lattice if for each pair x,y \in \mathcal X, there is a well-defined state x\wedge y \in \mathcal X such that

\displaystyle x\wedge y = y \wedge x

If \mathcal X has such a semi-lattice structure, we can define Arrhenius-type rates of the form

\displaystyle r(x,y) = R_0 e^{-\beta(E(x) - E(x\wedge y))},

for which we are guaranteed the detailed balance property

\displaystyle \pi(x)r(x,y) = \pi(y)r(y,x),

holds with respect to Boltzmann equilibrium distribution

\displaystyle \pi(x) = \frac{1}{Z} e^{-\beta E(x)}.

The proof of this is straightforward calculation from the definition of detailed balance.

For example, consider an adatom diffusing on the surface of a single-species crystal. Below, we indicate the initial and final state of a hop x, y, respectively. We also show the standard choice for intermediate state x\wedge y, namely the state with the atom hopping removed.

If bonding energies do not depend on coordination number (e.g. simple SOS KMC without reactions) and energies are calculated via bond counting, the difference in energies E(x) - E(x\wedge y) is exactly the sum of the bonding energies involving the diffusing adatom:

Above, we have indicated those bonds that are in x but not in x \wedge y. Those atoms whose local energies do not change between x and x\wedge y are in lighter blue.

However, when energies depend on coordination number, the expression E(x) - E(x\wedge y) may no longer be obtained solely from the bonds of the diffusion adatom. Here we have indicated the atoms whose coordination number changes in the transition x \rightarrow y and in x \rightarrow x\wedge y.

Therefore, if the bond energies depend on coordination number, the energy difference E(x) - E(x\wedge y) depends on the bonds involving the diffusing adatom (green below) as well as those bonds between atoms whose coordination number changes (red and purple).

Explicitly, E(x) - E(x\wedge y) = \Gamma_g + (\Gamma_p - \Gamma_r), where \Gamma_g is the sum of the green bonds, those between the diffusing adatom and its neighbors and \Gamma_r, \Gamma_p are the sum of the red and purple bonds, respectively. We may interpret this transition as a diffusion (whose energy barrier is \Gamma_g) and a reaction (whose latent heat is given by \Gamma_p - \Gamma_r).

In the simulations, we do not use E(x) - E(x\wedge y) in our transition rates. This leads to a violation of detailed balance. Specifically, we only use those bonds between the diffusing adatom and its neighbors. That is, the correct (in the above sense, i.e. to get detailed balance and hence the Boltzmann equilibrium distribution) rates should be of the form

\displaystyle R_0 e^{\beta(E(x) - E(x\wedge y)} = R_0 e^{\beta(\Gamma_g + (\Gamma_p - \Gamma_r))}

while the rates in our simulations are

\displaystyle  R_0 e^{\beta\Gamma_g}.

That is, we get the reactions at no cost.

The energy \Gamma_p - \Gamma_r can be zero, in which case no reaction has taken place. To understand when this happens, it is important to note that the difference in energies occurs because an atom reduces its coordination number by one. For example, suppose atom A is of species a(3) in state x but is of species a(2) in state x\wedge y and suppose B is of species b(i) in both states. Consider the bond between atoms A and B as indicated below:

This bond contributes to the latent heat of transformation \Gamma_p - \Gamma_r only if \gamma(b(i), a3) - \gamma(b(i), a2) \neq 0. If we think of the energies \gamma(b(i), a(j)) as tabulated in a matrix, we see that this bond contributes to the heat of transformation if the energy at entry (i,3) differs from that one column over at entry (i,2). More generally, there is a potential for some bond to contribute to the heat of transformation whenever there are two unequal adjacent entries in the matrix \gamma(i,j).

Consider the form of the pairwise energies in our simulation

\displaystyle \gamma(G(i), A(j)) = \begin{pmatrix} \gamma_u & \gamma_u & \gamma_u & \gamma_u \\ \gamma_u & \gamma_u & \gamma_u & \gamma_p & \\ \gamma_u & \gamma_u & \gamma_p & \gamma_p & \\ \gamma_u & \gamma_p & \gamma_p & \gamma_b \end{pmatrix}

\displaystyle \gamma(G(i), G(j)) = \begin{pmatrix} \gamma_G & \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime \\ \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime & \\ \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime & \\ \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime & \gamma_G^\prime \end{pmatrix}

\displaystyle \gamma(A(i), A(j)) = \begin{pmatrix} \gamma_A & \gamma_A & \gamma_A & \gamma_A \\ \gamma_A & \gamma_A & \gamma_A & \gamma_A & \\ \gamma_A & \gamma_A & \gamma_A & \gamma_A & \\ \gamma_A & \gamma_A & \gamma_A & \gamma_A \end{pmatrix}

We can identify those places in the above matrices where adjacent entries are unequal, and hence transitions that lead to a reaction. For example, for \gamma(G(i), A(j)), the a reaction occurs when an atom crosses the indicated blue lines:

For \gamma(G(i), G(j)):

There are no such boundaries in \gamma(A(i), A(j)). The length of such boundaries is indicative of how often such reactions occur. However, not all transitions are equally likely to occur in the simulation, and so line segments must be weighted by the probability of finding such a transition in the simulation.

We may run the simulations in order to determine empirically how often we cross these reaction boundaries. That is, we can determine how often \Gamma_r - \Gamma_p is non-zero, and what its empirical average is. This was done by rerunning the homoepitaxy trials with the following parameters:

\gamma_G = 0.3 eV, \gamma_G^\prime = 0.25 eV, \gamma_A = 0.10 eV

\gamma_u = 0.30 eV, \gamma_p = 0.60 eV, \gamma_b = 1.0 eV

\mu_{As} = 1.38 eV, E_d = 1.35 eV

T = 734 K

\displaystyle \frac{R_{As}}{R_{Ga}} \in \left\{ 0, 1, 2, 4, 8, 10 \right\}

The results are tabulated below:

\frac{R_{As}}{R_{Ga}} \#\left\{\Gamma_p - \Gamma_r \neq 0 \right\} # total swaps Percentage \mathbb E \left[\Gamma_p - \Gamma_r\right] (eV)
0 10958917 41175412 26.62% 0.229039
1 4284101 49533378 8.65% 0.530630
2 4252524 49712888 8.55% 0.864026
4 3977355 49376893 8.06% 0.867867
6 4413732 49080222 8.99% 0.855852
8 3355452 47858852 7.01% 0.865336
10 3471576 48235325 7.20% 0.864649

Note: The number \#\left\{\Gamma_p - \Gamma_r\right\} includes those exchanges within the droplet and at the droplet/substrate interface in addition to those on the surface.


1 Comment »

  1. […] the previous post, we described how detailed balance can be violated for certain events. For such an event, we had […]

    Pingback by More on the new energy barriers « Kris's Research Notes — May 19, 2011 @ 6:51 am

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: