TFOx: A versatile kinetic Monte Carlo program for simulations of island growth in three dimensions

Computational Materials Science 91 (2014) 292–302

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

TFOx: A versatile kinetic Monte Carlo program for simulations of island growth in three dimensions Qing Zhu a,⇑, Chris Fleck a, Wissam A. Saidi a, Alan McGaughey b, Judith C. Yang a,c a

Department of Chemical and Petroleum Engineering, University of Pittsburgh, Pittsburgh, PA 15261, United States Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, United States c Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, United States b

a r t i c l e

i n f o

Article history: Received 26 November 2013 Received in revised form 21 March 2014 Accepted 23 April 2014 Available online 27 May 2014 Keywords: Kinetic Monte Carlo Film oxidation Three dimensional Potential gradient Ehrlich–Schwöbel barrier

a b s t r a c t A three dimensional (3D) kinetic Monte Carlo (KMC) code has been developed that simulates the general behavior of the 3D irreversible nucleation and growth of epitaxial islands, as motivated by experimental observations of oxide nuclei formation and growth during the early stages of copper oxidation. This package was originally a versatile two dimensional (2D) KMC code [Thin Film Oxidation (TFOx)] that considered a variety of elementary steps, including deposition, adsorption, surface diffusion, aggregation, desorption, and substrate-mediated indirect interactions between static adatoms. We extended TFOx to describe 3D island growth. This new version of TFOx is composed of a C++ console program and Python graphical user interface (GUI), such that parameterized simulation, parallel execution, and 3D growth capabilities are feasible. We examined the effects of the potential gradient and the Ehrlich–Schwöbel barrier and found that the 3D island morphology is significantly influenced by the incorporation of these two factors. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Metal alloys and oxides have been at the foundation of modern industry for several decades, as they provide the building materials for components that operate under harsh environments such as blades in the turbine engine, and are also widely used in catalysis industry [1–3]. Understanding the process of oxidation on metal surfaces is crucial for both corrosion control and the use of metals as catalysts in reactions involving oxides, especially due to evidence that metal oxides can also act as catalysts [4]. In the case of metallic alloys, the formation and maintenance of a thermally grown oxide scale is often needed for extended surface protection. Zhou and Yang demonstrated that the transient early stages of metal oxidation bear a striking resemblance to heteroepitaxy [5]. This discovery is not only important with respect to the control of metal corrosion, but also promises advanced engineering methods towards the fabrication of novel metal oxide nanostructures. However, the oxidation reaction of metal surfaces is highly complex. Different processes, many of them coupled, are involved from the onset of reaction: that is, from the absorption and dissociation of gaseous molecules at the surface through to the competitive nucleation and growth processes that ultimately ⇑ Corresponding author. Tel.: +1 8147532677. E-mail address: [email protected] (Q. Zhu). http://dx.doi.org/10.1016/j.commatsci.2014.04.053 0927-0256/Ó 2014 Elsevier B.V. All rights reserved.

lead to the establishment of a steady-state scale. To attain sustained surface stability in corrosion control or successful growth of desired metal oxide nanostructures, it is necessary to understand the energetics and dynamics of the initial- and early-stage oxidation processes that occur at the nano- and micro-scales. To date, such an understanding has been limited due largely to difficulty in taking measurements, which in turn is due to the transient nature of these processes. Therefore, there exists a significant gap between the experimental techniques that operate in large temporal scale and the rapid-occurring reactions at the initial stage of the surface oxidation. Benefiting from recent advances in computation (software and hardware) capabilities, we are now poised to fill this gap through the assistance of theoretical simulations. Theoretical studies at various levels, from atomistic quantum mechanical calculations to classical mechanical molecular dynamic (MD) and Monte Carlo (MC) simulations and ultimately to continuum scale microkinetic modeling have been widely explored in surface sciences [6]. Theoretical calculations can provide a deeper understanding of solid surfaces, supplying critical information such as molecular adsorption energies, chemical reaction energy barriers, and the dynamical process of surface reactions. Through the combination of different computational methods, a multi-scale modeling framework of surface oxidation may be established that can be used to better interpret experimental measurements.

293

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

An important bridge that connects atomistic scale simulations with continuum models is kinetic Monte Carlo (KMC) [7]. A KMC simulation is built upon a stochastic model of the desired reaction system [8], where the behavior of the system over time is determined from the rates of relevant physical or chemical events. Unlike classical MC methods, KMC deals with systems that are not in equilibrium where dynamical phenomena can be simulated. Thus sometimes the term Dynamic Monte Carlo (DMC) is used. Although KMC and DMC are usually interchangeable, recently it has come to a common understanding that KMC refers to the case where the rate constants for the system evolution are known, while DMC is used when such a rate table is unavailable [9]. As KMC allows for simulation of the time and length scales of surface science experiments, there has been increasing interest in modeling surface reactions with KMC [10–18]. It has been applied successfully to study the initial nucleation process for homoepitaxial and heteroepitaxial thin film growth [19–24]. However, KMC modeling of 3D nanostructure formation during surface reactions has rarely been explored, even though this is an important part of film growth [22]. Such KMC modeling is more complicated, as it must account for lattice strain effects and the surface conditions during the reaction. Yet oxidation has additional complexities, such as the reaction of oxygen with the metal surface, that are not taken into account in heteroeptiaxial modeling. In this paper, we describe the KMC simulation of heteroepitaxial thin film formation to include the additional complexities of nano-scale oxidation. A versatile KMC simulation package – the Thin Film Oxidation simulator (TFOx), was previously developed to study the effects of elastic strain and other physical parameters on the morphology of oxide islands in two dimensions [8,25]. TFOx is constructed on top of a Diffusion Limited Aggregation (DLA) model [26]. Atoms deposited on the film surface are subject to a random walk due to Brownian motion and can aggregate to form nuclei. In TFOx, the nucleation reaction is irreversible and is only limited by the diffusion process as we will explain later. Through the inclusion of potential energy gradients [27,28], the elastic strain effect during oxide growth is accounted for which led to a good surface morphology agreement with experimental observations for the Cu(1 0 0) surface oxidation [8]. Despite the success of the 2D TFOx simulator, the growth of 3D nanostructures is yet to be addressed. Herein, to improve modeling of nano-oxidation, we expand TFOx to 3D, and provide an implementation of the Ehrlich–Schwöbel barrier effect [29,30] on surface diffusion. 2. Theoretical background of TFOx In the canonical Metropolis algorithm, a Markovian chain representing a sequence of random samples can be efficiently generated from a probability distribution even if its direct sampling is difficult [31,32]. The Metropolis algorithm is purely governed by thermodynamic equilibrium restrains with no dynamic information included, and thus its direct application to KMC simulations completely fails to achieve the dynamical hierarchy, e.g., all transitions to a state with lower or equal energy are accepted with probably of unity [7]. Instead, the correct description of a dynamic process is provided by the master equation [7,33]:

