Multifractality of ab initio wave functions in doped semiconductors

TL;DR

This paper demonstrates how combining linear-scaling DFT with tight-binding Hamiltonians enables the computation of multifractal wave functions, providing insights into the critical properties of the Anderson MIT in doped semiconductors, and resolving previous discrepancies in critical exponent values.

Contribution

It introduces an ab initio approach to analyze multifractality at the MIT in doped semiconductors, improving understanding of critical exponents and wave function properties.

Findings

Critical multifractal spectrum at the MIT for Si:S obtained.

Resolution of the 'exponent puzzle' with ab initio methods.

Large, realistic samples of sulfur-doped silicon analyzed.

Abstract

In Refs. [1,2] we have shown how a combination of modern linear-scaling DFT, together with a subsequent use of large, effective tight-binding Hamiltonians, allows to compute multifractal wave functions yielding the critical properties of the Anderson metal-insulator transition (MIT) in doped semiconductors. This combination allowed us to construct large and atomistically realistic samples of sulfur-doped silicon (Si:S). The critical properties of such systems and the existence of the MIT are well known, but experimentally determined values of the critical exponent $\nu$ close to the transition have remained different from those obtained by the standard tight-binding Anderson model. In Ref. [1], we found that this ``exponent puzzle'' can be resolved when using our novel \emph{ab initio} approach based on scaling of multifractal exponents in the realistic impurity band for Si:S. Here, after a short review of multifractality, we give details of the multifractal analysis as used in [1] and show the obtained \emph{critical} multifractal spectrum at the MIT for Si:S.

Figures

Equations

$$g(\lambda )={\lambda }^{\kappa }$$

$$⟨{\mu }^q{⟩}_L\sim {\lambda }^{D+{\tau }_q},$$

$$R_q(\lambda )=\sum _i{\mu }_i^q\sim {\lambda }^{{\tau }_q}.$$

$$a\ll l<L\ll \xi$$

$$N_{\lambda }(\alpha )\sim {\lambda }^{-f(\alpha )}\text{and}{\alpha }_i=\frac{\log {\mu }_i}{\log \lambda }.$$

$$R_q(\lambda )=\sum _i{\mu }_i^q=\sum _i{\lambda }^{q{\alpha }_i}=\int N(\alpha ){\lambda }^{q\alpha }d\alpha \sim \int {\lambda }^{q\alpha -f(\alpha )}d\alpha .$$

$$f_q=q{\alpha }_q-{\tau }_q\text{and}{\alpha }_q=\frac{d{\tau }_q}{dq}.$$

$${\alpha }_q+{\alpha }_{1-q}=2d{f}_{1-q}=f_q+d-{\alpha }_q.$$

$$S_q(\lambda )=\frac{d{R}_q(\lambda )}{dq}=\sum _i{\mu }_i^q\log {\mu }_i.$$

$${\tau }_q=\underset{\lambda \to 0}{\lim }\frac{\log R_q(\lambda )}{\log \lambda },{\Delta }_q={\tau }_q-q(d-1),{\alpha }_q=\underset{\lambda \to 0}{\lim }\frac{S_q(\lambda )}{R_q(\lambda )\log \lambda }.$$

$$f_q^{\text{ens}}=q{\alpha }_q^{\text{ens}}-{\tau }_q^{\text{ens}}\text{and}{\alpha }_q^{\text{ens}}=\frac{d{\tau }_q^{\text{ens}}}{dq}=\underset{\lambda \to 0}{\lim }\frac{⟨S_q(\lambda )⟩}{⟨R_q(\lambda )⟩\log \lambda }.$$

$${\sigma }_{{\alpha }_q}={\alpha }_q^{\text{ens}}\sqrt{\frac{{\sigma }_{S_q}^2}{⟨S_q{⟩}^2}+\frac{{\sigma }_{R_q}^2}{⟨R_q{⟩}^2}-2\frac{\mathrm{cov}(S_q,R_q)}{⟨S_q⟩⟨R_q⟩}}\text{ and }{\sigma }_{{\tau }_q}=\frac{{\sigma }_{R_q}}{⟨R_q⟩\log \lambda },$$

