Structures Built by Steps Motion during Sublimation from Annealed

Article pubs.acs.org/crystal

Structures Built by Steps Motion during Sublimation from Annealed GaN(0001) Surface Magdalena A. Załuska-Kotur,*,†,‡ Filip Krzyzė wski,† Stanisław Krukowski,#,§ Robert Czernecki,#,∥ and Michał Leszczyński#,∥ †

Institute of Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland Faculty of Mathematics and Natural Sciences, Cardinal Stefan Wyszynski University, ul Dewajtis 5, 01-815 Warsaw, Poland # Institute of High Pressure Physics, Polish Academy of Sciences, ul. Sokołowska 29/37, 01-142 Warsaw, Poland § Interdisciplinary Centre for Materials Modeling, Warsaw University, Pawińskiego 5a, 02-106 Warsaw, Poland ∥ TopGaN, Sokolowska 29/37, 01-142 Warsaw, Poland ‡

ABSTRACT: Monte Carlo (MC) data, obtained during simulation of an annealed surface GaN(0001), are compared with those obtained from experiment. It is shown that simulations reproduce the majority of patterns observed in experiments provided that a suitable selection of the simulation parameters is made. MC simulations show the development of the characteristic step patterns that depend on the annealing temperature and the time. The pattern, observed for all processes after a very short evaporation time, consists of rough steps. The long time evolution develops into two different patterns. The first, composed of single curly steps, is typical for lower temperatures. In the second, after a sufficiently long time, the steps merge together creating bunches which, at higher temperature, bend into wavy-like structures.

■

INTRODUCTION Epitaxial growth of nitride layers, both pure GaN or solid AlGaN and GaInN solutions, either by metal organic vapor phase epitaxy (MOVPE) or by molecular beam epitaxy (MBE), is the core of nitride based device technology.1,2 Device performance depends on the quality of the synthesized layers, both in terms of crystallographic structure and of their chemical composition. It is well-known that the results of the growth depend critically on the structural and chemical state of the surface, prior to epitaxy. Therefore, the substrates are mechanochemically polished to attain an atomically smooth surface, denoted as an epi-ready surface.3−5 However, when introduced into the epitaxy growth chamber (MOVPE or MBE), the substrates contain a variety of molecules attached when they are exposed to air.6 Therefore, prior to the epitaxial growth, the substrates are annealed to clean the surface.7,8 The temperature and the time of the annealing, as well as the atmosphere, has to be optimized for every kind of substrate, because the inappropriate choice of these parameters may result in deterioration of the surface crystallographic quality, causing considerable degradation of the quality of the grown epitaxial layers. In the work below, we present an atomic scale model which is used in Monte Carlo simulations of the sublimation of the uppermost GaN substrate layer and compare the obtained simulation results with experimentally observed structures. In the experiment, the initial arbitrary configuration, created by mechanochemical polishing,3,7−9 undergoes a transition to a © 2013 American Chemical Society

well developed step pattern. This is achieved in the initial stage of MOVPE growth when the GaN substrate is annealed for 5 min at 1040 °C in a mixture of ammonia, nitrogen, and hydrogen, which remove foreign chemical species, and also some portion of gallium nitride. During this stage, the atomistic surface structure undergoes transformation which could be controlled to some extent by the thermodynamic conditions of the annealing and by the duration of the process. On using Monte Carlo simulations, we study basic features of the step motion that lead to morphological instability and the emergence of double step surface morphology, frequently observed after the annealing in MOVPE experiments. We consider GaN, the wide-band gap semiconductor used for fabrication of white LEDs, BluRay lasers, and many other optoelectronic and electronic devices. GaN physical properties differ significantly from other III−V semiconductors (AlGaIn)(AsP). One of these differences, explored in this paper, is that GaN crystallizes in hexagonal wurtzite, different from cubic zinc blende structure, typical for the other III−V compounds. In wurtzite structure, the consecutive atomic steps of the (0001) surface are different, the feature having immense consequences not only for the growth morphology, but also for the surface changes during annealing prior the epi-growth. As a result, the difference in the bonding energy at steps rather than etching Received: April 2, 2012 Revised: January 4, 2013 Published: February 6, 2013 1006

dx.doi.org/10.1021/cg301521a | Cryst. Growth Des. 2013, 13, 1006−1013

Crystal Growth & Design

Article

anisotropy10,11 leads to step bunching instability of the GaN(0001) surface. The paper reports the results of the Monte Carlo simulations of the GaN surface evolution during annealing. The simulations are compared with the experimental data obtained for the GaN substrates annealed in the MOVPE chamber in a flow of ammonia, nitrogen, and hydrogen with a ratio of 2:1:1.

