Comment on “Diffusion of n-type dopants in germanium” [Appl. Phys. Rev. 1, 011301 (2014)]

A recent Applied Physics Review has discussed self-diffusion and B diffusion during irradiation at high temperature.¹ Unfortunately, erroneous assumptions in B diffusion data analysis led the authors mistakenly to critique recent work by us that identified two forms of self-interstitial in Ge.² Here, we show that their work when correctly interpreted confirms our conclusions. In the following discussion, C_X and D_X represent the concentration and diffusivity of species X, E_f^X and E_m^X (S_f^X and S_m^X) its formation and migration energies (entropies), respectively, and D_SD^X= D_XC_X^eq/C₀, where C₀ is the lattice density, is the contribution of X to equilibrium self-diffusion with activation energy E_SD^X and activation entropy S_SD^X. The species described are the vacancy, V, self-interstitial, I, and B-interstitial pair, BI, and we consider two distinct forms of I; I and I. The first is the compact I, well known from the literature, the second is the i-morph—an extended self-interstitial with properties of a small amorphous pocket.² This entity, in some ways reminiscent of the high-temperature “liquid drop” proposed by Seeger,^3,4 however, needs to be understood from a fundamentally different perspective; the key feature is an amorphous-like property with corresponding energetic and entropic behavior.² 

Our evidence has been disputed by Bracht and coworkers.^1,5 Their objection appears to be based on (a) a misunderstanding of our analysis of long-range BI migration in Ref. 2, (b) an erroneous analysis of BI mediated B diffusion in Ref. 1 and a precursor paper.⁶ To clarify the issues, we first briefly review the disputed analysis of Ref. 2. B diffuses in Ge, as in Si,⁷ via a fast migrating BI pair formed by the reaction B_s + I ↔ BI.⁸ A parallel reaction, B_s ↔ BI + V, also occurs but has no significant influence under the conditions of Refs. 1 and 2. BI in Ge, as in Si, has a large migration length, λ—a quantity closely connected to the difference in Gibbs free energy between BI and I.² This leads to exponential diffusion tails after anneal times short enough that only a fraction of C_Bs experiences a reaction with I to form BI. This behavior occurs under both equilibrium and irradiation conditions, with λ independent of I supersaturation.^7,8

To analyse this behavior, the diffusion of I, BI, and V can be modeled by numerical solution of the coupled equation system as in Refs. 1 and 6 and elsewhere. Under certain conditions, the full system can be reduced to one equation which has an analytical solution involving just g, λ, position, and time—the g-λ solution.^7,9 This is a mathematical approximation to the full equation system that describes the detailed properties of dopant diffusion and is applicable under equilibrium and non-equilibrium conditions. The necessary and sufficient conditions for accuracy of this approximation are that C_BI ≪ C_Bs, the Fermi level at the diffusion temperature, T, is slowly varying in the local region of interest, and there are no significant gradients in C_I or C_V. The latter condition prevails if D_BC_B ≪ D_IC_I^eq, and no significant gradients are generated by external processes. It is not necessary to assume point-defect equilibrium. When the preceding inequality is relaxed towards D_BC_B ∼ D_IC_I^eq, the rate of the reaction I + B_s → BI is slightly modified by “chemical pump” effects but the g-λ solution still accurately describes λ.¹⁰ 

In the experiments of Ref. 2, all these conditions were satisfied, so the g-λ approach could be used to extract accurate λ values from our experimental secondary-ion mass spectrometry (SIMS) profiles, thus avoiding the costly use of a general diffusion solver as kernel in least-squares minimization. Following this analysis, we deduced a T-dependent free-energy difference between BI and I, indicating that the latter has two distinct forms. The first, dominating I-mediated transport at low T, is the simple I. The second, dominating D_SD^I at high T, is an extended I; the i-morph.² A transitional region, where both defect forms contribute, exists over a ∼100 °C range centered on 475 °C. Reference 1 cites this analysis, claiming in error that the g/λ approach is inapplicable under non-equilibrium conditions.

