Calculated surface X-ray scattering with alternating step types

We present calculations showing how the intensity distribution along the CTRs is sensitive to the fraction of the surface covered by α or β terraces. Our calculations include the effect of surface reconstruction, using relaxed atomic coordinates that have been obtained previously²⁸. For a vicinal surface, the CTRs are tilted away from the crystal axes, so that the CTRs from different Bragg peaks do not overlap. The X-ray reflectivity along the CTRs can be calculated by adding the complex amplitudes from the substrate crystal and the reconstructed overlayers, with proper phase relationships^29–31. Details of our calculations are given in a separate paper³².

Figure 3a shows calculated intensity distributions along $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$ CTRs for the GaN (0001) surface, demonstrating how their L-dependences vary with f_α. Bragg peak locations have integer indices H₀K₀L₀ in reciprocal lattice units; these indices identify the CTR associated with each peak. For this comparison, we use the 3H(T1) surface reconstruction²⁸ and a fixed surface roughness, independent of f_α, as discussed in the “Methods” section below and Supplementary Discussion 1. The same qualitative behavior is obtained using other surface reconstructions. For f_α = 0 and f_α = 1, the regions between the Bragg peaks have alternating stronger and weaker intensities, with the alternation being opposite for $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$. For f_α = 0.5, the intensities between all Bragg peaks are about the same, and there is no difference between the $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$ CTRs. As required by symmetry, the $$ \left(01\overline{1}L\right)$$ CTRs with f_α = X are identical to the $$ \left(10\overline{1}L\right)$$ CTRs with f_α = 1 − X, for any value X. Figure 3b shows calculations of the reflectivity as a function of f_α at positions near L = 1.6 on the $$ \left(01\overline{1}2\right)$$ and $$ \left(10\overline{1}2\right)$$ CTRs. The variation in reflectivity is almost monotonic in f_α at these positions. These curves are used below to extract f_α(t) during dynamic transitions.

Fig. 3

Calculated CTR intensities for a vicinal GaN (0001) surface.

a Top and bottom rows show the families of CTRs along $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$, respectively. Black, red, green, blue, cyan, and magenta curves are for CTRs from the L₀ = −1 to 4 Bragg peaks, respectively. Values of f_α for each column are given at the top. b Reflectivity of selected CTRs as a function of terrace fraction f_α for fixed values of L = 1.6.

In situ X-ray scattering measurements during growth

We studied four OMVPE conditions having different net growth rates at the same temperature 1076 ± 5 K, summarized in Table 1 (see “Methods” for further details). The substrate used was a GaN single crystal. Figure 4a shows its initial surface morphology determined by ex situ atomic force microscopy (AFM). One can see straight steps almost perpendicular to the y or $$ \left[01\overline{1}0\right]$$ direction over large areas. An analysis of the step spacing shows a slight tendency towards pairing, with one of the two alternating terrace types having an area fraction of 0.47. AFM is insensitive to whether this fraction corresponds to the α or β terraces. We also characterized the miscut angle by measuring the splitting of the CTRs. Figure 4b shows a transverse cut through the CTRs in the Q_y direction near (000L) at L = 0.9. Both the AFM and X-ray measurements give a double-step spacing of w = 573 Å corresponding to a miscut angle of 0.52°. To relate the α terrace fraction to the behavior of A and B steps, it is critical to determine the sign of the step azimuth. By making measurements as a function of L, we verified that the peak at high Q_y is the CTR coming from (0000), while the peak at low Q_y is the (0002) CTR. This confirms that the downstairs direction of the vicinal surface is in the +y direction, as drawn in Fig. 1. It is also useful to know the precise angle of the step azimuth with respect to the crystal planes, which determines the minimum kink density and thus the predicted values of some kinetic coefficients. X-ray measurements found this to be 5° off of the $$ \left[01\overline{1}0\right]$$ direction towards $$ \left[10\overline{1}0\right]$$, which gives a maximum kink spacing of 33 Å. The kink spacing could be smaller due to thermally generated kinks¹, as discussed in Supplementary Discussion 2. With this low-dislocation-density substrate and the low growth rates used, the previously reported instability to step bunching during growth³³ was not observed.