X X @Pðr; tÞ ¼  Wðrjr0 ÞPðr; tÞ þ Wðr0 jrÞPðr0 ; tÞ @t r0 r0

ð1Þ

where P(r, t) is the probability that the system is in state r at time t, and W(r|r0 ) is the probability per unit time for the transition from state r to r0 . To maintain the dynamical hierarchy, W(r|r0 ) is closely related to the kinetic rate constant C described in transition state theory (TST):

  Q  kB T DE C¼ exp  QR h kB T

ð2Þ

Here h is the Plank constant, Q and QR are respectively the partition functions of the transition state complex and the reactant, and DE is the energy barrier between them. In practice, W(r|r0 ) is usually scaled by the largest C among all possible transition events:

Wr ¼

Cr maxfCr g

1

ð3Þ

so that the dynamical hierarchy can be maintained in a Monte Carlo implementation. In many systems, there are large separations in the time scales between processes of interest. For example, the surface diffusion is typically slower than the thermal vibrations. Thus, if this is not accounted for using clever KMC approaches, large amounts of computer time can be wasted sampling uninterested events. This is known as the stiffness problem. An on-lattice KMC model can overcome the stiffness problem by integrating out the uninterested events such as thermal vibrations, and can only account for slow events such as adatom hopping [33]. Such on-lattice KMC model is widely used in the field of surface science [7,14,18,28,34–37], and is also the foundation of the TFOx program discussed in this paper. For on-lattice models, the master equation can be rewritten as a difference–differential equation [33]:

dri ¼

X

X

p

p

Cþip ðrÞdt 

Cip ðrÞdt

ð4Þ

which describes the change of the occupation (ri) of site i in terms of the transition rates that lead to addition (C+ip) and deletion (C ip) due to different possible processes p. Here the rate for process p (i.e., Cp) can also be represented as:

Cp ¼ v 0 f p

ð5Þ

where m0 is the event attempt frequency, and fp is the successful ratio (e.g., the jumping probability for each jumping attempt). m0 is usually several magnitudes smaller than the prefactor in Eq. (4) for stiff systems, and fp is proportional to the exponential of the event energy barrier. In the current implementation of TFOx, different approximations are utilized. First, molecular dissociation processes are neglected so that only monoatomic adatoms are deposited onto the lattice. Second, the desorption process is only governed by a pre-set lifetime of the adatom. Third, only the hopping to the first neighbor sites are considered in the diffusion process, long jumps that span more than one lattice-site distance are omitted. Therefore, the most important elementary events included in the probability evaluation are the adatom surface diffusion and adatom nucleation/aggregation. 3. Details of TFOx program The TFOx program consists of a main simulation console that is written in C++ and a Python GUI that controls the input parameters for the simulation and provides a live animation tracking of the simulation process. The simulation console is also enabled for parallel execution; thus large-scale calculations are feasible. TFOx simulator is constructed using the lattice-gas model in which adsorption process is unactivated [38]. Currently, TFOx simulations are restricted to a square lattice, which is a good fit for surfaces configured with four nearest neighbors around each atom, such as the (1 0 0) surface of face centered cubic (fcc) crystals. The surface morphology of a simulation layer that is built from such a square lattice is shown in Fig. 1. Other types of lattices may be needed for different types of surfaces, e.g., the (1 1 1) surface of fcc crystals are better represented by hexagonal lattices. Flexible

294

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

Fig. 1. An example of a 2D lattice. Different adsorption sites are shown, red sites are adatom occupied adsorption sites, hollow sites are the adsorption sites available for adatoms, and the brown balls represent substrate crystal atoms. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

choices of lattice types will be implemented in future TFOx releases as needed. Previous 2D TFOx implementation was limited to the study of adatom motion and island growth in a 2D lattice [8]. Here, we have extended TFOx to 3D growth where oxidation in the direction perpendicular to the interface can be simulated. This feature is a step forward towards a complete film growth simulation. In the rest of this section, we discuss the methods implemented in TFOx in details. Also, to illustrate the new features of the 3D TFOx code, we use the oxidation of Cu(1 0 0) as an example and compare to previous results obtained using 2D TFOx and also to experimental results. Thus, for the remainder of the paper, adatom refers to monoatomic oxygen atoms exclusively. 3.1. Surface diffusion Adatoms are allowed to diffuse or ‘‘jump’’ along the surface sites until they become static (i.e., form a nucleus). Before nucleation, these active adatoms are referred to as free adatoms. A MC step in TFOx is defined as performing a jump attempt for all free adatoms in the system. Thus, m0 in Eq. (5) is the same for all free adatoms whereas Cp can be different due to the differences in the jumping probability per jumping attempt fp. In the lattice representation of KMC as in the TFOx simulator, the diffusion on a square lattice is either orthogonal along the h1 0i and h0 1i directions, or diagonal along the h1 1i direction. To select a jumping direction for site i, a jumping scale Sj based on the eight surrounding sites j is created, then its value is compared to the product (r) of a [0,1) random number 11 and the maximum of surrounding jumping scales:

r ¼ 11 maxfSj g

ð6Þ

A surface diffusion to site j is made if Sj1 < r 6 Sj, where Sj1 and Sj corresponds to the jumping scale of two different lattice sites. If r is less than S0, i.e., the smallest jumping scale among all eight sites, the adatom will stay at its current position. Thus for each jump attempt, there are nine possible outcomes, including the jumping to the eight different directions and the no-jumping event. Double occupancy is not allowed, i.e., jump attempt to the site already occupied by another diffusing adatom will be rejected automatically. To better illustrate this algorithm, we show in Fig. 2(A) how an adatom diffuses on the surface in detail. In this example, the jumping probability per jumping attempt f, which we will simply refer

