Crystal Truncation Rods from Miscut Surfaces with Alternating Terminations

A long-standing experimental challenge has been to identify the orientation of the {\alpha} and \{beta} terraces on basal plane surfaces of crystals with hexagonal close-packed and related structures. To demonstrate how surface X-ray scattering can be sensitive to such {\alpha} vs. \{beta} terminations, we develop a general theory for the intensity distributions along crystal truncation rods (CTRs) for miscut surfaces with a combination of two terminations. We consider fractional-unit-cell-height steps, and variation of the coverages of the terraces above each step. Example calculations are presented for the GaN (0001) surface with various reconstructions. These show which CTR positions are most sensitive to the fractional coverage of the two terminations. We compare the CTR profiles for exactly oriented surfaces to those for vicinal surfaces having a small miscut angle, and investigate the circumstances under which the CTR profile for an exactly oriented surface is equal to the sum of the intensities of the corresponding family of CTRs for a miscut surface.


I. INTRODUCTION
Crystal truncation rods (CTRs) are features in the scattering patterns of crystals that arise from the truncation of the bulk crystal lattice at a surface 1 . They consist of streaks of scattering intensity extending away from all Bragg peaks in the direction normal to the crystal surface. The intensity distributions along the CTRs provide a sensitive measure of the atomic structure of the surface. In particular, they provide information about the arrangement of atomic-scale steps on the surface, and the reconstruction of surface layers. In situ X-ray measurements of CTRs have been used to determine the atomistic mechanisms of crystal growth 2-8 , such as the classical homoepitaxial growth mechanisms of 1-dimensional step flow, 2-dimensional island nucleation and coalescence, and 3-dimensional roughening 9,10 . Because of their penetrating nature and atomic-scale sensitivity, X-ray methods are uniquely suitable for in situ studies in non-vacuum environments. Here we calculate the behavior of CTR intensities to address a longstanding issue in crystal growth, the variation in the morphology of surfaces having alternating terminations and step types, driven by differences in step kinetics during step-flow growth.
The standard model of a miscut (vicinal) crystal surface consists of a stair-like sequence of identical terraces, separated by full-unit-cell-height steps, as shown in Fig. 1(a). The average step spacing is determined by the miscut angle of the surface (i.e. the angle of the surface away from the closest low-index crystallographic orientation). During stable step flow growth, the terrace widths are typically all nearly equal, since the steps have identical kinetic properties. However, when the space group of the crystal includes screw axes or glide planes, the terrace and step arrangement and resulting growth behavior on surfaces perpendicular to one of these symmetry elements can have fundamentally different characteristics 11 . In this case, each succeeding terrace has the same atomic arrangement, but a different in-plane orientation. These terraces are separated by fractional-unit-cell-height steps, and the step structures and properties can vary from step to step, as shown in Fig. 1(b). A commonly encountered version of this effect occurs on the basal-plane {0001}-type surfaces of crystals having hexagonal close-packed (HCP) or related structures, which are normal to a 6 3 screw axis. The close-packed layers in HCP crystals have 3-fold symmetry alternating between 180 • -rotated orientations from layer to layer, as shown by the α and β terrace structures in Fig. 2. The αβαβ stacking sequence typically results in half-unit-cellheight steps on miscut surfaces. Often the lowest energy steps are normal to [0110]-type directions. The alternating structures of these steps are conventionally labelled A and B 10,12 as shown on Fig. 3. When the in-plane azimuth of an A step changes by 60 • , e.g. from [0110] to   12,19,22,[26][27][28][29][30][31] . In particular, different attachment kinetics at A and B steps can produce a tendency to step pairing during growth and thus to different local fractions of α and β terraces. In limiting cases, the α terrace fraction f α can approach zero or unity, when pairs of half-unit-cell-height steps join to form full-unitcell-height steps. However, for {0001} surfaces of HCPrelated systems it has been difficult to distinguish experimentally the terrace orientation, and thus to determine whether a given set of steps is of A or B type. Thus there is a clear need for a method for experimental determination of the difference between adatom attachment kinetics at A and B steps, especially an in situ measurement in the relevant growth environment.
Streaks of scattering intensity extending away from Bragg peaks normal to the surfaces of finite-sized crystals are inherent in exact treatments of X-ray scattering, extending back to early work 32,33 . As synchrotron sources enabled detailed study of these surface-sensitive features, theoretical treatments specific to CTRs from truncated perfect crystals were developed 1,34 . Subsequent work included the effect of overlayers or reconstructions 35 and miscut surfaces [36][37][38][39][40] . Previous general treatments of miscut surfaces typically considered only full unit cell steps; half-unit cell steps on Si (001) have also been considered 37,39 .
Here we develop expressions for the CTR intensity distributions for surfaces with two terrace types, α and β. As in previous calculations we include the effect of surface reconstruction and surface roughness. New aspects of our calculations include a general treatment of fractionalunit-cell-height steps, with fractional terrace coverages; generalization from a cubic to an orthorhombic bulk unit cell, which allows us to consider HCP type crystals using orthohexagonal coordinates; and the effects of absorption and extinction. We use a mathematical formalism that emphasizes how exact results can be obtained through the summation of geometric sequences. We compare the CTR profiles for exactly-oriented and miscut surfaces, and consider the relationships between the two cases.
To provide a concrete example, we make calculations for the GaN (0001) surface, for which we have recently carried out a surface X-ray scattering study 41 . We show how the CTRs vary as a function of the surface fraction of α vs β terraces. This demonstrates that surface X-ray scattering can be used to distinguish the fraction of the surface covered by α or β terraces during growth, and thus unambiguously determine differences in the attachment kinetics at A and B steps.

