The structure of heavily doped impurity band in crystalline host

The effect of disorder has been extensively studied since Anderson’s seminal paperAnderson (1958). Diluted magnetic semiconductors (DMS) doped with a small concentration of charged impurities constitute an interesting magnetic system that has a number of novel features for study by numerical simulationAvérous and Balkanski (1991). Much of the research has been focused on II-VI (such as CdTe or ZnSe) and III-V (such as GaAs) compound semiconductors doped with a low concentration ($x\sim 1-8\%$) of Manganese (Mn) impurities. Of particular interest in this field is Ga_1-xMn_xAs which has been shown to exhibit ferromagnetic behavior above 100KOhno (1998). In these samples, the Manganese is substitutions with the Gallium and acts as an acceptor (donating one hole to the crystal), so that the material is p-type. The holes bind to the impurities with an energy of around 130 meV around $x\sim 10\%$Beschoten et al. (1999). Since $x\ll 1$, the overlap between different impurity states can be ignored, thus the interaction between the charge carriers can be neglected. The system can be simply described by a noninteracting tight-binding model. When the system contains only one impurity, and the binding energy is large enough, an impurity state appears below the conductance band (we assume the impurity potential is attractive). It is locally distributed in space near the impurity potential within a localization length $\zeta$. As increasing the concentration $x$, the overlap between different impurity states extends the single impurity energy to an impurity band in the density of state (DOS) and eventually merges with the conductance band. Simultaneously, the states in the impurity band are expected to become more and more extended and ultimately regain their bandlike character Cook and Berciu (2012). However, the details inside the impurity band are rarely studied.

One reason for lacking such a study is the computation difficulty even in the non-interacting case. Generally, the percentage of the state in the impurity band in the total number of states is about $10\%$ at the concentration we are interested in. Taking a 3-dimensional Anderson model with lattice size $30\times 30\times 30$ as an example, the number of states which we need to know in the impurity band is about 3000. The exact diagonalizationWeiße and Fehske (2008) for such a system is very difficult due to the large dimension. On the other hand, we have to do a large number of sample averages. The sparse matrix diagonalization, such as the Lanczos methodOJALVO and NEWMAN (1970), can be adapted to obtain a few lowest-lying states or a few states nearby special energy (the simplest way is diagonalizing $(H-\epsilon I)^{2}$ by using the original Lanczos diagonalization method).

Machine learning methods have recently emerged as a valuable tool to study the quantum many-body physics problemsCarleo and Troyer (2017); Carrasquilla and Melko (2017); Ch’Ng et al. (2017); Van Nieuwenburg et al. (2017); Venderley et al. (2018); Wetzel (2017); Rodriguez-Nieva and Scheurer (2019); Lidiak and Gong (2020); Hsu et al. (2018); Hendry et al. (2021); Choo et al. (2019); Pfau et al. (2020); Sharir et al. (2020); Hendry et al. (2022); Chen et al. (2022). Its ability to process high dimensional data and recognize complex patterns have been utilized to determine phase diagrams and phase transitionsWang (2016); Ohtsuki and Ohtsuki (2017); Tanaka and Tomiya (2017); Mano and Ohtsuki (2017); Broecker et al. (2017); Schindler et al. (2017); Li et al. (2019); Dong et al. (2019); Kotthoff et al. (2021); Zhang et al. (2019a, b); Käming et al. (2021). In particular, Convolutional Neural Network(CNN)Krizhevsky et al. (2012) model, which initially is designed for image recognition, was widely used to study different kinds of phase transition problems including the Bose-Hubbard modelBohrdt et al. (2021), spin 1/2 Heisenberg modelThéveniaut and Alet (2019), quantum transverse-field Ising modelZhang et al. (2019a) and etc. The power of using machine learning to recognize quantum states lies in their ability to finish tasks without the knowledge of physics background or the Hamiltonian of the system. Even if the neural network is trained in a small energy region of the system, it can be used to obtain the whole phase diagramOhtsuki and Ohtsuki (2017); Mano and Ohtsuki (2017). Also, it can discriminate quantum states with high accuracy even if they are trained from a totally different Hamiltonian. This special feature of machine learning inspires us to try to identify the delocalized states in the “impurity band”.