as the jumping probably from now on, along h1 0i and h0 1i is 0.25, and 0.10 along h1 1i. The different settings of the jumping probabilities along the orthogonal and diagonal directions reflect the different energy barriers along the diffusion pathways. For the eight sites surrounding the adatom, the jumping scale of each site is the summation of the jumping scale from all sites with smaller labels, starting from the top left corner running clockwise:

Sn ¼

n X Sj

ð7Þ

j¼1

where ‘‘n’’ refers to the site index of the 8 surrounding neighbors. For example, site [2,2] has a jumping scale S1 = 0.10, site [3,2] has S2 = 0.35 [=0.10 + 0.25], site [4,2] has S3 = 0.45 [=0.10 + 0.25 + 0.10], site [4,3] has S4 = 0.70 [=0.10 + 0.25 + 0.10 + 0.25] and so on. If the random number 11 generated was 0.74, then the product of this random number and the max jumping scale 1.40 (S8 from site [2,3]) is r = 1.04, which is larger than the 0.8 jumping probability at site [4,4] but smaller than the 1.1 jumping probability at site [3,4]. Hence the atom would jump in the down direction towards site [3,4,25]. The diffusion algorithm that is implemented in TFOx is very similar to the rejection free direct KMC method introduced by Gillespie [39], in which the jumping probability for all nine possible jumping choices, for every adatom is calculated prior the jump attempt. In TFOx, all adatoms that have not nucleated yet are included in the diffusion evaluation process simultaneously at each MC step, i.e., multiple events (multiple adatom jumps) take place in each MC step. This is different from common KMC implementations that allow for only one event at each MC step. The main advantage of the multi-event method is that this allows for faster simulation of the surface diffusion process. However, this approach breaks the mathematical connection between the physical time and the MC time, and is more suitable in the studies where the equilibrium structures are of importance. If the dynamical process is of interest, TFOx allows also the setting of only one free diffusing adatom on the surface, so no other adatoms are deposited until the current adatom becomes static or re-evaporates. In this case, only one event takes place at each MC step, and the time interval s for each MC step follows Poisson process and is given by [7]:

lnð12 Þ

s ¼  P8

j¼0

Cj

ð8Þ

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

295

Fig. 2. An example of jumping, jumping probability gradient and ES effect. (A) As an example of how the jumping probability scale works. (B) The same jumping example using an ES multiplier of 2.0.

where 12 is another random number between [0,1), and Cj is the jumping rate for each possible jumping choice. 3.2. Nucleation and island growth Diffusing adatoms can stick to a lattice site and become static, which corresponds to a nucleation process and island growth. The evaluation for the adatom sticking probability is performed simultaneously with the diffusion process at every MC step. In TFOx, three different types of nucleation attempts can take place: type I), nucleation on a clean surface can occur spontaneously for random diffusing adatoms, or type II) at the event when multiple free adatoms collide with each other, and most commonly type III) nucleation or aggregation can be induced by a pre-existing static adatom. There is no general formula to calculate the sticking probability and experimental data are often absent as well. Thus, different approximations are often considered in KMC programs. In TFOx, the sticking probabilities for type I and type II nucleation processes are usually assigned with small values justified by the agreement between the simulation results and experimental observations. As for the type III nucleation, an implementation scheme that is similar to the bond-counting KMC method [16,17,34] is implemented in TFOx. In this case, the nucleation probability for a candidate lattice site is equal to the sum of the probabilities of first neighbor sites. This scheme is illustrated in Fig. 3(A). The sticking probability for a random nucleation when no static adatom is present is zero (type I nucleation is blocked), and the sticking probability near a static adatom is 0.25 along the h1 0i orthogonal direction and 0.10 along the h1 1i diagonal direction. Thus, the probability for sticking at site [3,3] in this example is 0.60 [= 0.25 + 0.25 + 0.10]. The adatom is considered to stick if the sticking probability is larger than a random number between 0 and 1, or otherwise the nucleation attempt is rejected and the adatom is free for further diffusion. Notice that the nucleation process is irreversible in TFOx, i.e., once an adatom sticks, it will no longer be subject to further evaluation for diffusion again. In 3D simulations, the adatoms are allowed to diffuse and nucleate on top of the oxide islands from the interface layers, which lead to multilayer 3D film growth. For an adatom to be able to stick as part of a second or higher layer of oxidation, it must have four oxide sites that form a square directly beneath it to form a supporting base. The sticking probability on top of an oxide island

can be set differently than the interface layer if desired. The nucleation and island growth activities are independent among the different layers. An example of an adatom sticking on the second layer is shown in Fig. 3(B). The ambient sticking probability is also zero for this layer. Notice that the empty sites directly above the bottom nuclei sites have zero sticking probabilities, due to the sticking probabilities of each layer that is independent from the other layers. Thus, the probability for sticking at the candidate site is 0.35 in this example. 3.3. Ehrlich–Schwöbel (ES) barrier effect on surface diffusion The Ehrlich–Schwöbel (ES) barrier is an essential feature for describing 3D island (e.g., mound) growth [40]. The basic idea is that the 3D island is higher than the ground interface layer, and thus in the vicinity of the island edge, the diffusion activity for adatoms moving towards or away from the island will differ from the ones on the flat surface. An intuitive way for seeing this is that an adatom that approaches the lower side of the step edge, i.e., the step foot, sees a larger number of nearest neighbors compared to an adatom on the terrace. This typically translates into a larger stability or stronger binding energy Eb as schematically illustrated in Fig. 4. On the other hand, an adatom that approaches the upper side of the step, i.e., step top, sees fewer nearest neighbors and this usually leads to a higher diffusion barrier Es [29,30]. While Es hinders the diffusion process in both ascending and descending directions, the additional binding energy Eb gained during the ascending movement can help to offset this elevated diffusion barrier. Both density functional theory (DFT) calculations [41–43], molecular dynamics (MD) simulations [44], and AFM experiments [45] have indicated that the ascending movement across a step edge can play an important role in epitaxial film growth. Thus the overall ES barrier enhances adatom diffusion from terrace to the step. This will be referred to as normal ES barrier. The ES barrier for homoepitaxy film growth is usually of the normal type, whereas this may not be the case in heteropexitaxy process. For example, DFT calculations showed that the binding of O atoms to the lower side of step is weaker than on the upper side of step on Pt(1 1 1) [46]. Similar observations for Xe atoms on Pt(1 1 1) are inferred from STM experiments [47,48]. On the other hand, Xe atoms on Cu(1 1 0) show opposite behavior where Xe is more stable on the lower side of the step than the upper side [49]. The ES barrier is of the inverse type in systems where the