II. SURFACE X-RAY SCATTERING THEORY
In this section we develop expressions for the intensity distributions along CTRs and demonstrate their sensitivity to the surface termination arrangements. The X-ray reflectivity along the CTRs can be calculated by adding the complex amplitudes from the bulk crystal and the terminating overlayers, with proper phase relationships. We start with the simple case of an exactly-oriented surface, and then extend this to the case of a vicinal surface having an array of straight steps. After developing general expressions, we show calculations for the GaN (0001) surface for both cases.
Although our motivation for this work is based on hexagonal crystals, it is more convenient to use the Cartesian coordinates of an orthorhombic unit cell. As shown in Figs. 2 and 4, the hexagonal lattice can be mapped to an orthorhombic lattice using a doubled unit cell and an orthohexagonal coordinate system 42 . The lattice parameters of the orthorhombic unit cell are a, b, and c in the x, y, and z directions, respectively. The corresponding Cartesian components of the scattering wavevector Q are given by (Q x , Q y , Q z ), which are related to reciprocal lattice units HKL by Q x = (2π/a)H, Q y = (2π/b)K, For a hexagonal crystal, the orthohexagonal coordinate system has a lattice parameter b ≡ a √ 3, where a is the in-plane lattice parameter of the conventional hexagonal unit cell, as shown in Fig. 2. The out-of-plane lattice parameter c is the same in both coordinate systems. This gives Cartesian x, y, and z axes parallel to the hexagonal [2110], [0110], and [0001] directions, respectively. In reciprocal space, the orthohexagonal coordinates HKL are related to the hexagonal Miller-Bravais reciprocal space coordinates hki by H = h, K = h + 2k, L = , as shown in Fig. 4. In particular, the hexagonal (1010) and (0110) correspond to the orthohexagonal (110) and (020), respectively. Because of the doubled unit cell, orthohexagonal indices H and K must both be even or odd integers at allowed Bragg peaks. For the remainder of this paper, we will use 3-index orthohexagonal indices HKL rather than 4-index hexagonal Miller-Bravais indices hki .