Table 1

Growth conditions studied.

Growth condition index	TEGa flow (μmol min−1)	H2 frac. in carrier (%)	Net growth rate G (monolayer s−1)		Measured values from CTR fits	Best fit from BCF model
1	0.000	50	−0.0018	$$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$	0.111 ± 0.013	0.136
2	0.000	0	0.0000	$$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$	0.461 ± 0.018	0.440
3	0.033	50	0.0109	$$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$	0.811 ± 0.014	0.836
4	0.033	0	0.0127	$$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$	0.868 ± 0.011	0.847
1–2				trel	2200 ± 200 s	2478
2–4				trel	340 ± 30 s	331

Fig. 4

Imaging and scattering from steps.

a AFM image of steps of height c/2 (2.6 Å) on the vicinal GaN substrate used in X-ray measurements. To emphasize the positions of steps, we plot the amplitude error signal, which is proportional to the height gradient in the scan direction (y). Image was obtained ex situ at room T after an anneal for 300 s at 1118 K in zero-growth conditions (0% H₂, 0 TEGa). Average fraction over a 2 × 2 μm² area of terrace family marked black at top is 0.47. The average double-step spacing of w = 573 Å corresponds to a miscut angle of $$ {\tan }^{-1}\left(c/w\right)=0.52^{\circ }$$. b Profile of split CTRs from the vicinal surface, measured at T = 1170 K during growth at 0.053 μmol min⁻¹ TEGa and 50% H₂. The splitting of the CTRs ΔQ_y = 0.0110 Å⁻¹ corresponds to a miscut angle of $$ {\tan }^{-1}\left[\mathrm{\Delta }{Q}_y/\left(2\pi /c\right)\right]=0.52^{\circ }$$, in agreement with the AFM result.

Figure 5 shows the measured steady state CTR intensities as a function of L, for both the $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$ CTRs and at all four conditions. The qualitative behavior agrees with that expected from a variation in f_α shown in Fig. 3a, with alternating higher and lower intensities between the Bragg peaks under some conditions, and opposite behavior of the two CTRs.

Fig. 5

Measured and calculated crystal truncation rods.

Symbols show measured net intensities of the $$ \left(01\overline{1}{L}_0\right)$$ CTRs and the $$ \left(10\overline{1}{L}_0\right)$$ CTRs families (left and right) for CTRs from the L₀ = 0, 1, 2, and 3 Bragg peaks at each of four conditions. Curves show fits of all CTRs to obtain steady state α terrace fraction $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ at each condition.

Steady state and dynamics of terrace fraction during growth

To obtain values of the steady state terrace fraction $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ for each of the four conditions, we fit calculated CTR intensities to the measured profiles. We performed fits using different possible surface reconstructions (see “Methods” for further details). While the 3H(T1) reconstruction gives the best fit to all conditions, similar $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ values are obtained using alternative reconstructions. For each condition, both the $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$ CTRs were simultaneously fit. The values of $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ obtained as a function of net growth rate G are given in Table 1. The marked increase in $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ as G is increased reveals the qualitative difference between the kinetics at A and B steps during OMVPE of GaN: adatom attachment coefficients for A steps are larger. Thus, a surface with initially balanced α and β terrace fractions at zero net growth rate will evolve to one with higher $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ during positive net growth, because of the initially higher adatom attachment rate at the A steps. Likewise, during evaporation the initially higher detachment rate at A steps will give a lower $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$.

We also observed the dynamics of the change in f_α by recording the intensity at a fixed detector position as a function of time before and after an abrupt change between conditions, as shown in Fig. 6a. We chose positions near L = 1.6 where the X-ray reflectivity R changes almost monotonically with f_α, as shown in Fig. 3b. It is thus straightforward to obtain f_α(t) from the intensity evolution using the calculated R(f_α), as shown in Fig. 6b. The characteristic 1/e relaxation times t_rel were 2200 ± 200 and 340 ± 30 s for the transitions from conditions 1 to 2 and 2 to 4, respectively.

Fig. 6