■

THE MODEL We construct the model taking into account overwhelming excess of active nitrogen in the system as the annealing is conducted in ammonia rich ambient conditions; therefore all bonds are saturated by N atoms, including those pointing vertically from the topmost Ga atoms. It is natural to assume that the whole surface is covered by nitrogen atoms and the process is controlled by desorption and diffusion of Ga atoms only. Thus lattice gas model12 of wurtzite structure of GaN is built by Ga atoms, which due to the presence of nitrogen between them interact via the layer and orientation dependent forces. Such many-body interactions determine the equilibrium and barrier energies, i.e., the diffusion and desorption rates, which define kinetics of the system. Due to the orientation dependent interactions, a 3-fold symmetry is enforced, natural for the GaN(0001) surface, rotated by 60° for even and odd terraces. We will demonstrate a consequence of these assumptions for the surface evolution during annealing at low and high temperatures. Surface evolution is presented using the results of Monte Carlo simulations attempting to grasp the nature of the changes and the physical forces responsible for the scenario observed. GaN crystal is a complicated system and many aspects of its structure and interaction energies can shape its dynamical behavior. However as we observe many details as the exact values of the energy, diffusional paths of individual particles or even height of the Schwoebel barrier (up to some level) has little influence on the collective system behavior. Hence we have oversimplified our model wanting to catch the most important parameters that decide the surface dynamics as an effect of many particle collective behavior. For this purpose, the KMC method is a very good one as it helps to study large systems for a long time. The critical assessment of the obtained results and the conclusions, pertinent to the ammonia based MOVPE homoepitaxial growth of the nitrides, are drawn. The GaN crystal is periodic with respect to translation by two double Ga−N atomic layers along the hexagonal axis, i.e., along the [0001] direction. Therefore, the atomic structure of the GaN(0001) surface has to be built considering two consecutive double atomic layers that differ by the location of the Ga sites and by Ga−N bond orientations. Due to this feature, the two successive atomic steps on the GaN(0001) surface have different atomic structures. As illustrated in Figure 1, three bonds connecting the Ga atom with the nearest N atoms within a single double atomic Ga−N layer are rotated by 60° with respect to the bonds in the consecutive Ga−N layer. Hence, at the parallel steps, depending on the layer, the particles are attached to the steps either saturating one or two dangling bonds. Thus, the step changes its properties either by stepping one terrace down or by rotation of 60°. Such crystallographic difference and consequently number of stepparticle bonds may lead to a different step attachment energy. However, a closer analysis of the GaN surface structure under heavy excess of active nitrogen proves that adatoms attached at both types of steps have the same number of first neighbors (atoms N) and second neighbors (atoms Ga). The number of

Figure 1. The model used for calculation of the total energy GaN crystal decompose the crystal into a set of tetrahedrons containing N atom inside, presented from the top (upper diagram) and the side (bottom diagram). Shaded triangles represent top and side view of tetrahedrons we consider in the model. The full tetrahedrons are shaded blue, and the ones with missing Ga atom - green. Ga atoms at two consecutive layers are plotted by red (upper layer) and gray circles (lower layer); nitrogen atoms are green. The terrace height decreases by one layer at a step (black line) having two different atomic structures, either rotated by 60° (top), or shifted one atomic layer down (bottom). The blue line limits the simulation region.

first neighbors is identical when we assume that each bond from a Ga is saturated by nitrogen, as we have full occupation of N lattice by adatoms, either bonding Ga adatom to the step or to the layer below. Thus, the model that describes different properties of the steps cannot be limited only to additive twobody interactions, but, as we have shown in previous works,13−15 the character of the step has to be described by the orientation dependent many-body interactions. The main difference between steps is in the step structure and in the number of bonds by which the particle is attached. The simplest way to model such difference is to introduce a fourbody interaction between Ga atoms. An assumption of overwhelming presence of nitrogen simplifies the lattice model, reducing the role of N atoms as bonds between Ga atoms. Accordingly, the dynamics of the surface is described by modeling the energy and the dynamics of the Ga atoms only. According to the GaN crystallographic structure, each Ga atom belongs potentially to four tetrahedrons, the three in the single Ga−N layer and the one linking this layer and the layer above (thus the tetrahedrons point down). Naturally the gallium atoms of the flat surface topmost layer belong to three tetrahedrons only. The total energy of the system could be expressed summing over the Ga atoms, using the abovedescribed decomposition into N-atom centered tetrahedrons, labeled by i, as 4

E(J ) = J ∑ ni i=1

(1)

The parameter J scales the energy of the bonds, ni is a number describing occupancy of the ith tetrahedron (i.e., number of gallium sites occupied), and the sums run over four tetrahedrons, surrounding the Ga atom considered. By the 1007