In this paper, we develop a method to obtain the correct density of states (DOS) and other localization properties, such as inverse participation ratio (IPR)Brndiar and Markoš (2006) and thouless numberEdwards and Thouless (1972), by using Jacobi-Davidson sparse matrix diagonalizationBollhöfer and Notay (2007) with an importance sampling statistics method. Meanwhile, we train a 3-dimensional CNN model using the data generated from the Anderson model, and then the trained model is used to identify the existence of extended states in the impurity band. This manuscript is organized as follows: In sec.II we describe the tight-binding model on the cubic lattice and numerical methods; Sec. III demonstrates the effect of heavy doping studied by studying the IPR and Thouless number; Sec.IV demonstrates the implementation of the deep learning approach and the results from the trained neural network model; finally, we close with a conclusion.

We consider a tight-binding model on a D-dimensional hypercubic lattice with the nearest neighbor hopping t, and on-site energies $\epsilon_{i}$:

The hopping term simulates the iterative electrons and the on-site energy has a bimodal distribution $\epsilon_{i}=-W$ with probability $x$, and $\epsilon_{i}=0$ with probability ($1-x$). This model a host lattice with a single relevant band, with a fraction $x$ of substitutional impurities. For one-dimensional (d = 1) free electrons, the energy-momentum dispersion relation is $E(k)=2t\cos(k)$, it is easy to get the DOS with the formula

The result for 1D is:

There is no analytic solution for higher dimensional systems, however, an approximation that is accurate to roughly $2\%$ was given by Andres et al Andres et al. (1981). Instead, the DOS can be calculated numerically by exact diagonalization as shown in Fig. 1 where $t$ has been set to the unit. After introducing the impurities, all states become localized in 1D and 2D based on the scaling theory of localization Abrahams et al. (1979). Part of the states become to be localized and develop into an impurity band at the edge of the conducting band. To determine the localized/extended state, namely the location of the mobility gap, we calculate the inverse participation ratio (IPR)Brndiar and Markoš (2006)

for each state, where the $\psi_{i}$ is the weight of an eigen wave function on the $i$’th site. Heuristically, if we compare two trivial states with wave functions for a $N$-site system:

and

where $\Psi_{extended}$ is an extended state which has equal weight on each site and $\Psi_{localized}(j)$ is a localized state which only has weight on the $j$’th site. It is easy to see that the IPR of $\Psi_{extended}$ decreased with the order of $\frac{1}{N}$ and a constant for $\Psi_{localized}(j)$. On the other hand, the Thouless numberEdwards and Thouless (1972) is defined as:

where $\delta E$ is the energy difference while the boundary condition changes from periodic boundary condition (PBC) to anti-periodic boundary condition (APBC) and the $|\Delta E|$ is the average energy distance around $E$. Since only the extended states are sensitive to the change of boundary condition, $g(E)$ grows linearly as a function of the system size for the extended state, and conversely, it reduces for the localized state. In this work, we determine the localization properties by systematically studying the IPR and Thouless number for different system sizes, and the crossover points of the Thouless number give us a hint of the mobility edge.

For three dimensional cubic lattice of size $L$, the Hamiltonian matrix has a dimension of $L^{3}$. General full exact diagonalization methods, such as Lapack library Anderson et al. (1999), can only deal with small system sizes. The computation time of diagonalizing one matrix with size $L^{3}$ grows dramatically as a function of the system size. As shown in Fig. 2, we deal with a system with size $19\times 20\times 21$ with doping concentration $x=5\%$, after averaging thousands of samples, we obtained the DOS for different doping energies. It is shown that a peak emerges gradually near the band edge as increasing the doping energy $W$. This peak becomes more prominent around $W\sim 4.5$ at which an obvious depletion is developed at the junction between the impurity band and the conduction band.

