Resolution of the exponent puzzle for the Anderson transition in doped semiconductors

The Anderson metal-insulator transition (MIT) is central to our understanding of the quantum mechanical nature of disordered materials. Despite extensive efforts by theory and experiment, there is still no agreement on the value of the critical exponent $\nu$ describing the universality of the transition --- the so-called"exponent puzzle". In this work, going beyond the standard Anderson model, we employ ab initio methods to study the MIT in a realistic model of a doped semiconductor. We use linear-scaling DFT to simulate prototypes of sulfur-doped silicon (Si:S). From these we build larger tight-binding models close to the critical concentration of the MIT. When the dopant concentration is increased, an impurity band forms and eventually delocalizes. We characterize the MIT via multifractal finite-size scaling, obtaining the phase diagram and estimates of $\nu$. Our results suggest an explanation of the long-standing exponent puzzle, which we link to the hybridization of conduction and impurity bands.

The Anderson metal-insulator transition (MIT) is the paradigmatic quantum phase transition, resulting from spatial localization of the electronic wave function due to increasing disorder [1]. As for any such transition, universal critical exponents capture its underlying fundamental symmetries. This universality allows to disregard microscopic detail and the Anderson MIT is expected to share a single set of exponents. The last decade has witnessed many ground-breaking experiments designed to observe Anderson localization directly: with light [2][3][4][5][6][7][8][9][10][11], photonic crystals [9,12], ultrasound [13,14], matter waves [15], Bose-Einstein condensates [16] and ultracold matter [17,18]. The mobility edge [19], separating extended from localized states, was only measured for the first time in 2015 [20]. The hallmark of these experiments is the tunability of the experimental parameters and the ability to study systems where many-body interactions are absent or can be neglected. Under such controlled conditions, the observed exponential wave function decay, the existence of mobility edges and the critical properties of the transition [21,22] are in excellent agreement with the non-interacting Anderson model [1]. Furthermore, scaling at the transition [23] leads to highprecision estimates of the universal critical exponent ν from transport simulations (ν = 1.57 (1.55, 1.59) [24]) and wave function statistics (ν = 1.590(1.579, 1.602) [25]).
Anderson's original challenge was to describe localization in doped semiconductors. For these ubiquitous materials, the existence of the MIT was confirmed indirectly by measuring the scaling of the conductance σ ∼ (n − n c ) ν when increasing the dopant concentration n beyond its critical value n c . However, a puzzling discrepancy remains: a careful analysis by Itoh et al. [26] highlights that the value of ν can change significantly with the control of dopant concentration around the transition point, the homogeneity of the doping, and the purity of the sample. Following Stupp et al. [27], they suggest that the intrinsic behaviour of an uncompensated semiconductor gives ν ≈ 0.5 [28], while any degree of compensation results in ν ≈ 1 [29]. Evidently, these values disagree with the aforementioned theoretical and experimental studies. The inability to characterize the Anderson transition in terms of a single, universal value for ν is known as the "exponent puzzle" [27,30].
Most theoretical models that have been applied to this problem lack the ability to capture the full complexity of a semiconductor. The Anderson model, for example, ignores the detail of the crystal lattice and the electronic structure, and also simplifies the physics by ignoring many-body interactions and interactions between dopant and host material. These factors are known to change the universal behavior [31,32] and the value of ν, as shown in studies on correlated disorder [33][34][35] and hydrogenic impurities in an effective medium, where ν ≈ 1.3 [36,37]. Here we propose a fundamental shift from studying localization using highly-simplified tightbinding Anderson models, to atomistically correct ab initio simulations [38,39] of a doped semiconductor. We illustrate the power of our approach for sulfur-doped silicon, Si:S, where the MIT occurs for concentrations in the range 1.8-4.3 × 10 20 cm −3 [40]. We model the donor distribution in Si:S by randomly placing the impurities in the lattice [41]. While we concentrate on Si:S here, our method is straightforwardly applicable to other types of impurities (Si:P; Si:As; Ge:Sb), hole doping (Si:B) and co-doping (Si:P,B; Ge:Ga,As).
With this approach we observe the formation of the impurity band (IB), upon increasing n, and its eventual merger with the conduction band (CB). States in the IB become delocalized, as measured directly via multifractal statistics of wave functions [42], and we observe and characterize the MIT. In Fig. 1 we plot how n c and ν vary for energies ε in the IB below the Fermi energy ε F . For ε ∼ ε F , the values are ν ∼ 0.5, while deeper in the IB the exponents increase to about ν ∼ 1, reaching values around 1.5. As we will show below, our simulations of an uncompensated semiconductor suggest that the reduction in ν at ε F is due to the hybridization of IB and CB. Deep in the IB the physics of the Anderson transition reemerges with ν reaching the range of its proposed universal value [24,37,43]. Experiments can readily access these higher values by moving ε F via compensation [26] -intentional or otherwise. ization [38] and discovery [39]. With the choice of Si:S, we can observe the transition in systems of up to 11 × 11 × 11 unit cells, i.e., 10648 atoms. These large system sizes can in principle be reached by linear-scaling DFT [44], but despite this, the necessity to average over many hundreds of disorder realizations, makes repeated DFT calculations impractical for our purposes [45]. We therefore devise a hybrid approach: linearscaling DFT calculations are performed, using the ONETEP code [46], on prototype systems of 8 × 8 × 8 diamond-cubic unit cells (4096 atoms), employing geometry optimization to allow for the lattice to accommodate single or multiple S impurities. We include nine in-situ-optimised non-orthogonal local orbitals for each site (in atomic Si, atomic orbitals are occupied up to level 3p; for better convergence we additionally consider the five 3d orbitals). When embedded in silicon, sulfur, like the other chalcogens, acts as a deep donor. Such defects have highly localized potentials that are well-described in a local orbital basis [47]. The resulting Hamiltonians and overlap matrices, represented in terms of the nonorthogonal local orbital basis, are used to construct three catalogs of local Hamiltonian blocks (cf. Fig. 2). For each system size L, concentration of impurities n and disorder realization, we build the effective tight-binding Hamiltonians H and overlap matrices O from these catalogs (cf. Fig. 2) and solve the large generalized eigenvalue problem [48,49] Hψ j = ε j Oψ j , j = 1, . . . , 9L 3 (1) for eigenenergies ε j and normalized eigenvectors ψ j . In binding matrices of size 95832 × 95832). We average up to 1000 different disorder realizations for each L and n (cf. Table  I).
Characterizing the IB and its DOS is interesting for its spin and charge transport properties [50,51]. We compute the density of states (DOS) of the IB from the ε j 's while changing the number of impurities N S . We define ε F as the midpoint between the highest occupied IB state at energy ε IB and the lowest unoccupied CB state at ε CB . To obtain the average DOS for given N S and L, we shift the spectrum of each realization such that ε F = 0. The DOS shown in Fig. 3 is calculated by summing over Gaussian distributions of standard deviation σ = 0.05 mHa = 1.36 meV centered on ε j − ε F . We find that the IB has a peak at ε − ε F ∼ −0.1 eV and a tail extending towards the VB with increasing n. This agrees with known features of the IB in doped semiconductors [52]. We emphasize that Si:S is particularly interesting for intermediate-band photovoltaic devices, where the efficiency increases when deep IB states can capture low-energy photons [50]. In order to avoid electron-photon recombination, the IB states should be delocalized such that they can contribute to the photocurrent. The determination of n c and the pronounced tail of the IB as presented in Fig. 3 therefore provide essential information for future device applications.