Article

interchange their roles at the next step. The situation is illustrated in Figure 1, where the difference between [11̅00] and [2̅110] is clearly seen. The total detachment energy from every step of the latter orientation is the same. When the step disappears, on coming up, the upper layer is converted into the lower one, and a new layer is opened on the bottom of the terrace. In such a way, a continuity of particle−particle interactions at the step is guaranteed. We start our simulations with an even number of straight steps at the surface, equally spaced by d lattice constants. Heights of the neighboring steps differ by one Ga atomic layer. Ga atoms detach from the step, diffuse along the terrace and desorb from the surface. Periodic boundary conditions are applied in the lateral direction, and in the direction in which the crystal grows they are corrected by a constant height difference between both ends of the system. In this work we consider the evaporation from the surface, differently from ref 16, where crystal growth was investigated. In the present model the individual particle dynamics is essentially identically determined on the microscopic level by the interactions between individual admolecules at the surface.

construction the three Ga atoms of any tetrahedron are located within the layer and the fourth atom is located below, in the layer underneath. The model defines effective four-body interactions by the determination of ni in the following way: ⎧1, when tetrahedron has all atoms ⎪ ⎨ ni = 1 ⎪ rη , when tetrahedron has empty sites ⎩3

(2)

where η is a number of Ga-occupied neighboring sites of tetrahedron i and r describes the relative strength of the fourbody and the two-body interactions in the system. When r = 1 two body Ga−Ga interactions (or more precisely triple Ga− N−Ga interactions) sum up to the value characteristic for a fully occupied tetrahedron; i.e., no additional four-body Ga interactions are present in the system. When r < 1 Ga−Ga bonds do not sum up to a value of one multiparticle bond. In such a case, the tetrahedron energy is not a simple sum of standard Ga−Ga interactions only and this is the case studied here. The relative strength of the interactions, given by the parameter r are not known from the experiment, nor from the DFT calculations; however its exact value does not influence the collective behavior of the system so much. We set r = 0.4 across the whole paper. We checked that such quantity provides that anisotropy of the every second step movement is similar to that observed in the experimental systems. We study GaN crystal surface dynamics by controlling kinetics of Ga atoms in the lattice gas model. At a given temperature, particles can jump between neighboring adsorption sites which are identical to the GaN lattice sites. The target jump sites are limited to the sites only that belong to the same tetrahedron. A probability of a jump from the initial to the target site is given by

D = D0e−β ΔE

■

STEP PATTERN EVOLUTION DURING ANNEALING In the experiments we used bulk GaN substrate obtained by the HVPE method17 in Unipress. Before annealing process the substrate was prepared as follows: first, it was cut into pieces of about 1 cm2; second, orientation of the substrate surface was set at about 0.8 degrees to the surface (0001) using X-ray diffractometer; third, the Ga-side surface was mechanochemically polished,18 and then samples were etched in a mixture of H2O2/H2SO4, and rinsed in DI-water. All samples were annealed in a vertical hot-wall quartz homemade reactor (diameter 60 mm) with the following parameters: temperature rise time 5 min, annealing time in constant temperature 5 min, temperature falling time 5 min, pressure 900 mbar, the flow of process gases: nitrogen 2 s l m supplied continuously and ammonia 2 s l m delivered when the substrate temperature was higher than 300°C. The substrate temperature was controlled by a thermocouple placed inside the 1 in. graphite susceptor inductively heated by an RF generator. Atomic force microscopy (AFM) images were taken with a microscope Veeco Digital Instrument DI-300 running under tapping-mode. We used a standard cantilever model RTESPW with a diameter of about 10 nm. The GaN crystal with exposed GaN(0001) surface was warmed up and kept at a high temperature for strictly defined annealing time. In the simulations time scale is set by the diffusion of the particles over the surface. Thus, the diffusion barrier, which was estimated to be about 0.4 eV,19 adds to the value of the desorption barrier given by the potential v. The value of this last parameter, that in our calculations gives results closest to that observed experimentally, is 1.6J which translates into the barrier for desorption of the single particle from the surface equal to 1.12 eV. In order to compare this quantity with the real surface energy, one should add the diffusion barrier of 0.4 eV, and subsequently the resulting sum should be rescaled according to the change of the system scale. In order to simulate large systems in the periods of time corresponding to the annealing time, one has to account that the simulated surface consisted of the 10a wide terraces (a is lattice constant) considerably smaller than the terrace width of 60a, resulting from the experimental miscut of around 0.8°. When terraces of the simulated system are wider the obtained results are similar; nevertheless, due to computer limited resources, the system