$$f(\alpha )\simeq d-\frac{(\alpha -{\alpha }_0{)}^2}{4({\alpha }_0-d)}\text{and}{\alpha }_0=d+\gamma .$$

Full text

††journal: Physica E

1 Introduction

The Anderson metal-insulator transition (MIT) [3] is one of the fundamental manifestations of the quantum mechanical nature of disordered materials [4, 5, 6]. In his 1958 publication [3], Anderson studies the localization of electrons in doped semiconductors. The existence of the MIT in these materials was later confirmed by measuring the scaling of the conductance when increasing the dopant concentration beyond a critical value [7, 8, 9, 10]. The critical properties of the MIT, such as the exponent $\nu$ of the conductivity, should be universal quantities [11, 12]. For classical waves [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] and cold atom systems [24, 25, 26, 27, 28, 29, 12] results agree well with many of the non-interacting Anderson model estimates [30, 31]. However, in experiments with semiconductors, $\nu$ is found to vary with sample-specific properties, namely the dopant concentration around the transition point, the homogeneity of the doping, and the purity of the sample itself [32]. The term “exponent puzzle” [8, 9] was hence coined to describe this inability to characterize the Anderson transition in terms of a single, universal value for $\nu$.

Recently, we presented a study [1] that moves beyond the paradigmatic, highly-simplified, tight-binding Anderson model, and employs atomistically correct ab initio simulations [33, 34] of the doped semiconductor Si:S [35]. With this approach we observe how the impurity band (IB) forms and eventually merges with the conduction band upon increasing the dopant concentration. We then exploit the multifractal nature of the (near-)critical electronic wave functions as a basis for a finite-size scaling analysis that aims to retrieve the critical properties of the transition in the thermodynamic limit. In order to reach sufficiently large system sizes, we devised a hybrid approach: linear-scaling DFT calculations using the ONETEP code [36], on prototype systems of $8\times 8\times 8$ diamond-cubic unit cells ($4096$ atoms) to construct catalogs of local Hamiltonian blocks to describe the no/single/double dopant situation. For each concentration of impurities and disorder realization, we then built much larger, effective tight-binding Hamiltonians $H$ and overlap matrices $O$ from these catalogs for system size $L$, and solved the large generalized eigenvalue problem $H\psi_{j}=\epsilon_{j}O\psi_{j}$ for eigenenergies $\epsilon_{j}$ and eigenvectors $\psi_{j}(\mathbf{r})$, with $\mathbf{r}=(x,y,z)$ for coordinates $x,y,z=1,\ldots,L$.

In this paper, we show in detail how to perform the multifractal analysis of the $\psi_{j}$ and present the resulting critical multifractal spectrum. Our results suggest that it is different from the spectrum in the Anderson model [2]. In addition, we give further details in the finite-size scaling analysis needed to ascertain the existence and the properties of the MIT when $L\rightarrow\infty$.

2 Some of the basics of multifractals

For the Anderson transition from a metal to an insulator upon increasing the disorder [3], the absence of length scales at criticality means that the wave function intensity $|\psi(\boldsymbol{r})|^{2}$ at the critical disorder is self-similar [37, 4, 6]. It needs to have a “filamentary” structure [38], i.e. it needs to be extended throughout the volume, a property of the metal phase, but also to occupy only an infinitesimal fraction of it, a property of the localised phase. This structure allows the critical phase to be continuously connected to both the extended and localised phases. In conjunction with its self-similar property, the critical wave function qualifies as a fractal [39], at least as long as we can disregard the lower limit imposed by the lattice spacing $a$. Castellani et al. [40] realised, based on the earlier work of Wegner [41], that the critical wave function is not a simple fractal, but rather an “interwoven family” of fractals, each with its own dimension and distribution. Such an object is a multifractal [42, 43].

2.1 Self-similarity