A. Exactly oriented surface with two terminations
For an exactly oriented surface normal to the z direction, the CTRs extend continuously in the Q z direction at fixed Q x and Q y through each Bragg peak. The CTRs from the Bragg peaks of different L 0 at the same H 0 K 0 all overlap, as shown in Fig. 5(a).
The contribution to the complex amplitude r b of the reflectivity from the truncated bulk crystal below the terminating overlayers can be obtained by coherently summing the scattering from each unit cell in a column normal to the surface, taking into account a factor Z difference between the scattering amplitude from each unit cell 35 . The sum of this geometrical sequence is given by where r f ≡ 4πir 0 /AQ, r 0 = 2.817 × 10 −13 cm is the Thomson radius of the electron, A is the in-plane area of the unit cell, and Q is the magnitude of the wavevector. The bulk structure factor F b is Here the first sum is over the k chemical elements present in the crystal, is the atomic form factor f k (Q) modified by the Debye-Waller thermal vibration length u k for element k, the second sum is over the n bulk atoms of type k in a unit cell, and r b kn is the position of bulk atom n of type k. The quantity Z in Eq. (1) is the ratio of the scattering amplitude from one unit cell to that from the unit cell at ∆z = −c below it. It consists of a phase factor and an absorption factor, where = 4π/(λ abs ) is related to the photon wavelength λ and absorption length abs . Typically the absorption factor is very close to unity and only becomes important in Z − 1 when the phase factor is unity, e.g. at Bragg peaks. One can see from Eq. (1) that the reflectivity amplitude is built up by summing the scattering from each layer of the semi-infinite crystal in the z direction from m = −∞ to m = 0.
To provide a direct comparison with the results for miscut surfaces obtained below, we allow the surface to have a mixture of two terminations, α and β. For an exactly oriented surface, this consists of islands of α in a matrix of β, or vice versa, as shown in Fig. 6. We also allow reconstructions in which the atoms in the top layer of unit cells at the surface are relaxed from their bulk crystal positions, and there can be extra atoms bonded to the surface. The reconstructed reflectivity amplitude per unit area of the x = α or β overlayer is where the structure factor of the reconstruction F x r is Here the first sum is over the j possible domain orientations of the reconstruction, θ xj is the fraction of domain j for the x = α or β component, the second sum is over the k chemical elements present in the reconstruction, the third sum is over the n atoms of type k in a unit cell, and r jkn is the position of atom n of type k in domain orientation j. The total reflectivity amplitude is the sum of the complex amplitudes from the bulk and both reconstructed components,  where f α is the fraction of the surface covered by the α termination. The total reflectivity amplitudes calculated above are for the kinematic limit in which they are much smaller than unity. Near the Bragg peaks, where the reflectivity amplitude approaches unity, the kinematic expressions can be corrected using which insures that the amplitude of the reflectivity does not exceed unity. The intensity reflectivity is the square of the modulus of the amplitude reflectivity, where the final factor S has been introduced to account for surface roughness. Here we use a form for the roughness factor as a function of L for the exactly oriented surface adopted from that obtained for a miscut surface 40 , where σ R is the surface roughness (related to the total step roughness given in 40 by σ R = cσ tot ), and the final approximation holds for σ R < 0.3c. Figure 7 shows plots of this function for various values of σ R /c. Away from the Bragg peaks, so that the dynamical scattering correction can be neglected, the intensity reflectivity can be written as Here we have neglected an overall factor of |Z| 2 , which can be set to unity when absorption is small. To illustrate the behavior for the basal plane of an HCP-type system, we consider Ga-face GaN (001) surfaces with a Ga termination for the bulk, so that the top-layer atoms shown in Fig. 3 are Ga. Table I in Appendix A lists the atomic coordinates used for the bulk. We use atomic coordinates for the reconstructed overlayers calculated previously 43 . Since these were calculated using a 2 × 2 unit cell, for consistency the unit cell sums used in calculating the structure factors here are carried out over two adjacent orthohexagonal unit cells having an area A = 2ab, which is normalized out in the denominator of r f . To compute the scattering from surfaces terminated at α and β terraces, we terminate the bulk at a β terrace, and incorporate an extra half unit cell of atoms in their bulk positions into the bottom of the reconstructed overlayer for the α terrace regions. We also reverse the relaxation amounts in the y direction for the α terraces, relative to those for the β terraces. Figure 6 illustrates these arrangements. Appendix A gives tables of typical atomic coordinates r x jkn used for the x = α and β structure factors. For all calculations presented in this paper, we use a photon energy of 25.78 keV (λ = 0.4809Å) to correspond with recent experiments 41 . We use atomic form factors for each type of atom 44 with resonant corrections for this energy 45 , an absorption length of abs = 103 µm for GaN at this energy, a Debye-Waller length of u k = 0.16Å for all atoms, and a surface roughness of σ R = 1Å. Figure 8 shows the calculated reflectivity as a function of L for different integer H 0 K 0 values and three α fractions, f α = 0, 0.5 and 1. Here we show calculations for an unreconstructed surface, with all atoms in their bulk positions. The (00L) CTR is insensitive to the difference between the α and β terminations (f α = 1 and f α = 0, respectively); both give the same intensity distribution. Furthermore, the (02L) scattering from the α terrace is identical to the (11L) scattering from the β terrace, and vice versa, as required by symmetry. For surfaces with mixed terrace coverages, the (02L) CTR for f α = X and the (11L) CTR for f α = 1 − X are similar but not identical. This can be seen for the f α = 0.5 case shown, where the minima between the Bragg peaks all have about the same depth. The small difference between the two CTRs is determined by whether the α terraces or β terraces are the top layer. The calculations shown here are for α terraces on top. We have performed calculations using atomic coordinates for all of the GaN (001) reconstructions found previously 43 . On this surface, there are 6 reconstruction domain orientations, related by 3-fold rotation about the 6 3 axis, and/or reflection about a plane passing through the axis at x = 0. Figure 9 shows calculations for f α = 0 with equal fractions θ xj = 1/6 of all six domains. All show the same qualitative behavior as the unreconstructed surface, with small quantitative differences. Furthermore, because the X-ray scattering is dominated by the Ga atoms, which occupy an HCP lattice, the same qualitative behavior found for GaN is also obtained for an elemental HCP crystal.