(3)

where D0 = 1 is the diffusion time scale, β = 1/kBT is the inverse temperature, ΔE depends on the initial Ei (J) and the final Ef (J) bonding energy of the jumping atom, ⎧ if Ei(J ) < Ef (J ) ⎪ 0, ΔE = ⎨ ⎪ ⎩ Ei(J ) − Ef (J ), otherwise

(4)

The particles attached at the step can detach from the step with the probability defined by eq 3, which depends on the energy of a particle bonding to the step via eq 4. In addition the particles could migrate along the step. When the particle is attached at the flat terrace only, it migrates over the surface executing the jumps to the target sites. At each jump the particle could leave the system by desorption with the probability given by pd = ad e−βν

(5)

where v is desorption potential and prefactor ad sets the desorption time scale. In our simulations we use v = 1. Thus the desorption rate is controlled via temperature β and potential v. No readsorption from the vapor is allowed. A crystal surface microstate is defined by setting the two uppermost layers of atoms. The surface is originally misoriented along one direction. Below, we will present results for [2̅110] step orientation. The particles detaching from the step of this specific orientation have to break two kinds of bonds, different at the left side and at the right side. Left and right side 1008



Article

sizes that can be studied are much smaller. We can compare results of different miscuts, when the ratio of the diffusion time from one step to another and the mean particle lifetime at the surface is set to be the same. In order to get corresponding results the particles should desorb 36 times more quickly from the surface having 6 times narrower terraces than in the annealing experiment, taking into account that the mean diffusion time is proportional to the second power of the distance. As defined in eqs 3 and 5 time scales of the atom diffusion and desorption are set by two parameters D0 and ad. In all our simulations we keep these parameters the same and equal to one. Such assumption means that both attempted frequencies of the diffusion and desorption processes are very close in their values. To make the desorption process faster, prefactor ad should be enlarged. Instead, we reduced desorption energy of the particles by about 0.5 eV. In such a way we have only one control parameter  the energy. The correction of the energy barriers by ΔE reduces the time scale for desorption by a factor exp[ΔE/(kbT)], and this relation is then used when we set the relative time scales at different temperatures. In summary the rescaling gives the barrier for desorption in the modeled system on the order of 2 eV, close to 2.1 eV discussed in ref 20 and lower than 2.8 eV from ref 21. By proper rescaling of the desorption energy, we can analyze results for a higher miscut than in the experiment, having comparable simulation times and the step density at the surface. MC simulation steps are related to the experimental real times by the time scale of the adatom diffusion over the surface. We set the time scale by comparing the number of desorbed layers of the crystal in the experiment and for the simulated systems. The temperature dependence of the time scale is given by τ = τ0 exp[(Ediff + ΔE)/(k bT )]

Figure 2. AFM topography image of the crystal surface (left plot) is compared with the MC simulated system (right plot). Left panel presents substrate prepared during mechanochemical preparation of the crystal. Cutting angle is 0.8°, with the step oriented along the direction [112̅0]. The sample is of 2 μm × 2 μm size. Steps in a long, average scale are straight and parallel, but in a short scale they are rough. The same results were obtained in the simulated system in the initial stage of the evaporation process as shown in the right panel. The size of the simulated system was 400a × 600a which corresponds to 124 nm × 186 nm, kBT/J = 0.2 and initial terrace width was 10a; the other parameters of the simulation are r = 0.4 and v/J = 1.6.

(6)

where Ediff is the energy barrier for the diffusion over the surface, and τ0 is the prefactor independent of the temperature . This factor is used further in order to rescale the number of MC steps for the different temperatures. During the simulation process step patterns transform from one form to the other, starting from the regular pattern of equally spaced steps. In the left panel of Figure 2 we show an AFM topographical scan of the surface as it is prepared by mechanochemical polishing before sublimation. At the beginning the steps are straight and rough. This can be compared with the initial phase of simulations in the right panel of Figure 2, when primary smooth steps become rough, but they are still straight on a larger scale. Such structure of the steps, similar in the experimental observation and in the simulated results, is achieved by different atomistic processes. During the mechanochemical preparation stage, the atoms that are weakly bonded are detached from the steps by chemical action. During the initial stage of the simulation various processes occur at the beginning of the annealing of straight steps. Particles are randomly detached from steps by thermal agitation. In the short time scale steps become rough. The only step evolution that takes place occurs at very short distances developing the pattern similar to the substrate surface before the annealing procedure. In order to analyze pattern similarity more precisely we compare in Figure 3 the surface height autocorrelation function (ACF) for experimental and simulated data in perpendicular and parallel to the steps directions. The height ACF function in the x direction is calculated as

