You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: |
    The pyramid-to-dome transition in Ge$_{x}$Si$_{1-x}$ on Si(100) initiated by step formation on pyramidal quantum dots is atomistically simulated using a multistate lattice model in two-dimensions incorporating effective surface reconstructions. Under quasi-equilibrium growth conditions associated with low deposition rates, the transition occurs at island size $n_c$ following $\sqrt{n_c}
      \sim x^{-1.69}$ independent of temperature and deposition rate. The shape transition is found to be an activated process. Results are explained by a theory based on simple forms of facet energies and elastic energies estimated using a shallow island approximation. An asymptotic scaling relation $n_c^{1/d} \sim
      x^{-2}$ for $x\rightarrow 0$ applicable to $d$ = 2 or 3 dimensions is derived. The shape transition energy barrier can be dominated by the interface energy between steep and shallow facets.
author:
- 'Chi-Hang Lam'
bibliography:
- 'island.bib'
title: Kinetic Monte Carlo simulation of shape transition of strained quantum dots
---

introduction
============

The self assembly of quantum dots in heteroepitaxy exhibits very interesting physics and has possible application to device fabrication [@Politi2000; @Shchukin2003; @Berbezier2009]. Growth of Ge or GeSi alloy nanostructures on Si(100) is the prototype example most widely studied. Under typical deposition conditions, pyramidal islands bounded by shallow (105) facets form spontaneously on a wetting layer. They can then undergo transitions into multi-faceted dome islands dominated by much steeper (113) facets and bounded also by other facets [@Ross1999]. This shape transition gives rise to a bimodal island size distribution with enhanced dome size uniformity [@Ross1998]. An atomic pathway based on step bunching on the pyramids has been identified [@Montalenti2004; @Montalenti2007].

In this work, we report atomistic dynamic simulations of the pyramid-to-dome transition using a fast kinetic Monte Carlo (KMC) approach based on a multistate solid-on-solid model in two-dimensions (2D). Extensive simulations under a wide range of conditions are performed and a simple analytical description is presented. A scaling relation for the transition island size generalizable to three dimensions (3D) is investigated.

Kinetic simulation of a strained film is much more challenging computationally than the unstrained case because of the long-range nature of elastic interactions. First principles calculations [@Fujikawa2002; @Montalenti2005; @Lu2005; @Shklyaev2005; @Montalenti2007] and molecular-dynamics simulations [@Maxim2003; @Retford2007] have provided important information on the energetics of the relevant surfaces and steps, but are in general limited to the studies of static properties at small system sizes. Continuum simulations in contrast are instrumental for investigating large scale and long time behaviors [@Tu2004; @Chiu2007], but lacking atomic discreteness, nucleation events associated with island formation and shape transition cannot be naturally simulated. KMC simulations based on lattice models are hence unique in allowing large scale atomistic studies on the dynamics of strained heteroepitaxy.

Using a ball and spring lattice model for elastic solids, Orr et. al [@Orr1992] conducted early KMC simulations of strained layers in 2D. In the simulations, the elastic energy of the system has to be computed repeatedly in order to simulate the atomic hopping events responsible for the morphological evolution. Applying more advanced algorithms for the solution of the elastic problem and the sampling of atomic hopping events, simulations using large and moderate system sizes in 2D [@Lam2002; @Gray2005; @Schulze2009] and 3D [@Lung2005; @Smereka2006; @Lam2007; @Lam2008; @Smereka2009] respectively became possible. The model was extended recently to model (105) facets [@Lam2010] and atomic intermixing with substrate atoms [@Smereka2010]. Alternatively, KMC simulations can also be performed efficiently using more approximate forms of elastic interactions [@Ratsch1996; @Meixner2001; @Zhu2007].

This paper is organized as follows. Section II explains our multistate model for elastic solids which can account for both a shallow and a steep facet. The KMC simulation results are presented in Sec. III. In Sec. IV, island energies and the island transition rate are calculated theoretically. In Sec. VI, a scaling relation between the transition island size and the Ge concentration is derived. We conclude in Sec. VI with some further discussions.

multistate solid on solid model
===============================

Our model is based on a ball and spring square lattice model of elastic solids for Ge$_x$Si$_{1-x}$ on Si [@Lam2002]. The substrate lattice constant is $a_s=2.72$ Å  while the film material admits a lattice misfit $\epsilon = 0.04 x$. Each node on the lattice represents an atom and it is connected to its nearest and next nearest neighbors by elastic springs with force constants $k_N=13.85$ eV/$a_s^2$ and $k_{NN}=k_N/2$ respectively. This choice gives the correct modulus $c_{11}$ of silicon and a shear modulus constant along tangential and diagonal directions.