Let us consider a system occupying a finite region of space $\Sigma\subset\mathbb{R}^{D}$ with a local density $\rho(\boldsymbol{r})$. Following Ref. [44], we define the pair correlation function, $g(\boldsymbol{r})=\langle\rho(\boldsymbol{r}+\boldsymbol{r}^{\prime})\rho(\boldsymbol{r}^{\prime})\rangle_{\boldsymbol{r}^{\prime}}$, which gives the probability that two points separated by $\boldsymbol{r}$ both belong to the region $\Sigma$. For simplicity, we now assume that the correlation function is isotropic, $g(\boldsymbol{r})=g(r)$. In the absence of length scales, $g$ obeys homogeneity laws (or scale-invariance) with respect to a resolution or coarse-graining $\lambda$. More specifically, if we rescale lengths as $r\rightarrow r^{\prime}=\lambda r$ we have that $g(r^{\prime})=\lambda^{\kappa}g(r)$, where $\kappa$ is a homogeneity exponent. The solution to this equation is given by a power-law behaviour, $g(r)\propto r^{\kappa}$. If we then fix $r$ as the reference length scale, and $g(r)=1$, we see that

[TABLE]

translates the concept of self-similarity into a mathematical relation. We can hence describe fractal objects as self-similar structures whose observed spatial extent (e.g. volume) depends, with a power-law behaviour, on the resolution at which we look at it. For fractals originating from a mathematical relation, the dependence on the resolution can extend over an infinite range. For fractals appearing in physical systems, instead, the range of $\lambda$ is usually limited by macro- and/or microscopic scales. A very comprehensive list of examples can be found, e.g., in [45].

2.2 Measures, fractals and multifractals

Let $\psi(\boldsymbol{r})$ be the wave function of an electron in a $L\times L\times L$ volume. The modulus square $|\psi(\boldsymbol{r})|^{2}$ defines a normalised measure on this volume, which we partition in boxes of linear size $l=\lambda L$. The number of boxes will then be $\lambda^{-d}$, where $d=3$ is the Euclidean dimension of the support of the system. The probability of finding the electron in box $\mathcal{B}_{i}$ is the box-probability $\mu_{i}=\int_{\mathcal{B}_{i}}|\psi(\boldsymbol{r})|^{2}\operatorname{d}\boldsymbol{r}$. We can then compute the fractal dimension $D$ of the system by counting the number of boxes where the box-probability does not vanish:111It is customary to use $\sim$ to indicate that the proportionality constant is independent of the resolution and can thus be ignored. This constant might appear, for instance, when the boxes, whichever their shape, do not perfectly cover the system. Since in our analysis we are covering boxes with boxes, most of the relations in this section are actually equalities. $N(\lambda)\sim\lambda^{-D}$. Because the electron can access any portion of the volume, i.e. there is no region of space with vanishing probability, we trivially conclude that the fractal dimension is $D=d=3$.

Compared to the fractal dimension, more insightful is actually the study of the powers of the box-probability $\mu^{q}_{i}$, which is the idea behind multifractal analysis. If the wave function is a multifractal, we expect to see the power-law behaviour of (1):

[TABLE]

where $\langle\ldots\rangle_{L}$ denotes the average over all boxes in the volume. Equivalently, we can introduce the partition sum $R_{q}(\lambda)$ (also the generalised inverse participation ratio) as $\langle\mu^{q}\rangle_{L}=\lambda^{D}R_{q}(\lambda)$ and write

[TABLE]

The mass exponents $\tau_{q}$ describe the scaling behaviour of the moments and do not depend on $\lambda$.

Let us stress again that multifractality holds if, in the power-law relation of Eq. (2), we have $\tau_{q}\neq 0$ for a finite range of $\lambda$: the box size $l$ should be smaller than the system size, but also larger than the microscopic scale $a$, usually equal to the lattice spacing in the Anderson case. At the same time, for critical states at the Anderson transition, the system size is much smaller than the correlation length $\xi$, such that

[TABLE]

Additionally, the wave function is truly critical (and hence multifractal) only in the thermodynamic limit, where both $L$ and $\xi$ diverge. Thus $\tau_{q}$ is defined in the limit $\lambda\rightarrow 0$. For finite systems, instead, we choose states and coarse-grainings that satisfy (4). In this case, we can estimate $\tau_{q}$ by fitting the slope of $\log R_{q}(\lambda)$ versus $\log\lambda$. We are assuming here that multifractality survives in finite systems [46], and postpone the discussion of this non-trivial assumption to section 5.