C. Miscut surface with alternating terminations
We now consider a miscut surface, with an array of straight steps parallel to the x axis, distributed along the y axis, as shown in Fig. 1. We assume that the sur-face height decreases by a full unit cell c every M unit cells in y, so that the period of the step array is M b. The surface miscut angle γ relative to (001) is given by tan γ = c/(M b), and the surface is parallel to (01M ) planes. As shown in Fig. 5(b), the CTRs from this surface are tilted in the Q y direction at an angle γ from (001). Because of the tilt, there are M times as many CTRs as in the exactly oriented case, indexed not just by H 0 K 0 but also by values of L 0 from 0 to M − 1. The Q y value varies with L along the CTR according to Q y = (2π/b)[K 0 + (L − L 0 )/M ], where H 0 K 0 L 0 is the primary Bragg peak associated with the CTR. The spacing in L along a given CTR between Bragg peak positions is M , rather than unity as in the exactly oriented surface. (This neglects systematic absences in the Bragg peaks; for the orthohexagonal case, the spacing between allowed peaks is 2M .) We assume that the surface has alternating terminations of α and β terraces along the y axis. Figure 10 shows the bulk and reconstructed unit cells used to calculate the CTRs for the vicinal surface. The width of the α terraces is N unit cells, and the width of the β terraces is M − N unit cells. The α terrace fraction is given by For a miscut crystal, the bulk reflectivity amplitude can be built up by summing the scattering from each unit cell in the y direction from m = −∞ to m = 0, according to where the quantity Y is the ratio of the scattering amplitude from one unit cell to that from the unit cell at ∆y = −b beside it, made up of a phase factor and an absorption factor This is similar to the geometric sequence for the exactly oriented surface, Eq. (1), except that the summation is carried out in the y direction rather than the z direction, as indicated in Fig. 10. Substituting the above expression for Q y along the CTR into the phase factor, the term involving K 0 drops out. This gives Comparison with Eq. (3) shows that Y M = Z. The reflectivity amplitude per unit area of the reconstructed overlayer on the α terraces can be written as A similar expression applies to the β terraces, In these equations the sums run over the fraction of the surface periodicity covered by each termination. Like the exactly oriented surface, the total reflectivity amplitude r t is the sum of the complex amplitudes from the bulk and the reconstructed layers on the α and β terraces. Expressions similar to Eqs. (6)(7)(8) give the intensity reflectivity R L0 for the CTR from the H 0 K 0 L 0 Bragg peak, with the roughness factor for each CTR on a vicinal surface now given by a simple Gaussian form 34,40,46 , Away from Bragg peaks, so that the dynamical scattering correction can be neglected, the intensity reflectivity can be written as where we have neglected an overall factor of |Y | 2 .  gives reflectivities that agree within ∼ 10% at all L. As in the case of an exactly oriented surface, the (00L 0 ) CTRs are identical for f α = 0 and f α = 1, and they are very different for f α = 0.5, with the CTRs for even L 0 becoming stronger and the CTRs for odd L 0 becoming very weak. The (02L 0 ) and (11L 0 ) CTRs have a more monotonic dependence on f α . For f α = 0 and f α = 1, there are alternating stronger and weaker intensities between the Bragg peaks, with the alternation being opposite for (02L 0 ) and (11L 0 ). For f α = 0.5, the intensities between the Bragg peaks are about the same, and there is no difference between the (02L 0 ) and (11L 0 ) CTRs. The (02L 0 ) CTRs with f α = X are identical to the (11L 0 ) CTRs with f α = 1 − X, for any value X. As with the exactly oriented surface, other GaN reconstructions have the same qualitative behavior, as shown in Fig. 12, as do elemental HCP structures.