Figure 3. Autocorrelation functions (ACF) calculated along vertical (top panels) and horizontal directions (bottom panels) of surfaces shown in Figure 2. ACF functions for experimental data at the lefthand side can be compared with corresponding functions for simulated systems at right-hand side.

G (x ) =

1 Ly(Lx − x)

Ly Lx − x

∑∑ l=1 k=1

hk + x , lhk , l

(7)

where hk,l means the height of the surface at the position k,l, calculated as the distance from the mean surface level. The size of the system is Lx × Ly and G is averaged over the y axis, when calculated along x. The ACF function along the y direction is given by an analogous formula. The patterns shown in Figure 2, described by ACF functions, display sharp oscillations in the direction perpendicular to the steps, having one peak for one step. The experimental curve is less regular due to the surface irregularities. Steps are almost straight, and hence the ACF function along the steps does not vary much for both surfaces. Weak and short-range correlations can be seen in both curves 1009



Article

in the case of experimental data emerging from the background of surface irregularities. In the next stage of the step evolution, different surface patterns are created, depending on the temperature and the evaporation time. The obtained structure, typical for lower temperatures, consists of steps which are bent in a characteristic way creating a specific curly pattern. Observed from the top, the pattern looks like the system of parallel lines, similar to the initial one, but more sharp and locally wider. In Figure 4

Figure 5. Autocorrelation functions (ACF) calculated along vertical (top panels) and horizontal directions (bottom panels) of surfaces plotted in Figure 4. ACF functions for experimental data on the lefthand side can be compared with the corresponding functions for simulated systems on right-hand side.

steps are at larger distances than in the previous plot and they are higher. By analyzing the minima of the plots we can verify that the steps are doubled or tripled compared to the structure of the separate, single steps shown in Figure 2. The ACF function in both cases  experimental and simulated  is not so regular and sharp that reflects the fact that steps are no longer so well-defined and straight. When we look at correlations along steps, i.e., in the vertical direction, we can see that some correlations are built at a distance of 0.7 μm in the experimental situation and of 50 nm in the simulated systems. This corresponds with the period of the delicate wavy pattern visible in Figure 4. The same system was annealed at a higher temperature. In Figure 6, we can see a GaN surface of the same orientation and miscut as in Figure 3, but the structure is measured after annealing at 1050 °C. The annealing time is the same as in the previously discussed cases, i.e., 5 min; nevertheless the resulting surface pattern is completely different. In the place of the structure of parallel alternatively bent lines, characteristic for lower temperatures, at higher temperatures we can see step bunches that form large waves. The same pattern was created in simulations at kBT/J = 0.22, which corresponds to T = 1050 °C and for the simulation times comparable to these in the experiment. In the experimental system crystal decayed by 50 nm and during simulations 180 layers (47 nm) desorbed. The relation between experimental and simulated system was corrected in agreement with eq 6. At a higher temperature the time of the single particle jump is shorter, and thus more MC steps are needed to obtain the same real time value. In the bottom panel of Figure 6 we see waves of step bunches. For a more detailed analysis let us compare ACF functions plotted in Figure 7. The upper curves are calculated for the direction perpendicular to the initial step orientation; thus the distinct step bunches are clearly visible in the form of strong maxima and minima of the correlation function. In the experimental situation we see 10 main bunches formed in the system, which gives five in Figure 7 at the distance 1 μm equal to half of the total system width. In the simulated system we observe five bunches and all of them are represented by consecutive maxima in Figure 7. The horizontal correlation function illustrates a

Figure 4. AFM topography image of the crystal surface (left plot) is compared to MC simulated system (right plot). At left-hand side GaN(0001) surface of miscut 0.8° of steps oriented along the direction [112̅0] after heating at T = 900°C shows the characteristic curly structure. The size of the sample is 2 μm × 2 μm. Steps bend alternatively creating impression that terraces are broadened. The step pattern of the simulated annealed system at KBT/J = 0.2 with steps oriented perpendicularly to [12̅10] and separated by 10a is shown in the right panel. Simulation was carried out for system of the size 400 × 600 lattice constants what corresponds to 124 nm × 186 nm, r = 0.4 and v/J = 1.6. The step evolution, after 2.5 × 106 MC steps, is presented at the right side. The experimental crystal decayed by a few nanometers and in the simulation 12 layers evaporated, which equals 3.1 nm.