We now consider the analysis of B diffusion in Refs. 1 and 6. Fig. 17 of Ref. 1 presents B profiles in Ge measured by SIMS after proton irradiation of a B-doped Ge superlattice at 550 and 630 °C. To show clearly the detailed B profile shape evolution during diffusion, we have selected and plotted the data for a single marker layer in Fig. 1. The profile shows characteristic exponential-like tails (showing up as almost straight lines on the logarithmic plot of Fig. 1) on each side of the B-doped marker layers. The curves turn up at the edges of the plot owing to overlap of diffusion from neighboring markers. The data at 550 °C show significantly more diffusion than at 630 °C, because at lower T both λ ( = ${DBI/(k1+C0)}1/2$ in Bracht's notation) and the forward reaction rate g ( = $k1−CI$⁠) are larger. The larger λ reflects the increased number of BI diffusion jumps per migration event as the thermal energy available for dissociation, BI → B + I, is reduced. The larger g reflects the increased number of lattice sites each beam-generated I visits before recombining with V. The static peak represents those B_s which have not yet undergone reaction (1)—a statistical effect due to the finite arrival rate of I at B_s.

The approach taken in Ref. 1, following Ref. 6, fails to recognize and model these key effects. This seems to be caused by unrealistic assumptions (a) on B clustering during annealing of initially substitutional B, (b) that C_BI ≫ C_Bs. In relation to point (a), Ref. 6 assumed a priori that, at all considered anneal times, t, each B marker had a large clustered component, adjusted for each T/t combination to keep C_Bs ≤ 5 × 10¹⁸cm⁻². This ignores the transient dynamics of B clustering in the MBE-grown doping structure as BI migrates and traps on other B atoms. In the simulations of Ref. 6 (Fig. 17), clustering, unrealistically, actually decreases with time. Point (b) is a result of assuming S_f^BI≈ 30 k—a problematic choice as the entropy of D_B, S_f^BI+ S_m^BI≈ 20 k¹¹ and negative S_m^BI is highly unlikely. It is unclear why such a large S_f^BI has been used, unless it is to prevent D_B from varying as (p/n_i)² as the model assumes BI is in a singly positive charge state.⁶ The result of these several choices is that all the simulated profiles have Gaussian shapes at low B concentration (blue curves in Fig. 1). This is a poor fit to the data, which show a clear exponential-like trend, thus directly demonstrating that C_BI ≪ C_Bs, refuting assumption (b) above and rendering equation (20) and Fig. 6 of Ref. 6 invalid. This key point is further underscored by the fact that proton irradiation experiments with almost identical Frenkel-pair production rates¹² to those used in Ref. 1 explicitly show $g∝ϕ$⁠, where $ϕ$ is the beam flux⁸ (a test not reported in Refs. 1 and 6). The failure of the assumption C_BI ≫ C_Bs is most graphically evident in the lower panel in Fig. 1, where we present data from an earlier study⁸ using very similar processing conditions. The data show essentially the same exponential tails as in Ref. 1, although in this case clustering is entirely absent, all B is available to diffuse, and the static peak represents those B atoms which have escaped interaction with I during the short annealing time. The imposition of C_BI ≫ C_Bs, however, identifies essentially all unclustered B as continuously diffusing BI, leading to a Gaussian diffusion profile.

Fully coupled models as in Refs. 1 and 6 easily reproduce observed exponential tails if model parameters are correct. A first step towards this goal is to eliminate the unrealistic saturation of C_BI. This can be done by reducing S_f^BI from 30 k to below 20 k. This then allows extraction of other key parameters, inaccessible with the assumption C_BI ≫ C_Bs, such as the charge states of BI and I (from data on the Fermi-level dependence of B diffusion), and E_f^I, S_f^I (from exponential tails, since E_f^I, S_f^I determine the parameter k₁⁺ in $λ={DBI/(k1+C0)}1/2$ (Refs. 2 and 6)).

In Refs. 1 and 6, the peak C_B is ∼10× higher than in earlier experiments.² In this situation, clustering, chemical-pump, and Fermi-level effects may all influence diffusion, so the data in Ref. 1 are a more complex resource for parameter determination than those in Ref. 2. Nevertheless, to illustrate the robustness of λ extraction with the g-λ approximation, we apply it, outside its strictly applicable range, to the “short-time” (1 h) B data of Refs. 1 and 6. The results, using the same λ values as in Ref. 2, are shown in Fig. 1 (red curves).¹³ The fits are essentially perfect—far better than those obtained in Refs. 1 and 6. Moreover, unlike the fits in Ref. 1 they respect the data from Ref. 8, which explicitly showed C_BI ≪ C_Bs. It should now be clear that our analysis in Ref. 2, where conditions were optimal for λ extraction, is extremely robust. Moreover, the B data discussed in Ref. 1 support the analysis in Ref. 2, not refuting it as claimed.¹ The B model parameters used in Ref. 1, however, are far from correct. We now restate the established position:² our data (further supported by high-T data in Refs. 1 and 6) show that I-mediated self-diffusion in Ge involves a simple I at low T (significantly below 475 °C), an i-morph, I, with S_SD^I≈ 30 k, E_SD^I≈ 6.1 eV at high T (significantly above 475 °C), and a transitional region around 475 °C where both are significant. Finally, it should be noted that the open triangle^14,15 in Fig. 2 shows that B also diffuses via BI under equilibrium conditions. An alternative model based on vacancy exchange¹ would imply jump lengths of only 0.25 nm.

Having dealt with B diffusion analysis in some detail in this comment, we would like to emphasize that Section V of Ref. 1 also references an elegant analysis of experiments by the Munster group and coworkers on the diffusivity of I in irradiation experiments on Ge isotope superlattices.⁶ That analysis revealed E_m^I= (1.84 ± 0.26) eV. This value far exceeds estimates of 0.5–0.6 eV obtained from perturbed angular correlation measurements at low T,^16,17 0.6 eV obtained for simple I configurations from density functional theory using accurate LDA + U functionals,¹⁸ and <1 eV indicated by jump rates exceeding ∼1 s⁻¹ at RT for I directly observed in aberration-corrected TEM.¹⁹ Thus in retrospect one can see that the 1.84 eV value rules out the simple I assumed in Ref. 1 and strongly favors the i-morph mechanism we proposed in Ref. 2. This has vast implications for defect physics which remain to be explored. Finally, taken together with our observed E_SD^I= 6.1 eV at high T,² E_m^I= 1.84 eV implies E_f^I≈ 4.3 eV, in the range of recent atomistic calculations in course of publication.²⁰ In conclusion, discussion prompted by conflicting analyses of experiments in Refs. 1, 2, and 6 has significantly progressed understanding of the complex behavior of self-interstitials in Ge.