From (2) and the normalisation of the wave function, it is possible to show that $\tau_{0}=-D$ and $\tau_{1}=0$. This implies that we can generalise the definition of the fractal dimension to the anomalous dimensions $D_{q}$ such that $D_{0}=D$ and mass strengths $\tau_{q}=D_{q}\,(q-1)$. In the case of a simple fractal $D_{q}\equiv D$, while for a multifractal $D_{q}$ has a non-trivial dependence on $q$. The deviation from the simple-fractal case is captured by the anomalous scaling exponent $\Delta_{q}=(D_{q}-d)(q-1)=\tau_{q}-d(q-1)$.

2.3 The multifractal spectrum

The scaling of the moments $R_{q}$, yielding $\tau_{q}$, is enough to fully characterise the multifractal nature of the wave function. Now we present an equivalent description of the multifractal that will be useful, in the sections that follow, to validate our results and compare them to the 3D Anderson model. This description is founded on a multifractal measure [47], a distribution such that, around each box $\mathcal{B}_{i}$, $\mu_{i}=\lambda^{\alpha_{i}}$. The set of boxes with $\alpha_{i}\in[\alpha,\alpha+\operatorname{d}\alpha]$, then, constitutes a simple fractal with dimension $f(\alpha)$, such that the number of said boxes is

[TABLE]

This is the formalisation of the idea of Castellani [40] that the multifractal is composed of different simple fractals. We re-express the partition sum of Eq. (3) as

[TABLE]

For small $\lambda$, we can use the saddle point approximation and find that the biggest contribution in the integral (6) comes from the value of $\alpha$ that maximises (since $\lambda<1$) the argument of the exponential, i.e. the $\alpha_{q}$ such that $f^{\prime}(\alpha_{q})=q$. We can then write, from (3), $\tau_{q}=q\alpha_{q}-f(\alpha_{q})$. If we identify $f_{q}=f(\alpha_{q})$ we can see that $(q,\tau_{q})$ and $(\alpha_{q},f_{q})$ are related by a Legendre transformation

[TABLE]

It can be proven, e.g. in [48], that $\tau_{q}$ is a monotonically increasing function in $q$, which implies that $\alpha_{q}>0,\forall q$.

We can combine singularity strengths $\alpha_{q}$ and the singularity spectrum $f_{q}$ to obtain the multifractal spectrum $f(\alpha)$. This function is equivalent to the generalised dimensions $D_{q}$ in characterising the multifractal, and in the case of a simple fractal analogously reduces to the point $(D,D)$ in a $(\alpha,f(\alpha))$ plot. As shown in the example of Fig. 2, $f(\alpha)$ is a convex function reaching its maximum at $\alpha_{0}$ with a value $f_{0}=\tau_{0}=D$. From (7) we further notice that $f_{1}=\alpha_{1}$, since $\tau_{1}=0$. The spectrum is therefore tangential to the functions $f_{0}(\alpha)\equiv D$ and $f_{1}(\alpha)=\alpha$.

2.4 Symmetry of the multifractal spectrum

Using the nonlinear $\sigma$ model, Mirlin et al. [49] have analytically proven that at criticality the multifractal exponents (7) satisfy the exact symmetry relation

[TABLE]

Assuming the universality of the critical properties at the Anderson transition, this result is expected to generally hold for the Wigner-Dyson symmetry classes [6]. Indeed, Eq. (8) was confirmed numerically for different systems, including the 3D Anderson model [50, 51] and in experiments [20]

3 Multifractal analysis of the wave function

We are mainly interested in the singularity strengths $\alpha_{q}$, which, together with $\tau_{q}$ and $\Delta_{q}$ are called multifractal exponents (MFE). In this section we recast the exponents derived in Sec. 2 in a form that is more convenient for numerical calculations, mostly by reducing the loss of precision. We then extend our definitions to include a disorder ensemble average.