CHARACTERIZING THE METAL-INSULATOR TRANSITION
In the last decade, multifractal analysis [42,53,54] has become the method of choice to reliably and accurately extract the localization properties from wave functions [14,25,55]. In its essence, it describes the scaling of various moments of the spatial distribution of |ψ j | 2 , which is encoded in the singularity strengths α q . At criticality, the universality class of the transition determines the scaling of α q with n and L. We capture this behavior using the well-established framework of finite-size scaling. Following [25], we assume that the data for each L meet at the critical point w = 0 with a value α crit q , and scale polynomially with ρL 1/ν , where ρ(w) = w + m ρ m=2 b m w m includes higher-order dependencies on the dimensionless concentration w = (n − n c )/n c . We hence fit the data against the function [56] α q (n, L) with n c , ν, α crit q , the a i 's, and the b i 's as fitted parameters, and m L and m ρ as expansion orders [57]. We illustrate the localization and scaling properties of the wave functions using the moments α 0 and α 1 . Figure 1 shows the results of the fits as ε is varied, obtained from (2) (see tables 1 and 2 [57]). They include systems with up to L 3 = 10648 atoms, varying n and number of realizations for each L as given in Table I. Crucially, we only accept estimates of n c and ν after consistently and rigorously checking their robustness against perturbations in n and stability when increasing m L and m ρ [24,25,57].
Following this recipe, we identify the Anderson MIT and reconstruct the energy dependence of the mobility edge n c (ε) in Si:S. It exhibits (i) a maximum close to ε F and a decrease until ε − ε F ≈ −0.09 eV. (ii) For lower energies, n c increases again and the mobility edge moves towards the tail of the IB (cf. Fig. 3). These findings suggest a natural split into two different regimes, as also seen in the energy dependence of ν. Values of ν in regime (i) increase continuously from ν ≈ 0.5 at ε F to about ν ∼ 1. In regime (ii), we find a larger spread of values 1 ν 1.5. This spread is consistent with the statistical uncertainty of each estimated ν in Fig. 1, which is dominated by the range of L and the ensemble size N (cf. Tab. I). However, the trend in ν observed in regime (i) requires a different explanation.