surfaces ordered in such a way are shown. The left panel presents an experimental image of the surface annealed at 900°C, and the right panel presents a similar pattern of the simulated surface. In these simulations the temperature was T = 0.2J/kB which corresponds to the experimental 900 °C, assuming that the mean activation energy for particle detachment from the step is equal to 0.7 eV. This is the energy value which sums the Ga bonding to two other Ga atoms within the same layer, while the third one is in the layer below and does not affect the barrier for detachment from the step as the adatoms are stuck identically to the layer below during the whole diffusion process over the terrace. When particle detaches from the step it breaks two additional Ga−Ga bonds as compared to the single particle at the terrace. Thus, taking also into account the presence of N atoms in between and the orientation dependence of interaction strength, the energy of the adatom attached to the step should be slightly higher than the double of the Ga adatom energy adsorbed at the surface. This last value was estimated to be about 0.3 eV,20 and hence 0.7 eV for the mean activation energy for particle detachment from the step is a reasonable value. The characteristic features of the surface can be more precisely studied in Figure 5, where ACF functions calculated for surfaces from Figure 4 are plotted. Top plots show correlation across steps. It can be easily seen that consecutive 1010



Article

Figure 7. Autocorrelation functions (ACF) calculated along vertical (top panels) and horizontal directions (bottom panels) of surfaces plotted in Figure 6. ACF functions for experimental data on the lefthand side can be compared with corresponding functions for simulated systems on the right-hand side.

the value of the lattice constant c = 0.52 nm is 0.07 nm, 0.1 nm, and 0.4 nm. The simulated values are very close to the experimental ones. At higher temperature the annealing processes is faster. Let us analyze the time dependence of the annealing process at higher temperature. When the surface is annealed longer, the steps have sufficient time to group into step bunches. This can be seen in Figure 8, where an intermediate stage of the surface evolution at higher temperature (Figure 6) is shown. Its characteristics resemble the surface at lower temperature illustrated in Figure 4. Roughness parameter changes on coming from the top panel down as the following 0.13c, 0.34c, and 0.77c. Roughness increases with time, similar to what it does with temperature. Time dependence of the roughness parameter at both temperatures can be seen in Figure 9. Circles represent the system at 1050 °C and crosses at 900 °C. As it can be seen for longer times roughness coefficient scales as t1/2 at both temperatures. Such coefficient is characteristic for bunching and meandering processes happening together.28,29 We see that at a higher temperature the curve bends at around 5 × 106 MC steps. It means that the surface roughening stops when steps start to create bunches, and after that the roughness parameter increases further with the same scaling factor 1/2. For lower temperatures we do not observe a bunch stage of the evolution either in experimental nor in simulated systems. At higher temperatures step bunching is observed in systems having various misorientations and step directions. These parameters determine how many MC time steps are needed to obtain such structures. Moreover, when step trains are created easily, then after they are formed they start to bend, developing wavy-like structures. We presented such structures in the top panel of Figure 6, emerging in the system annealed at 1050 °C, and the comparable structures obtained in the MC simulated system are shown in the bottom panel. Both surface structures  step trains and meanders result from the step movement backward during the evaporation process. The step movement in an obvious way is a source of the surface particle desorption flux, which in the annealed system leads to the bunching phenomenon.22 On the other side, in some conditions the

Figure 6. AFM topography image of the crystal surface (top panel) is compared with MC simulated system (bottom panel). Top panel presents GaN (0001)surface of miscut 0.8° steps were initially oriented along [112̅0]. The pattern after heating at T = 1050 °C shows a wavy structure. Size of the sample is 2 μm × 2 μm. In the bottom panel step pattern which evolves toward a meandered step train structure is shown in the system evaporated at kBT/J = 0.22 with steps initially oriented perpendicularly to [12̅10] and separated by 10a. Simulation pattern obtained, for system of size 600 × 400 lattice constants which corresponds to 186 nm × 124 nm, r = 0.4, and v/J = 1.6 after 107 MC steps. At this temperature 180 atomic layers evaporated, which means 47 nm, close to the experimental crystal decay of 50 nm.

wavy structure along the x axis, and so we can see oscillations along the ACF curve, in the case of both experimental and simulated structures. The period length in both cases is equal to the half of the system size. We analyzed the simulation of the surface structure at various annealing temperatures. These results can be summarized as follows: For short evaporation times only a part of the surface layer is removed and steps move by a relatively small distance only and become rough, as it was illustrated in Figure 2. In the longer evolution, the steps bend over longer distances forming characteristic curly structures like those in Figure 4 at lower temperature and in Figure 6 at higher temperatures. We have compared these surface patterns with experimental surfaces annealed at two different temperatures quantitatively calculating their ACF functions in the directions across and along the steps. They confirm similarities observed visually in the corresponding surface structures. The roughness coefficient (rms) of the analyzed surfaces increases from 0.1 nm for substrate in Figure 2, 0.13 nm for surface at 900 °C in Figure 4 up to 0.51 nm for 1050 °C (Figure 6). This can be compared with the roughness coefficient of simulated surfaces, which are consecutively: 0.13c, 0.19c, and 0.77c, which take into account 1011