Before moving to the determination of these exponents, let us briefly discuss their expected values [51, 52]. In the case of very low disorder, the wave function intensities $|\psi(\mathbf{r})|^{2}$ will be nearly plane-wave like, i.e. extended throughout the whole of the volume $L^{3}$. Hence, as shown in section 2.2, $\alpha_{q}=D=d=3$ for all $q$. This implies that $\tau_{q}=d(q-1)$, $\Delta_{q}=0$ and $f_{q}=d$. Hence $f(\alpha)$ is seen to contract, and eventually converge, to one point $f(d)=d$. On the other hand, for very large disorder, well into the insulating/localized regime, $|\psi(\mathbf{r})|^{2}$ can be seen as localized on a few (single) sites only. In this case $\alpha_{q}\rightarrow\infty$ for $q\leq 0$, and $\alpha_{q}\rightarrow 0$ for $q>0$. Similarly, $\tau_{q}\rightarrow-\infty$, $\Delta_{q}\rightarrow-\infty$ for $q\leq 0$, and $\tau_{q}=0$, $\Delta_{q}=d(1-q)$ for $q>0$. The $f(\alpha)$ spectrum broadens and in the limit of strong localization, converges to the points $f(0)=0$ and $f(\infty)=d$. At criticality, the behaviour is even richer [52], with a non-parabolic $f(\alpha)$ [53] and $\alpha_{0}=4.048(4.045,4.050)$, $\alpha_{1}=1.958(1.953,1.953)$ the current best estimates for the non-interacting 3D Anderson model [52]. As a rough guide for the following sections, the $\alpha_{q}$ tend towards $d=3$ for extended states as weak disorder, while $\alpha_{0}$ increases without bounds in the insulating regime.

3.1 Numerical calculation

Following [54], it is convenient to define, from (3) and (7), the auxiliary quantity

[TABLE]

This ratio can be interpreted as an average with respect to the measure defined by $\mu^{q}$. The latter is also called $q$-microscope, because it increases the large (small) box-probabilities for $q>0$ ($q<0$). A computationally-friendly formulation of the MFE reads

[TABLE]

To comply with (4), we choose $\lambda\leq 1/2$, namely we consider boxes of linear size up to $l\leq L/2$. We coarse-grain the wave function using the partitioning scheme proposed by [55]. Here, the box size $l$ can take any integer value (up to $L/2$), so that $\lambda^{-1}=L/l$ can take non-integer values. This is achieved by first periodically replicating the original system, such that it can be exactly covered by an integer number of boxes, and then by averaging over the possible equivalent box origins. The increased number of available box sizes translates, in the linear fits, in reduced uncertainties in the estimated slopes.

3.2 Ensemble averaging

So far we have computed the multifractal properties of a single wave function. The multifractal analysis of the Anderson transition is usually performed by taking an average over the disorder realisations. The definitions of the MFE can be extended by defining the ensemble average of the partition sum as $\langle R_{q}(\lambda)\rangle\sim\lambda^{\tau_{q}^{\text{ens}}}$, such that $\tau_{q}^{\text{ens}}=\lim_{\lambda\rightarrow 0}{\log\langle R_{q}(\lambda)\rangle}/{\log\lambda}$. We then proceed to take the Legendre transform and define

[TABLE]

Notice that, in the ensemble average of $\alpha_{q}$, the $q$-microscope $\mu_{i}^{q}$ is normalised by $\langle R_{q}\rangle$, namely the averaged partition sum of the wave function. If we normalised the $\mu_{i}^{q}$ terms for every wave function we would obtain the typical average $\alpha_{q}^{\text{typ}}=\langle S_{q}/R_{q}\rangle/\log\lambda$ ($\lambda\rightarrow 0$). While in the ensemble average all wave functions, including rare events, are equally weighted, the typical average is dominated by the behaviour of “typical” wave functions. The presence of rare events translates in the appearance of negative fractal dimensions (see Sec. 4.1), a feature of the $f(\alpha)$ that is best captured by ensemble averaging [6].

4 Results for eigenstates of the effective Hamiltonians