296

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

Fig. 3. (A) An example of an adatom sticking on the interface layer. (B) An example of an adatom sticking on the second layer.

3.4. Potential gradient effect on surface diffusion

Fig. 4. An illustration showing the Ehrlich–Schwöbel effect. Ea is the activation barrier for diffusion on terrace, Eb the additional binding energy at the step foot, and Es is the additional activation barrier or diffusion across the step top cliff. The diffusion described in (A) is more favorable, while that in (B) is less favorable. This phenomenon can be accounted for in TFOx using the ES multiplier.

adatom is more stable on the top rather than the foot of the terrace [50]. In this case, the adatom diffusion trends are expected to be opposite. Most KMC simulations can only deal with the normal ES barrier effect, and typically a biased method of adding Es only to the adatom descending movements is used while leaving the effect from Eb out [36,40]. In the 3D TFOx code, both the normal and inverse ES barriers can be simulated. This is done using a multiplier method, in which the ES barrier is represented by a constant multiplier kES that rescales the jumping probability at the step when adatoms move across ledges. Normal ES barrier is simulated using kES > 1, while the inverse ES effect is simulated with kES < 1. kES = 1 corresponds to the case where the ES barrier is not activated. This idea is illustrated in Fig. 2(B) for a normal ES barrier effect. The same jumping example from Fig. 2(A) is now altered by using a kES = 2.0 where the jumping probability up onto the existing island is a factor of two larger. Using the same random number of 0.74, the product of the random number and the largest jumping scale 1.85 (from site [2,3]) is r = 1.37, which is larger than the jumping probability of 1.35 at site [2,4] and smaller than the jumping probability of 1.85 at site [2,3]. Thus, it is likely that the adatom will now jump to the left direction towards the oxide island site [2,3].

When an adatom ‘‘sticks’’ to a vacant site and forms a nucleus, strain on the surrounding substrate atoms starts to build, and the movement of other diffusing adatoms in its vicinity can be altered. Thus, thin film growth morphology may be influenced significantly from the substrate-mediated interactions [51]. The magnitude of this effect varies depending on the system and the temperature [51]. There are two common ways of simulating the elastic lattice strain effect. The first approach is called the ‘‘ball-and-spring’’ model, in which the short-ranged interactions are included using the bond-counting approach and the long-ranged interactions are evaluated from imaginary Hookean springs [15]. The second approach is the ‘‘off-lattice’’ method where effects of the strain are incorporated using a new potential term that describes the substrate-adsorbate interactions [27,28]. TOFx adapts the second approach, where strain effects are accounted for using a potential gradient function that describes the strength of the elastic interaction effect on different sites around a nucleus. The jumping probability at each site surrounding a nucleus is rescaled by multiplying the corresponding potential gradient parameter with the original jumping probability. At different temperatures, along with different event attempt frequency m0, different forms of potential gradient function are used. Thus the potential gradient function plays a key role to bridge the relationship between temperature and the jumping probability f. The choice of potential gradient function is critical for describing the correct island morphology during film growth. With various oxide island shapes observed experimentally [5], there are different potential gradients representing different energetic or kinetic mechanisms under different temperature that may dominate the morphology variation during the simulation. We previously developed two sets of 2D potential gradient for the Cu(1 0 0) surface oxidation process, namely the 4-fold and the 2-fold potential gradient [8]. The development of the 4-fold potential gradient function is inspired by the experimental observation that the Cu2O(1 1 0) surface has the slowest growth rate during oxidation at moderate to high temperatures [52] and thus determines the final shape of the expanding oxide island [53,54]. The 4-fold potential gradient function has the form:

G4fold ¼

1þ

1 jija þ jjja

!c b

ð9Þ

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

where [i, j] represents the lattice site coordinate, a is the ‘‘orientation sensitivity’’ parameter, b is the ‘‘distance sensitivity’’ parameter, and c is the exponent of potential gradient function. For the ‘‘deep well’’ 4-fold potential gradient function used in this paper, a, b, and c are set to 0.6, 0.9 and 2.0, as suggested before [8]. The 4-fold potential gradient function enhances the jumping probability along the h1 0i orthogonal direction and leads to the growth of a square shape island. The 2-fold potential gradient function is applicable when the strain region introduced during island growth is inhomogeneous and calls for an asymmetric strain potential gradient. These conditions occur at relatively low temperatures during the Cu(1 0 0) surface oxidation, in which the lattice parameter difference between copper substrate and the Cu2O island is relatively large. The 2-fold potential gradient function has the form:

G2fold ¼ xðe½i; jÞ

2lþ1

þv

ð10Þ

where x is the potential gradient ‘‘scaling constant’’, l is the exponent of the potential gradient, and v is the ‘‘potential gradient offset’’. In this study, we used x, l, v equal to 1.25, 2.0 and 0.8, as used before [8]. e is the strain associated with the displacement of the substrate atoms from their equilibrium lattice positions due to the presence of the adatom at site [i, j], which has been discussed in more detail elsewhere [8]. The 2-fold potential gradient function enhances the jumping probability along the h1 1i diagonal direction and leads to the growth of a rod-shape island. Lattice strains under different temperatures would call for the development of other variations of potential gradient, this is however not included in this work. Higher-orders of surface reconstruction are not explored indepth in this study. Currently, two extensions to the existing 2D potential gradient are implemented and tested in 3D TFOx, namely cylindrical and spherical extensions of the 2D potential gradients. The difference between the different extension methods is presented in Fig. 5. Fig. 5(A) shows the original 2D gradient, Fig. 5(B) shows the cylindrically extended gradient, and Fig. 5(C) shows the spherically extended gradient. The cylindrically extended gradient applies the same set of potential gradient parameters to all layers:

Gði;j;kÞ ¼ Gði;j;0Þ

ð11Þ

G(i,j,k) is the potential gradient at site [i, j, k] defined by the [i, j] coordinate on XY plan and the k coordinate on the Z axis corresponding to the plane height. The spherically extended gradient gradually decreases the magnitude of the potential gradient and reduces its effect when moving away from the interface layer. The general formula for implementing the spherically extended potential is:

Gði;j;kÞ ¼ Gðiþ1;j;k1Þ cosðhÞ2 þ Gði;jþ1;k1Þ sin ðhÞ2

ð12Þ

where h is the polar angle on the XY plan determined by the vector from the nucleus center to site [i, j] with respect to the X axis. The potential gradient on each site is built upon the previous layer’s outer site, thus the potential gradient effect shrinks with increasing layer height. We will discuss the performance of these two different extended potentials in Section 5. 4. Calculation details We previously demonstrated using the 2D TFOx simulator that the potential gradient that affects the atom jumping and attachment probabilities along different lattice directions is of key importance in the growth of the island morphology [8,25]. Here, we investigate the effects of 3D potential gradients and step barriers, and compare the results to similar simulations that are confined

297

to 2D. We choose the oxidation of the Cu(1 0 0) surface as an example model for this paper. A preset nucleation site is chosen at the center of a 200  200 lattice where each site represents a single adsorption site. Multiple adsorption sites are not taken into account, this is a good approximation in the low oxygen concentration limit where only the lowest energy adsorption sites are occupied. For Cu(1 0 0), DFT calculations showed that the hollow sites are the most stable adsorption site for oxygen adatoms [55]. Periodic boundary conditions in the surface plane are used for the simulation on the X and Y directions. The spontaneous nucleation probability for oxygen adatom at the interface layer is set to zero, thus oxide island growth is only feasible near the preset nucleation site. A delicate choice of simulation parameters require the rate tables for all elementary steps to be obtained first, since such a table is not available yet, here the simulation parameters are set to be the same as in the 2D TFOx simulation which shows satisfactory film growth behavior compared to experiment [8]. This choice additionally makes the comparison between results obtained using the 2D and 3D TFOx simulations straightforward, and servers the main purpose of demonstrating the capability and features of 3D TFOx program. More specifically, the jumping probability along the h1 0i and h0 1i orthogonal directions is set to 1.0 and that along the h1 1i diagonal directions is 0.7. The sticking probability is set at 0.05 and 0.1 along h1 0i and h0 1i, and h1 1i directions for both the interface and island layers. The jump attempt frequency m0 is 1011/s and the ratio of the deposition flux to the jump attempt frequency is 2.5  104 ML/loop (monolayers per simulation loop), which represents a real-time deposition flux of 2.5  107 ML/s. By the end of the simulation 3000 oxygen atoms are nucleated, corresponding to 0.075 ML oxide island coverage. This low-coverage limit justifies the inclusion of only one adsorption site. 5. Results and discussion 5.1. Different 3D extensions of the potential gradient The morphology differences between the results obtained using the traditional 2D simulation and the new 3D simulations employing three variants of potential gradient are shown in Fig. 6. In the second column of Fig. 6 labeled with ‘‘None’’, the potential gradient corresponds to the case where the traditional 2D potential gradient is only applied on the interface layer while leaving the extra layers unaffected. For all the cases, we note that the islands in the 3D simulation occupy a smaller base area that can be explained because we nucleated a fixed number of oxygen atoms in all the simulations, and thus the extra dimension that allows the stacking of oxygen atoms would naturally reduce the base area. Additionally, for the 4-fold potential gradient shown in the first row of Fig. 6 where no 3D extension is applied, the resulting island has a similar square shape compared to the 2D simulation result. This is to be expected considering that after the Cu2O island has been formed on the interface layer, further growth on the layers above this island would simply covers up the island surface evenly, as there is no more potential gradient to induce the island growth behavior. The large portion of red color coverage indicates that the growth of the second layer Cu2O film is preferable. Higher layer film growth is less probable, as indicated by the limited green, blue, etc. color coverage. This result is better demonstrated in the 3D view of the lattice in Fig. 7(B), in which many Cu2O islands are seen on the second layer. The third layer also contains an appreciable amount of Cu2O, sites but not as many as the second layer. The number of Cu2O islands on higher layers is much smaller. When the potential gradient is extended to 3D space, more refined details can be observed upon multilayer oxidation.

298

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

Fig. 5. A 3D visualization of the different options available for extension of the potential gradient. Color represents the magnitude and the arrow represents the direction in gradient space. (A) Original 2D gradient. (B) Cylindrically extended gradient. (C) Spherically extended gradient. (D) Top-down view of a cylindrically extended gradient. (E) Slightly elevated view of a spherically extended gradient. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Top view of the potential gradient effects. The color representation from zero layer of island growth to higher layers follows the order of: white, black, red, green, blue, yellow, purple, cyan, dark grey, light grey, dark red, dark green, dark blue, dark yellow, dark purple, dark cyan. (Not all colors are presenting in this figure.) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Both the extended 2D and the spherical 3D potential gradients lead to the growth of a diamond shaped island with a concave square base. The film growth towards higher layers beyond the second layer is enhanced, as can be seen from the larger portion of green and blue color coverage, which can also be seen in the 3D view of the lattice in Fig. 7(C) and (D). This result is consistent with the fact that the 4-fold potential gradient applied here increases the diffusion activity around the static adatoms, thus promoting island growth with an increased height. For the extended 2D potential gradient, the same form of the traditional 2D potential gradient is applied to the multiple layers uniformly [see Fig. 5(B) and (D)], thus each layer is affected by the potential gradient in the same fashion. It should be noted that the 4-fold potential gradient enhances the diffusion activity along the orthogonal directions. When the growth is limited to 2D, the preset nucleus site at the center of the lattice promotes the adatom diffusion around it along the orthogonal axis. As a result, nucleation along the orthogonal axis is enhanced, while