Article

moving steps become unstable with respect to the fluctuations of selected wavelength and they start to meander.14,22 Interestingly, in the case discussed above we observe both kinds of instability: the step bunching happens as the first and the meandering is the second following process. Bunching instability is caused by the particle flux driven by natural step movement of the evaporated surface. The origin of the meandering instability is not so clear, but it is due to both the combined step movement and the presence of step bunches. The wavy structure is constructed out of the meandering steps bunches. It seems that meandering instability starts simultaneously with the bunching instability, but the meanders build up more slowly and become visible for longer simulation times than bunches. This is to some extent analogous to the step roughening in a shorter time scale observed at the beginning of the simulated annealing. It has to be stressed that some of the structures, described above, are different from the structures built during the growth.13,14 The main discrepancy lies in the tendency of forming bunches. The MC simulated surfaces of growing crystals build single step patterns rather than create step bunches. During evaporation step trains are formed much more easily. It has to be added that in our simulations the bunches form out during evaporation without barrier anisotropy at the step edge. Typically a Schwoebel barrier23 is taken as the main source of such anisotropy as a bunching mechanism. In the simulations reported above no Schwoebel barrier was used. The relative particle movement directed down the terraces occurs due to a natural parallel step flow at the surface during evaporation. Such flow is sufficient to induce step bunching.15 If the crystal grows, the diffusing particles are driven in the opposite direction. Opposite flow direction leads rather to meandering of steps than to bunching. Similar tendency is observed experimentally during evaporation and growth of GaN surfaces. Easily built bunches during surface evaporation spread out during further crystal growth process. At the same time, in the experimentally grown GaN crystals we observe various structures: bunches and meanders, as well as straight, parallel steps. However, the main factor controlling development of these structures is not elucidated. Step bunching is observed during electromigration24 or often can be attributed to the presence of stress25 or impurities at the surface.26 Any surface anisotropy in the particle attachment or diffusion may also cause bunching during growth.27 The presence of hydrogen adatoms also changes the kinetics of surface growth. We have shown that in the model without any additional anisotropies and the natural surface structure is composed of the step bunches.

Figure 8. Step pattern built in the evolution of surface evaporated at temperature kBT/J = 0.22. From top to bottom patterns obtained after 104, 2.5 × 106, and 107 MC steps are shown. Other parameters are the same as in Figure 6.

■

CONCLUSIONS Step patterns created in the process of temperature annealing of the GaN surface were modeled by a Monte Carlo procedure and compared with the experiment. For the simulated systems similar surface patterns were found as those emerging at real, experimentally investigated GaN(0001) surfaces. The main issue of this study is to show that several simple assumptions are enough to catch some basic behavior of the sublimated GaN surface. We also have shown that the Schwoebel barrier has not so a large an impact on the basic dynamical behavior, because if not too large it plays exactly the same role as the step movement only. We have also analyzed the bunching process during surface sublimation and show that roughness parameter scales as t1/2 for long evolution times. We study the simplest

Figure 9. Surface roughness evolution in time measured as rms roughness coefficient for simulated systems. Upper curve is plotted for the system at temperature 1050 °C and lower at 900 °C. Dashed lines represent t1/2 scaling.

1012



Article

(16) Xie, M. H.; Gong, M.; Pang, Y.; Wu, H. S.; Tong, S. Y. Phys. Rev. 2006, 74, 085314−1−6. (17) Lucznik, B.; Pastuszka, B.; Weyher, J. L.; Kamler, G.; Grzegory, I.; Porowski, S. Phys. Status Solidi C 2009, 1610−1642. (18) Weyher, J. L. Cryst. Res. Technol. 2012, 333−340. (19) Zywietz, T.; Neugebauer, J.; Scheffler, M. Appl. Phys. Lett. 1998, 73, 487−489. (20) Adelmann C.; Lymperakis L.; Brault J.; Mula G.; Neugebauer J.; Daudin B., GaN and Related Alloys-2002. Wetzel, C., Yu, E. T., Speck, J. S., Arakawa, Y., Eds. Materials Research Society Symposium Proceedings 2003; Vol. 743, pp 91−96. (21) Adelmann, C.; Brault, J.; Mula, G.; Daudin, B.; Lymperakis, L.; Neugebauer, J. Phys. Rev. 2003, 67, 165419−1−9. (22) Misbah, C.; Pierre-Louis, O.; Saito, Y. Rev. Mod. Phys. 2010, 82, 981−1040. (23) Schwoebel, R. L. J. Appl. Phys. 1969, 40, 614−8. (24) Xie, M. H.; Cheung, S. H.; Zheng, L. X.; Ng, Y. F.; Wu, H.; Ohtani, N.; Tong, S. Y. Phys. Rev. 2000, 61, 9983−5. (25) Tersoff, J.; Phang, Y. H.; Zhang, Z.; Lagally, M. G. Phys. Rev. Lett. 1995, 75, 2730−3. (26) Van der Eerden, J. P.; Mïller-Krumbhaar, H. Phys. Rev. Lett. 1986, 57, 2431−3. (27) Pascale, A.; Berbezier, I.; Ronda, A.; Videcoq, A.; Pimpinelli, A. Appl. Phys. Lett. 2006, 89, 104108−1−3. (28) Yu, Y-M; Voigt, A.; Guo, X.; Liu, Y. Appl. Phys. Lett. 2011, 99, 263106−1−3. (29) Verga, A. arXiv:1207.4354v1 [cond-mat-stat-mech].