For finite systems, the multifractal exponents are computed, as explained in Sec. 2.2, by estimating the slope of a $\log\langle R_{q}(\lambda)\rangle$ vs. $\log\lambda$ plot, in the case of $\tau_{q}^{\text{ens}}$. Accordingly, the statistical uncertainties at fixed $\lambda$ have to be multiplied by a factor $\log\lambda$.222The standard deviation $\sigma_{\alpha_{q}}$ associated to $\alpha_{q}^{\text{ens}}$ (at a fixed $\lambda$) is related to the standard deviations $\sigma_{S_{q}}$ and $\sigma_{R_{q}}$ and the covariance $\operatorname{cov}(S_{q},R_{q})$ via propagation of the variance [52]

$\sigma_{\alpha_{q}}=\alpha_{q}^{\text{ens}}\sqrt{\frac{\sigma^{2}_{S_{q}}}{\langle S_{q}\rangle^{2}}+\frac{\sigma^{2}_{R_{q}}}{\langle R_{q}\rangle^{2}}-2\frac{\operatorname{cov}(S_{q},R_{q})}{\langle S_{q}\rangle\langle R_{q}\rangle}}\quad\text{ and }\quad\sigma_{\tau_{q}}=\frac{\sigma_{R_{q}}}{\langle R_{q}\rangle\log\lambda}\,,$

(12)

so that $\sigma_{f_{q}}^{2}=\sigma_{\alpha_{q}}^{2}+\sigma_{\tau_{q}}^{2}$. The standard error of the mean is obtained by dividing every standard deviation by $\sqrt{\mathcal{N}}$, where $\mathcal{N}$ is the number of available realisations. In Fig. 3 we show the average singularity spectrum for the ensemble of $10\,648$ atoms with $140$ impurities, a system that is critical at energy close to $-0.249\text{\,}\mathrm{e}\mathrm{V}$ [1]. The increase in the ensemble size does not change the shape of the spectrum significantly, but has the effect of reducing the error bars on the data points. In particular, it is the calculation of the average $\langle R_{q}\rangle$ and $\langle S_{q}\rangle$ in Eq. (11) that benefits from larger ensembles, since smaller error bars in the data used for the linear fits results in smaller uncertainties on the fit parameters, see Fig. 3(b). We quantitatively report the quality of said fits in the lower panel of Fig. 3(a), where we show the linear correlation coefficient $r^{2}$ and the $p$ value. As noted in [50], while $r^{2}\approx 1$ indicates a good linear behaviour, small $p$ values suggest that the uncertainties on the data point are too small to support the deviation from the linear behaviour we are fitting. This is likely due to the limited number of realisations available for the ensemble averaging. For comparison, at the end of the two branches, i.e. for large $|q|$ values, error bars are larger and hence the quality-of-fit increases again.

4.1 Negative fractal dimensions

Error bars increase on the two ends of the spectrum, for negative $q$ (right) and positive (left). For $q<0$, the $q$-microscope increases the weight of small values of the wave function, which are more sensitive to numerical fluctuations from the diagonalization. The other end of the spectrum ($q>0$, left) describes instead the presence of rare critical wave functions with small values of $\alpha$ and hence large $|\psi|^{2}\sim L^{-\alpha}$ [53]. The set of these values scales with a negative fractal dimension $f(\alpha)$, which means that their occurrence frequency vanishes in the $L\rightarrow\infty$ limit. This is a known effect arising from ensemble averaging [56] and known since the pioneering work of Mandelbrot [39]. This finite-size effect is further observed and commented in [50], where larger systems are accessible and studied. In comparison, the single state used to produce Fig. 2 does not show said rare boxes with large probability amplitudes, as indicated by $f_{q}>0$.

4.2 Width of the multifractal spectrum