the nucleation along the diagonal direction occurs at a reduced rate, resulting in a square-shaped island. Once the 2D growth limitation is removed, the adatoms can diffuse on top of the Cu2O island layers. In this case, less adatoms are available for diffusion at the interface layer, further weakening the island growth along the diagonal direction, resulting in a concave square island base. The images in Fig. 6 show that the higher layers of the island (represented by the green and blue colors) adapt a cross shape, which indicates that the enhancement of the orthogonal direction diffusion is amplified as the island layers grow higher. This result can be understood because the enhanced orthogonal direction pattern is applied to all layers, and thus the potential gradient effect that promotes the orthogonal directional island growth is accumulated towards higher layers. For the spherical 3D potential gradient, the effect of the potential gradient is also applied to higher layers, however, in a non-uniform fashion. When the height of the layer increases, the strength of the potential gradient effect decreases [Fig. 5(C) and (E)].

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

299

Fig. 7. 3D top view of the 4-fold potential gradient effects. (A) 2D simulation with a planar oxide island. (B) ‘‘None’’ 3D potential gradient where oxide island grows vertically uniformly. (C) Extended 2D potential gradient that results in an uneven oxide island height and forms a cross-shape pattern. (D) Spherical 3D potential gradient that also results with a cross-shape pattern but more homogenously in terms of the oxide island height. Yellow balls represent the copper atoms and red balls represent the Cu2O nuclei. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

From Fig. 6, we see that the result from the spherical 3D potential gradient is quite similar to that of the extended 2D potential gradient, while the cross pattern of the higher layers is more regular, illustrated by the more evenly distributed Cu2O sites at the higher layer in Fig. 7(D) comparing to Fig. 7(C). This is because the reduced potential gradient strength can offset the layer accumulation effect when growing towards higher layers. Thus, the island growth at higher layers is more uniform along the orthogonal direction, instead of creating ‘‘hot spots’’ at which the Cu2O island is much higher than the neighbors. For the 2-fold potential that enhances the diffusion activity along the diagonal direction, the behavior of the three different variants of the potential gradient is quite similar to the ones from the 4-fold potential. When no 3D potential gradient expansion is used, the rod island grows into a similar shape compared to the 2D simulation. Multi-layer growth is mostly limited to the second and third layers, as seem by the large red and green color coverage in Fig. 6 and the 3D view of the lattice in Fig. 8(B). Upon applying the 3D potential gradient, the growth of the higher layers is enhanced. The extended 2D potential gradient leads to an enhanced diagonal rod growth on the higher layers with several ‘‘hot spots’’ at which the Cu2O island is much taller than the neighbor [Fig. 8(C)], while the spherical 3D potential gradient leads to multi-layer growth more evenly [Fig. 8(D)]. Both the 3D extended 2-fold and 4-fold potential gradients produce Cu2O island morphologies that are in good agreement with the experimental results shown in Fig. 9 [5]. The improvements realized upon using the 3D extension form of the potential gradient is particularly apparent for the 4-fold potential gradient. Compared to Fig. 9(B), the cross pattern towards higher layer Cu2O island growth was missing from the 2D simulation, while it can be reproduced in the 3D simulation. This finding shows that the spherical 3D potential gradient is a better choice as it agrees with the fact that the potential gradient effect from a nucleus site decreases with distance. Additionally, the simulation results obtained with the spherical 3D potential gradient shows smoother growth towards higher oxide island layers, which is consistent with experimental observation [5].

5.2. ES barrier effect As introduced in Section 3.3, the ES barrier effect is incorporated in TFOx through the ES multiplier method. The effect of the ES barrier on the 4-fold and 2-fold potential gradient induced film growth is shown in Fig. 10. The middle right column corresponds to no ES barrier as kES = 1. When kES > 1 is used for normal ES barrier effect, the island grows higher as can be seen from the increased maximum height of the island with increasing kES. This result is consistent with the enhanced jump probability towards the oxide island at kES > 1. On the other hand, when kES < 1 is used for inverse ES barrier effect, the island grows lower as can be seen from the decreased maximum height of the island with decreasing kES. This is because of the demoted jump probability towards the oxide island for kES < 1. The influence of the ES multiplier is better represented when combined with the 2-fold potential gradient, as shown in the second row in Fig. 10. Because the base area of the rod island is smaller than that of the diamond island, the rod island height increment is more apparent when a larger kES is applied. Also, the length of the rod island decreases at the same time because the number of the nucleated oxygen atoms is held constant. With a small kES (103), the rod island is mainly composed of single and double oxide layers as represented by the large portion of black and red colors. This result indicates that the diffusion towards the island is stronger, thus higher layer growth is greatly depreciated. Similar trends can be seen for the diamond island too. Given that the Cu2O island observed in TEM experiments shown in Fig. 9 actually consists of multiple layers, we can predict from our KMC simulations that oxygen atoms on Cu(1 0 0) surface should comply with a normal ES barrier effect. 5.3. Discussion Previously-developed KMC simulators used algorithms confined within the scope of very specific models, and thus are not versatile [11,13,15–17]. The ‘‘off-lattice’’ approach of including elastic lattice stain effect was only explored in a one dimensional

300

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

Fig. 8. Similar to Fig. 7 except we show 3D side views of the island growth from the 2-fold potential gradient effects.

600 °C

150 nm

(A)

750 °C

300 nm

(B)

Fig. 9. The morphology of Cu2O islands formed during in situ oxidation of Cu(1 0 0) at a oxidation partial pressure of 5  104 Torr and oxidation temperatures of: (A) 600 °C; (B) 750 °C [5].

Fig. 10. ES Barrier effects for sample 2-fold and 4-fold potential gradients. Color scheme is the same as in Fig. 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

model previously [28], and our effort here has extended this approach in TFOx to 3D with ‘‘on-lattice’’ model. Additionally, the ES barrier effect has been overlooked by earlier KMC implementations, with the recent exceptions of Wu et al. [36] and Leal et al. [40], where the ES barrier were both applied by introducing an additional activation barrier Es for adatom descending movement while the effect from Eb is ignored.