E. Relation between CTR profiles for miscut and exactly oriented surfaces
For a miscut surface with full-unit-cell-height steps and a single termination, e.g. N = 0 or N = M , one can show that the intensity reflectivity R for an exactly oriented surface is simply the sum of the intensity reflectivities R L0 for all the CTRs with the same H 0 K 0 , in the limit of small miscut angle and neglecting absorption, For N = 0 or N = M (x = β or α, respectively), using Y M = Z, Eq. (17) becomes Neglecting the absorption factors in Z and Y , one obtains |Z −1| 2 = 4 sin 2 (πL) and |Y −1| 2 = 4 sin 2 [π(L−L 0 )/M ]. For small miscut angle (large M ), the second expression becomes |Y −1| 2 = 4π 2 (L−L 0 ) 2 /M 2 . Substituting these into Eq. (19) gives (20) While the first two factors depend upon Q and thus vary slightly with L 0 at fixed L due to the spacing of the CTRs in K, for small miscut angle this is negligible, and these factors are the same as for an exactly oriented surface. Comparison with Eqs. (8), (9), and (16) shows that the summation Eq. (18) holds. This is illustrated in Fig. 13(a), where the summation over L 0 of the miscut CTR intensities R gives the same result as the calculation for the exactly oriented surface for f α = 0. Fig. 13(b) shows that the summation does not hold for the mixed termination, f α = 0.5.
It is interesting that, for the single termination surface, the exactly oriented CTR is equal to the sum of the intensities of the miscut CTRs, not the square of the sum of the amplitudes of the miscut CTRs. This is demonstrated in Fig. 13(a). From an experimental perspective, this means that there is no difference between the CTR profile measured from an exactly oriented surface and that from a surface with a miscut so small that the splitting of the CTRs cannot be resolved. For the mixed termination surface, Fig. 13(b), the exactly oriented CTR is not equal to the miscut CTRs summed in either intensity or amplitude.