Let us comment on the width of the distribution in Fig. 3, as compared to the Anderson model studied in [50]. A narrow $f(\alpha)$ spectrum implies that extreme values (either large or small) occur less frequently. This means that, in our case, the average state near criticality in our model looks more homogeneous or extended than in the Anderson model. In Fig. 4, instead, we show the singularity spectrum for a system of $10\,648$ atoms, which, for $230$ impurities, is close to criticality at the Fermi energy $\varepsilon_{\text{F}}=0$ and deeper in the impurity band at $-0.320\text{\,}\mathrm{e}\mathrm{V}$ (estimated from [1]). While both spectra are narrower than the Anderson model, the critical wave function at $\varepsilon_{\text{F}}$ appears on average more extended than deeper in the impurity band.

Our results are reminiscent of those found by [57] for the power-law random banded matrix (PRBM) model, which describes a 1D chain with random long-range hopping decaying as $r^{-\alpha}$ over distances larger than a band width $b$. For the critical value $\alpha=1$, the model undergoes an Anderson transition for any value of $b$, which parametrises a family of critical models that can be studied from the weak- ($b\gg 1$) to the strong-coupling ($b\ll 1$) regime. For $b\gg 1$ the model shows a “quasi-metallic” behaviour, where the critical wave functions has statistical properties similar to the delocalised phase. The singularity spectrum becomes correspondingly narrower with a parabolic shape, a regime called weak multifractality. In this case the multifractal spectrum follows the parabolic approximation [48]:

[TABLE]

In Fig. 3 we fit our full-ensemble data to (13) to find the estimate value $\alpha_{0}\approx 3.55$. We only report an approximate value without uncertainty because the parabolic behaviour is an approximation and does not necessarily hold for the whole spectrum, since $f(\alpha)$ is defined only for positive $\alpha$. In fact, for $q_{c}=(d+\gamma)/2\gamma$ [6] we have $\alpha_{q_{c}}=0$ and $f(0)$ is finite (a termination point), whereas our results in Fig. 3 and those in [50] seems to suggest that $f(\alpha)\rightarrow-\infty$ in the limit $\alpha\rightarrow 0$ (hence no termination point).

4.3 Symmetry of the multifractal spectrum

In Fig. 3 we also show the symmetrised spectrum obtained by computing and plotting $\alpha_{1-q}$ and $f_{1-q}$ from Eq. (8). As expected from the previous paragraphs, the uncertainty on the data points increase at the extremities. Within these error bars, the spectra are in good agreement with each other. We verify the same symmetry relation also for the spectra in Fig. 4. While at $-0.320\text{\,}\mathrm{e}\mathrm{V}$ there is excellent agreement between the spectra, at the Fermi energy there is a slightly higher discrepancy, especially at extreme values of $q$. Since this discrepancy would be resolved by taking two standard deviations as confidence intervals, instead of one, we cannot attribute this discrepancy to any specific underlying physical factor or systematic error.

5 Validity of the scaling assumption

The question that arises when dealing with finite systems is whether the wave function is still a multifractal. Formally speaking, the wave function is a true multifractal only at the critical point. For a finite $L$, however, the effective critical point shifts away as $L^{-1/\nu}$ from its thermodynamic limit [11]. Luckily this is not a problem, since, as shown by [46], states on the two sides of the transition still show multifractal features characteristic of a critical wave function.

Now that we can construct a multifractal measure from the wave functions away from the critical point, we can actually check the most important assumption we have taken so far, namely that a localisation-delocalisation transition occurs in our model. The histogram distribution $N_{\lambda}(\alpha)$ of the measure $\alpha$ depends, at the critical point, only on the coarse-graining $\lambda=l/L$, rather than separately on the system size $L$ and the box size $l$. At the critical point, then, $N_{\lambda}(\alpha)$ has the same shape for any $L$, provided that the wave functions are coarse-grained with the matching $l$ box size. The dependence of $N_{\lambda}(\alpha)$ on $L$ gradually reappears away from the critical point, where in the strong (weak) disorder regime, larger systems become more localized (delocalized). This is shown in Fig. 5, where we plot the ensemble histogram $\operatorname{PDF}(\alpha)=N_{\lambda}(\alpha)\lambda^{d}/\mathcal{N}$ at $\lambda=1/2$ for three values of $n$: the lowest in the localized regime, the intermediate close to the critical point and the highest in the delocalised regime. The ensemble $\operatorname{PDF}(\alpha)$ is built by filling one single histogram with the $\mathcal{N}$ available wave functions.