Since only the developed impurity band is the interesting part we are focusing on. The number of states in the impurity band is about the lowest $10\%$ of states in the whole band, thus we do not have to fully diagonalize the Hamiltonian. On the other side, we just need to calculate the DOS, IPR, and Thouless number for these lowest $10\%$ states after averaging thousands of samples. According to our demand, we use the sparse matrix diagonalization with Jacobi-Division (JADA) method Bollhöfer and Notay (2007) which can search a few (10 to 20) states efficiently near specific points. For a given sample at fixed doping strength, we randomly distribute the reference points (30-50 points) in the impurity band, taking $W=-4.5$ as an example, the reference points are picked randomly in the region $[-8:-4]$, about 10-20 states can be obtained by JADA around each reference point. The reference points could also be picked by importance sampling based on the DOS for a small system from the full diagonalization. We collect all these energies for each reference point in one sample. After thousands of sample averaging, we obtain the same DOS as that from the full diagonalization for a small system. It is obvious that the JADA method can easily go beyond the limit of the full exact diagonalization. At least on the same price of the computation time, we can nearly double the system size compared to the Lapack method. In this work, we calculate the properties for system sizes up to $40^{3}$ sites by using the JADA method.

As analyzed in the previous section, with typical doping concentration $x=5\%$, we find that a clear impurity band in the DOS is developed at about $W=-4.5$. We plot the DOS and IPR together for different system sizes as shown in Fig.3. The line of the DOS for different system sizes collapses to a single curve and it is the same as that from ED as shown in Fig. 2, which tells us that we have already obtained the essential information of the impurity band. As increasing the system size, the IPR does not change on the edge of the band which means the states on the edge of the whole band are localized. The IPR in the bulk decreases as enlarging the system, especially at the center of the impurity band ($E\sim-6.7$), the IPR drops to zero, which is the same as in the system bulk ($E\sim-4.0$). However, there is a small peak near $E\sim-5.5$ which is at the right edge of the impurity band. The IPR in the vicinity of this point tends to saturate to a fixed value as increasing the system size. The nonzero saturation of the IPR at this energy means another possible mobility edge existing near the junction between the conduction band and the impurity band.

In order to justify our conjecture, we systematically study the value of IPR for several system sizes. As shown in Fig. 4(a), we choose four points from the knowledge of the DOS and IPR. (1) $E=-4.2$ is in the bulk of the conduction band, at which the state is extended. (2) $E=-5.4$ is at the right edge of the impurity band. The state here is localized according to our conjecture. (3) $E=-6.2$ is in the bulk of the impurity band, which is extended according to its zero IPR value in large $L$ limit. (4) $E=-6.8$ is on the left edge of the impurity band and thus at the edge of the whole energy band. The state at the band edge is supposed to be localized. In Fig. 4(b) we again compare the DOS from JADA with that from Lapack which shows a convergence in large system size. According to the way of choosing these four points, (1) and (3) should have similar behavior as increasing the system size, and vice versa for (2) and (4). Fig. 4(c) shows the IPR for these four energies in different system sizes. We plot the data in log scale and fit it by function

The sign of the curvature $C$ tells us whether the state is localized or not. For (1) and (3), $C<0$ means they are extended states, and oppositely $C>0$ for localized states at points (2) and (4).

As another criterion, we calculate the Thouless number $g(E)$ for different system sizes. The results are shown in Fig.5 in which we plot the DOS together with the same horizontal axis. The impurity band has been divided into several regions at the crossover of $g(E)$ for different sizes. We label these regions by “L” (localized) and “E” (extended) to demonstrate different behavior $g(E)$. As increasing the system size, it is obvious that the $g(E)$ increases in the “E” region and decreases in the “L” region. The energies with vertical lines are the locations of the mobility edges, or the boundaries between the localized states and extended states.

Convolutional neural network(CNN), which is originally designed for 2D image recognition, has been widely adopted in studying phase transition and achieves high accuracy in recognition. A standard image recognition model can be used for a 3D electron system by integrating the 3D electron density in one direction. But the drawback of this approach is that the information of the electron density along one direction is lost during integration. So, we design a 3D CNN model for our 3D lattice model. To distinguish the localized and delocalized state, the CNN model will return two real numbers to represent the probability of the extended state $P$ and localized state ($1-P$) for the given wave function. If the probability of the extended state is larger than 0.5, we think the eigenstate is delocalized, and vice localized. Due to the limitation of the graphics memory (8GB) of our graphics card (NVIDIA GTX 1080), we consider a 3D $20\times 20\times 20$ lattice. The hidden layers in the CNN model consist of convolutional layers, max-pooling layers, and fully connected layers. The loss function is defined by the cross entropy $H(x)=-\sum_{x}p(x)\log q(x)$. During the training, we use the RMSPropOptimizer solver defined in TensorflowAbadi et al. (2015) as the stochastic gradient descent solver to minimize the loss function. The details of the neural network model are in Appendix A.

