Double step structure and meandering due to the many body interaction at GaN(0001) surface in N-rich conditions
Magdalena A. Załuska-Kotur, Filip Krzyżewski, Stanisław Krukowski
aa r X i v : . [ c ond - m a t . m t r l - s c i ] N ov , Double step structure and meandering due to the many bodyinteraction at GaN(0001) surface in N-rich conditions.
Magdalena A. Za luska–Kotur ∗ Institute of Physics, Polish Academy of Sciences,Al. Lotnik´ow 32/46, 02-668 Warsaw,Poland and Faculty of Mathematics and Natural Sciences,Card. Stefan Wyszynski University,ul Dewajtis 5, 01-815 Warsaw, Poland
Filip Krzy˙zewski † Institute of Physics, Polish Academy of Sciences,Al. Lotnik´ow 32/46, 02-668 Warsaw, Poland
Stanis law Krukowski ‡ High Pressure Research Center, Polish Academy of Sciences,ul. Soko lowska 29/37,01-142 Warsaw,Poland and Interdisciplinary Centre for Materials Modeling,Warsaw University, Pawi´nskiego 5a, 02-106 Warsaw, Poland bstract Growth of gallium nitride on GaN(0001) surface is modeled by Monte Carlo method. Simulatedgrowth is conducted in N-rich conditions, hence it is controlled by Ga atoms surface diffusion. It isshown that dominating four-body interactions of Ga atoms can cause step flow anisotropy. KineticMonte Carlo simulations show that parallel steps with periodic boundary conditions form doubleterrace structures, whereas initially V -shaped parallel step train initially bends and then everysecond step moves forward, building regular, stationary ordering as observed during MOVPE orHVPE growth of GaN layers. These two phenomena recover surface meandered pair step patternobserved, since 1953, on many semiconductor surfaces, such as SiC, Si or GaN. Change of terracewidth or step orientation particle diffusion jump barriers leads either to step meandering or surfaceroughening. Additionally it is shown that step behavior changes with the Schwoebel barrier height.Furthermore, simulations under conditions corresponding to very high external particle flux resultin triangular islands grown at the terraces. All structures, emerging in the simulations, have theircorresponding cases in the experimental results. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] . INTRODUCTION Crystal growth remains subject of intensive study, connecting research on fundamentalunderstanding of the pattern selection, or the kinetic properties and the technical aspectsof synthesizing useful materials and structures for present day advanced information tech-nologies, optoelectronics, medicine, or molecular sensing, etc. The shape selection growthphenomena are intimately related to the morphological instability of the flat or curved sur-faces of growing crystals [1–3]. In the macroscale, the variety of the growth shapes extendsfrom flat surfaces, cellular patterns [4] to dendritic patterns [5], and finally to astonishing va-riety of snow crystal shapes [6, 7]. The rich variety of the shapes is characterized by differentsizes, starting from the above mentioned macroscale patters, to diversity of microcrystals,e.g. spherullitic crystals, and ultimately to nanocrystals, such as fullerenes [8], nanotubes[9]or graphene [10]. The shape selection phenomenon encompasses different dimensionality,extending from these three dimensional structures to the various 2-dimensional phenomena,observed on flat crystal surfaces, such as step repulsion, step meandering, and also stepbunching. These phenomena, though observed for very long time, still escape from basicunderstanding. Notably, step meandering, in which the pairs of close steps are bend, in-terchanging with the other pairs for different direction, observed more than half centuryago on silicon carbide surfaces [11], was subsequently reported to exist on many other semi-conductor surfaces[12], such as silicon[13] or gallium nitride [14]. It still remains unsolvedpuzzle, despite many efforts to understand its dynamics. This work is devoted to elucidationof this phenomenon in the specific case of the growth of gallium nitride by metal organicvapor phase epitaxy (MOVPE) or hydride vapor phase epitaxy (HVPE) on polar GaN(0001)surface [15].A model of GaN(0001) surface growing in nitrogen rich condition is proposed. Since,at the standard growth conditions for both MOVPE and HVPE, the system is kept isnitrogen-rich state, Ga atoms motion is considered explicitly as they control the crystalgrowth process. It is assumed that N atoms are at such abundance at the surface, that theyfill lattice sites instantaneously after the incorporation of Ga atoms in the neighboring sites.The basic assumption, used in the model, is that complete tetrahedral Ga structure aroundN atom has lower energy than mere sum of two body Ga–Ga interactions connected via singleN atom. It is shown that the presence many-body forces between Ga atoms plays dominant3ole in the crystal growth dynamics. Strongly connected tetrahedral structures lead to thesubstantial difference in the speed of the step flow. Growth of wurtzite GaN crystal ismodeled by kinetic Monte Carlo (MC) simulations in which stationary step patterns areidentified. Using results of the MC simulations we demonstrate that selection of variousstationary surface step patterns is realized, depending on the assumed interaction energiesand the kinetic growth parameters.In the proposed model bulk gallium nitride crystal is built by subsequent deposition ofconsecutive Ga-N atomic layers on GaN(0001) surface, forming stationary step patternsduring growth simulations. Initially straight parallel steps, which are oriented along one ofthe main crystallographic directions, form double terrace structure when tightly fixed byperiodic boundary conditions. In the case when the parallel steps are initially bent at 120 o forming the V -shape form, they change their orientation to more sharp one first and thenthey build double step structure. In the case of initial, larger scale, hexagonal structures,the tendency of step reorientation leads to the characteristic surface wavy pattern (stepmeandering). In addition to such evolution it is also shown that for high Schwoebel barrierthe steps become unstable and such unstable step pattern ultimately evolves to random,rough surface. II. THE MODEL
In the MOVPE and HVPE growth of GaN, the fluxes of the reactive species at GaN(0001)surface, are not equivalent or even comparable. Ammonia flux is typically two (HVPE), orthree (MOVPE) orders of magnitude larger than flux of gallium containing molecules. Atthe typical growth temperature which, for both methods, is close to 1300K, Ga transportingagents, TMG or GaCl, are highly unstable. They decompose at the GaN surface easily,leaving Ga atoms. Nitrogen transporting agent, ammonia is also unstable. Differently fromgallium transporting agents, ammonia decays in the gas phase that creates highly nonreactivemolecular nitrogen and hydrogen, reducing availability of the active nitrogen at the surface.Thus larger ammonia flux is needed to flush out these nonreactive decomposition products.It was shown recently that ammonia is attached at GaN (0001) surface in molecularform without energy barrier [23, 29, 30]. The fast adsorption process of NH molecules isaccompanied by desorption of hydrogen in molecular form (H ). Thus the entire surface is4 IG. 1: (color online) Model of GaN crystal. Two types of large circles correspond to Ga atoms attwo terrace types, differing by step velocity. From bottom to top terrace heights increase by onelayer at each step. Positions of N atoms are marked by dots. b + E b + E + E ( J ) E b + E + E ( J ) Eb FIG. 2: (color online) Potential felt by a particle close to the two step types. Steps go up fromright to left side. Particle that attach step from upper step has to overcome Schwoebel barrier.Steps differ by bonding energy, what decides about the ability of a particle to detach. covered by NH radicals, at which, any Ga atom, arriving at the surface, can be attached[31].Thus the creation of GaN layer proceeds via accumulation of Ga atoms on top of layers ofthese radicals.Due to crystallographic structure of wurtzite, two consecutive dia-atomic layers differby the location of sites, which are occupied by Ga atoms and by Ga-N bond orientations.Third layer returns to the original position as it is reflected by stacking ABABAB sequenceof wurtzite. Due to this the two consecutive surface steps have different atomic structure.As illustrated in Fig 1 three bonds connecting Ga atom with the nearest N atoms within5ne layer are rotated by 60 o with respect to the bonds in the consecutive Ga layer. As aconsequence of this crystallographic buildup, every second step has different atomic struc-ture. Rotation of the step orientation by 60 o exchanges that two kinds of step. Such crys-tallographic arrangement can result in different step motion under given supersaturation,depending on the step precedence and also on its orientation.More detailed analysis of the structure of GaN surface indicates that Ga adatoms, at bothtypes of the step, have the same number of first neighbors (atoms N) and second neighbors(atoms Ga). It follows then that the step difference cannot be demonstrated within themodel, limited to additive two-body interactions only. It turns out that the main differencebetween both steps lies in the relative orientation of the second neighbors of Ga atom, i.e.the closest Ga adatoms, those from the layer below and those within the layer where Gaatom is located when attached at the surface. The simplest way to model such difference isto introduce four-body interaction between Ga atoms. GaN wurtzite lattice can be build oftetrahedrons of Ga atoms, centered on N atoms. Three atoms of such Ga tetrahedron areshown in Fig 1, surrounding any displayed N atom there and the fourth Ga atom is locatedin the layer underneath (not shown there). Alternatively, the N based tetrahedrons withGa atom inside can be also created, recovering GaN lattice, but this will be not used inthe present paper as we are interested in the growth model in which motion of Ga atoms isgrowth controlling factor.In four-body interaction, each full Ga tetrahedron, bonded by single N atom inside, hasdifferent energy than the simple sum of two-body Ga-Ga bonds, connected by N atom (i.e.when the third Ga atom is missing). Note that by assumption of overwhelming presenceof nitrogen, we simplified our model, reducing the role of N atoms to the links connectingGa atoms in the lattice. Also we assumed that the Ga atoms underneath is always presentwhile these in Ga layer above are always empty. Accordingly our system is described bymodeling the energy and dynamics of the single layer Ga atoms only.Accordingly to the described buildup of GaN lattice from the Ga tetrahedrons centeredon the N atoms, we calculate the given tetrahedron contribution to the energy of Ga surfaceatom, that depends on the number Ga neighbors, connected by N atom in the tetrahedroncenter, in the following way: ǫ i = , when tetrahedron has all atoms; rη, when tetrahedron has empty sites, (1)6here η is a number of Ga occupied neighboring sites, belonging to a given tetrahedronand r describes the relative strength of the four-body and the two-body interactions in thesystem. When r = 1 two body Ga-Ga interactions sum up to the value characteristic forfully occupied tetrahedron i.e. no additional four-body Ga interactions are present in thesystem. When r < α ( J ) = J X i =1 ǫ i . (2)where parameter J scales energy of bonds, and the sum runs over four tetrahedrons, sur-rounding every Ga atom.It is assumed that, in the surrounding vapor, at the distances comparable to distancebetween steps, the differences of the concentration of Ga transporting agent are negligible,therefore the Ga atoms are adsorbed at the surface uniformly. The Ga adsorption is ac-counted for by creation of an adatom at any empty adsorption site, at each MC step, witha probability p a = ν a e − βµ , (3)where µ is a chemical potential, ν a is an attempt frequency for the adsorption process and β =1 /k B T , T is a temperature and k B is Boltzmann constant. Each adsorbed particle diffusesover the terrace until it is attached at the steps. Thus, the possibility of reevaporation isneglected. Probability of a jump from the initial to the final site, in the diffusional movement,is given by p d = ν d e − β ( E + E D ) , (4)where E D is energy barrier for diffusion, E depends on the initial α i ( J ) and the final α f ( J )energies of the Ga diffusing atom, given by equation (2), in the following manner: E = α i ( J ) − α f ( J ) , if α i ( J ) > α f ( J );0 , otherwise. (5)7nd ν d is an attempt frequency for the diffusion process. We use ν d exp( − βE D ) = 1 and ν a = 10 − . The first of these parameters sets the timescale of the whole process. Results ofsimulations can be easily rescaled to any other timescale. The second value relates timescaleof adsorption process on the timescale of the diffusion process. Note that having ν a fixed,we control rate of the adsorption process by setting value of µ .Incorporation of Ga adatoms at the step is determined by the rate at which, the Gaadatoms, once adsorbed at the terrace jump to the sites at the step. This rate is additionallymodified by Schwoebel barrier [32] that sets up different probabilities for atoms jumping tothe step from upper and lower terraces. The Schwoebel barrier height is determined by theparameter B . Probability of the jump, over the barrier to the step, is given by: p Bd = e − βB p d . (6)so that the height of the barrier B modifies the standard jump rate, given by equation (4).The potential felt by a single particle at both step types is sketched in Fig 2. As typicallyassumed, we impose additional Schwoebel barrier for the jumps from upper terrace. In orderto keep the number of control parameters as low as possible we use the same barrier height B for both types of the steps. In general Schwoebel barrier can depend on the kind of thestep, but in order to have similar influence on the step kinetics it should depend on the steporientation also. Technically that means that the barrier should depend on the energy ofthe particle adsorbed at the step, but such dependence is in our model already accountedfor as probability of the jump out the step. We show below that such assumption, being themost natural one, is sufficient to reproduce various variants of the surface kinetics.Crystal surface microstate is given by setting the two uppermost layers of atoms. Whennew step appears, the upper layer is converted into the lower one, and a new layer is built ontop of the terrace. The lower terrace is shifted into the bulk of the crystal, not modeled here.In such a way continuity of particle-particle interaction at the step is guaranteed. Everysecond layer of Ga atoms has different bond orientation. Due to different positions of Natoms, orientation of tetrahedral bonds is rotated by 60 o in alternate layers. Such geometrycauses step flow anisotropy characteristic for GaN surface. Step flow anisotropy in the caseof Eq. (4) is realized by the particle jump in the potential that is illustrated in Fig. 2. Notethat it differs from the potential proposed in Ref.33. In our model it is potential minimumat the step, which is modified by the interactions with other particles, whereas in Ref. 338
0 50 100 150 200 0 50 100 150 200 9 12 15 18
FIG. 3: (color online) Stationary pattern of double terrace structure of different width and straightsteps. The initial terrace width was 20 lattice constants. Simulation was carried out for system ofsize 200 ×
240 lattice constants, r = 0 . B = 0, βJ = 5 and µ = 15. Steps are perpendicular to[01¯10] orientation. the barrier height B for a jump at the step is modified in order to realize different stepvelocities. As a result when r = 1 we get different bonding not only at every second stepfor also for the step orientation at the misoriented surface (0001) of GaN crystal. For thesame density of adatoms, the step velocity depends on the particle interaction and thereforeit changes with its location and orientation.We start our simulations with an even number of equally spaced, parallel steps on thesurface. Heights of the neighboring steps differ by one Ga atomic layer. New Ga atoms areadsorbed at the surface. Periodic boundary conditions are applied in the lateral directionwhile in the direction in which the crystal grows they are corrected by constant heightdifference between both sides of the simulated area.9 IG. 4: (color online) Step pattern evolution for system with steps initially oriented perpendicularlyto [10¯10] and [01¯10]. Surface patterns of a system for following evolution times 0, 1 . , and 3 . MC steps are presented from top to bottom. Step anisotropy is given by r = 0 . βµ = 15, βJ = 5, B = 0 and at it the system of size 200 ×
160 lattice constants, eight V shaped terraces wereprepared. For comparison, initial step shape is drawn in thin black line. II. STEP DOUBLING
We begin our studies from the case starting from the initial configuration, having straightterraces of equal widths and the steps descending along [01¯10] direction, i.e. the steps areoriented in [11¯20] direction. After very short initial time, around 10 MC steps, the stepsystem evolves to the stationary double terrace pattern, shown in Fig. 3. Faster steps moveforward, and in the result, every second terrace shrinks, being close to disappearance. Suchpattern is very stable and it is preserved, essentially in the same form, during simulationsbuilding up to eight consecutive layers, as illustrated in Fig 3. The height of the lowestterrace of the modeled area was initially set to zero, thus its displayed height indicates thenumber of newly deposited layers during simulation. Let us stress that the stationary statelike the presented one was attained only for the step anisotropy caused by four-body Gainteractions. For the isotropic case (r = 1), the width of the terraces does not change,indicating the same kinetics of both steps. The parameter range for which system evolvesinto such pattern is determined by the complex interplay of the temperature, particle flow,Schwoebel barrier height and the anisotropy. In the case when the terraces are wider or forhigher particle flux, for given temperature and anisotropy, the new phenomenon of the stepmeandering emerges.As proposed initially by Burton, Carbera and Frank (BCF), dislocations of the screw ormixed type, serve as continuous step sources on otherwise atomically flat crystal surfaces[34]. Thus it could be assumed that steps are created at the location where such dislocationline emerges at the surface. Under supersaturation, the steps attain the spiral shapes, ex-tending over large part of the surface. The exact spiral shape depends on the anisotropy andthe temperature, being rounded, close to Archimedean, for high temperatures, and havingstraight line fragments for low temperatures. GaN growth corresponds to low tempera-ture case in such fashion that the steps, limiting flat terraces, encircle the dislocation line,creating hexagonal pyramid. Thus these steps are directed along six sides of hexagon.In order to model such case we prepared array of V-shaped steps, broken in the center,with its parts oriented along two neighboring sides of the hexagon, i.e. along the followingdirections: [10¯10] and [01¯10]. Such design could be treated as a segment of a step pyramid,corresponding to two sides of the hexagon. Such initial configuration of the system is depictedin top panel in Fig.4. Both halves of this system create regular pattern of parallel steps,11
IG. 5: (color online) Step pattern evolution for systems with anisotropic and isotropic steps.Initially steps were bend 19 o to the crystallographic directions. Consecutive evolution plots arepresented from top to bottom for the system with anisotropic steps described by parameter r = 0 . r = 1 at the right hand side. System was 200 ×
240 latticeconstants with 12 V shaped terraces and other parameters: βµ = 15, βJ = 5, B = 0 Number ofMC steps changes from top down as follows: 0, 5 . , and 10 V -shape becomes more sharp. Stepsmove closer, but their motion is not the uniform, some moves faster, the other slower, inthe result they start to form meanders. Slowly they bend to different orientation, evolvinginto one which seems to be the stationary state for the system. When proper orientation isreached the step system evolves to the double terrace structure, preserving the selected steporientation, at the approximate angle 19 o to main, initially used crystallographic direction.It is worth mentioning that the terraces widths are not the same, green and black terracesbeing wider on the left and right side, respectively.In order to find out whether such curved step line is dynamically stable, V -shaped con-figuration, we started our simulations with the initial pattern of steps, oriented at angle19 o to the main crystallographic direction, as shown in Fig.5. With steps initially bent thesystem evolved into stationary state presented at the left side of the Fig.5. System withanisotropic steps forms patterns of wide and narrow terraces following each other. Whenthe step changes its orientation, wide terraces smoothly evolve into narrow ones buildingcharacteristic pattern, which is often observed experimentally at the surface of MOVPEgrown GaN layers [14]. The step pattern, obtained in the case for which the steps areanisotropic, plotted at left hand side of Fig. 5, could be compared with the one obtainedfor the isotropic case, plotted at the right hand side of Fig. 5. In the latter case, the ter-races are equally spaced and therefore, for isotropic case, dynamically stable step pattern iscompletely different.It is interesting to find out which factor is responsible for such difference in the stepevolution from a set of parallel steps, shown Fig 3 and from V -shaped pattern, demonstratedin Fig 4. It seems that the difference between these two cases stems from initial conditionscoupling to boundary conditions. In the first case, the steps are fixed by periodic boundaryconditions, making the whole structure very rigid. Each step, attached at both sides, couldmove forward only, and finally characteristic periodic step pair structure is formed. It ispossible that for larger size, the system will evolve into the V shaped pattern, observed forthe second case, where periodic boundary conditions are also applied, but V shaped structureleaves freedom for step bending, and it appears that the steps change their orientation to13 ( k - κ ) / λ D (k - κ )/ λ D d=0 d=0d
The mechanism of a double step pattern creation, in the system having two types ofanisotropic steps can be analyzed employing Burton-Carbera-Franck (BCF) model [34] ofthe surface diffusion controlled growth. In the BCF model it is assumed that due to theincorporation of the adsorbed molecules at the steps, its local density ρ changes, being thesmallest at the steps. Accordingly, the particles diffuse towards the steps. In the simplestformulation one can consider a train of parallel steps, i.e. the problem is reduced to a onedimensional diffusion system. In the case of a stationary step motion with constant velocity V , we the mass balance equation, in the step coordinate system, can be written in the time14ndependent form: D d ( dz ) ρ + F − ρτ + V ddz ρ = 0 (7)where D is the diffusion coefficient, F describes the flux of incoming atoms, adsorbed atthe terrace and τ is an average lifetime of an atom in a state of mobile adsorption on thecrystal surface. In the MC model presented above, adatoms do not desorb, but doubleoccupancy is forbidden that leads to the relation τ = 1 /F . Equation (7) is formulatedin the frame attached to the step that moves with constant velocity V , hence the latter,velocity dependent term arises from Galilean transformation from crystal lattice to stepattached coordinate system. We express local density of the adsorbed particles ρ in unitsof a number of particles per one lattice site, and similarly all distances are measured in anumber of intersite distances (i.e. lattice constants a ). Thus z is dimensionless and the unitof D , F and V is 1 /s . Solution of this equation (7), has the form ρ ( z ) = F τ − A cosh( λz ) − B sinh( λz ) (8)where, for small V : λ = s Dτ (9)The choice of the boundary conditions completely specifies problem. With two differentstep types we have to assume two density profiles (8) with different constants A , B . Foreach of density profiles we have two edges, what gives four boundary conditions. Let thelength of the first of two consecutive terraces be equal to d , and of the second to l − d . Thuswe assume that the double step structure repeats with period length l , i.e. the surface iscovered by an infinite train of parallel steps, at which the following boundary conditions areimposed on the density of adatoms at the terraces: D dρdz | ( − d )+ = k ( ρ − ρ +1 ) | ( − d )+ − D dρdz | − = κ ( ρ − ρ − ) | − D dρdz | = k ( ρ − ρ +2 ) | (10) − D dρdz | ( l − d ) − = κ ( ρ − ρ − ) | ( l − d ) − In order to make our calculations simpler we position one step at z = 0, and the otherone at − d , which is periodically equivalent to l − d . Parameters k , κ and k , κ describe15he rates with which particles attach from the right and left side to the first and to thesecond step and ρ ± i is equilibrium density at step i . The velocity of the step motion shouldbe equal for both steps, and is given by V = D dρdz | + D dρdz | − = D dρdz | ( − d )+ + D dρdz | ( l − d ) − (11)When conditions (10) and (11) are imposed, the solution (8) it turns out that the width d can be determined from the relation k + κ k − κ coth( λd ) − k + κ k − κ coth[ λ ( l − d )] = (12) κ k + λ D λD ( k − κ ) − κ k + λ D λD ( k − κ )The above equations describe complete solution of the stationary dynamics of the doublestep structure, obtained from Eq.(7). Step permeability was neglected here. It general,it can be accounted for as in Ref.35, slightly complicating the above expressions but notchanging the conclusions. The solution of Eq. (12) can be obtained providing that thecondition ( k − κ )( k − κ ) > d and l − d of two consecutiveterraces in a stationary solution. However when κ becomes close to k or κ ≈ k the firstor the second type of terrace disappears, the condition (13) is not obeyed, and Eq.(12) doesnot describe relative terrace widths. Double steps form and they propagate differently froma pairs of steps described above.In Fig. 6 we plot regions of parameters for which different step behavior is observed.Double step structure is observed when the terrace width d = 0 or d = l . This was the case,described in the previous Section. When the Schwoebel barrier is zero i.e. when the ratesof the particle attachment to the step from upper and lower terrace are the same, k = κ and k = κ , the condition (13) is not fulfilled. In such case, the faster step overtakes theslower one, and double step structure forms. In the special case, when all rates are thesame, the steps move with the same speed, thus the terrace widths do not change. Theystay essentially identical to the initial state. Note that in simulations presented in Figs. 3, 4and 5, every second terrace is very narrow, but it does not disappear totally. This can be dueto the discrete character of the system, interactions between diffusing particles and the fact16 IG. 7: (color online) Step pattern evolution for two different Schwoebel barriers. Steps initiallywere oriented as in Fig.5. Consecutive evolution plots are presented from top to bottom for thesystem with low Schwoebel barrier βB = 0 . βB = 1 at right handside, other parameters of the systems were βµ = 15, βJ = 5,initial step width equal to 20 latticeconstants and number of MC steps changes from top down as follows: 0, 5 . , and 10 κ = κ = 0 equation (12) simplifies to the formula[36] coth( λd ) − coth[ λ ( l − d )] = Dλ ( 1 k − k ) (14)This equation could be solved for any value of the pair of parameters k , k , hence thissystem always attains stationary state, with two terraces having different, finite widths. V. STEP MEANDERING
Step anisotropy in the system leads to the stable double terrace structure for sufficientlyhigh disorientation, low particle flow and low Schwoebel barrier. For wider terraces the stepsstart to meander and cannot overtake each other. We do not observe double step structurefor wider terraces. Instead, patterns of meandering steps change during surface evolution.The value of terrace width at which meandering starts, depends on the Schwoebel barrierand on the flux of incoming particles. In Fig. 7 step evolution for the two systems, havingdifferent Schwoebel barrier only, is shown. Plotted, at the left hand side of Fig. 7, theevolution of the pattern, obtained for the barrier B = 0 . J , does not differ much from thatfor B = 0, shown in Fig.5. When r = 0 .
4, the barrier B = 0 . J corresponds also to the caseof the condition (13) not fulfilled, i.e. the same as for B = 0 at any flux. For the barrierincreased to above 0 . J , the growth, in accordance to the fulfilled condition (13), tendsto attain stable pattern of the two terraces having different widths. At very high barrier,i.e. increased up to B = J , the steps are destabilized by fragmentation so that the systemeventually evolves into the rough surface, the mechanism known as Ehrlich-Schwoebel effect[32, 37, 38]. In Fig. 7, the patterns obtained for two different Schwoebel barriers at thesame evolution times, are compared. Destabilization, due to the particle flux increase, lookssimilarly to this one. When particles attach step so frequently that the diffusion alongstep is too slow for step smoothening, meandering is overcome by the step roughening andfragmentation. Note, that beside difference in the step detachment rates we assume nodiffusion, or potential anisotropy, which were present in other MC models [39, 40].18 IG. 8: (color online) Domains for large particle flux, µ = 13. Upper row shows system with r = 0 . r = 1 at righthand side and anisotropic system after longer evolution. Three first cases are shown after 10 MCsteps and the last one after 10 MC steps. In these simulations βJ = 5 and βB = 0 (no Schwoebelbarrier). In the case when the chemical potential of the vapor phase is set very high, i.e. whenthe external particle flux is very intense, the particles tend to form islands at the terraces(nucleation). The phenomenon occurs because incoming particles reach surface so often thatthey rather find others, stick together and create islands (nuclei of new layer), before theycould reach steps. In Fig. 8 in two top plots, we see triangular domains oriented antiparallelyat every second terraces, obtained for the two different, perpendicular orientations of thesteps. Left hand side top panel shows steps perpendicular to [11¯20] direction and righthand side top panel steps perpendicular to [10¯10] direction. Both were grown for r = 0 . r = 1 at the same growth conditions as in the upper, right diagram. In this case, thenuclei of the new layer are more regular, close to hexagonal. The latter picture, shown inthe right bottom, presents rough surface, obtained during longer evolution of an anisotropicsystem shown the left, upper diagram. From these diagram we conclude that the particledensity at the terrace, or equivalently, the flux of incoming atoms, should be low enough(low supersaturation) in order to attain stationary step flow. VI. STABILITY ANALYSIS
In order to check when parallel steps start to meander we carry out stability analysis [41–43]. In such analysis we perturb the shape of each step with a finite harmonic oscillationsof wavevector k and frequency ω . Position of i -th step ( i = 1 ,
2) becomes z i ( x ) = z i + δ i sin( kx + ωt ) (15)where δ i is the amplitude of small perturbation of the step position. The particle density ateach terrace up to the first order is given by ρ ( x, z, t ) = A sinh( λz ( x, t )) + B sinh( λz ( x, t )) (16)+ ǫ [ A sinh(Λ z ) + B sinh(Λ z )] sin( kx + ωt )with small parameter ǫ and Λ = √ k + λ In the boundary conditions (10) full expressionfor the equilibrium density should be used[41, 42] ρ ± i = ρ ± i − Γ z xx ( √ z x ) (17)where Γ is capillary length of the step [3]. The velocity of bent step should be corrected asfollows[41, 42] V = V + ˙ z √ z x (18)We look for such value of k which fulfills the following two conditions [41, 42] ω ( k c ) = 0 ∂ω∂k k = k c = 0 (19)The conditions determine such relation between parameters that separates stable straightmoving steps from that for which unstable modes emerge.20et us first analyze the case where k ≈ κ and k ≈ κ , like the one presented in Figs3,4, 5, and in Fig. 7 at left hand side. These step attachment rates obviously do not fulfillcondition (13), so double step structure is formed. Such double step pattern is a uniformsystem with terraces of width l and the steps having growth rate equal to the slower of k , k . In such a case the solution of (19) in the uniform step system is k c = 0, and stepinstability is found when the following relation between the system parameters is fulfilled λ Γ F τ [ coth( λl ) k + κk − κ + ( Dλ ) + kκDλ ( k − κ ) ] < , (20)where we assumed that k = k < k and κ = κ . According to the Eq. (20) the steps arenot destabilized when κ → k . In the classical relations used in Refs 41, 42, the value of λl isassumed to be very large and k infinite. In the limit of large λl and k , with κ = 0, a standardcondition ( λ Γ) / ( F τ ) < / F , the steps are more unstable. It follows also from a full version of the relation(20) that large terrace width l or high step adsorption rate k also lead to step meandering.According to Eq. (20), large diffusion coefficient has no influence on step stability, howeverwhen diffusion is very low, its increase can destabilize steps.In the case when the condition (13) is fulfilled, the terraces of two different widths exist,terminated by steps having different rates. We assume that at both steps parameters δ and ǫ are different. Two different steps are assumed and then two different values of δ = δ , δ and ǫ = ǫ , ǫ are present. Moreover as in Eqs. (10) z = − d, z = 0. Finally the followingrelation is obtained λ Γ F τ [ coth( λd ) k + κ k − κ + ( Dλ ) + k κ Dλ ( k − κ ) ] < , (21)According to Eq. (12) the above condition can be formulated for the second terrace width l − d replacing d and for k , κ substituting k , κ . When equilibrium conditions (12) arefulfilled on the base of (21) both step destabilize at once. In fact such behavior is visible inFig. 7. VII. CONCLUSIONS
At GaN(0001) misoriented surface two crystallographically different kinds of steps arepresent. The difference stems for the way how particle are bonded at the step. Ga adatoms21ttach step via one N atom, closing tetrahedral Ga structure or via two independent Natoms leaving two open structures. We assumed that closed tetrahedral Ga structure hasmuch lower energy than sum of simple Ga-N-Ga bond. As an effect particle step detach-ment rates change at every second steps. Rotation of the step orientation by 60 o exchangesthat two kinds of steps. We show that assumption of strong bonding of closed tetrahedralstructures has its consequences in the specific stationary surface patters. We have demon-strated that equally distanced parallel steps evolve towards structure of pairs of terraces ofdifferent width. For some parameters values, one of terraces reduces to zero and double stepstructure is build. The condition for such behavior was calculated using a one-dimensionalBCF formulation. In the other considered case, the steps, which initially form V -shapedstructure change their orientation first, after which every second step moves forward formingpairs of step structure. When V -shaped steps are initially prepared along lines of highercrystallographic numbers, they strictly evolve into the double step structures. Both thesetypes of system behavior together explain step patterns observed in the growth on manysemiconductor surfaces, such as SiC[11], Si[12, 13] and GaN[14]. For small structures, thatare fixed at both ends, like in spiral structure two opposite step bending, we observe pairof terraces terminated by parallel steps. For longer hexagon edges, the steps bend and weobserve wave-like step pair structures. When external particle flow increases or for higherSchwoebel barrier, the step starts meandering and after that surface becomes rough. Sta-bility analysis of double step structure allows to find conditions for which step destabilizes.It appears that both steps destabilize for the same conditions. Therefore the main distinc-tive feature of the proposed model i.e. the existence of many-body interactions between Gaatoms allowed us to recover various surface GaN(0001) phenomena, observed experimentally[14]. Acknowledgments
Research supported by the European Union within European Regional DevelopmentFund, through grant Innovative Economy (POIG.01.01.02-00-008/08) [1] W.W. Mullins, R.F. Sekerka, J. Appl. Phys. , 323(1963).
2] W.W. Mullins, R.F. Sekerka, J. Appl. Phys. ,444(1964).[3] J.S. Langer, Rev. Mod. Phys. , 1(1980).[4] D. Wollkind, L. Segel, Phil. Trans. Roy. Soc.(London) , 45(1984).[6] U. Nakaya, Snow Crystals, Harvard University, Cambridge, Mass. 1954.[7] W.A. Bentley and W.J. Humphreys, Snow Crystals, Dover, London 1962.[8] H. W. Kroto, J. R. Heath, S.C. O’Brien, R.F. Curl, R.E. Smalley, Nature , 162(1985).[9] S. Iijima J. Cryst. Growth , 675(1980).[10] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zang, S. V. Dubonos, I. V.Grigorieva, and A. A. Firsov, Science , 666(2004).[11] A. R. Verma, Crystal Growth and Dislocations, Butterworth Scientific Publications, London1953.[12] I. Sunagawa and P. Bennema, in Preparation and Properties of Solid State Materials, vol. 7.Ed. W. R. Wilcox, Marcel Dekker New York 1982, p. 1.[13] A. D. Verga, Phys. Rev. B , 174115 (2009).[14] S. Krukowski, P. Kempisty, P. Strak, G. Nowak, R. Czernecki, M. Leszczynski, T. Suski, M.Bockowski, and I. Grzegory, Crys. Res. Technol. , 1281(2007).[15] S. Nakamura, M. Senoh, and T. Mukai, Jpn. J. Appl. Phys. Part-2-Letters. (10A), L1708(1991).[16] T. Suski, E. Litwin-Staszewska, R. Piotrzkowski, R. Czernecki, M. Krysko, S. Grzanka, G.Nowak, G. Franssen, L.H. Dmowski, M. Leszczynski, P. Perlin, B. Lucznik, I. Grzegory, R.Jakiela, Appl. Phys. Lett. ,172117 (2008).[17] A.E. Romanov, T.J. Baker, S. Nakamura, and J.S. Speck, J. Appl. Phys. ,23522(2006).[18] E. Frayssinet, W. Knap, P. Lorenzini, N. Grandjean, J. Massies, C. Skierbiszewski, T. Suski,I. Grzegory, S. Porowski, G. Simin, X. Hu, M.A. Khan, M.S. Shur, R. Gaska, D. Maude, Appl.Phys. Lett. ,2551 (2000).[19] R. Langer, J. Simon, V. Ortiz, N.T. Pelekanos, A. Barski, R. Ander and M. Godlewski, Appl.Phys. Lett. (1999) 3827.[20] J. S. Im, H. Kollmer, J. Off, A. Sohmer, F. Scholz and A. Hangleiter, Phys. Rev. B , R9435(1998).[21] D.A. B. Miller, D.S. Chemla, T.C. Damen, A.C. Gossard, W. Wiegmann, T.H. Wood, and .A. Burrus, Phys. Rev. Lett. 53, 2173 (1984).[22] M.P. Halsall, J.E. Nicholls, J.J. Davies, B. Cockayne, and P.J. Wright, J. Appl. Phys. ,907(1992).[23] S. Krukowski, P. Kempisty, P. Strak, G. Nowak, R. Czernecki, M. Leszczy`nski, T. Suski, M.Bockowski, I. Grzegory, Cyst. Res. Technol. , 1281(2007).[24] C. Skierbiszewski; ZR Wasilewski, I. Grzegory, S Porowski, J. Cyst. Growth ,1632 (2009).[25] I. Grzegory, J. Phys.: Condens. Matter 14, 11055 (2002).[26] Yu.A. Vodakov, and E.N. Mokhov, in: Vacuum Science and Technology: Nitrides as Seen bythe Technology, Eds T. Paskova and B. Monemar, Research Signpost, Trivandrum 2002, pp.59-78.[27] F. Dwikusuma, S.E. Babcock, and T. F. Kuech, in: Vacuum Science and Technology: Nitridesas Seen by the Technology, Eds T. Paskova and B. Monemar, Research Signpost, Trivandrum2002, pp. 79-105[28] R. Dwilinski, R. Doradzinski, J. Garczynski, L.P. Sierzputowski, A. Puchalski, Y. Kanbara,K. Yagi, H. Minakuchi , and H. Hayashi, J. Cryst. Growth , 3911 (2008).[29] S. Krukowski, P. Kempisty, J. Cyst. Growth , 900 (2008).[30] S. Krukowski, P. Kempisty, A.F. Jalbout, J. Chem. Phys. , 234705 (2008).[31] J. Neugebauer and C.G. Van de Walle, Phys. Rev. Lett. , 4452(1995).[32] R. L. Schwoebel, J. Appl. Phys. (1969)614.[33] M. H. Xie, M.Gong., E.K.Y. Pang, H. S. Wu, and S. Y. Tong, Phys. Rev. B , 085314 (2006)[34] W. K. Burton, N. Carbera, and F. C. Frank, Philos. Trans. R. Soc. London, Ser. A , 299(1951).[35] B. Ranguelov, M.S. Altman , and I. Markov, Phys. Rev. B , 245419 (2007).[36] M.A. Za luska-Kotur, F. Krzy˙zewski and S. Krukowski J Non. Solids 2010,doi:10.1016/j.jnoncrysol.2010.05.029.[37] M. Rusanen, I. T. Koponen, J. Heinonen, and T. Ala-Nissila1, Phys. Rev. Lett. 86, 5317(2001).[38] M. Rusanen, I. T. Koponen, T. Ala-Nissila, C. Ghosh, and T. S. Rahman Phys. Rev. B 65,041404(R) (2002).[39] R.Kato, M. Uwaha, Y. Saito, Surf. Sci (2003) 64.[40] R.Kato, M. Uwaha, Y. Saito, Surf. Sci (2004) 149.
41] G. S. Bales and A. Zangwill, Phys. Rev. B, (1990) 5500.[42] I. Bena, C Misbah, and A. Valence, Phys. Rev. (1993) 7408.[43] C Misbah, O Pierre-Louis and Y. Saito, Rev. Mod. Phys. , 981(2010)., 981(2010).