Dynamics following a change of condition at t = 0.

a Measured CTR intensities. b Calculated terrace fractions f_α. Blue curves: conditions 1–2 on the $$ \left(10\overline{1}2\right)$$ CTR at L = 1.627. Red curves: conditions 2–4 on the $$ \left(01\overline{1}2\right)$$ CTR at L = 1.603. Circles show 1/e relaxation times t_rel of 2200 ± 200 and 340 ± 30 s, respectively.

BCF model for surface with alternating step types

To quantitatively relate the behavior of the terrace fraction to the kinetic properties of A and B steps, we have developed a model³⁴ based on BCF theory. Such models have been used extensively to understand growth behavior such as the step-bunching instability³⁵, pairing of steps³⁶, and competitive adsorption³⁷, typically where all steps in a sequence are equivalent. In our model, we consider an alternating sequence of two types of terraces, α and β, and two types of steps, A and B, with properties that can differ, as shown in Fig. 7. Related BCF models with an alternating step or terrace properties have appeared previously^{12,16,17,38–40}. We include the effects of step transparency⁴¹ (i.e., adatom transmission across steps) and step–step repulsion².

Fig. 7

Schematic of alternating terraces and steps for the BCF model.

Vicinal {0001} surfaces of HCP crystals have alternating α and β terraces separated by A and B steps. Notations are indicated for the kinetic coefficients for adatom attachment from below and above and for adatom transmission. Schematic of alternating terraces and steps of BCF model for HCP basal-plane surfaces, showing kinetic coefficients for the A and B steps.

The rate of change in the adatom density per unit area ρ_i on terrace type i = α or β is written as

To solve this model, we develop a quasi-steady state expression for the dynamics of the terrace fraction f_α. Under fairly general assumptions³⁴, the behavior of f_α can be written as a function of the net growth rate,

The full steady state df_α/dt = 0 is obtained at a growth rate of

To calculate BCF model results to compare with the experimental conditions, we make three assumptions: (1) the only parameter affected by the TEGa supply rate is the deposition flux F; (2) the only parameter affected by the carrier gas composition (0 or 50% H₂) is the adatom lifetime τ; and (3) F and τ enter only through the net growth rate G given by Eq. (8), listed in Table 1 for each condition. We use the known values $$ {\rho }_0=2a^{-2}/\sqrt{3}=1.13\times 10^{19}$$ m⁻² and $$ w=c/\sin \left(0.52^{\circ }\right)=5.73\times 10^{-8}$$ m, where a = 3.20 × 10⁻¹⁰ m and c = 5.20 × 10⁻¹⁰ m are the lattice parameters of GaN at the growth temperature⁴². We performed fits to the measured quantities (four $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ and two t_rel) using the general expressions for K^ss(f_α) and K^dyn(f_α)³⁴. The relaxation times for the model were calculated by integrating Eq. (9) to obtain f_α(t), and then extracting the relaxation time with the same normalization procedure used for the experimental data. The best fit, shown in Fig. 8, was obtained in the limit in which the $$ {\kappa }_{+}^A$$ coefficient is large, while the $$ {\kappa }_{-}^A$$, $$ {\kappa }_{-}^B$$, and $$ {\kappa }_0^A$$ coefficients approach zero. In this case, K^ss(f_α) and K^dyn(f_α) can be written as³⁴

Fig. 8

Steady state α terrace fraction $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ as a function of growth rate G.

Circles are experimental values given in Table 1, showing monotonic increase of $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ with increasing G. Curve is best-fit BCF model calculation; also fit simultaneously with these four points were the two relaxation times t_rel given in Table 1.

Calculated CTR intensities

Derivation of the CTR calculations for miscut surfaces with alternating terrace terminations is provided in a separate paper³². As described in Supplementary Discussion 1, in fits shown in Table 2, and in recent predictions⁴⁹, we expect that the GaN surface under OMVPE conditions has a 3H(T1) reconstruction, in which 3 of every 4 Ga atoms in top-layer sites shown in Fig. 2 are bonded to an adsorbed hydrogen. The calculations in Fig. 3a include the effect of this reconstruction, with equal fractions of all reconstruction domains on both terraces. Since no fractional-order diffraction peaks from long-range-ordered reconstructions are observed in experiments, we expect that the domain structure has a short correlation length and all domains are present. We use a surface roughness of σ_R = 0.74 Å to match the experimental fits described below.