In the study of Leal et al., a height-dependent ES barrier model is constructed, where the adatom movement across multilayer steps is modeled through a one-dimensional random walk algorithm. This study found that the introduction of ES barrier can significantly increases the surface roughness during film growth and that the aspect ratio (height/width) of the 3D island is found to increase with the ES barrier [40]. Interestingly, these findings are

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

comparable to our results where we have observed that the oxide island height increases with a larger ES multiplier as well. Compared to the implementation of Leal et al., the ES multiplier method in TFOx is more robust since it does not require the counting of lateral bonds on the fly, nor it requires extra computational effort for the one-dimensional random walk along the ledge. Even though the current ES multiplier method in TFOx does not account for the effect of edge height on the ES barrier, an extension for a height-dependent ES multiplier can be incorporated in the future with minimal computational cost, similar to the extension of the potential gradient from 2D to 3D. To our knowledge, TFOx is the first program that includes both effects of elastic lattice stain and ES barrier in a coherent 3D KMC framework. Extending the 2-fold and 4-fold potential gradient function to 3D and the inclusion of ES barrier have reproduced a surface morphology in agreement with experimental observations for the Cu(1 0 0) surface oxidation. This agreement gives confidence that TFOx can be applied to a broad range of nano science fields, such as the modeling of 3D quantum dots and nano wires. Previous experimental observations [23,24,56] and first-principles DFT calculations [55,57–59] have shown that oxygen atom can diffuse into the subsurface layers of the Cu film and that the growth of the Cu2O island takes place both below and above the metal surface. However, the current TFOx program only accounts for the growth in the upward direction. It would be necessary to include the 3D diffusion of the adatom into layers below the interface layer and allow the nucleation and island growth towards the lower layers. One possible approach to alleviate the difficulty of interlayer adatom diffusion is to apply an artificial inverse ES barrier effect using kES < 1 for the subsurface layers to limit the diffusion rate. Thus, two kES values can be incorporated, where a large value is used to enhance the upwards island growth and a smaller value for the downwards island growth. As an important bridge between microscopic and macroscopic modeling, KMC simulations need information provided by more accurate atomistic scale calculations such as DFT, such as a rate table for the diffusion process of adatoms. Thus, for a more accurate understanding of oxidation, atomistic calculations that improve the accuracy of the TFOx input parameters are needed. Also, physical underpinnings of the symmetric potential strain call for improvement in interface energy calculations, from which more accurate potential gradient functions for TFOx can be developed. Finally, as oxidation often takes place on ill-defined, rough surfaces, with either added dopants or existing impurities, incorporation of these surface inhomogeneities into 3D TFOx is also needed to truly achieve a first principles multi-scale simulation of oxidation to compare with experiments.

