Fluctuation corrections to Lifshitz tails in disordered systems

Quenched disorder in semiconductors induces localized electronic states at the band edge, which manifest as an exponential tail in the density of states. For large impurity densities, this tail takes a universal Lifshitz form that is characterized by short-ranged potential fluctuations. We provide both analytical expressions and numerical values for the Lifshitz tail of a parabolic conduction band, including its exact fluctuation prefactor. Our analysis is based on a replica field integral approach, where the leading exponential scaling of the tail is determined by an instanton profile, and fluctuations around the instanton determine the subleading pre-exponential factor. This factor contains the determinant of a fluctuation operator, and we avoid a full computation of its spectrum by using a Gel'fand-Yaglom formalism, which provides a concise general derivation of fluctuation corrections in disorder problems. We provide a revised result for the disorder band tail in two dimensions.

Disorder is an inherent property of physical systems, arising, for example, as defects or impurities in doped semiconductors [1], frozen non-equilibrium degrees of freedom in quenched alloys and glasses [2,3], or thermal fluctuations [4].It can also be engineered in photonic structures [5,6] or quantum gases in an optical speckle potential [7,8].For the electronic density of states in a semiconductor, disorder has two effects: First, it induces a narrowing of the band gap, where valence-and conduction-band levels are raised and lowered, respectively, by a characteristic impurity energy.Second, it causes band tailing, where fluctuations in the disorder potential give rise to localized states inside the band gap.Describing the disorder-induced band narrowing and tailing is important because it affects semiconductor devices such as transistors [1,9,10] or limits the efficiency of solar cells [11,12], and it sets the variable-range hopping conduction that dominates transport in localized systems at low temperatures [13].
In general, potential fluctuations that are large enough to generate states with energies far away from the band edge are rare, and the tail of localized levels takes an exponential form (written here for a single conduction band) [1,[14][15][16], where E is the energy below the band edge and ⟨•⟩ denotes a disorder average.Of particular interest is the limit of high impurity density, which is universal in the sense that the band tail is insensitive to details of the disorder correlations over a large energy range [1].This is the Lifshitz region or Lifshitz tail [17,18], where (1) takes a stretched exponential form [i.e., the argument B(E) of the exponential has a power-law dependence on E].Remarkably, even for the seemingly simple case of nonin-teracting particles in a Gaussian-correlated disorder potential, the disorder-averaged density of states can only be computed exactly in one space dimension (1D) [19].In higher dimensions, dating back to works by Halperin and Lax [19][20][21] and Zittartz and Langer [22], the tail exponential B(E) is obtained from a variational argument that determines the most likely disorder potential that gives rise to a bound state at large negative energy E. Such variational arguments are quite general [23,24] and may be extended to include the effect of correlated disorder [25], magnetic fields [26], interactions [27], and also to describe the spectra of random operators [28][29][30].However, they do not capture the prefactor A(E) in the density of states (1), which is set by fluctuations of the disorder potential and which requires the evaluation of a functional determinant.The established result for A(E) derived by Brézin and Parisi [31] is based on an analysis of large-order perturbations in ϕ 4 theory discussed in a classical paper by the same authors [32], which evaluates the fluctuation determinant by exact diagonalization [32,33].At the same time, the computation of fluctuation corrections is a recurring problem in other areas, for example to determine the decay rate of metastable states [34].Indeed, one of the first studies on the subject by Gel'fand and Yaglom [35] and Levit and Smilansky [36] in the context of one-dimensional path integrals in quantum mechanics avoids the direct evaluation of the fluctuation spectrum by mapping the functional determinant to a differential equation, a more tractable problem compared to exact diagonalization.
In this Letter, we show that the fundamental problem of computing the band tail prefactor A(E) of the Lifshitz tail (1) has an efficient solution that does not rely on evaluating the full fluctuation determinant.This is achieved using a Gel'fand-Yaglom approach, which in particular provides a revised result for the band tail in two dimensions.While the calculation is very general, we focus on the universal Lifshitz regime of infinitely dense but infinitesimally weak point scatterers, which is described by a Gaussian white noise disorder potential, and consider a parabolic conduction band described by the one-electron Hamiltonian Throughout this Letter, we set ℏ = 2m = 1.We base our calculation on a functional integral representation for the disorder-averaged density of states, which we evaluate in a saddle-point approximation including fluctuations.The central results of our work are the following revised expressions for the density of states in the deeptail regime (i.e., for expressed here in relation to the density of states of the free Hamiltonian, ).The 1D and 3D cases agree with previous results derived using other means [19,20,31], but the 2D result corrects an existing literature value [31].
The starting point of our analysis is a functional integral representation of the disorder-averaged density of states Here, the probability measure for the random field V (x) is normalized to unity, L d is the d-dimensional volume, and E V n is the nth energy eigenvalue for a realization of the Hamiltonian (3) with potential V (x).In the second line, we introduce the partition function which we express as a path integral over a scalar field ϕ(x).Written in this form, Eq. ( 5) is equivalent to the Edwards-Jones formula for the eigenvalue distribution of the random matrix H V [37,38].The factors (E V n − E) −1/2 in the partition function have a branch cut for E > E V n , resulting in a discontinuity between complex energies E ± i0, and thus a nonzero contribution to the density of states (5).The analytic continuation of the path integral to energies E ± i0 is performed by rotating the fields ϕ in Eq. ( 6) into the complex plane by an angle e ±iπ/4 , i.e., on a contour C ± along the positive and negative complex diagonals [34,[39][40][41].Note that with these contour choices, the integral representation ( 6) is convergent for all E with nonzero imaginary part.
The disorder average ( 5) is difficult to compute directly due to the presence of the logarithm; therefore, we apply the replica trick [42], ln Z = lim N →0 (Z N − 1)/N , which avoids the logarithm in favor of introducing N copies (replicas) of the original system.The disorder average is now performed exactly as a Gaussian integral over the potential V (r), which gives where the subscript indicates that we consider the difference over both contours, Φ(r) = −w 2 /2E (ϕ 1 , . . ., ϕ N ) is a vector of dimensionless replica fields that depend on a dimensionless position r = x √ −E, and we introduce the scaling parameter The effective action corresponds to a ϕ 4 Ginzburg-Landau theory [40], where the disorder field is removed at the expense of introducing a nonlinear term.Crucially, for d < 4, the prefactor g −1 of the effective action (7) is large in the deep-tail regime E → −∞.This allows us to evaluate the integral over the replica field Φ in the saddle-point approximation [43], that is, by expanding around the dominant contribution to the action.Before evaluating the saddle point, however, note that the theory (9) requires renormalization.The divergent one-loop self-energy corrections are shown in Fig. 1, where a continuous line denotes the free Green's function and the vertex is split to connect two replica indices.Physically, the divergences are a consequence of the white-noise approximation (2) to the disorder potential, which breaks down as the field Φ varies over length scales smaller than the separation between scatterers.The divergences are subtracted by adding a counterterm to the action (9), which is linear in g and will thus give a subleading correction in the saddle-point approximation.Here, is the free Green's function at coinciding points, which diverges in d ≥ 2, and the factor N + 2 in Eq. ( 10) arises from the different contractions of the replica fields.The self-energy Σ(E) computed including the counterterm is finite and represents the band narrowing that shifts the reference point E 0 of the energy scale.Including the finite part of G 0 in the counterterm corresponds to a renormalization choice Σ(E)| 0 = 0 that puts this reference at the origin [44].Another possible choice is ∂Σ(E)/∂E| 0 = 1, which is motivated from a coherent potential approximation [31,45].Results computed in different schemes are identical and linked by a shift in the energy parameter (with possible logarithmic corrections in 2D).
The leading saddle point of the action ( 9) is the trivial configuration Φ = 0; however, this contribution is the same for both contours C ± and cancels when taking the difference in Eq. (7).The parameter B(E) in the exponential band tail (1) is thus determined by the action at the subleading saddle point of Eq. ( 7), which is of order S[Φ] = O(g −1 ).The field configuration that minimizes the effective action ( 9) is of the form Φ(r + r 0 ) = f (r) u [41,46], where r 0 is an arbitrary center, u an arbitrary unit vector in replica space, and f (r) is a nodeless radial function known as the instanton.The saddle-point equation is with boundary conditions f ′ (0) = f (r → ∞) = 0. We can interpret it as a Schrödinger equation for the potential −f 2 (r) that describes the most likely shape of the disorder potential with a bound state at energy E [27].
Figure 2 shows the instanton profile in d = 1, 2, and 3.In 1D, the solution f (r) = √ 2 sech(r) is known analytically, and for d ≥ 2, we employ a spectral renormalization method to solve Eq. ( 12) numerically [47].The exponen- tial band tail is then given by the saddle-point action where The corresponding values for I 4 are well established [33,41] and tabulated in Table I.
Our main task in this Letter is to determine the leading contribution to the prefactor A(E) of the exponential band tail (1) by evaluating the next order in the saddlepoint approximation.To this end, we expand the action to second order in new fields that represent longitudinal (δf ∥ ) and transverse (δf k ⊥ ) fluctuations, where the set {u, v 1 , . . ., v N −1 } forms an orthonormal basis of replica space.However, certain fluctuations do not represent small corrections to the saddle-point action.These so-called zero modes correspond to changes of the orientation u of the instanton in replica space and of its center r 0 .These changes leave the action invariant and thus cannot be treated perturbatively but must be integrated exactly.This is done using the method of collective coordinates [48,49], resulting in an overall factor Here, terms in parenthesis are the Jacobian of the transformation to the collective position r 0 and angle coordinate u.Note that this contribution introduces additional factors of g and thus contributes to the energy dependence of the prefactor A(E).The remaining terms correspond to the measure of the degenerate symmetry spaces, which is the N -dimensional surface element in replica space and the real-space volume L d .Formally, the full fluctuation correction reads where the factor I 2 is the saddle-point value of the field Φ 2 in Eq. ( 7), the exponential prefactor is the counterterm contribution (10), and we define the fluctuation kernel which describes the decoupled longitudinal (z = 1) and transverse (z = 1/3) fluctuations.The excluded zero modes that correspond to an infinitesimal shift along one coordinate axis α = 1, . . ., d are linked to d zero modes of the longitudinal fluctuation kernel ∆ 1 given by ∂ α f (r).
Likewise, the operator ∆ 1/3 has a zero mode f (r) corresponding to rotations in replica space.The path integral in Eq. ( 17) over the fluctuations δf ∥ and δf ⊥,k excludes the zero modes.It is then a Gaussian integral and results in factors involving the functional determinant of ∆ z , where for the shifted operator ∆ z +ε the zero modes have eigenvalue ε, and are thus removed in the limit ε → 0 (m is the degeneracy of zero modes, with m = d for z = 1 and m = 1 for z = 1/3).Taken together, the longitudinal and transverse fluctuation contribution to the path integral (17) reads where the exponent in the third factor accounts for the N − 1 transverse replica field directions, and we split up the counterterm contribution.Note that the longitudinal kernel D(1) is negative since it contains a single eigenfunction f with negative eigenvalue −2.The square root is still well defined after analytic continuation, and the resulting imaginary unit cancels with that in the prefactor of Eq. (17).Moreover, in d > 1, the fluctuation determinant D(z) diverges.The mathematical origin of this divergence is identified by expanding the determinant in powers of z, where G 0 is introduced in Eq. ( 11) and higher-order terms in z are finite.The divergent term in Eq. ( 21), however, is precisely canceled by the leading-order counterterm (10), rendering the fluctuation determinant (20) manifestly finite.
To evaluate the finite renormalized fluctuation determinant, one could in principle compute the eigenvalues of ∆ z and evaluate their product, which was done in a similar context by Brézin and Parisi [32] to determine the large-order behavior of perturbation series in ϕ 4 theories.Here, we make use of a more direct method using spectral functions.This method was originally proposed

Re[λ] Im[λ]
Initial contour γ enclosing the eigenvalues F i n a l c o n t o u r γ ′ e n c l o s i n g t h e b r a n c h c u t θ FIG. 3. Complex structure of the integral (24).Crosses indicate eigenvalues (shown here for longitudinal fluctuations with one negative eigenvalue) and the wavy line denotes the branch cut.The dashed line indicates the initial integration contour γ that encloses all eigenvalues, and the solid line shows the deformed contour γ ′ that runs around the branch cut.
by Gel'fand and Yaglom [35], and applied to higherdimensional problems with radial symmetry in Ref. [50] (for a review, see Refs.[51][52][53]).Crucially, the Gel'fand-Yaglom method does not require knowledge of the eigenvalues, but instead expresses the functional determinant in terms of the solution to an ordinary differential equation, which greatly simplifies the overall calculation.
The Gel'fand-Yaglom method is based on the formal identity [36] det where λ n are the eigenvalues of the operator ∆ z , defined in Eq. ( 18), and We introduce an analytic function F z (λ) with zeros that coincide in location and multiplicity with the eigenvalues λ n , in terms of which Here, the contour γ encircles all eigenvalues, and we place the branch cut of the integrand at an angle θ with the real axis to avoid overlapping with the eigenvalues [51]; see Fig. 3 for an illustration.To obtain the last identity, we deform the contour to enclose the branch cut, which is shown as a continuous line in Fig. 3.We may now take the derivative with respect to s to obtain where we use the rotational symmetry of ∆ z to factorize the determinant and by extension F into angular components F (l) [50], and we denote the angular degeneracy of the eigenvalues by d(l; d).Note that in 1D the rotational symmetry reduces to a parity symmetry; thus, l labels only nondegenerate even (l = 0) and odd (l = 1) contributions.Intuitively, the function F z (λ) generalizes the characteristic polynomial of finite-dimensional matrices to the operator ∆ z .Equation ( 25) then corresponds to the result that the constant term in the characteristic polynomial is the determinant.
To obtain F (l) z (λ), we follow the same heuristic as used to find eigenvalues numerically with the shooting method [51]: Consider the function ϕ (l) z (λ, r) obtained by solving the initial value problem where we include the factor ε ≪ 1 in Eq. ( 19), and If λ is such that the function vanishes at infinity, ϕ z (λ, ∞) = 0, ϕ will also be an eigenfunction of ∆ (l) z with eigenvalue λ.But this is exactly the defining property of F z , and hence we set contributes to Eq. ( 25), so it is useful to derive an equation for R with boundary conditions R z,ε (0) = 1 and R (l)′ z,ε (0) = 0, inherited from Eq. (26).Here, I is the modified Bessel function of the first kind.
The renormalized determinant, finite to all orders in z and without zero modes, reads where δ l is equal to unity if ∆ z has a zero mode, and zero otherwise.The second term in the brackets is the counterterm contribution, which in our discussion of Eq. ( 21) we identified as the term of linear order in z, and which is thus written as where ν = l + d/2 − 1.In the last identity we used the connection between R and the established asymptotic expansion of the Jost function in scattering theory [50].The values of the renormalized functional determinant e 3zI2G0 D(z), Eq. ( 29), are listed in Table I.In the 1D case, the results agree with the known exact values derived from the eigenvalues of the Pöschl-Teller potential [34].The calculations for two and three dimensions were previously carried out by Brézin and Parisi [32] by a direct computation of the spectrum of ∆ z .We agree with their results in 3D, but not in the 2D case.However, repeating their calculation in 2D, which also illustrates the comparative efficiency of the Gel'fand-Yaglom approach, we confirm the result derived in this Letter [47].
Taking everything together, our final result for the density of states is (taking the replica limit N → 0) The numerical results for these band tails are stated in Eq. ( 4).Again, our results agree with the explicit result in 1D [1].In 2D, we correct the values for the prefactor reported in Refs.[31,32], and in 3D we agree with Ref. [31] when adjusting for a factor e −I4/16π that accounts for a different renormalization choice (i.e., a shifted band origin).
In summary, we have quantified disorder effects on the electronic density of states in the universal Lifshitz regime including fluctuations, which set the prefactor of the Lifshitz tail.Using a replica path-integral approach, we derive the leading exponential form of the tail from an instanton saddle-point profile, which describes the most likely disorder potential that gives rise to a particular deep bound state.Importantly, we also obtain the fluctuation correction to this result, which sets the magnitude of this tail.Here, the main technical advance of our work is a generalization of the Gel'fand-Yaglom approach to evaluate the fluctuation determinant around the saddlepoint solution, which avoids a direct evaluation of the fluctuation spectrum.We thus provide an and efficient calculation of disorder tails that confirms results in 1D and 3D obtained using other means, and which corrects a result for the important case of two-dimensional systems.The derivation presented in this Letter has the dual advantage of simplicity and adaptability to more complex scenarios.For instance, extending our calculations to describe a broader class of band structures with, for example, fourth-order corrections, spin-orbit coupling, or band structure asymmetries, or to include more general disorder correlations are avenues for future work.
The saddle-point equation for the instanton configuration, Eq. ( 12) of the main text, is commonly solved using the shooting method.Here, we discuss the spectral renormalization method proposed by Ablowitz and Musslimani [S1] as a rapidly converging and robust alternative.It is based on the mapping between functions implied by the Fourier transform of ( 12) Fixed points of this map are solutions of Eq. ( 12).Unfortunately, a direct iteration of this equation does not usually converge to the instanton profile f (r), but to different and undesirable fixed points (e.g., zero or infinity).To improve the convergence, we follow [S1] and impose an additional integral condition characteristic of the fixed points of (S1), ) This is done by rescaling the function f n → γ n f n at each iteration, where the value of γ n is chosen to fulfill (S2) The rescaled map is given by fn+1 For general instanton profiles, the fixed-point map (S4) is evaluated using a fast Fourier transform, which converges to the desired solution within a small number of iterations.An additional simplification in our case comes from the observation that f (r) has rotational symmetry.This property is preserved by both the map (S4) and the Fourier transform.For rotationally symmetric functions, the Fourier transform simplifies to the one-dimensional integral where J is the Bessel function of the first kind.This integral is numerically evaluated using a Gauss-Legendre quadrature rule.The results shown in Fig. 1 and Tab.I of the main text are obtained using the spectral renormalization iteration (S4) with an initial Gaussian profile.The lowest eigenvalues (in any dimension) are 1/3 and 1, corresponding to the zero modes f (r) and f ′ (r) respectively.

II. DIRECT DIAGONALIZATION OF THE FLUCTUATION DETERMINANT
The results for the fluctuation determinant that we derive in the main paper using the Gel'fand-Yaglom method agree with the exact result in d = 1 and also with the result in d = 3 obtained by Brézin and Parisi [S2].However, Ref. [S2] also states results in d = 2, and our result differs from this work.To resolve this discrepancy and to confirm our result, this supplemental section presents an alternative derivation in d = 2 following [S2], which is based on an evaluation of the complete spectrum of the kernels of longitudinal and transverse fluctuations.The results obtained in this section agree with our calculation in the main text and correct the values reported in [S2].Furthermore, it highlights the use of the Gel'fand-Yaglom approach, which avoids the explicit computation of the fluctuation spectrum.
To evaluate the fluctuation determinant, we follow Brézin and Parisi and compute the generalized eigenvalues µ n (not to be confused with the λ n used in the main text) of the operator In d = 2, the eigenvalues depend on a radial angular momentum number l and are doubly degenerate for l ̸ = 0. Figure S1 shows all eigenvalues up to a cutoff of μ = 1000, computed using the shooting method.
In order to express D(z) in terms of the generalized eigenvalues µ n , we introduce the Weierstrass product where the numerical values are for d = 2.
In practice, we evaluate Eq. (S7) by explicitly computing a subset of eigenvalues, µ n ≤ μ, and estimating the contribution of the remaining large-µ terms using an asymptotic approximation for the spectrum of eigenvalues:

3 FIG. 2 .
FIG.2.Instanton solution of the saddle-point equation(12), which determines the leading scaling of the band tail B(E).
where z = 1, 1/3 for longitudinal and transversal fluctuations respectively, and the index n runs over all multiples of degenerate eigenvalues.The exponential factor in the the Weierstrass product is only present for d ≥ 2 and removes the divergent linear-in-z term.This procedure agrees with the counterterm renormalization applied in the main text.We then haveΩ(z) ≡ lim z→µ0 W (z) (1 − z/µ 0 ) d(µ0) = D(µ 0 ) lim z→µ0 det ∆ z (1 − z/µ 0 ) d(µ0) ,(S8)where d(µ 0 ) denotes the degeneracy of the zero mode with generalized eigenvalue µ 0 , and the term in brackets is to be evaluated in the subspace of zero modes.For µ 0 = 1, 1/3 the normalized zero modes were stated in the main text.

TABLE I .
(14)rical values for the moments of the instanton profile(14)and for the normalized longitudinal and transverse fluctuation determinant e 3zI 2 G 0 D(z) in d = 1, 2, and 3.