Table 2

Values of fit parameters for each of four OMVPE conditions.

Growth condition index	Reconstruction	$$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$	σR	χ2
1	3H(T1)	0.111	0.75	106
	Ga(T4)	0.144	1.26	130
	NH(H3) + H(T1)	0.098	0.94	187
	NH(H3) + NH2(T1)	0.106	0.88	200
	NH(H3)	0.095	0.91	167
2	3H(T1)	0.461	0.93	57
	Ga(T4)	0.476	1.26	81
	NH(H3) + H(T1)	0.460	1.15	76
	NH(H3) + NH2(T1)	0.460	1.11	67
	NH(H3)	0.459	1.13	99
3	3H(T1)	0.811	0.85	118
	Ga(T4)	0.670	1.46	218
	NH(H3) + H(T1)	0.876	1.19	205
	NH(H3) + NH2(T1)	0.869	1.15	248
	NH(H3)	0.869	1.16	168
4	3H(T1)	0.868	0.47	80
	Ga(T4)	0.942	1.06	112
	NH(H3) + H(T1)	0.892	0.90	174
	NH(H3) + NH2(T1)	0.879	0.85	220
	NH(H3)	0.891	0.87	135

In situ microbeam X-ray experiments

We performed in situ measurements of the CTRs in the OMVPE environment at the Advanced Photon Source beamline 12ID-D⁵⁰. At an incidence angle of 2°, the 10 μm X-ray beam illuminated an area of 10 × 300 μm. To obtain sufficient signal, we used a wide-bandwidth pink beam setup^51,52. Details of the measurements and data analysis are given in Supplementary Methods 1. Two types of measurements were performed. We determined the steady state terrace fractions $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ under four different growth/evaporation conditions by scanning the detector along the $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$ CTRs while continuously maintaining steady state growth or evaporation. We also observed the dynamics of f_α after an abrupt change in conditions.

Under the conditions studied, deposition is transport limited, with the deposition rate proportional to the supply of the Ga precursor (triethylgallium, TEGa), with a large excess of the N precursor (NH₃) constantly supplied. We investigated conditions of zero deposition (no supply of TEGa) as well as deposition at a TEGa supply of 0.033 μmol min⁻¹. The NH₃ flow in both cases was 2.7 s.l.p.m. or 0.12 mol min⁻¹, and the total pressure was 267 mbar. The V/III ratio during deposition was thus 3.6 × 10⁶. For both of these conditions, we studied two carrier gas compositions: 50% H₂ + 50% N₂ and 0% H₂ + 100% N₂. The addition of H₂ to the carrier gas enhances evaporation of GaN, so that the net growth rate (deposition rate minus evaporation rate) is slightly lower; at zero deposition rate, the net growth rate is negative. We determined the net growth rate for all four conditions as described in Supplementary Methods 2. These values are given in Table 1. Substrate temperatures were calibrated to within ±5 K using laser interferometry from a standard sapphire substrate⁵⁰. While we used the same heater temperature for all conditions, the calibration indicates that the substrate temperature was slightly higher in 50% H₂ (1080 K) than in 0% H₂ (1073 K). Based on an expected activation energy of  ~ 0.4 eV for surface kinetics, this temperature difference produces a negligible (few percent) change in kinetic coefficients.

Fits to steady state CTRs