found that the ES multiplier kES can strongly influence the oxide island height. When kES > 1 is used, normal ES barrier effect is realized and the diffusion towards the oxide island is enhanced resulting in higher oxide island growth. On the hand, kES < 1 which corresponds to inverse ES barrier can inhibit the height of the oxide island. By comparing the results obtained by TFOx to TEM measurements, it is concluded that a normal ES barrier effect is suggested for oxygen atom on Cu(1 0 0) surface. Overall, we developed a versatile 3D KMC simulator that is capable of capturing many important features during a thin film oxidation process, including the lattice mismatch strain effect, the Ehrlich–Schwöbel barrier, and orientation-biased epitaxial growth. The incorporation of these key elements into a single program makes TFOx the first KMC simulator that can successfully describe the unique 3D oxide island morphologies under different film oxidation conditions, Furthermore, TFOx can describe many different systems by choosing different lattice parameters, jumping probabilities, sticking probabilities, and potential gradient functions. Though the motivation for this work is based on the in situ TEM work of Cu oxidation, we anticipate that the results of this body of work will have major impact in other fields, such as nanosynthesis, catalysis, extreme fuel reactions, as well as oxidation and corrosion. Acknowledgements This work is supported by the United States Department of Energy (#DOE BES-ER46446 and #DOE DE-FG02-03ER1547), and by the internal funding through Swanson School of Engineering, University of Pittsburgh. We are grateful for computing time provided in part by the University of Pittsburgh Center of Simulations and Modeling, and by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (#OCI-1053575). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

6. Conclusion Through enabling the diffusion and sticking of adatom on top of the interface layer, we have extended the TFOx code from two to three dimensions. Two major aspects are accounted for in this expansion. First, the extension of the potential gradient field in the 3D space, where three different methods are implemented, namely ‘‘None’’, extended 2D, and spherical 3D. Simulation results show the choice of the potential gradient has a strong influence on the morphology of the resulting oxide island. When a 3D potential gradient is applied, the growth towards higher oxide island is promoted. Also, the potential gradient effect is accumulated from lower to higher layers, using a spherical 3D potential gradient that can offset this accumulation effect to an extent that results in oxide island morphologies agreeing better with experimental observations. The ES barrier effect near the edge of an oxide island is implemented through a relatively simple and robust scheme. We

301

[13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

N.P. Padture, M. Gell, E.H. Jordan, Science 296 (5566) (2002) 280. C.L. Hill, C.M. Prossermccartha, Chem. Rev. 143 (1995) 407. C.T. Campbell, Surf. Sci. Rep. 27 (1–3) (1997) 1. H. Over, Y.D. Kim, A.P. Seitsonen, S. Wendt, E. Lundgren, M. Schmid, P. Varga, A. Morgante, G. Ertl, Science 287 (5457) (2000) 1474. G.W. Zhou, J.C. Yang, Appl. Surf. Sci. 210 (3–4) (2003) 165. L.J. Broadbelt, R.Q. Snurr, Appl. Catal. A- Gen. 200 (1–2) (2000) 23. K.A. Fichthorn, W.H. Weinberg, J. Chem. Phys. 95 (2) (1991) 1090. X.T. Han, R. McAfee, J.C. Yang, J. Comput. Theor. Nanosci. 5 (1) (2008) 117. M.L. Guerra, M.A. Novotny, H. Watanabe, N. Ito, Phys. Rev. E 79 (2) (2009) 026706. J.J. Lukkien, J.P.L. Segers, P.A.J. Hilbers, R.J. Gelten, A.P.J. Jansen, Phys. Rev. E 58 (2) (1998) 2598. A.S. McLeod, Catal. Today 53 (2) (1999) 289. H. Persson, P. Thormahlen, V.P. Zhdanov, B. Kasemo, Catal. Today 53 (2) (1999) 273. C. Ghosh, A. Kara, T.S. Rahman, Surf. Sci. 601 (15) (2007) 3159. M. Stamatakis, D.G. Vlachos, J. Chem. Phys. 134 (21) (2011) 214115. P.P. Petrov, W. Miller, Comp. Mater. Sci. 60 (2012) 176. P.F. Zhang, X.P. Zheng, S.P. Wu, D.Y. He, Comp. Mater. Sci. 30 (3–4) (2004) 331. X.P. Zheng, P.F. Zhang, Comp. Mater. Sci. 50 (1) (2010) 6. M. Stamatakis, D.G. Vlachos, ACS Catal. 2 (12) (2012) 2648. J.A. Venables, Philos. Mag. 27 (3) (1973) 697. J.A. Venables, G.D.T. Spiller, M. Hanbucken, Rep. Prog. Phys. 47 (4) (1984) 399. G.W. Jones, J.M. Marcano, J.K. Norskov, J.A. Venables, Phys. Rev. Lett. 65 (26) (1990) 3317. J.A. Venables, Surf. Sci. 299 (1–3) (1994) 798. J.C. Yang, M. Yeadon, B. Kolasa, J.M. Gibson, Appl. Phys. Lett. 70 (26) (1997) 3522. J.C. Yang, M. Yeadon, B. Kolasa, J.M. Gibson, Scr. Mater. 38 (8) (1998) 1237. X. Han, R. McAfee, J.C. Yang, Multidiscip. Model. Mater. Struct. 3 (1) (2007) 12. T.A. Witten Jr., L.M. Sander, Phys. Rev. Lett. 47 (19) (1981) 1400. M. Schroeder, D.E. Wolf, Surf. Sci. 375 (1) (1997) 129. F. Much, M. Biehl, Europhys. Lett. 63 (1) (2003) 14. G. Ehrlich, F.G. Hudda, J. Chem. Phys. 44 (3) (1966) 1039. R.L. Schwoebe, E.J. Shipsey, J. Appl. Phys. 37 (10) (1966) 3682. N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (6) (1953) 1087.

302

Q. Zhu et al. / Computational Materials Science 91 (2014) 292–302

[32] A.R. Leach, Molecular Modelling: Principles and Applications, second ed., vol. 414, Prentice Hall, Harlow, England: New York, 2001. [33] A. Chatterjee, D.G. Vlachos, J. Comput-Aided Mater. 14 (2) (2007) 253. [34] L.G. Wang, P. Clancy, Surf. Sci. 473 (1–2) (2001) 25. [35] L. Xu, C.T. Campbell, H. Jonsson, G. Henkelman, Surf. Sci. 601 (14) (2007) 3133. [36] F.M. Wu, H.J. Lu, Y.Z. Fang, S.H. Huang, Commun. Theor. Phys. 50 (1) (2008) 231. [37] L.A. Li, G.W. Zhou, Surf. Sci. 605 (1–2) (2011) 54. [38] E.V. Albano, K. Binder, D.W. Heermann, W. Paul, Surf. Sci. 223 (1–2) (1989) 151. [39] D.T. Gillespie, J. Comput. Phys. 22 (4) (1976) 403. [40] F.F. Leal, S.C. Ferreira, S.O. Ferreira, J. Phys.: Condens. Matter. 23 (29) (2011) 292201. [41] K. Fichthorn, M. Scheffler, Nature 429 (6992) (2004) 617. [42] W.G. Zhu, F.B. de Mongeot, U. Valbusa, E.G. Wang, Z.Y. Zhang, Phys. Rev. Lett. 92 (10) (2004) 106102. [43] Y. Tiwary, K.A. Fichthorn, Phys. Rev. B 81 (19) (2010) 195421. [44] H.L. Yang, Q. Sun, Z.Y. Zhang, Y. Jia, Phys. Rev. B 76 (11) (2007) 115417. [45] F.B. de Mongeot, W.G. Zhu, A. Molle, R. Buzio, C. Boragno, U. Valbusa, E.G. Wang, Z.Y. Zhang, Phys. Rev. Lett. 91 (1) (2003) 016102.

[46] T. Ogawa, A. Kuwabara, C.A.J. Fisher, H. Moriwake, T. Miwa, J. Phys. Chem. C 117 (19) (2013) 9772. [47] P. Zeppenfeld, S. Horch, G. Comsa, Phys. Rev. Lett. 73 (9) (1994) 1259. [48] S. Horch, P. Zeppenfeld, G. Comsa, Appl. Phys. A— Mater. 60 (2) (1995) 147. [49] M. Dienwiebel, P. Zeppenfeld, J. Einfeld, G. Comsa, F. Picaud, C. Ramseyer, C. Girardet, Surf. Sci. 446 (1–2) (2000) L113. [50] F. Slanina, J. Krug, M. Kotrla, Phys. Rev. E 71 (4) (2005) 041605. [51] M.L. Merrick, W.W. Luo, K.A. Fichthorn, Prog. Surf. Sci. 72 (5–8) (2003) 117. [52] G. Zhou, Dynamics of Copper Oxidation Investigated by in situ UHV-TEM, University of Pittsburgh: PhD Thesis, 2004. [53] R.H. Doremus, Rates of phase transformations. (1985), Orlando: Academic Press. 176. [54] Z.Y. Zhang, M.G. Lagally, Science 276 (5311) (1997) 377. [55] W.A. Saidi, M. Lee, L. Li, G.W. Zhou, A.J.H. McGaughey, Phys. Rev. B 86 (24) (2012) 245429. [56] J.C. Yang, G.W. Zhou, Micron 43 (11) (2012) 1195. [57] M. Lee, A.J.H. McGaughey, Surf. Sci. 603 (24) (2009) 3404. [58] M. Lee, A.J.H. McGaughey, Surf. Sci. 604 (17–18) (2010) 1425. [59] M. Lee, A.J.H. McGaughey, Phys. Rev. B. 83 (16) (2011) 165447.