The training data for different phases are sampled from the 3-dimensional Anderson model using different disorder parameters. It’s well known that the critical disorder at $E=0$ for the 3D Anderson model is 16.54 ±0.01 MacKinnon and Kramer (1981, 1983); Kramer and MacKinnon (1993). When the disorder strength $W$ is larger than the critical value, the wave functions are exponentially localized and the system behaves as an insulator. Otherwise, the wave functions are delocalized and the system behaves as a metal. This phenomenon is known as Metal-Insulator Transition(MIT)Anderson (1958). We get 4000 eigenstates from $W\in[14.0,16.0)$ as the delocalized phase and 4000 eigenstates from $W\in[17.0,19.0)$ as the localized phase by steps of 0.1. For each W, we prepare 40 different realizations of randomness and for each realization, we take five eigenstates around $E=0$. For the validation data set, we get another 600 eigenstates from $W\in[10.0,16.0)$ and 600 eigenstates from $W\in[17.0,23.0)$ in steps of 0.1. During each step of the training, we randomly select 256 eigenstates from the training data set as the input and calculate the gradient of the loss function with respect to the parameters in the CNN model and update them. After every 50 steps, we test the prediction accuracy on the validation data set and save the model with the highest prediction accuracy.

To show the prediction accuracy for different disorder parameters $W$, we generate another 16000 eigenstates sampled from the Anderson model using $W\in[0.1,16.0]$ and $W\in[17.0,33.0)$. The prediction accuracy for different disorder strengths $W$ is shown in Fig.6(a), and the overall accuracy is $99.0\%$. The lowest prediction accuracy around the critical disorder $0.8W_{c}<W<1.2W_{c}$ is about $83\%$. We also test our trained model by producing the phase transition diagram of the 3D Anderson model. The testing data are sampled from $W\in[8.0,25.0]$ by steps of 0.1. In each realization of the same disorder parameter $W$, we pick 5 eigenstates around the band center($E=0$) as input data and use the averaged delocalized probability of the five eigenstates as the delocalized probability of this realization. We prepare 5 random realizations for each $W$ and average the delocalized probability. The phase diagram calculated using our trained CNN model is shown in Fig. 6(b). From Fig. 6(b), we see that the trained CNN model successfully captures the Metal-Insulator Transition(MIT).

Owing to its excellent classification accuracy, the trained neural network model is ready to find the extended state in the impurity band. We generate 1000 random realizations for the Hamiltonian in Eq.1 with doping probability $x=5\%$ and disorder parameter $W=-4.5$, and obtain all eigenstates using the exact diagonalization method in Lapack. These quantum states are used as the input data for our trained CNN model to calculate the delocalized probability. We average the probability over 1000 realizations and the result is shown in Fig.7. We can see that the CNN model confirms that delocalized states exist in the impurity band, which is in good agreement with the results obtained by IPR or Thouless number.

In this work, we numerically investigate the properties of the states in the “impurity band” of heavily-doped non-magnetic semiconductors. By using general full exact diagonalization and sparse matrix diagonalization with Jacobi-Division (JADA) method, we find that with a typical doping probability $x=5\%$, the impurity band in the DOS is developed at about $W=-4.5$. We calculate the IPR, Thouless number, and DOS together for different system sizes and study the relationship between them. The data fitting of IPR and system size on four points suggests the existence of the extended states in the impurity band. The Thouless number calculation supports the same conclusion and gives the exact location of mobility edges.

Besides, we also utilize the supervised deep learning method, which is the state-of-the-art method in pattern recognition, to distinguish the extended and localized states in the impurity band. We train a 3D CNN model using the data generated from the Anderson model and then apply the trained neural network model to classify the states in the “impurity band”. Our trained neural network model achieves high accuracy ($99.0\%$) in classifying different states in the Anderson model. The prediction of our trained model on “impurity band” also supports the finding from the relationship between IPR, Thouless number and system size though the predicted locations of mobility edges have small discrepancies. Our calculation gives direct evidence that there are three mobility edges in the impurity band for a specific on-site impurity potential in heavily-doped non-magnetic semiconductors.