For each growth condition, fits were performed of the calculated CTR intensities to the measured profiles of the $$ \left(01\overline{1}L\right)$$ and $$ \left(10\overline{1}L\right)$$ CTRs. In addition to a single value of f_α, parameters varied in the fit included a single surface roughness σ_R and two intensity scale factors, one for each CTR. We fit to $$ \log \left(I\right)$$ with equal weighting of all points. Fits were done using different possible surface reconstructions. Supplementary Fig. 10 shows the calculated reconstruction phase diagram for the GaN (0001) surface in the OMVPE environment²⁸, as a function of Ga and NH₃ chemical potentials. Based on the chemical potential values that correspond to our experimental conditions estimated in Supplementary Discussion 1, shown by the green rectangle, we considered the five reconstructions highlighted in Supplementary Fig. 10. (The estimate for Δμ_Ga has a large uncertainty because it depends on the nitrogen potential produced by decomposition of NH₃.) Table 2 shows the values of steady state terrace fraction $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$, surface roughness σ_R, and the goodness-of-fit parameter χ² from fits to reflectivity at four conditions for each of five reconstructions. The qualitative results for the variation of $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ with growth condition are independent of which reconstruction is assumed: $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ increases monotonically as the net growth rate G increases. The 3H(T1) reconstruction gives the best fit (minimum χ²) of the five potential reconstructions, for all four conditions. This is consistent with recent results on GaN (0001) reconstructions in the OMVPE environment⁴⁹, which found a phase stability region for the 3H(T1) structure even larger than that shown in Supplementary Fig. 10. Figure 5 compares the fits with the 3H(T1) reconstruction to the measured CTR intensities.

Extracting terrace fraction dynamics

We first convert the intensity evolution near L = 1.6 to reflectivity R(t) by normalizing it to match the predicted change in reflectivity for the transition in $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$. We then invert the R(f_α) relation calculated for the experimental L value, shown in Fig. 3b, to obtain f_α(t). For simplicity, we assume the surface roughness is not a function of growth condition, and use the average value of σ_R = 0.74 Å from the 3H(T1) fits to calculate R(f_α). If the surface roughness were to vary with condition by the amounts shown in Table 2, this would explain only a small fraction of the observed changes (decrease of 7% out of 40% for the transition from conditions 1 to 2, and increase of 16% out of 350% for the transition from conditions 2 to 4). The resulting f_α(t) are shown in Fig. 6b. To extract a characteristic relaxation time for each transition, we interpolate the normalized the change [f_α(t) − f_α(∞)]/[f_α(0) − f_α(∞)] to obtain the time at which it equals 1/e.

Fits of BCF theory to experimental results

The general expressions for K^ss(f_α) and K^dyn(f_α)³⁴ involve nine unknown quantities (D, $$ {\rho }_{\mathrm{eq}}^0{\ell }^3$$, $$ f_{\alpha }^0$$, and the six $$ {\kappa }_x^j$$) to be determined or constrained by the measurements. This is a challenge because there are only six measured quantities (four steady state α terrace fractions $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$ at different growth rates G, and two relaxation times for transitions in G.) The best fit allowing all nine parameters to vary was first determined by a numerical Levenberg–Marquardt nonlinear regression search algorithm. This minimized the goodness-of-fit parameter $$ {\chi }^2\equiv \sum {\left[\left(y_i-y_i^{\mathrm{calc}}\right)/{\sigma }_i\right]}^2$$, where the y_i and σ_i are the six measured quantities and their uncertainties. To estimate the uncertainties in the $$ f_{\mathrm{\alpha }}^{\mathrm{ss}}$$, we multiplied those obtained in the fits to the 3H(T1) reconstruction by a factor of 4, to account for the uncertainties in the atomic coordinates used. We estimated the uncertainty in the t_rel to be 10%.

This initial fit found that the minimum χ² occurs in the region of parameter space in which the $$ {\kappa }_{+}^A$$ coefficient is large, while the $$ {\kappa }_{-}^A$$, $$ {\kappa }_{-}^B$$, and $$ {\kappa }_0^A$$ coefficients approach zero. In this case, the limiting expressions (Eqs. (11) and (12)) for K^ss(f_α) and K^dyn(f_α) involve only the four unknown quantities $$ D/{\kappa }_{+}^B$$, $$ D/{\kappa }_0^B$$, $$ D{\rho }_{\mathrm{eq}}^0{\ell }^3$$, and $$ f_{\alpha }^0$$. A final fit using Eqs. (11) and (12)) allowing only these four parameters to vary gave the same solution as the nine parameter fit.