III. DISCUSSION AND CONCLUSIONS
The analysis presented above shows that CTR profiles are very sensitive to the difference between α and β terraces on vicinal basal plane surfaces of HCP-type systems. This enables in situ X-ray measurements during growth to determine the surface fraction covered by each type of terrace, and thus to distinguish the dynamics of A and B steps. Measurement of the CTR profiles as a function of L can be carried out by running scans during steady-state growth conditions. In addition, the dynamics of the terrace fraction f α (t) can be observed by monitoring CTR intensities at particular L values during changes in growth conditions. Figure 14 shows calculations of the reflectivity as a function of f α at fixed These curves can be used to convert observed intensity changes into terrace fraction dynamics f α (t), e.g. as growth or evaporation rate is changed during step-flow growth. Figure 15 shows calculations of the reflectivity as a function of f α at fixed L positions on three CTRs for an exactly oriented surface. The (00L) with L = 0.99 position shows behavior qualitatively similar to the (001) CTR of a miscut surface, with a minimum at f α = 0.5. This position (or the equivalent (13L) position) has been used to monitor growth under island nucleation and coalescence conditions, where oscillations in intensity occur with a period corresponding to growth of half-unit-cell monolayers 5,8 . Figure 15 also shows how the sum of the CTR intensities for a miscut surface, which agree with the exactly oriented surface at f α = 0 and 1, deviate at intermediate values of f α . This analysis extends previous calculations of CTR profiles for miscut surfaces to include the case of fractional unit cell steps and terrace coverages. While the example CTR calculations presented here are for wurtzitestructure GaN, the ability to distinguish α and β terrace fractions applies directly to many other HCP-type systems with a 6 3 screw axis, including other compound semiconductors, as well as one third of the crystalline elements and many more complex crystals.

ACKNOWLEDGMENTS
Work supported by the U.S Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Materials Science and Engineering Division.

Appendix A: Atomic coordinates
To provide a detailed example of how we calculate the CTR intensities including the effects of reconstruction, we here provide an example of the atomic coordinates for a particular reconstruction. The qualitative behavior we observe does not depend upon the reconstruction chosen or the exact values of the atomic coordinates used. Tables I, II, and III give the atomic coordinates for the 3H(T1) reconstruction obtained in 43 . The fractional coordinates x, y, and z given in the tables are the components of the positions r kn , r α jkn , and r β jkn used to calculate the structure factors, normalized to the respective orthohexagonal lattice parameters a, b, and c, i.e. r = (ax, by, cz). A 2 × 2 surface unit cell is used, equivalent to two orthohexagonal unit cells, so there are 8 Ga and 8 N sites in each. These coordinates place a bulk Ga site on a β layer at the origin. We use u = 0.3768 for the internal lattice parameter of bulk GaN, i.e. the frac-tional distance between Ga and N sites, which deviates slightly from the ideal 3/8 value, in agreement with ab initio calculations 43,47 and experiments 48 . Relaxed positions were calculated for a one-unit-cell thick layer at the surface. For the α terrace, and extra half unit cell of bulk (unrelaxed) atoms is attached to the bottom to account for the difference in height of the α and β terraces, as shown in Fig. 10. Coordinates for only one domain are given. Those for other 5 domains are obtained by 3-fold rotation about the 6 3 axis and/or reflection of the y coordinate. One can see that the Ga atoms bonded to the three adsorbed hydrogens of the 3H(T1) reconstruction relax to higher z positions.  Fractional coordinates x, y, z of atoms in domain j = 1 of the 3H(T1) reconstruction used to calculate the α terrace contribution to the CTRs, as well as their differences ∆x, ∆y, ∆z relative to bulk lattice positions. The differences for H atoms are relative to N sites. The lowest four Ga and N sites are an extra half unit cell of bulk lattice to account for the difference in height of the α and β terraces.
Atom Site