Because of the limited spread in system sizes and the large common coarse-graining, the difference in the PDFs outside the critical point is not very well pronounced in Fig. 5. An alternative check consists in fixing the system size and study how the PDF renormalises with the coarse graining [58]. The “scaling” variable $\xi/L$ [59] scales like $\lambda^{-1}$, which implies that, with increasing $\lambda$, $\xi/L$ becomes smaller. Physically this means that, upon coarse-graining, localised (delocalised) states become more localised (delocalised), or, equivalently, that the renormalisation flow rescales the disorder away from its critical value,if a phase transition, and hence a critical point, exists [2]. We verify this in Fig. 6 where the PDF’s move in opposite directions upon increasing the box size $l$ in a system of $L^{3}=$4096$$ atoms.

6 Conclusions

After reviewing the foundations of multifractal analysis in the study of a disordered system at criticality, we have shown the multifractal nature, captured by the $f(\alpha)$ spectrum, of the Kohn-Sham wave function. This result is in line with previous experimental [60], theoretical [61], and numerical [62, 63, 64] studies, where the critical fluctuations of the wave function at criticality are expected to survive in the presence of the Coulomb interaction.

We have shown that multifractality persists also in the wave functions near the critical concentrations computed in the effective tight-binding model presented in Refs. [1, 2]. The multifractal spectrum follows the symmetry relations derived from field-theoretical models [49] and the ensemble average over hundreds of realizations results in negative fractal dimensions ascribed to rare events. These are known features of (finite) multifractal critical states that have been studied also in the non-interacting 3D Anderson model [31]. A difference with this case, however, lies in the nearly-parabolic shape of the $f(\alpha)$ spectrum observed here. This behaviour, referred to as weak multifractality, is a common trait with the Anderson transition in $2+\epsilon$ deminsions, with $\epsilon\ll 1$, and with the power-law random banded matrix model with $b\gg 1$ [6].

Finally, we compute the $f(\alpha)$ spectrum near the two mobility edges, respectively, far from and at the Fermi energy. In the latter case we notice that the spectrum is narrower, indicating that the critical wave function at this energy is, on average, more delocalized than deeper in the band. In Refs. [1, 2] we propose the origin of this quasi-metallic behaviour to be the hybridization of the impurity states near the Fermi energy with states from the conduction band. The presence of this second extended band might alter the physics of the metal-insulator transition close to the Fermi energy, leading to a varying critical exponent across the impurity band.

A similar phenomenon has very recently been reported in a Hartree-Fock study of the Anderson transition in the presence of tunable Coulomb interaction [64]. There the authors show that a second mobility edge appears near the Fermi energy, where the Coulomb gap forms. At this energy, the critical state shows a narrower multifractal spectrum compared to the mobility edge at higher energy. While in this model quasi-metallic behaviour appears where the Coulomb gap forms in the centre of the band of the Anderson model, it is present in our model when the impurity band forms and then merges with the conduction band.

To conclude, the numerical analysis of the Anderson transition in the presence of interactions, and, in particular, in real materials, is still very challenging. Progress relies on a large amount of resources needed for the simulations, with recent first-principles or self-consistent calculations reaching system sizes of the order of $10^{2}$–$10^{4}$ sites, and of hundreds of realisations [1, 64, 65] — while studies of the non-interacting Anderson model can nowadays easily achieve more then $10^{6}$ sites and $10^{4}$ samples [52, 58]. Still, ab initio studies of real materials have observed new phenomena and traced these back to the electronic interaction. Future investigations with increasing system and sample sizes will undoubtedly clarify the universal properties of the transitions and critical points.

Acknowledgements

RAR thanks Presidency University, Kolkata, where part of this work was written, for their warm hospitality. We thank EPSRC for support via the ARCHER RAP project e420 and the MidPlus Regional HPC Centre (EP/K000128/1). UK research data statement: data is available at http://wrap.warwick.ac.uk/id/eprint/92910.