possible version of the MC model of the GaN system. It allows us to analyze the behavior of a large, two-dimensional system, and we show that the main temperature dependence of the system behavior is controlled by a few model parameters only. The work under additional model elements is under study. Thus, the main types of experimentally observed surface morphologies were reproduced by Monte Carlo simulations employing a lattice model of GaN crystal with many body interactions. The tendency of the formation of particular surface patterns depends strongly on the temperature. When the heating temperature is too high, steps tend to group in bunches and form wavy patterns. There is only a small window of annealing temperatures, for which a crystal surface attains a parallel, equally distanced, somewhat bent, characteristic curly step structure. Such a structure is the best initial state for growing regular, microscopically flat crystalline layers for application in electronic devices. In this work we compare surface patterns created during sublimation of the model and the real system.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Research supported by the National Science Centre (NCN) of Poland (Grant NCN No. 2011/01/B/ST3/00526). REFERENCES

(1) Nakamura, S.; Fasol, G.; Pearton, S. J. The Blue Laser Diodes; Springer: New York, 2000. (2) Skierbiszewski, C.; Wasilewski, Z. R.; Grzegory, I.; Porowski, S. J. Cryst. Growth 2009, 311, 1632−1639. (3) Weyher, J. L.; Mueller, S.; Grzegory, I.; Porowski, S. J. Cryst. Growth 1997, 182, 17−22. (4) Pearton, S. J.; Shul, R. J.; McLane, G. F.; Constantine, C. In Gallium Nitride and Related Materials First International Symposium; Ponce, F. A.; Dupuis, R. D.; Nakamura, S.; Edmond, J. A., Eds.; Materials Research Society: Pittsburgh, 1996; pp 717−722. (5) Sato, H.; Sugahara, T.; Maosheng, H.; Naoi, Y.; Kurai, S.; Yamashita, K.; Nishino, K.; Sakai, S. Jpn. J. Appl. Phys. 1998, 37, 626− 631. (6) Mönch, W. Semiconductor Surfaces and Interfaces: Springer-Verlag: Berlin, 1993. (7) Weyher, J. L.; Tichelaar, F. D.; van Dorp, D. H.; Kelly, J. J.; Khachapuridze, A. J. Cryst. Growth 2010, 312, 2607−2610. (8) Weyher, J. L.; Lazar, S.; Macht, L.; Liliental-Weber, Z.; Molnar, R. J.; Muller, S.; Sivel, V. G. M.; Nowak, G.; Grzegory, I. J. Cryst. Growth 2007, 305, 384−392. (9) Nowak, G.; Xia, X. H.; Kelly, J. J.; Weyher, J. L.; Porowski, S. J. Cryst. Growth 2001, 222, 735−740. (10) Garcia, S. P.; Bao, H.; Hines, M. Phys. Rev. Lett. 2004, 93, 166102−1−4. (11) Garcia, S. P.; Bao, H.; Hines, M. J. Phys. Chem. B 2002, 106, 8258−8264. (12) Saito, Y. Statistical Physics of Crystal Growth; World Scientific: Singapore, 1996. (13) Załuska-Kotur, M. A.; Krzyżewski, F.; Krukowski, S. J. Appl. Phys. 2011, 109, 023515−1−9. (14) Załuska-Kotur, M. A.; Krzyzewski, F.; Krukowski, S. J. Cryst. Growth 2012, 343, 138−144. (15) Załuska-Kotur, M. A.; Krzyżewski, F. J. Appl. Phys. 2012, 111, 114311−1−8. 1013


Structures Built by Steps Motion during Sublimation from Annealed

Recommend Documents