In this work, (100), (105) and (113) surfaces must be effectively simulated. A leveled surface in the model naturally accounts for a (100) surface. However, lattice models generally lead to islands with a single type of sidewalls at 45$^\circ$ inclination or steeper [@Orr1992; @Lam2002; @Gray2005; @Schulze2009; @Lung2005; @Smereka2006; @Lam2007; @Lam2008; @Smereka2009]. A multi-state extension has been introduced recently in Ref. [@Lam2010] to effectively model the much shallower (105) facets of a pyramid in 2D. We now further generalize it to simulate both shallow and steep facets with slopes $$\label{eq:s}
s_1=1/5 \mbox{~~~ and ~~~~} s_2=1/2$$ comparable to those of realistic facets in pyramids and domes. Specifically, atoms are normally represented by squares. To effectively model surface reconstructions leading to specific facets, we allow local deformation of the topmost atoms in the film or substrate into trapezoids each characterized by a tilt variable $\sigma_i$ and an extension variable $\kappa_i$. Here, $\sigma_i$ gives the slope of the upper surface of an atom and equals $$\sigma_i = 0,  ~ \pm 1/5 , \mbox{~ or ~} ~\pm 1/2$$ at a locally undeformed region, a shallow facet or a steep facet respectively. Allowing atomistically flat shallow and steep facets further requires additional freedoms of vertical stretching or compression of the topmost atoms by $$\kappa_i = 
    \left\{ 
        \begin{array}{ll}
            0 & \mbox{for $\sigma_i = 0$}\\
        0 ,  ~ \pm 1/5,  \mbox{~ or ~} \pm 2/5 ~~~ & \mbox{for $\mid \sigma_i\mid=1/5$}\\
        \pm   1/4 ~~~ & \mbox{for $\mid\sigma_i\mid= 1/2$}
        \end{array}
    \right.$$ This characterizes a total of $15$ possible local deformation states. All lengths are measured in unit of $a_s$ throughout this paper. Atomic column $i$ with $h_i$ atoms is thus trapezoidal in general with the left and right edges of heights $h_i^a$ and $h_i^b$ given by $$\begin{aligned}
  h_i^a &=& h_i + \kappa_i - \frac{\sigma_i}2\\
  h_i^b &=& h_i + \kappa_i + \frac{\sigma_i}2.\end{aligned}$$ A surface step in between column $i$ and $i+1$ has a height $$\delta_i = \mid h_{i+1}^a - h_{i}^b \mid$$ projected along the vertical direction. Figures \[F:surf\](a)-(b) show examples of atomic configurations.

[0.75]{}

Misfit induced elastic strain is assumed to be completely independent of the local deformations associated with faceting introduced above. The elastic relaxation energy $E_s$ of the system is defined as the total energy storied in all springs at mechanical equilibrium compared with that in the homogeneously strained state. The bond energy of the system is defined relative to that of a flat surface by $$\begin{aligned}
\label{Eb}
E_b  &=& \sum_i \left[ 
\phi_{\alpha_i} +  \psi(i,i+1) +\omega_{\alpha_i \alpha_{i+1}}( {\delta_i})
\right]+ E_R ~\end{aligned}$$ where the facet-type label $\alpha_i$, depending on $\mid\sigma_i\mid$, indicates if column $i$ is locally undeformed ($\alpha_i=0$) or corresponds to a shallow ($\alpha_i=1$) or steep facet ($\alpha_i=2$). The facet formation energy per site $\phi_{\alpha_i}$ equals $\phi_0=0$, $\phi_1 = 5$ meV, or $\phi_2=50$ meV. These values control the relative stability of the facets in our simulations. They are chosen empirically so that shallow and steep facets start to emerge on islands of appropriate sizes. Also, the facet interface energy $\psi(i,i+1)$ is non-zero only at the boundary between either different facet types or different facet orientations (i.e. $\sigma_i \neq \sigma_{i+1}$) where it equals $\psi_{\alpha_i\alpha_{i+1}}$ with $\psi_{01}=\psi_{11}=0.35$ eV, $\psi_{12}=\psi_{22}=0.5$ eV, and $\psi_{02}=\psi_{01}+\psi_{12}$, assuming $\psi_{\alpha\alpha'}=\psi_{\alpha'\alpha}$. The $\omega$ term represents surface step energy. It equals $\gamma\delta_i/2$ on a locally undeformed region where $\gamma=0.5$ eV is the nearest neighboring bond energy. If the site $i$ or $i+1$ belongs to a shallow or a steep facet, it equals ${\beta_\alpha} [ 
    1+\chi - \chi \exp({ 1 - { {\delta_i}/s_\alpha }}) ] +
    {\gamma} (  {\delta_i}- s_\alpha) /2$ where $\alpha$ is the larger of $\alpha_i$ and $\alpha_{i+1}$. Here, $\beta_1=0.3$ eV and $\beta_2=0.2$ eV represent the single step energies on shallow and steep facets respectively and $\chi=0.3$ dictates the tendency of step bunching. To discourage very steep regions, the intrinsic step height defined by $\delta_i' = \mid h_{i+1}-h_i\mid$ disregarding local deformation is constrained to $\delta_i'\le 1$ and furthermore each pair of consecutive upward or downward intrinsic steps with $\delta_i'=1$ contributes 0.15 eV to the repulsion energy $E_R$. The constraint also limits the step heights $\delta_i$ to bounded values, although double steps in particular, which have heights $2/5$ and $1$ on shallow and steep facets respectively are still possible.

The KMC approach simulates every hopping event of a topmost film atom $m$ according to the rate $$\label{rate2}
\Gamma(m) = 
{R_0}\exp \left[ \frac{\Delta E_b(m)  + \Delta E_s(m) +
    E_0'}{k_{B}T}\right].$$ Atoms then lands random on any other site at most 8 columns away. Here, $\Delta E_b(m)$ and $\Delta E_s(m)$ denote the change in the bond energy $E_b$ and the strain energy $E_s$ of the system when the site is occupied versus unoccupied. We put $E_o'=-\gamma-0.67$ eV and $R_o = 4.1 \times 10^{11} s^{-1}$. This gives the appropriate adatom diffusion coefficient for silicon (100). Due to the long-range nature of elastic interactions, the repeated calculations of $\Delta E_s(m)$ dominates the simulation run time and we handle it using a Green’s function method together with a super-particle approach [@Lam2008]. Exposed substrate atoms are not allowed to hop. Elastic couplings of adatoms with the rest of the system are weak and are neglected for better computational efficiency. Atomic hoppings are assumed to preserve the local deformation states. After every period $\tau$, a set of deformation states will be updated. We put $\tau = 2/\Gamma_{ad}$ where $\Gamma_{ad}$ is the adatom hopping rate on an locally undeformed region easily calculable from Eq. (\[rate2\]). At an odd (even) numbered updating event, all odd (even) lattice sites will be updated. The variables $\sigma_i$ and $\kappa_i$ at those sites are re-sampled from the allowed set of 15 possible combinations using a heat bath algorithm based on the relative probability $\exp(-E_{b}/kT)$. Our model obeys detailed balance. The dynamical rules described above reduces back to those used in Refs. [@Lam2002] at locally undeformed regions. We have critically checked our software implementation, in particular using a Boltzmann’s distribution test [@Lam2008], which is found to be indispensable in verifying that practically all, but not only the dominating hopping pathways can be correctly simulated. We also have checked in small scale simulations that restricting hoppings to only nearest neighboring sites rather than allowing long jumps gives similar results except for an insignificant shift in the time scale. Wetting layers on the substrates are believed to be relatively immobile and are not simulated for simplicity.

Simulation results
==================

[0.79]{} [0.6]{}

Figures \[F:surf\](a) show a shallow facet with steps from a small scale simulation. The steps are subsequently smoothed out by the formation of a steep facet as shown in Fig. \[F:surf\](b). Figure \[F:profile\](a) shows snapshots of a surface from a typical simulation of deposition at 4 ML/s, 700$^\circ$C and $x=1$ on a substrate of width 2048. Successive profiles are displaced vertically for clarity. Stepped mounds first develop and some of them matures into pyramids bounded by shallow facets as explained in Ref. [@Lam2010]. Some of the pyramids further turns into domes bounded mainly by steep facets often with regions of shallow facets at the top. Figure \[F:profile\](b) shows the detailed evolution of one of the domes. Steep facets on either side of a pyramid form independently. The transition hence often goes through a meta-stable half-dome stage. The formation of most domes is preceded by steps appearing on the shallow facets as shown in Fig. \[F:surf\] and Fig. \[F:profile\](b) and as proposed in Ref. [@Montalenti2004]. A close examination reveal that these steps are highly dynamic and continuously bunches, separate and diffuses around. After accumulating a considerable total step height, they transform highly reversibly into a steep facet. As the total height of the steps increases, the resulting steep facet is more stable and eventually become fully stabilized. A smaller number of domes are initiated instead by the formation of steep facets at the base of the pyramid when shallow facets temporarily decay into unfaceted regions due to thermal excitations.

0.99

0.99

Large scale deposition simulations have been performed at temperature $T$ from 450$^\circ$C to 850$^\circ$C at $x=1$ on lattices of width 2048. The deposition rate $R$ varies from 0.006 to 8 ML/s chosen empirically to generate typically 3 to 5 pyramids or domes on each substrate. The low island density minimizes elastic interactions among islands which are known to alter the dynamics [@Floro2000; @Capellini2003]. An island is defined as one in which all constituent columns must be at least 2 atoms tall. We measure island size in number of atoms $n$ so that $\sqrt{n}$ is proportional to the linear size of the island. Also, island aspect ratio is defined by $r= h / 2l$ where $h$ is the height of the highest point of the island and $2l$ is its basewidth. Figure \[F:aspect\](a) shows a scatter plot of the aspect ratio $r$ against $\sqrt{n}$ for all islands from 3 independent runs at each temperature. Measurements are conducted throughout the evolution. Time averaging of values associated with individual islands over short periods are performed, but no ensemble averaging is done as each island develops in general at a different pace. We observe that $r$ first converges towards 0.1 as islands transform from stepped mounds into pyramids. It then rises again to around 0.2 characterizing the dome transition similar to experiment findings in Ref. [@Montalenti2004]. The much lower density of the data points at $0.1 < r < 0.2$ corresponding to highly unstable intermediate states was also observed in Ref. [@Montalenti2004]. The morphologies of these intermediate states have been shown in Fig. 2(b). Results in Fig. 3 reveals two distinct trends. For the runs at $T \agt 650 ^ \circ C$, all islands follows an identical evolution path with the dome transition occurring at size ${n_c}\simeq 900$. In contrast, at lower temperature $T\alt
600 ^\circ C$, the transitions are delayed randomly to increasingly larger sizes. At $T=450^\circ$ for instance, $n_c$ ranges from about $900$ to $2000$. We will explain these distinct trends at the end of this section. Similar deposition simulations are also performed at $T=700^\circ$C and $R=1$ ML/s for Ge concentration $x$ from 0.6 to 1 with 3 independent runs in each case. We find that the dome transitions occur at increasingly larger island sizes as $x$ decreases. The precise dependence is easily illustrated by a rescaled plot of $r$ against $\sqrt{n} /
x^{-\zeta}$ with $\zeta=1.69$ as shown in Fig. \[F:aspect\](b). Data for various values of $x$ collapse reasonably well into a single curve, implying $$r = f(n/x^{-\zeta})$$ where $f$ is a rescaled function. To calculate $\zeta$ used above, we have measured the transition size $n_c$ by averaging the sizes of all islands right at the transitions with $ 0.12 \le r \le 0.16$. The resulting plot of $\sqrt{n_c}$ against $x$ in log-log scales is shown in the inset of Fig. \[F:aspect\](b). A linear relation observed in the log-log plot implies $$\label{nc_x}
\sqrt{n_c} \sim x^{-\zeta}$$ and a linear fit gives $\zeta=1.69$. This scaling relation will be explained in Sec. IV.

0.99

To study the energies of individual islands, we have performed simulations on annealing of single pyramids directly constructed on smaller substrates each of width $256$. The annealing temperature is 550$^\circ$C. The pyramid is initially bounded by shallow facets and sits on an otherwise empty substrate surface. It contains 1230 atoms and has a basewidth slightly less than the lattice width. This number is chosen empirically to be just sufficient to ensure an irreversible dome transition. Figure \[F:theory\] shows a scatter plot of the system bond energy $E_b$, strain energy $E_s$ and total energy $E=E_b+E_s$ against $r$ measured during the annealing from 16 independent runs. Since a single island dominates, these energy of the whole system approximates those of an island. Only time averaging of the values over short periods but no ensemble averaging has been carried out. The symbols used indicate if the data points correspond to pyramids, half-domes, or domes. The geometries are identified reliably by the number of steep facets present. The result indicates that there is an energy barrier for the transition. Moreover, a number of data points associated with half-domes cluster around $r \simeq 0.12-0.14$ showing that the geometry characterizes a meta-stable state. We next show that the dome transition is an activated process. We have repeated the above simulations on the annealing of pyramids at $T$ from 450$^\circ$C to 850$^\circ$C. A dome transition time $\tau$ defined as the average annealing duration required to reach an aspect ratio $r \ge
0.12$ is measured. Values of $\tau$ each averaged over 16 independent runs are plotted against $1/T$ in the inset in Fig. \[F:theory\]. The data fits well to $$\label{tau}
\tau = \tau_o \exp( \Omega_0/kT)$$ with $\Omega_0=1.97$ eV and $\tau_o
= 2.9 \times 10^{-12}$ s. The Arrhenius temperature dependence is typical of activated processes. The value of $\Omega_o$ will be explained in Sec. IV. We further repeat the simulation 300 times at $700^\circ$C. The values of $\tau$ measured are histogrammed. The result is well fitted by an exponential distribution. This further supports that the dome transition is an activated process.

With the shape transition time $\tau$ known, the distinct trends followed by the data in Fig. \[F:aspect\](a) can now be explained. The deposition rate $R$ has been empirically chosen to produce a constant and sparse island density. The choice hence ultimately depends on pyramid nucleation and coarsening dynamics and is not directly related to the dome transition dynamics. For the $T \ge 650^\circ$C runs, we find that $1/\tau \gg R$. The shape transition is thus fast compared with deposition and hence also the island growth. The dome transition process is only limited by the availability of atoms. The geometry as characterized by the aspect ratio $r$ therefore only depends on $\sqrt{n}$ and is independent of $T$ and $R$ as shown in Fig. \[F:aspect\](a). In contrast, for the $T < 650^\circ$C runs, we get $1/\tau \alt R$. Island growth can then be fast enough to out-run the dome transition, which becomes rate-limited. There is a significant random waiting time for the dome transition following an exponential distribution during which the island may already have grown to a larger size. The transition thus occurs at a more broadly distributed island size $n_c$ as observable from Fig. \[F:aspect\](a). Note that if we consider for example a higher island density by increasing the values of $R$ used, the characteristic temperature separating the two trends, which is found to be around $650^\circ C$ here, will increase.

Theory of shape transition
==========================

We now present a detailed theoretical analysis based on generic forms of elastic and facet energy for the transition of a pyramid into a half-dome in 2D, which is applied to interpret our KMC simulation results. Half-domes are meta-stable and they quickly transform into domes. Our formulation is consistent with that of Montalenti [*et al.*]{} who have shown using energy parameters for Ge/Si from first-principle calculations that the dome transition is energetically favorable for sufficiently large pyramids [@Montalenti2004].

Consider an island of size $n$ initially in the form of a pyramid with a half-basewidth $l_0$. We have $$\label{n}
n=s_1l_0^2.$$ Geometrical rearrangements can lead to the formation new atomic layers of total vertical thickness $v$ on one of the facets as shown in the schematic diagram in Fig. \[F:theory\]. The new half-basewidth $l$ is then related to $v$ by $$\label{n2}
n= s_1 l^2 + (u_B-u_A) v$$ where $u_A$ and $u_B$ denote the positions of the midpoints A and B on the edges of the new layers measured from the apex of the base pyramid. A single atomic step on a shallow facet has a height $s_1$. We assume for simplicity that there are $v/s_1$ single steps at point B. Using Eq. (\[Eb\]), the bond energy of the pyramid is $$\label{py}
  E_b^{py}= 2 \phi_1 l + 2 \psi_{01} + \psi_{11}  + \beta_1 v/s_1.$$ The pyramid becomes a half-dome when all the steps at point B turn into a steep facet of $\lambda = v /
(s_2-s_1)$ columns wide. The bond energy of the resulting half-dome also follows from Eq. (\[Eb\]) and we get $$\label{dome}
{E_b}^{hdome} = 2 \phi_1 l + 2 \psi_{01} + \psi_{11}  +   \frac{
  \phi_2-\phi_1}{s_2-s_1} v + 2\psi_{12}$$ Neglecting the small difference in the elastic energies of the two geometries, the island takes the form with the lowest bond energy $E_b ={\mbox minimum}\{E_b^{py},E_b^{hdome}\}$.

The elastic energy of the pyramid and the half-dome is assumed to be identical and is calculated by approximating both edges of the new layers as vertical walls located at A and B. A shallow island approximation [@Tersoff1994; @Tersoff1999] gives $$\label{eq:Es}
E_s = C \epsilon^2  \int\int dx dx' s(x)
s(x') \ln \frac{\mid x - x' \mid}{a_c},$$ where $s(x)$ denotes the local surface slope of the island at position $x$ and $a_c$ is a spatial cutoff. We put $s(x) = - sgn(x) s_1 + v
\delta(x-u_A) - v \delta(x-u_B)$ where $sgn(x)$ and $\delta(x)$ represent the sign function and the Dirac delta function. In 2D, $C=\sigma_b^2a_s^2/\pi\epsilon^2 Y$ where $\sigma_b \propto \epsilon$ is the $xx$ component of the bulk misfit stress and $Y$ is the Young’s modulus. For our model, a simple calculation based on lattice elasticity gives $C = 4 k_n
a_s^2/3\pi$.

Performing the integrations in Eq. (\[eq:Es\]), we have

$$\label{Es}
  E_s = -2 C \epsilon^2\left\{ 2 \ln(2) s_1^2 l^2 +  s_1 v \sum_{p=A,B}\xi_P 
  \left[ l \ln \left(\frac{l+u_p}{l-u_p}\right) +
  u_p \ln \left(\frac{l^2}{u_p^2}-1 \right)\right] 
  + { v^2} \ln  \left(\frac{u_B-u_A}{a_c} \right)\right\}$$

where $\xi_A=-1$ and $\xi_B=1$. In the following, we put $a_c =
e^{-3/2}\bar{\lambda}$ where $\bar{\lambda}=(\lambda+2u_A)/2$ is the average spatial extent of the misfit force monopoles applied over the edges at A and B. It can be shown that this choice gives the correct energy when approaching the point force limit. From simple geometry, $u_A=v/4s_1$. We calculate $u_B$ by minimizing the total island energy $E$ at small $v$. Linearizing $E$ from Eqs. (\[n2\])-(\[dome\]) and (\[Es\]) w.r.t. $v$, it can be shown after some algebra that for both pyramids and half-domes, $E$ is minimized at $$\label{uB}
{u}_B = l_0 \left [ 1 + 4 \exp\left (- \frac{\phi_1}{2 C
      \epsilon^2 s_1^2 l_0}\right) \right ]^{-1/2}.$$ The total energy cost $\Delta E$ of an island compared with the initial pyramid can then be calculated. For a half-dome, we get, up to linear order in $v$, $$\label{DElinear}
\Delta E =  \left[ \frac{
   \phi_2-\phi_1}{s_2-s_1}
- 2 C \epsilon^2  s_1 l_0 \ln \left ( \frac{l_o+{u_B}}{l_o-{u_B}} \right )\right]
v 
+ 2 \psi_{12}.$$ Equation (\[uB\]) gives the energetically most favorable position for the initial formation of a steep facet. Using the KMC model parameters for $x=1$, it gives $u_B/l_o
\simeq 0.53 $ for islands around the transition size. The result is in general close to a limiting value $u_B/l_0=5 ^{-1/2}
\simeq 0.447$ obtained by neglecting the shallow facet formation energy $\phi_1$. The energies $E_b$, $E_s$ and $E$ are numerically calculated for various layer thickness $v$ adopting the KMC model parameters for $x=1$ using Eqs. (\[n2\])-(\[dome\]) and (\[Es\])-(\[uB\]). The island aspect ratio $r$ is also calculated as a function of $v$ using $r = s_1/2 + v/4l$ and Eq. (\[n2\]). In Fig. \[F:theory\], the energies are plotted as solid lines against $r$. We have assumed an island size ${n}= 1183$ atoms measured during the dome transition in the KMC annealing simulations responsible for the data points in Fig. \[F:theory\]. The only tunable parameter is a fitted additive constant $3.4$ eV for $E_b$, which accounts for the bond energies of all other excitations in the system. It nevertheless plays no role in further calculations. The theoretical estimates of the energies generally show reasonable quantitative agreement with simulation results as observed in Fig. \[F:theory\]. The main discrepancies are due to errors in $E_s$, since the shallow island approximation is known to overestimate the elastic relaxation at large $r$. Nevertheless, important features including a shape transition energy barrier and a meta-stable half-dome state are correctly reproduced and these will be further studied.

From Fig. \[F:theory\], both theory and KMC simulation show an energy barrier for the dome transition associated with a maximum in the total energy $E$. Its location follows theoretically from $E^{py}=E^{hdome}$. We get a barrier height $\Delta
E_{max}=0.88$ eV which occurs at $r= 1.03$ or $v = 0.81$. It corresponds to $v/s_1\simeq 4$ new atomic layers on the shallow facet. For $r$ below and above 1.03 respectively, pyramid and half-dome are the energetically preferred states. Due to the small value of $v$ at the barrier, $\Delta E_{max}$ is dominated by the steep facet interface energy term as can be deduced from Eq. (\[DElinear\]), i.e. $\Delta E_{max} \simeq 2\psi_{12} = 1$ eV. The dominance of the steep facet interface energy on the transition energy barrier may be a general feature applicable also in 3D.

The existence of an energy barrier confirms that the dome transition is an activated process as suggested in Sec. III. The transition rate $R$ hence follows the Arrhenius form $ R = \nu \exp(- \Delta
E_{max}/kT)$, where $\nu$ denotes an attempt frequency. Assuming that the transition is limited by the diffusion of adatoms on the shallow facet, one expect $\nu \propto \rho D$, where $\rho=\exp(-E_{ad}/kT)$ and $D \propto \exp(-\Omega_{ad}/kT)$ are the adatom density and diffusion coefficient on the shallow facet. Here, $E_{ad}=0.6$ eV and $\Omega_{ad}=0.57$ eV are the adatom formation energy and hopping energy barrier on the shallow facet calculated using Eqs. (\[Eb\]) and (\[rate2\]). In particular, $\Omega_{ad}$ is not far from a previous estimate from first-principles calculations [@Montalenti2004b]. The dome transition time $\tau \propto 1/R$ is hence given by $$\label{tau2}
\tau
\propto \exp\left(\frac{\Delta E_{max} + E_{ad} + \Omega_{ad}}{kT}\right)$$ A comparison with Eq. (\[tau\]) leads to $\Omega_0 = \Delta E_{max} + E_{ad} + \Omega_{ad}$. It gives $\Omega_o= 2.05$ eV in agreement with $\Omega_0=1.97$ eV obtained above from KMC simulations.

scaling of shape transition size
================================

We first assume quasi-equilibrium conditions in which the dome transition is fast compared with island growth. It applies to the cases of slow deposition and small transition barrier and is valid for our KMC simulations at $T \agt 650^\circ C$ (see Fig. \[F:aspect\](a)). The island energy $E$ from theory as shown in Fig. \[F:theory\] exhibits a local minimum representing the meta-stable half-dome state. The energy rises again at larger $r$ because the base pyramid then becomes too small to relief the elastic energy efficiently. The dome transition occurs only if the half-dome is sufficiently stable, say of energy $kT$ below that of the initial pyramid. For island at the transition size $n_c$, the minimum of $\Delta E$ hence follows $\Delta E_{min}=-kT$. We can then calculate $n_c$ numerically using Eqs. (\[n2\]), (\[dome\]), (\[Es\]) and (\[uB\]) and the values are plotted against $x$ as a solid line in the inset of Fig. \[F:aspect\](b). Note that no tunable parameter is involved. The values are in reasonable agreement with the KMC results and supports the scaling relation in Eq. (\[nc\_x\]) with $\zeta \simeq 1.49$ consistent with $1.69$ found in simulations.

In addition to the numerical estimate of the exponent above, better insights are obtained by deriving an exact exponent $\zeta=2$ valid asymptotically in the small misfit limit, i.e. $\epsilon \propto x
\rightarrow 0$. Assume that the relative position of the steep facet is independent of $\epsilon$ so that $u_B \propto l_o$, which will be justified later. Simple scaling properties in 2D elasticity imply that the change in the elastic energy of a half-dome compared with the initial pyramid follows $\Delta E_s = - \epsilon^2 l_0^2 {g_2}(v/l_o)$ for some function ${g_2}$. This is also explicitly derivable from Eq. (\[Es\]). The total energy cost of a half-dome is hence $$\label{DE2}
\Delta E = A_2  v - \epsilon^2 l_0 ^2 {g_2}(v/l_o) + B_2$$ where $A_2 = ({\phi_2-\phi_1})/({s_2-s_1})$ and $B_2=2 \psi_{12}$. At island transition size $n_c$ and considering a layer height $v$ minimizing $\Delta E$ to $\Delta E_{min}=-kT$, we have $$\label{DE3}
\epsilon^2 l_o^2 {g_2}(v/l_o) - A_2 v =  B_2 + kT$$ It means that the elastic energy gain must overcompensate the facet formation energy cost by an excess amount $B_2 + kT$. As $\epsilon
\rightarrow 0$, we will see in the following that the shape transition occurs at a larger island size. Both energy terms on the l.h.s. of Eq. (\[DE3\]) increase unboundedly and must balance each other, while the constant energy excess becomes negligible. Therefore, we have $\epsilon^2
l_o^2 {g_2}(v/l_o) = A v$. The meta-stable half-dome state at transition is thus characterized by the scaling solution $l_o \sim v \sim
\epsilon^{-2}$. Using $\sqrt{n_c} \propto l_0$ and $\epsilon=0.04 x$, we get $$\label{scaling}
\sqrt{n_c} \sim x^{-2}$$ i.e. ${\zeta=2}$. This solution is consistent with the assumption that $u_B$ is independent of $\epsilon$ as deduced using Eq. (\[uB\]). It also justifies that both terms on the l.h.s. of Eq. (\[DE3\]) increase unboundedly as $\epsilon \rightarrow 0$.

For finite $\epsilon$ and $x$, the energy excess $B_2 + kT$ in Eq. (\[DE3\]) is not negligible. It gives a finite size correction to the exact scaling in Eq. (\[scaling\]). This results at an effective scaling $\sqrt{n_c} \sim x^{-\zeta}$ with $\zeta \alt 2$ consistent with $\zeta =1.69$ found in our KMC simulation.

In 3D, the dome transition is initiated by step bunching close to the mid-level on a pyramid [@Montalenti2004]. Generalizing our discussion above, we consider the formation of a square ring of steep facet of vertical thickness $v$ on a 3D pyramid of basewidth $2l_o$. Generalizing Eq. (\[DE2\]), the energy cost is $$\label{DE3d}
\Delta E = A_3  v l_o - \epsilon^2 l_0 ^3 {g_3}(v/l_o) + B_3 l_o$$ where the terms on the r.h.s. similarly denote the facet formation energy, the elastic energy gain, and the steep facet interface energy respectively, for some smooth function $g_3$ and constants $A_3$ and $B_3$ independent of $v$ and $l_o$. A similar calculation leads to $n_c^{1/3}
\propto l_o \sim x^{-2}$. This generalizes Eq. (\[scaling\]) to $$\label{scaling2}
{n_c}^{1/d} \sim  x^{-\zeta}$$ with the same exponent $\zeta=2$ in dimension $d=2$ or 3 for $x\rightarrow 0$. For finite $x$, an effective exponent $\zeta\alt 2$ is expected in both dimensions.

The asymptotic scalings derived above for quasi-equilibrium conditions essentially follows from the balance between the steep facet formation energy and the elastic energy gain which scales with the island size differently. It is analogous to the scaling predicted for island formation size based on the Asaro-Tiller-Grinfeld instability theory [@Politi2000]. Due to the simplicity, it is also rather robust as will be shown below. Instead of assuming quasi-equilibrium conditions, the dome transition can be limited by the kinetics. This may be appropriate in particular at 3D for small $x$ since the barrier predicted above can become too large to overcome. The transition is then delayed to a larger island size which lowers the barrier. Let us then assume a very simple transition criterion that the energy barrier $\Delta E_{max}$ must not exceed a given value, say a few times of $kT$. A similar calculation for $x \rightarrow 0$ again arrives at Eq. (\[scaling2\]) with the same exponent $\zeta=2$ in both 2D and 3D. The solution also requires $l_o \sim \epsilon^{-2}$ but $v \sim
\epsilon^0$ and we have used $g_d(z) \propto z$ for $z\rightarrow 0$ which readily follows in the 2D case from Eq. (\[DElinear\]) after neglecting a logarithmic factor.

discussions
===========

In the calculations above, $x$ denotes the actual Ge concentration in the film so that $\epsilon=0.04 x$. In experiments, it can differ greatly from the nominal concentration due to intermixing with substrate atoms and this complicates interpretation of experiment results. As a very rough estimate, experiments on deposition of pure Ge at e.g. 450$^\circ$C and 700$^\circ$C have found dome transition occurring at island volume $2800$ nm$^3$ [@Drucker2002] and $2
\times 10^5$ nm$^3$ [@Rastelli2006]. Neglecting compositional non-uniformity, it has been estimated that the actual Ge concentrate in the islands is $x = 0.82$ and 0.43 respectively at 450$^\circ$C and 700$^\circ$C [@Drucker2003]. This gives a very preliminary estimate for the scaling exponent $\zeta\simeq 2.2$ consistent with the asymptotic value $\zeta= 2$ derived above, although more experiments are required for a reliable conclusion.

Our simulations and theoretical calculations show that occurrence of well-defined dome transitions depends strongly in particular on the formation and interface energies of steep facets. The detailed dependences of these energies on the Ge concentration and temperature are neglected. We have also neglected the spatial non-uniformity of Ge concentration, surface stress and realistic crystal elastic anisotropy. They should have significant quantitative impacts on the shape transition, but are not expected to alter the dynamics described here qualitatively.

In summary, we have generalized a multi-state lattice model for elastic solids to account for both shallow and steep facets with tunable energy parameters. Using this model, we perform kinetic Monte Carlo simulations to study the pyramid-to-dome transition in the heteroepitaxy growth of Ge$_{x}$Si$_{1-x}$ on Si in 2D. For sufficiently slow deposition, the shape transition occurs at an island size independent of temperature and deposition rate. A scaling relation between the transition size and the Ge concentration is observed. For fast deposition, the transition can be delayed randomly to a larger island size. For annealing simulations, the shape transition time is found to follow an Arrhenius form. A theory based on elastic energy in the shallow island approximation and simple forms of facet formation energies is derived. Numerical solutions of the energetic equations give island energies, shape transition size, and shape transition rate in reasonable agreement with simulations. The shape transition energy barrier is dominated by the interface energy between the shallow and the steep facets. We have also derived analytically an exact scaling rule between the transition size and the Ge concentration applicable in the limit of small Ge concentration which is expected to be valid in both 2D and 3D. A finite size correction to the scaling at higher Ge concentration is explained.

This work was supported by HK GRF, Grant No. PolyU-5009/06P and PolyU Grant No. G-U354.