HYBRIDIZATION OF IMPURITY AND CONDUCTION BAND STATES
In Fig. 4, we present the distribution of states resolved in both energy ε and α 0 . Perfectly extended states correspond to α 0 = 3, while increasing localization results in α 0 → ∞. The data for N S = 40 (n = 4.9 × 10 20 cm −3 ) show metallic states of the CB with α 0 ≈ 3 at ε ≈ ε F . The IB is characterized by (i) a majority region of states with α 0 ≈ 3. For N S = 40 we show the density plot of the distribution (from blue for low to red for high density, see color scale) and the contour lines enclosing 68% (white) and 95% (black) of the α 0 's. For N S = 100 we indicate the same contours (red, dashed). As in Fig. 3, the shading denotes the delocalized region (in the L → ∞ limit) according to Fig. 1.
(the metallic limit) close to ε F . This observation is intriguing when tensioned against the simultaneous decrease in the value of ν at ε ∼ ε F (cf. Fig. 1). Apparently the localization of the IB states is substantially modified by the presence of the states from the CB. In Fig. 5, we show the α 0 data for N = 4096 as a function of ε and n. For small impurity concentrations, the IB consists of localized states with some of the largest values of α 0 ∼ 3.6, while the CB contains delocalized states with α 0 3. Upon increasing n, the IB develops and its states become more delocalized. Initially, this trend is most pronounced where the DOS of the IB is large (see Fig. 3), i.e. around ε − ε F ∼ −0.12 eV. Simultaneously, states at the top of the IB exhibit α 0 values close to those denoting extended states in the CB, even before the band gap has fully closed. When reliable scaling is possible, we eventually see how the two mobility edges emerge. At the lower mobility edge, we find values of ν ∼ 1 − 1.5. At the upper mobility edge, we observe lower estimates for ν coinciding with lower α 0 values at the transition due to the strong hybridization of IB and CB. Let us discuss how this observed hybridization and the resulting enhanced metallic behavior can affect the value of ν. The leading scaling behavior from (2) is α crit 0 − α 0 ∼ wL 1/ν for w > 0. A decrease in the effective α 0 yields an increase in α crit 0 − α 0 , which is consistent with a reduced exponent ν as observed in Fig. 1 for (ε−ε F ) −0.1 eV. An argument similar to the famous "Gang of Four" result [23] can be made directly for the transport experiments, where an increase in σ ∼ w ν for 0 ≤ w 1, i.e. close to the critical point, is also consistent with a reduced ν.  [26] that in experiments a change from 0.5 to ν ∼ 1 can be induced by compensation. Taken together, compensation and band hybridization provide two important pieces to complete the "exponent puzzle": modelling the Anderson transition in doped semiconductors needs to include the CB (VB for hole-doped materials) together with the IB provided by the Anderson modelthe experiments obviously include both and hence find ν values which, depending on their state of compensation, are occasionally different from the predictions based solely on the Anderson model of the IB.
How exactly the hybridization changes the value of ν, as well as whether the value of ν deep in the IB is different from the non-interacting predictions, remain challenges for future high-precision studies. Still, the approach we present here exploits and transfers the accuracy and versatility of modern ab initio simulations to the study of Anderson localization in doped semiconductors-at a fraction of the computational cost. Beyond bulk semiconductors, other disordered systems [59], 2D [60][61][62] and layered materials [63] are also well within reach, as is the investigation of the influence of many-body physics by, e.g., studying the interaction-enabled MIT in 2D Si:P [64,65]. We find that the critical concentration agrees quantitatively with a previous experiment in Si:S by Winkler et al. [40]. Our approach is hence capable of modeling fundamental physical phenomena while also mak-ing material-specific predictions.

METHODS
Our simulations use the ONETEP linear-scaling DFT package [46] with the PBE exchange-correlation functional [66]. We include nine orbitals of radius 10 Bohr radii on each site, and a psinc grid with an 800 eV plane-wave cutoff [67]. This gives an accuracy equivalent to plane-wave DFT for Si and other materials [68]. The first catalog of local Hamiltonian blocks describes the Si host material, i.e., a set of onsite energies and hopping terms, starting at a central Si atom and extending to 10 shells of Si neighbors. The second corresponds to the energies and hopping terms when the central atom is S, and the third catalog to pairs of neighboring S atoms. Here, we define a "neighbor" as being at most 4 shells apart. If two S atoms are 5 or more shells apart, each S atom is unaffected by the presence of the other [57]. The impurity distribution is generated by randomly substituting the impurity atoms onto lattice sites. This follows the experimental techniques used to achieve high S concentrations, combining ion implantation with nanosecond pulsed-laser melting and rapid resolidification [40]. The impurities are effectively trapped in the substitutional sites [41]. With φ α denoting the non-orthogonal local orbitals, we write the eigenvectors ψ j = α M α j φ α of Eq. (1) in a "site" basis by summing over the nine orbital coefficients of each site k, i.e. |Ψ j k | 2 = α∈k,β M α j O αβ M β j . We coarse grain Ψ by fixing a box size l < L and partitioning the domain in (L/l) 3 = λ −3 boxes. The amplitudes of the coarsegrained wave function µ are given by µ s = k∈B s |Ψ k | 2 , i.e. by summing all |Ψ k | 2 pertaining to the same box B s . After rescaling the amplitudes as log µ s / log λ, we compute their arithmetic mean α 0 = log µ s s / log λ and weighted mean α 1 = s µ s log µ s / log λ (proportional to the von Neumann entropy [69]). Finally, for each n and L we take the ensemble average α q (n, L), where q = 0 or 1. Further details on the DFT simulations and the scaling analysis can be found in [57].