following form for two- and three-dimensional systems: ..... continuous variable. An alternative way to look at the SchrÃ¶dinger equation. (often usef...

2 downloads 14 Views 308KB Size

arXiv:cond-mat/9606181v1 25 Jun 1996

a

Lyman Laboratory of Physics, Harvard University, Cambridge, MA 02138, and b NEC Research Institute, 4 Independence Way, Princeton, NJ 08540 (April 13, 2017)

Here d is the dimensionality of the system and t may have the meaning, e.g., of the current relaxation time, or the normalized local density of states at the Fermi energy ρ (EF , r) /νd , or the normalized square of the local wave 2 function amplitude V |ψE (r)| of an eigenstate with energy E. νd is the average d-dimensional density of states and V is the volume of the sample. One notable excep2 tion is the e− ln law obtained in Ref. [7] for the distribution of wavefunction amplitudes in three-dimensional samples. We will discuss a possible origin of this difference in Section IV. The coefficients Cd in general depend on the strength of disorder and, for d = 2, on the sample size L. In the two-dimensional case the calculations based on σ-models βπ 2 ν2 D give C2 = [2,6–9], where pF is the Fermi mo2 ln Ll 1 mentum, D = lvF is the electron diffusion constant, l is 2 the mean free path (which is assumed to be much larger than the electron wavelength p−1 F ), and vF is the Fermi velocity. β is a numerical coefficient which, depending on the symmetry of the ensemble of random potentials takes the following values: β = 1 in time-reversal-invariant systems (the orthogonal symmetry class), and β = 2 when symmetry with respect to time reversal is completely broken (the unitary symmetry class). In Ref. [9] only the unitary ensemble was considered. We will argue below that this value of C2 can only be regarded as an orderof-magnitude estimate. In the one-dimensional case all eigenstates are localized, and the distribution of wave function amplitudes has a simple Gaussian form (see below, Eq. (31)). Current relaxation times t˜, on the other hand, are characterized by a much broader logarithmically normal distribution P (t) ∼ exp −C1 ln2 t˜ , (2)

Asymptotic behavior of the distribution functions of local quantities in disordered conductors is studied in the weak disorder limit by means of an optimal fluctuation method. It is argued that this method is more appropriate for the study of seldom occurring events than the approaches based on nonlinear σ-models because it is capable of correctly handling fluctuations of the random potential with large amplitude as well as the short-scale structure of the corresponding solutions of the Schr¨ odinger equation. For two- and three-dimensional conductors new asymptotics of the distribution functions are obtained which in some cases differ significantly from previously established results.

I. INTRODUCTION.

It has been well known for at least a decade [1] that an adequate description of the Anderson localization transition in disordered mesoscopic conductors must necessarily be formulated in terms of the full distribution functions of conductance g and/or other quantities characterizing the sample. As a consequence, even well into the metallic regime it should be possible to observe the onset of localization by studying deviations of asymptotic tails of the distribution functions from their behavior in the infinite conductance limit. Such a study was first performed in Ref. [2] within the framework of the replica-based nonlinear σ-model [3]. It was demonstrated that the tails of the distributions of such quantities as conductance, local density of states, current relaxation times, etc. are all described in two spatial dimensions by rather similar logarithmically normal asymptotes whereas Gaussian distributions would be expected for g → ∞. Recently, it was discovered [4] that these asymptotes can be obtained in a more straightforward and elegant manner using the optimal fluctuation method in conjunction with the supersymmetric variant of the nonlinear σ-model [5]. In a series of subsequent publications the method was employed to study statistics of eigenfunction amplitudes in weakly localized twodimensional conductors [6] and then extended to systems in three spatial dimensions [7–9]. The common feature of all the results obtained so far with the use of various σ-models is that asymptotic behavior of a distribution function P (t) for large t has the following form for two- and three-dimensional systems: P (t) ∼ exp −Cd lnd t . (1)

where C1 = l/2L for a sample of length 2L [9,11]. In Refs. [7] and [8] saddle-point solutions of the supersymmetric nonlinear σ-models were obtained for the three-dimensional case. However, conventional nonlinear σ-models used in Refs. [2,4,6–8] are low-“energy” effective field theories in which the role of energy is played by the diffusion operator D∇2 . As such they are only applicable to the description of phenomena that can be characterized as diffusive – the mean free path l is the smallest length scale in these theories. Formally this is expressed as a requirement that the scale of spatial vari-

1

with those obtained by employing nonlinear σ-models. A short summary and a list of open questions can be found in Section V.

ations of the fundamental variables of the theory – the Q-matrices – must be much larger than l. The optimal fluctuations of the Q-matrices computed in Refs. [7,8] were found to vary rapidly over distances ∼ l in threedimensional systems, which made it impossible to obtain rigorous results. The coefficient C3 was estimated in Ref. 2 [8] to be of the order of (pF l) . A similar calculation in Ref. [7] led to the same estimate for C3 but a log-normal 3 rather than e−C3 ln t functional form of the distribution function. In an attempt to overcome the limitations of diffusive σ-models and account exactly for the spatial variation of the optimal fluctuation in 3D systems at the scale of the mean free path, a generalized version of the model (ballistic nonlinear σ-model) was introduced by Muzykantskii and Khmelnitskii in Ref. [10] and used in Ref. [9] to compute the distribution of current relaxation times. The functional form of Eq. (1) was reproduced in that calcuπ 2 lation and a numerical value √ (pF l) of C3 was found. 9 3 We will argue based on the calculations presented below that all the results for the three-dimensional systems quoted above contain serious errors. In this article we propose an alternative method for the investigation of distribution function asymptotics. Instead of integrating out the disorder degrees of freedom and then looking for a saddle point of the resulting effective field theory (the nonlinear σ-model), we suggest that large-t behavior of the distribution functions is governed by a saddle point of the original theory in which the distribution functions are expressed as functional integrals over both electronic and disorder degrees of freedom. In other words, such a saddle point corresponds to an optimal fluctuation, i.e. the highest-probability configuration of disorder that lets the electronic eigenstates at a given energy E have the desired property – for example, an anomalously large intensity V |ψE (r)|2 at some point r. In our view, this approach possesses the advantage of being applicable to systems of arbitrary dimensionality as well as to disorder models more general than the commonly used model with a Gaussian distribution of the potential fluctuations. It is also conceptually simpler and makes the physical origin of the results much more transparent. In many respects our approach is similar to the ideas utilized in Refs. [12,13]. We were able to reproduce the general ln P ∼ − lnd t behavior of the distribution functions for d = 2, 3 as well as the log-normal or Gaussian (depending on the precise definition of t) form of P for d = 1. In two- and three-dimensional cases we obtain new values for the coefficients Cd correcting what we believe to be the errors in their previously published estimations. The rest of the paper is organized as follows. In Section II we briefly describe the method and present the main results. A detailed derivation of the basic equations and the analysis of the saddle point solutions is deferred until Section III. In Section IV we will compare our results

II. ASYMPTOTES OF THE DISTRIBUTION OF WAVE FUNCTION AMPLITUDES.

A typical wave function in a metallic sample is spread more or less uniformly throughout the sample volume V ∼ Ld so that its amplitude does not differ much from √ the average value of 1/ V anywhere in the sample. In the metallic (or weakly localized) regime such states account for the bulk of the distribution of the local density of states ρ (r) or of the current relaxation times t˜. The chances of observing anomalously large values of these quantities are related to the probability of finding an “anomalously localized state”. Such a state would be characterized by an amplitude reaching a value much larger than the average at some point r inside the sample. In what follows we will concentrate on the distribution of wave function amplitudes. Other distributions, such as those of current relaxation times, local density of states, etc. can be more or less straightforwardly derived within the framework of the same formalism. The problem can be formulated in the following way. Let us consider a spherical (in 3D) or a disk-shaped (in 2D) conductor of radius L and compute the probability that an eigenstate ψ of energy E (which we take p to be equal to the Fermi energy EF ) has an amplitude t/V in the center of the sample (r = 0), with t ≫ 1 [14]. The distribution of the disorder potential is assumed to be Gaussian, Z πνd τ W [U (r)] = NU exp − U 2 (r) dr , (3) 2 where NU is the normalization constant and τ is the mean free time. Such a probability is then naturally expressed as D E 2 P (t) = δ V |ψ (0)| − t , (4) U

where h. . .iU denotes the averaging with the weight W over all possible configurations of U , and ψ (r) is the solution of the Schr¨odinger equation with the Hamiltonian ˆ 2 /(2m) is a H = H0 + U (r) and energy E. H0 = p ˆ is the Hamiltonian of free particles with a mass m, and p canonical momentum operator. The requirement that E must be an eigenvalue of H was enforced in Ref. [6] explicitly by introducing the corresponding δ-function into the definition of P. In the alternative approach proposed here – the direct optimal fluctuation method – this requirement is easier to impose, if necessary, at a later stage in the calculations through appropriate boundary conditions for the saddle point equations [15].

2

function of the resonant state in such a potential is characterized by exponential decay in the interval a < r < b. Assuming that the amplitude of the √ wave function for V we can r > b is close to its average value 1/ n p o estimate t ≡ V ψ 2 (0) ∼ exp 2 2m (U0 − EF ) (b − a) . Antici-

Rewriting Eq. (4) as a constrained functional integral and introducing Lagrange multipliers to enforce the constraints (see Section III for details) one can demonstrate that the resulting “action” A possesses a saddle point, and the leading contribution to ln P (t) is given by −A0 , the value of the “action” at that point. With exponential accuracy the results are: P (t) ∼ exp −κ (pF l) ln3 t , (3D) , (5a) P (t) ∼ exp −π 2 ν2 D

2

ln t ln (L/r0 )

, (2D) ,

pating the final result b ∼ p−1 F ln t we can neglect a in the exponent. We then have 1 U0 ≈ EF + 2m

(5b)

sin 2pF r rd−1

ln t 2b

2

,

(7)

πνd τ Vd U02 bd , 2

(8)

and for a ≪ b A=

where κ is a number which is approximately equal to 3 × 10−3 . The distance r0 in Eq. (5b) is of the order of the electron wavelength p−1 F . In both two- and three-dimensional cases, our answers have a different l-dependence compared to the results obtained in Refs. [2,6–9]. In the three-dimensional case ln P is proportional to the first power of pF l while all the results derived from the nonlinear σ-model in three dimensions lead to a (pF l)2 dependence. Similarly, in the two-dimensional case the logarithm in the denominator contains L/r0 instead of the L/l established in Ref. [2] and reproduced with minor modifications in Refs. [6–9]. While in the two-dimensional case this difference leads only to small corrections, the answer for a threedimensional sample indicates a substantially increased probability to observe an anomalously large eigenstate amplitude as compared to the previous results. These asymptotes correspond to the optimal configurations of the potential U (r) which are essentially Bragg mirrors, U (r) ∝

πνd τ 2

Z

U 2 (r) dr =

π d/2 is the volume of a d-dimensional Γ (d/2 + 1) sphere of unit radius. Substituting Eq. (7) into Eq. (8) and minimizing with respect to b we find r 4−d −1 ln t, (9) b = pF 4d where Vd =

and it follows then that U0 = 4EF / (4 − d). The resulting estimate for the asymptotic behavior of the distribution function reads n o P (t) ∼ exp −Kd (pF l) lnd t , (10a) 1 Kd = 2d

1 4π

d−1

d 4−d

4−d 2

Vd2 .

(10b)

In the three-dimensional case the coefficient in front of 1 (pF l) ≈ 0.032 (pF l). It differs only by a ln3 t is 2 · 35/2 number from a more rigorous estimate quoted in Eq. (5b) suggesting that local resonances play an important role in the three-dimensional case. As it should be, the coefficient K3 ≈ 0.032 is larger than κ since it would be na¨ıve to expect a box-shaped potential to be optimal. In the two-dimensional case, large values of t are not optimally produced by local resonances: according to Eq. (5a) the coefficient in front of ln2 t is proportional −1 L to ln , i.e. it depends on the sample size. r0 The hypothesis that log-normal (in one and two di3 mensions) and e−C3 ln t (in 3D) asymptotics of the distribution functions may be produced by Bragg reflection locking the states in (or out of) certain regions of the sample was put forward by Muzykantskii and Khmelnitskii in Ref. [9]. They have explicitly demonstrated that a potential of the form U (x) = pF ln t˜∆ /mL cos (2pF x), 1 where ∆ = is the mean level spacing, traps an elecνd V tron of energy EF in a one-dimensional sample of length 2L for time t˜. The probability of finding such a configuration is

(6)

as r → ∞. The r → 0 behavior is more complicated. In one- and two-dimensional systems the r → 0 region does not contribute to ln P in the leading in ln t order. In three dimensions the shape of the optimal fluctuation at small distances is a narrow potential well surrounded by a wider potential barrier (Fig. 1) such that the whole structure supports a narrow resonance in the s-wave channel. Anomalously large values of t are achieved by combining this local resonance with Bragg reflection at larger distances. The role of local resonances in producing large values of t for any d in the interval 2 ≤ d < 4 can be illustrated by the following simplified version of the direct optimal fluctuation method. Let us approximate the shape of the resonance-producing potential as U (r) = U0 θ (r − a) θ (b − r), where θ (r) is the step function, and U0 (> EF ) and b are optimization parameters. The constant a is determined from the condition that there is only one resonance in the potential well and it has a −1/2 given energy EF (a ∼ p−1 ). The wave F = (2mEF ) 3

l P t˜ ∼ exp − ln2 t˜∆ , 2L

large. The stationarity condition for the action Eq. (14b) leads to the following set of equations:

(11)

which is the same as the log-normal long-time asymptotics of the exact solution for the conductance [11]. Our results confirm this hypothesis in higher-dimensional samples. The most striking feature of our results is the fact that a single configuration of the potential – i.e. a single compact region of the configurations space – seems to be responsible for the large-t tails of the distribution functions. The alternative would be for a large number of more or less equal-weight configurations to be rable to prot , and the duce the required wave function amplitude V total probability P (t) would then be given by the sum of the weights of these configurations. A strong argument against such an alternative is the fact that the results obtained by simply singling out the most probable type of configurations can lead to higher probabilities (significantly higher in the 3D case) than the σ-model results which correspond to a sum over all configurations [17].

πνd τ U (r) − iχ (r) ψ (r) = 0, ˆ2 p + U (r) − E ψ (r) = 0, 2m

(15b)

√ ˆ2 p + U (r) − E χ (r) + 2λ V tδ (r) = 0, 2m

(15c)

(15a)

2

V |ψ (0)| − t = 0.

(15d)

In deriving the saddle-point equations we neglect the U -dependence of the normalization constant Nψ . The reasoning leading to this approximation is explained in detail in the Appendix. Below we will investigate only those solutions of Eqs. (15) which possess full rotational symmetry. While it is likely that even for rotationally invariant boundary conditions there exist solutions which break rotational symmetry we will assume that they are not optimal, i.e. that they correspond to extremum points other than the global minimum of A. The mean free time enters these equations only through the πνd τ factor in Eq. (15a). It can be absorbed into redefinitions of χ (r) and λ, and therefore the above assumption does not depend on the disorder strength. Eliminating χ in favor of U and ψ and switching to dimensionless variables v (r) = U(|r|) EF , r = pF |r|, we arrive (in the absence of magnetic field) at the following system ′ d (.) ≡ (.) ): of equations ( dr y (r) = ln′ r(d−1)/2 ψ (r) (16a)

III. THE DIRECT OPTIMAL FLUCTUATION METHOD.

In order to enforce the condition that ψ in Eq. (4) is a solution of the Schrodinger equation we rewrite Eq. (4) formally as Z Z P (t) = DU (r) W [U (r)] Nψ [U ] Dψ (r) Y × δ (H − E) ψ (r′ ) δ V |ψ (0)|2 − t , (12) r′

where Nψ [U ] is a U -dependent normalization constant defined by Z Y Nψ [U ] Dψ (r) δ (H − E) ψ (r′ ) = 1. (13) r′

y ′ (r) + y 2 (r) +

Introducing auxiliary variables χ (r) and λ associated respectively with the first and the second δ-functions in Eq. (12) enables us to represent the distribution function as an unconstrained functional integral: Z Z Z χ (r) P (t) = NU DU (r) Nψ Dψ (r) D 2π Z dλ −A[U,ψ,χ,λ] × e , (14a) 2π

where

1 δd,2 + 1 = v (r) , 4r2

′ 2y (r) rd−1 v (r) − rd−1 v (r) = λ, λ=

2iλt pdF , πνd τ EF2 Sd

(16b)

(16c)

(17)

δd,2 is the Kronecker symbol and Sd = dVd is the area of a d-dimensional sphere of unit radius. Equations (16a) and (16b) together are equivalent to the Schr¨odinger equation for the radial part of the wave function, while Eq. (16c) provides the nontrivial connection between the wave function and the potential that is necessary to achieve a given value of t at a minimal “cost” in terms of the weight function W. λ can be viewed as a parameter

Z πνd τ 2 dr U 2 (r) − iλ V |ψ (0)| − t A [U, ψ, χ, λ] = 2 2 Z ˆ p −i dr χ (r) + U (r) − E ψ (r) . (14b) 2m The next step is to make a saddle-point approximation in Eq. (14a) utilizing the fact that t is assumed to be 4

effecting a Legendre transformation of A, so (λ, t) – alternatively, λ, ln t – form a pair of conjugate variables. It follows from Eq. (16c) that λ must be real, so the saddle-point value of λ is purely imaginary. The saddlepoint action A0 is expressed in terms of the solution v (r) of Eqs. (16) as Z L πντ EF2 Sd A0 = rd−1 v 2 (r) dr. (18) 2 pdF 0

The normalization requirement Z L Sd rd−1 drψ 2 (r) = 1

provides the second condition. Another boundary condition is derived from regularity requirements on v (r). Since v (r) can have at most a square integrable singularity as r → 0, it follows that in three dimensions √ ′ (21c) lim rψ (r) = 0.

The choice of boundary conditions to Eqs. (16) is not entirely straightforward, and we would like to motivate it by appealing to an analogy with ordinary onedimensional Schr¨odinger equation. In the variational formulation of quantum mechanics the Schr¨odinger equation appears as a stationarity condition of a functional Z L ˆ − E ψ (x) , dxψ ∗ (x) H F= (19)

r→0

In two dimensions we simply have a condition that ψ ′ must be finite at r = 0, and for a one-dimensional system the corresponding requirement is ψ ′ (0) = 0. Finally, a fourth condition may be provided by another boundary condition on ψ, e.g.

0

ψ (L) = 0.

ˆ is the Hamiltonian, E is the particle enwhere H ergy, and ψ must satisfy two boundary conditions (e.g., ψ (0) = ψ (L) = 0) and also the normalization condiRL tion 0 dx|ψ|2 = 1. The overall number of conditions for this second-order equation is three, which makes it overdefined. As a result, solutions exist only for a discrete set of values of E. One of these values would correspond to the ground state and thus to an absolute minimum of F , while the rest would simply correspond to its stationary points. In the L → ∞ limit the set of allowed E values becomes infinitely dense so that E effectively becomes a continuous variable. An alternative way to look at the Schr¨odinger equation (often useful in numerical calculations [19]) is to consider E as an extra unknown function satisfying a trivial equation d E = 0. dx

ψ (0) =

t . V

(21d)

However the problem of minimizing A is not solved uniquely by the set of equations (16) together with the conditions (21). As in the ordinary quantum mechanical case, we are faced with a family of solutions, only one of which actually corresponds to the absolute minimum A0 of the action A. Finding this “ground state” solution requires a parametrization of the whole discrete family of solutions which presents a difficult task. It can nevertheless be somewhat simplified by employing the following observation. For large enough L the set of solutions of Eqs. (16) may be approximately considered continuous. It is then possible to replace the condition (21d) with a free boundary condition. The result will be that the family of allowed solutions will now admit a continuous parametrization, simplifying the task of determining the minimum of the functional A. In practice this parametrization depends on the dimensionality of the sample, and we will consider the cases of d = 1, 2 and 3 separately. Employing a continuous parametrization which effectively corresponds to relaxing the boundary condition Eq. (21d) means not enforcing the condition that EF is an eigenvalue of the Hamiltonian. As was already mentioned [15] such an approximation leads only to errors that are beyond the exponential accuracy with which asymptotes of P (t) are determined in the optimal fluctuation method and which can therefore be safely neglected.

(20)

One then has effectively three first-order equations with three extra conditions imposed on them. However, now these equations are nonlinear, and the three conditions do not specify the solution uniquely. Rather, there exists a family of solutions corresponding to various quantum states of the system, one of which minimizes the functional F . Again, in the L → ∞ limit the family of solutions becomes continuous. Returning to the problem of specifying a set of conditions for Eqs. (16) we see that equation (16c) is a generalization of Eq. (20). λ plays the role of an eigenvalue, so we can formally supplement Eqs. (16) with λ ′ = 0, increasing the number of equations to four. The following set of conditions must be imposed on this system of equations. First, a boundary condition is provided by Eq. (15d): r

(21b)

0

A. One-dimensional wire.

As a “toy” model, we consider first a purely onedimensional disordered wire of length 2L. All the eigenstates in this case are exponentially localized with a localization length being of the order of the mean free path l. A typical amplitude of √ such a state in its domain of √ localization is 1/ l ≫ 1/ 2L, so the optimal fluctuation method cannot be applied to the computation of

(21a) 5

the distribution function P (t) for t < ∼ L/l. Amplitudes larger than that are due to the states that are characterized by an anomalously small localization length and the probability of finding such states is determined by the probability of the corresponding optimal fluctuation. The system of Eqs. (16) takes the form v (x) = 1 + y 2 (x) + y ′ (x) , 1 2y (x) v (x) − v ′ (x) = λ Sgn x, 2

A = πν1 τ

′ 2

2

(y ) = y + 1

+ λ y + C.

(22)

From Eq. (25) and the fact that the change in the wave function amplitude over one period (increase for x < 0 or decrease for x > 0) does not depend on x, we deduce the exponential form of the wave function envelope ψ (x) ∝ e−

(23)

C

2|x| L

ln θ

.

(30)

Normalizing the wave function to 1 and using V = 2L we obtain t = 2Lψ 2 (0) = 2 (ln θ). We thus finally arrive at l 2 P (t) ∼ exp − t . (31) 2L

Here C is a constant of integration which parametrizes the stationary points of A as described at the end of the preceding subsection. Its value is fixed by a minimization procedure which is formally expressed as A0 = min A (C)

(28)

It is thus obvious that, at least to the lowest order in 1/L, the minimum of A is achieved when C = 0. Combining the results together we have for the asymptotic tail of the distribution function with exponential accuracy 2l (29) P (θ) ∼ exp − ln2 θ . L

where x is the dimensionless coordinate along the wire, and primes denote differentiation with respect to x. Note that when d = 1, y is simply the logarithmic derivative of the wave function. Eliminating v (x) and integrating once, we arrive (for x > 0) at 2

EF2 2πL 2 λ + 3C 2 . pF 32T

(24)

As explained at the beginning of this subsection the applicability of the above formula is restricted to the region L/l < t < pF L (the second inequality ensures the validity of the expansion in λ and C). As required, the corresponding value of the localization length L/ ln θ is much smaller than l. This model calculation can be used to illustrate the following points. First, to the lowest order in λ we have

Although Eq. (23) is exactly integrable in terms of elliptic functions, it is sufficient to make a perturbative expansion in λ and C in order to recover the leading order dependence on t. Rewriting Eq. (23) as q dy , R (y) = (y 2 + 1)2 + λ y + C, (25) dx = − R (y) we obtain the period T of the wave function oscillations Z ∞ dy T =2 (26) R (y) −∞

y (x) = cot (x + ϕ) +

λ sin2 (x + ϕ) , 4

(32)

where ϕ is a phase shift, and v (x) =

and the logarithm of the ratio of the wave function amplitudes at x = 0 and x = L: Z ψ (0) L ∞ ydy ln =2 . (27) ψ (L) T −∞ R (y)

1 λ sin 2 (x + ϕ) , 2

(33)

so the shape of the optimal configuration of the potential is indeed a Bragg mirror. Second, up to a constant, θ2 /∆ can be identified with t˜ – the delay time in the electric response of an open sample at EF [9], and the distribution in Eq. (29) is seen to coincide with the exact answer of Altshuler and Prigodin [11] quoted at the end of Section II. Third, the sign of λ is opposite to that of ln θ, so negative values of λ correspond to the wave function amplitude decreasing from the center of the sample outwards. The same will hold true in two and three dimensions. The amplitude of the oscillating potential is constant throughout the sample. A similar observation was made in Ref. [9] – in a one-dimensional wire the gradients of the supersymmetric density matrix corresponding to the optimal solution were found to be x-independent.

We assume that the length of the sample corresponds to an integer number of periods, so that if y (0) = 0 (from symmetry arguments) it follows also that y (−L) = y (L) = 0. The relative error introduced into Eq. (27) by neglecting the fractional part of L/T is O (1/L). It is convenient to use the logarithm in the l.h.s. of Eq. (27), which we will denote as ln θ, as a parameter of the distribution function. It can be easily seen that the contribution to P of the Jacobian of the transformation from θ to t is only a prefactor which does not change the leading order term in the exponent. Expanding the r.h.s. of Eqs. (26) and (27) in both λ and C we obtain −πλ/16 = (T /2L) ln θ and T = 2π + O (C). For the action A we then find to the lowest order (with Sd = 2) 6

√ where we have approximated r0 ψ (r0 ) ≈ ψ (0). For r0 ∼ 1 the error introduced by this approximation leads only to O (1) corrections to the large logarithm ln t. The normalization integral

B. Two-dimensional conductors

In dimensions higher than one the system (16) cannot be integrated exactly in terms of the standard elementary or special functions. Nevertheless, numerical methods combined with the asymptotic analysis allow us to fully investigate the behavior of its solutions in various regimes. In 2 dimensions the solution can be obtained as an asymptotic expansion in λ/r. The values of t for which this expansion breaks down turn out to lie close to the limiting value πL2 , and are thus not very interesting. The leading terms of the expansion of v (r) are v (r) ≈

η λ sin (2r + 2ϕ) + sin2 (r + ϕ) . 2r r

2

N = 2πψ (0)

λ sin2 (r + ϕ) . 4r

(34)

r 2 ln θ ln rL 0

r

0

(40)

This result has essentially the same form as the one obtained by Falko and Efetov [7] using the σ-model formalism. It differs, however, from their result in two important aspects. First, the logarithm in the denominator is cut off at distances ∼ 1 (p−1 F in conventional units) rather than at r ∼ l, as in Ref. [7]. We believe that this difference simply reflects the fact that both the approach employed here and the saddle-point solution of the σmodel represent the same saddle point of the underlying theory. However, in order to describe this saddle point exactly one has to include short-wavelength degrees of freedom which are lost in the derivation leading to the σ-model, and that accounts for its inability to obtain the correct cut-off scale. The second difference may turn out to be an indication of a deeper problem with the σ-model. Our calculation is performed explicitly for the case when there are no perturbations breaking the symmetry with respect to time reversal. Nevertheless, the numerical coefficient in the exponent coincides with the answer obtained in Ref. [7] for the unitary case (β = 2), when the symmetry with respect to time reversal is completely broken. In order to see what effect violation of time-reversal invariance would have in our approach, we explicitly introduce magnetic field into the system of equations (15). Choosing the direction of the field H along the zˆ axis (perpendicular to the two-dimensional sample) and writing A = 21 H × r we can immediately see that under the assumption of circular symmetry of the solution ψ, only the terms quadratic in A survive. For weak enough fields HL2 ≪ φ0 (where φ0 is the magnetic flux quantum) these terms can be neglected because they only lead to exponential decay of the wave function on a scale larger than the sample size. Thus the set of equations (16) is unchanged. On the other hand, the fields of this magnitude are sufficient to suppress the Cooperon contribution in the σ-model calculation leading to a crossover between the orthog-

(35)

(36)

√ in order to cancel out the overall 1/ r dependence of the wave function. θ is always large in the asymptotic region t ≫ 1. Integrating y (r) in the sense of the principal value we obtain Z L L |λ| + ln θ1 , (37) ln y (r) dr + ln θ1 ≈ ln θ = −P 8 r 0 r0 where the length scale r0 is defined as r0 ∼ max 1, |λ| , and ln θ1 represents the contribution to the integral from < distances r < ∼ r0 . When |λ| ∼ 1 this contribution can be neglected with logarithmic accuracy. Inverting Eq. L (37) we obtain |λ| = 8 ln θ/ ln . The dimensionless inr0 R L tegral in the saddle point action 0 v 2 rdr evaluates to L 8 ln2 θ/ ln , and we find r0 ln2 θ P (θ) ∼ exp −4π 2 ν2 D . (38) ln (L/r0 ) The envelope of the wave function corresponding to the solution described by Eq. (34) is ψ (0) ln θ/ ln rL 0 , ψ (r) ∼ √ (r0 /r) r

dr

has two distinct regimes. When θ2 is greater than L/r0 the integral is dominated by small distances, and the dependence of t on θ is weak. We will not be considering this regime here because it corresponds to wave function amplitudes that are close to the limiting value t = πL2 . In the opposite case we have N = 2πψ 2 (0) r0α L1−α / (1 − α), where α = ln θ2 / ln rL0 . Up to irrelevant constants we then have ln t ≈ ln θ2 , and the distribution function becomes ln2 t . (41) P (t) ∼ exp −πν2 D ln (L/r0 )

It is convenient to define the parameter θ as ψ (0) θ= √ Lψ (L)

L

r0

In this expression η is a constant of integration which plays a role analogous to the role of C in the onedimensional case – that of a minimization parameter. Just as in the one-dimensional case, the minimum of A is reached when η = 0. We will give a perturbative proof of this statement at the end of this subsection. With η = 0 the asymptotic expansion of y (r) has the form y (r) ≈ cot (r + ϕ) +

Z

(39)

7

onal and unitary symmetry classes. Therefore our results reproduce exactly the σ-model answer in the case of broken time-reversal invariance (the unitary symmetry class), but they predict a faster decay of the distribution function at large t in the orthogonal case. If the fault lies with our calculation, it would then mean that while in the unitary case a single optimal fluctuation of the potential is responsible for the high-amplitude tail of the distribution of the wave function intensities, an exponentially large number of non-circularly-symmetric configurations would have to contribute to the probability of finding large wave function amplitudes in the orthogonal case. (It is certainly very unlikely that accounting for the fluctuations around the saddle point defined by Eq. (15) can change the coefficient of a ln2 t term in the exponential). While this possibility cannot be completely discounted without a more detailed study of the system Eqs. (15), we would like to propose an alternative explanation. The expansion into the sum over quasiclassical paths that underlies the nonlinear σ-model, and the subsequent separation into the diffuson and Cooperon contributions implies that self-retracing paths constitute a set of measure zero in the overall sum. This is certainly a rigorous assumption when there is no “intrinsic” mechanism in the calculation that would strongly favor some types of quasiclassical paths over others – as is almost always the case in clean chaotic systems where the zero mode treatment of the σ-model is applicable. However, when a saddle point approximation is being made in order to compute the probability of some rarely occurring event, one effectively selects certain members of the ensemble of random potentials that are favorable for such an event. Within the σ-model formalism this selection occurrs as a restriction on the possible configurations of the effective order parameter (the Q-matrix, or the supersymmetric density matrix g in the formalism of Muzykantskii and Khmelnitskii [10,9]), and the information about the actual configurations of the potential that are being selected is lost. On the other hand, as indicated by the calculation presented here, these configurations are rather special. Indeed, one can interprete the appearance of a quasilocalized state possessing an anomalously large amplitude at some point r as an opening of a quasi-gap in the spectrum of s-waves with respect to r. In a weak periodic potential emergence of a gap can already be seen in the second order of the perturbation theory, and the second order terms in the expansion into a sum over quasiclassical paths by necessity contain only self-retracing trajectories. Such contributions are then overcounted by the σ-model when it takes into account the same self-retracing path twice – as a diffuson and as a Cooperon. While the above argument is rather speculative, it nevertheless presents a picture which, if confirmed, would restrict the applicability of saddle-point solutions of the nonlinear σ-model. We believe that this issue deserves further consideration. To complete the derivation of the distribution function presented in this subsection, we outline the proof

of the statement that η = 0 is the optimal choice. The solution ψ of equations (15) can be written as ψ (r) = A(r) √ sin (r + ϕ (r)). Unlike the η = 0 case, the phase ϕ r does not have a finite limit as r → ∞. The logarithmic derivative y˜ (r) of the amplitude function A (r) can be expanded in an asymptotic Fourier series of the form ( ∞ h ∞ X X 1 (c) y(n,m) cos m (r + ϕ (r)) y + y˜ (r) ∼ (n,0) n r m=1 n=1 io (s) (42) + y(n,m) sin m (r + ϕ (r)) ,

and ϕ′ (r) can be represented in a similar way. Substituting these expansions into Eqs. (15) we find that y(1,0) = λ/8. Thus the choice of the constant of integration η does not affect the relation between y(1,0) and λ. Integrating y (r) over r we obtain y(1,0) ln rL0 and therefore the relation between ln θ and λ established in Eq. (37) is also η-independent to the leading order in λ. On the other hand, adding the term ηr sin2 (r + ϕ) to v can RL only increase the value of the integral 0 v 2 rdr because the cross-term in the expansion of the square integrates to zero. It follows then that the minimum of the integral for a given value of θ is achieved by setting η = 0. Note, however, that, as in the one-dimensional case, this proof is perturbative – it relies on the possibility to expand the solutions in powers of λ/r. C. The three-dimensional case.

In the three-dimensional case, it can be self-consistently demonstrated that large values of t correspond to large negative values of λ. As a result, there exists a range of values of r where expansion in λ/rd−1 = λ/r2 is impossible. A typical solution which was obtained numerically using the so-called relaxation method [19] is shown in Fig. 1. One can distinguish three asymptotic regimes: (i) r ≪ r0 , (ii) r0 ≪ r ≪ r∗ , and (iii) q r ≫ r∗ ;

it will be shown below that r0 = 1/|λ| and r∗ = |λ|/2. The first region corresponds to a potential well, the second one – to a potential barrier. Taken together, these two regions support a resonance in the s-wave channel at the energy EF . However, besides a potential-well – potential-barrier combination, resonant scattering can be also caused by a weak periodic potential (Bragg reflection), and that is exactly what the third region corresponds to. An interesting consequence of the solution presented in Fig. 1 is that the optimal way to achieve large values of t in three dimensions is to combine the two effects in the “right” proportion. Analytically the three regions in Fig. 1 are described by the following asymptotic formulae. In the first region, v (r) diverges as 7 2 2 (43) v1 (r) ∼ λ/r + λ ln r + 1 + c + λ , 12 8

where c is an arbitrary constant which cannot be determined from the boundary conditions Eq. (21). Note that the singularity is weak enough so that it produces only a finite contribution to the saddle-point action A0 . The behavior of the wave function in this regime is given by y (r) ≈

1 2 1 1 + λ + λ r ln r + c r, r 2 3

r∗ ≈

(44)

.

The expansion in Eq. (43) breaks down for r ∼ 1/|λ| which gives the approximate value for r0 . In the second region the solution can be obtained with the help of the quasiclassical approximation, y 2 (r) ≈ v (r), and the result is v2 (r) ≈

|λ| 2r2

2/3

.

(47)

As c decreases further, η monotonically decreases as well, eventually covering the whole (+∞, −∞) interval. When η is positive or small negative (η > ∼ −|λ|) the behavior of v (r) and y (r) in the regions (i) and (ii) depends on η (or c) only very weakly so that this dependence can be ignored in computing the contribution of these regions to A. Larger negative values of η start affecting the length of region (ii), and may lead to an emergence of a hybrid regime in which v oscillates non-harmonically with an r−4/3 envelope. In what follows we will assume that similarly to the two-dimensional case the minimal value A0 is achieved by setting η = 0, even though the breakdown of the asympλ totic expansion in 2 at small distances makes it imposr sible to construct an analytical proof of this statement analogous to the proofs for one- and two-dimensional systems. This assumption is borne out by numerical analysis. It is certainly obvious that positive values of η can never be optimal because of a corresponding “costly” v ∼ 1 region. As for large negative values, they are ruled out by the fact that a pronounced “hybrid” regime never appears in numerical q solutions. Setting η = 0 reduces

which corresponds to ψ (r) ≈ ψ (0) 1 + λr/2 + O r2 ln r

1/4 1 2 λ + η2 . 2

(45)

Finally, in the third region, where the λ/r2 expansion works, v (r) has the asymptotic form λ sin (2r + 2ϕ) sin2 (r + ϕ) 1 v3 (r) ∼ . + η + O 2 r2 r2 r3

Eq. (47) to r∗ ≈ |λ|/2. With η ≈ 0, both v2 (r) and v3 (r) reach values ∼ 1 at r = r∗ . While approximation schemes devised for v ≫ 1 and v ≪ 1 break down around r∗ , we can calculate ψ(0) = θ1 θ2 from regions contributions θ1 and θ2 to θ = Lψ(L) ∗ ∗ r < r (i-ii) and r > r (iii) separately. Similar to the two-dimensional case, the L factor in the denominator of θ is introduced to cancel the overall 1/r dependence of the wave functions, which is an artefact of spherically symmetric boundary conditions. (i) - (ii) We combine together the first and the second regions because the wave function amplitude does not change appreciably in the very short region (i). Region (ii) corresponds to stretched exponential decay of p ψ. Integrating y (r) = − v (r) over r from 1/|λ| to 2 q

(46) The constants η and ϕ in this expression are analogous to their counterparts in the two-dimensional case. The phase variable ϕ again has the meaning of the wave function phase shift; it is finite due to the rapid decay of the potential at large r irrespective of the value of η. Constants η and ϕ are not independent: they both can be regarded as functions of c. Either η or c can be chosen as a minimization parameter. We have not been able to establish the analytical dependence η (c); however, based on the numerical analysis the qualitative features of this dependence can be described as follows. For a given λ there exists a critical value c0 λ such that if c > c0 , the third (oscillatory) region never develops. Instead, v (r) exhibits a singularity at some finite value of r, leading to a divergent integral in A0 . So the values c > c0 correspond to unphysical solutions. Exactly at the critical value c0 the oscillations are 2 also absent, and v (r → ∞) ≈ 1 + λ /(4r4 ). For slightly smaller c, η (c) is large and positive (≫ |λ|), and the oscillations appear only after a more or less protracted intermidiate regime in which v (r) ≈ 1. An estimate of the point r∗ at which the onset of the oscillatory behavior occurrs can be obtained by noting that the oscillations of v (r) are driven (through Eq. (16)) by oscillations of the wave function. Therefore the amplitude of the oscillations in Eq. (46) cannot significantly exceed 1 so as to preserve the oscillatory – rather than exponentially damped – behavior of the wave function. This requirement leads to

r∗ ≈

|λ|/2 we obtain ln θ1 ≈

Z

√

|λ|/2

1/|λ|

dr

|λ| 2r2

1/3

≈ 3r∗ .

(48)

(iii) Using Eq. (46) with η = 0 we find y (r) ∼ cot (r + ϕ) +

λ sin2 (r + ϕ) . 4r2

(49)

The first term in this expression stems from the oscillations of the wave function, while the second one describes the decrease (λ < 0) of the wave function envelope from the center of the sample outwards. The second contribuRtion ln θ2 is given by the principal value of the integral ydr: 9

Z L ln θ2 ≈ −P √

|λ|/2

y (r) dr ≈ r∗ /4.

Of course, the separation into regions r < r∗ and r > r∗ is approximate. It is possible that the crossover region r ∼ r∗ gives a contribution of the same order of magnitude as the two asymptotic regions (ii) and (iii). Thus the number 7/2197 ≈ 3.2 × 10−3 can only be considered as an order-of-magnitude estimate of the coeficient κ introduced in Section II. It must be mentioned that convergence to the asymptotic form of Eq. (56) is extremely slow because of the stringency of the requirement 4 ln θ ≫ 1. r∗ = 13 In the opposite regime L ≪ θ2 , ψ 2 (0) is close to its maximal value ∼ 1, and therefore its θ dependence is very weak. On the other hand, similarly to the onedimensional case, θ2 becomes proportional to the electric response time t˜, leading to the distribution of t˜ having the form identical to Eq. (56), P t˜ ∼ exp −κ (pF l) ln3 t˜∆ . (57)

(50)

Adding the two contributions we establish the relation between λ and θ: s 13 ∗ 13 |λ| ln θ = r = , (51) 4 4 2 verifying the self-consistency of the assumption |λ| ≫ 1 made at the beginning of this subsection. The dimensionless integral in the saddle point action is also evaluated separately in the regions (i)-(ii) and (iii): r∗

Z

r2 v 2 dr = 3r∗3 =

0

192 3 ln θ 133

(regions (i) and(ii)) (52a)

and Z

L

r2 v 2 dr = r∗3 /2 =

r∗

32 3 ln θ 133

Both Eq. (56) and Eq. (57) differ significantly from the corresponding σ-model results [7,9]. We discuss the possible origins of this difference in the next section.

(region (iii)). (52b) IV. DISCUSSION

Both contributions turn out to be proportional to the same power of ln θ indicating that local resonances and Bragg reflection play equally important roles in the formation of anomalously large wave function intensities. Combining the results we obtain 56 P (θ) ∼ exp − pF l ln3 θ . (53) 2197

A. General remarks.

The main result of the work presented in this paper is that statistics of seldom occurring events in disordered electronic systems can be successfully studied using a variant of the optimal fluctuation method. All previous approaches to the problem have been based on various formulations of the nonlinear σ-model, and they invariably seem to require an extension of the σ-model to, and sometimes beyond, its limits of validity. In retrospect, this situation seems rather natural. Indeed, the success of the σ-model in describing a wide variety of phenomena in chaotic and disordered systems can be traced to the fact that most such phenomena are quasiclassic in nature and therefore their most appropriate description is in terms of an expansion into sums over quasiclassical trajectories. Calculation schemes based on non-linear σ-models allow one to perform such summations in a very efficient manner by introducing collective variables (Q-matrices) that change slowly in space [3,5]. When low-order moments and correlation functions are described the integrals over Q-matrices get almost equal-weight contributions from the whole degenerate (or near-degenerate) manifold [3,5] of allowed Q-matrix configurations. This effectively corresponds to probing a large number of disorder configurations. An investigation of the statistics of rare events, on the other hand, presents quite a different type of problem. The configurations of the random potential that give rise to such events come from a small subset of all possible configurations. This circumstance was reflected in the approach developed in Ref. [4] – the large parameter t

The envelopes of the wave function in the two regions are ψ (r) ∼

1 ψ (0) exp{−3r∗2/3 r1/3 } (for ∗2 < r < r∗ ) r r (54a)

and ψ (r) ∼

ψ (0) exp{r∗2 /4r − (13/4) r∗ } (for r > r∗ ). r (54b)

The normalization integral is given by 1 r∗ − 13 2 2 . + 4Le N = πψ (0) 3r∗2

(55)

When L ≫ θ2 (with more realistic boundary conditions this inequality will become L3 ≫ θ2 ), the normalization integral is dominated by the contribution from region (iii). Then ψ 2 (0) ∝ θ2 , and we finally obtain with exponential accuracy 7 3 P (t) ∼ exp − (pF l) ln t . (56) 2197 10

dence between the saddle points of the theory defined by Eqs. (14) and the saddle points of the σ-model. Assuming the existence of such a correspondence we will try to elucidate the origins of the σ-model results by analyzing the corresponding optimal configurations of the potential U (r). Before proceeding with this analysis we would like briefly to address the question of the stability of the saddle point described by Eqs. (15). In the threedimensional case the dominant contribution to the saddle point action A0 comes from small distances where the optimal potential is much larger than its typical Gaussian fluctuation, and stability with respect to fluctuations does not pose a problem. In two dimensions, however, the outlying regions of the optimal configuration – where the magnitude of the potential vanishes as 1/r – must be taken into account. In order to argue that fluctuations do not destroy the saddle point we note that although we find it convenient to write the Gaussian distribution function W [U ] for the potential in the coordinate representation, it can be written in any orthonormal basis {fn (r)}. A typical amplitude of a dimensionless√basis function fn in a typical configuration U (r) is 1/ ν2 D. Let us now choose one of the basis functions, say f0 (r), to be proportional to the optimal solution v (r) R given by Eq. (34). Using the normalization condition drf02 (r) = 1 we obtain s sin 2r 1 , (58a) f0 (r) ≈ π ln (L/r0 ) r

served to strongly break the degeneracy of the realizations of the Q-matrices, and only the vicinity of a saddle point of the integral over Q’s was shown to give a significant contribution to the probability of observing large t. What was effectively encoded in such a scheme was the fact that the purpose of the calculation was to select those special potentials that are capable of producing a given value of t, and a corresponding selection of the relevant configurations of Q’s was a computational tool allowing to do that. The disadvantages of such an approach are obvious: (i) First, the assumption that the motion of electrons can be described entirely in quasiclassical terms imposes certain restrictions on the types of the potentials over which averaging is performed. As follows from the calculations presented here, in three spatial dimensions there exists a region r ∼ r∗ where the solution of the Schr¨odinger equation cannot be obtained in the quasiclassical approximation. As a result, the σ-model fails to recognize the existence of this region and it also misses entirely the local resonance formed at r < r∗ . (ii) Second, even when the σ-model calculations succeed in correctly – albeit implicitly – identifying the relevant disorder configurations (as in one- and twodimensional samples) they do not always produce correct answers because contributions from the short-wavelength degrees of freedom (massive modes) eliminated in the transformation from the fast to slow variables are missed. B. The direct optimal fluctuation method.

and It is not immediately evident, however, that a relatively na¨ıve approach based on the direct search for an optimal fluctuation should be more reliable. In order for it to work, the probability of observing a large value of t must be determined by a sum over disorder configurations which all come from a single compact region of the configurations space. The potentials forming this region differ only slightly from some optimal configuration which corresponds to the saddle point. Although we have not proved in this work that the saddle point identified here gives the dominant contribution to the functional integral, “the preponderance of evidence” based on the comparison of our results with those obtained using the σ-model would indicate that this is indeed the case. The starting point of our qualitative analysis is the observation that, apart from the difference in the cutoff scale, our variant of the optimal fluctuation method (the “direct” optimal fluctuation method) reproduces identically in the two-dimensional and one-dimensional cases the σ-model results for the unitary ensemble. Assuming that more than a chance coincidence is involved, it is reasonable to conclude that the optimal configurations of disorder found in this work are the same as the ones that are responsible for the saddle point of the σ-model. In other words, there must exist a one-to-one correspon-

ln θ v (r) ∼ p f0 (r) . ln (L/r0 )

(58b)

Thus the optimal fluctuation has a much larger ampliln (L/r0 ) . tude than the typical one as long as ln2 θ ≫ ν2 D This condition, of course, is just a natural requirement for the validity of the optimal fluctuation method | ln P| ≫ 1. It remains to be shown, however, that other components of a typical fluctuation U (r) (i.e. those orthogonal to f0 ) do not destroy the saddle point. A rigorous investigation of the fluctuations around the saddle point defined by Eqs. (15) will be the subject of a forthcoming publication. Nevertheless, a plausible argument in favor of the stability of this saddle point can be made based on the following observation. Since the appearance of anomalously localized states due to Bragg reflection is a phenomenon analogous to the emergence of a band structure in a periodic lattice, the suppression of this effect by fluctuations of the potential is equivalent to the localization transition which destroys the band structure in ordinary periodic lattices. Therefore L < Lc seems to be a sufficient condition for the stability of the saddle-point solution in the two-dimensional case.

11

tween “global” and local fluctuations of the potential gets blurred. It is interesting to note that a calculation based on the ballistic σ-model performed in Ref. [9] for quasione-dimensional conductors indicates the existence of a crossover from log-normal to power-law distribution at L ln t˜∆ ∼ . Although we have not investigated the quasil one-dimensional case here, it is likely that the physical picture of the interplay between the “global” and local fluctuations discussed above should not be much different. If that is the case, then it is probable that the crossover found in Ref. [9] has as its underlying cause the same mechanism of local fluctuations becoming comparable in size to the length of the sample. This hypothesis, however, leaves unexplained the difference in scales L at which the crossover occurrs – in Ref. [9] as opposed l to pF L in the argument presented above. It is possible that the quasi-one-dimensional case brings in some new features that are not recognized by the estimates based on the purely one-dimensional model. It should be mentioned, however, that none of the variants of the σmodel can provide an adequate description of the effects associated with local resonances, and it is possible that a σ-model estimation of the crossover scale may not be entirely reliable.

C. The asymptotics of the distribution functions. 1. The one-dimensional case.

We will now try to use the physical intuition afforded by the optimal fluctuation concept to compare the σmodel results of Eqs. (1) and (2) with the distribution functions Eqs. (29), (41), (56) and (57) derived in the preceding Section. The first question that can be easily answered is why the distribution function for the electric response times in the one-dimensional case is log-normal instead of a power law (i.e., exp−C1 ln t˜ ) obtained by a na¨ıve extrapolation of the exp −Cd lnd t dependence to d = 1. In one and two dimensions the optimal configurations found in the preceding Section are “global”, i.e. R the integration in the saddle-point action v 2 dr must be extended over the whole sample. In contrast, in the three-dimensional case the optimal fluctuation, even with the oscillating tail included, is local, so that the above integral converges at large distances. It is well known [21] d that distribution function tails of the e− ln t˜ type usually appear as probabilities of optimal fluctuations confined to a finite volume. Thus in order to explain the “anomaly” in the one-dimensional case we have to understand why a local fluctuation of the random potential necessary to achieve a given value of t˜ has a lower probability than the ”global” one proposed in Ref. [9] and rederived in Section III A. A local fluctuation of the potential leading to a large value of t˜ would have to be able to support a narrow resonance. Assuming as in Section II a rectangular shape for the potential barrier responsible for the formation of the resonance we can repeat the calculation presented there almost verbatim except that t must everywhere be replaced by t˜∆. We then obtain 2 ln P ≈ √ (pF l) ln t˜∆. In order for this behavior to 3 3 dominate we must have pF l ln t˜∆ ≪ (l/L) ln2 t˜∆, or ln t˜∆ ≫ pF L.

2. Two-dimensional conductors.

A rather comprehensive discussion of the results obtained by the direct optimal fluctuation method in the two-dimensional case and their counterparts established using the σ-model formalism has already been presented in Section III B. Here it seems appropriate to briefly reiterate the following two main points of that discussion. First, it comes as no surprise that when massive modes are taken into account, the short distance cut-off scale in the logarithm determining the system size dependence of ln P becomes of the order of the electron wavelength p−1 F rather than the mean free path l. This leads, however, only to a small relative change in the σ-model result. On the other hand, the apparent ensemble independence of the distribution function asymptotes obtained by our method stands in a striking contradiction with the letter and the spirit of the σ-model and is quite enigmatic.

(59)

However, from Eq. (9) we see that the corresponding values of b – which determine the size of this local fluctuation of the potential – become larger than the length of the sample 2L. This is clearly unphysical. Thus despite a slower t˜ dependence of the probability of local resonances, the same value of t˜ has a much larger probability to be produced by an accidentally formed Bragg mirror for all reasonable values of t˜. Note also that ln t˜∆ ∼ pF L corresponds to v ∼ 1 (or U ∼ EF ), which invalidates the perturbative expansion of Section III A. Essentially, the values of ln t˜∆ of the order of pF L or larger correspond to a trivial case: a sample is insulating because the potential is larger than EF almost everywhere except for a small island in the center where almost all the weight of an eigenstate at EF is concentrated. In this regime the distinction be-

3. The three-dimensional case.

The same as in the one-dimensional case issue of whether the σ-model calculations can correctly handle the role of fluctuations of the potential with the amplitude larger than the Fermi energy arises in connection with the ln2 t asymptotics of the distribution of wave

12

function amplitudes obtained in Ref. [7] for the threedimensional case. Derivation of the σ-model involves linearization of the spectrum near the Fermi energy. As a result, there is no intrinsic scale in the model that would relate the amplitude of the fluctuating potential to the electron energy. The only scale is provided by 1 . the dispersion of the fluctuations of the potential πνd τ This limitation of the model does not affect the computation of probabilities of typical events, or of the averages that are dominated by such events, because typical potentials are small compared to the Fermi energy. However, a problem arises when rare large-amplitude fluctuations of the potential U (r) become dominant. The σ-model cannot detect the existence of classical turning points around which the quasiclassical approximation breaks down. Thus a possible explanation of the result of Falko and Efetov [7] is that within the σ-model approach the classical turning point at r∗ is missed, and the Bragg mirror is effectively assumed to persist until distances of the order of the mean free path l. In the notation of Section III C that would correspond to RL |λ|2 |λ| ∼ l ln θ, and then A0 ∝ l v 2 r2 dr ∝ ∝ l ln2 θ, l leading to ln P ∼ − (pF l)2 ln2 t which coincides with the answer obtained by Falko and Efetov. The reason that such a cut-off scheme leads to a higher estimate for the probability P is that it corresponds to an incorrect solution of the Schr¨odinger equation in the region of large potentials U > EF which overestimates the rate of growth of the wave function amplitude towards the center of the sample. Another discrepancy between our results and those obtained with the help of nonlinear σ-models in the threedimensional case is the difference in powers of (pF l) in the exponents in Eqs. (1) and (57). We believe that this discrepancy has the same origin as the difference L L between ln in the two-dimensional case – and ln l r0 the error introduced by using the mean free path l to determine the cut-off scale. It is interesting to note that while the cut-off procedure used in the diffusive σ-model approach of Ref. [8] is capable of producing only an order-of-magnitude estimate C3 ∼ (pF l)2 , the ballistic σ-model of Muzykantskii and Khmelnitskii [10,9], while giving an illusion of computing the coefficient C3 exactly π 2 (C3 = √ (pF l) ), nevertheless leads to the same extra 9 3 power of (pF l). To explain this seemingly paradoxical situation we will first examine the cut-off procedure employed in Ref. [8]. It is based on the condition, pointed out in Ref. [4], that in order for the diffusive σ-model to be applicable, the spatial gradients of the Q-matrix components cannot exceed 1/l. Explicitly, for the calculation performed in Ref. [8] this condition reads d ln λ1 < 1 , (60) dr l

where λ1 parametrizes the non-compact bosonic sector of the Q-matrix [5] (not to be confused with the Lagrange multiplier λ used throughout this work). The distance l∗ where this condition is violated is used in Ref. [8] as a short-distance cutoff. It was conjectured in Ref. [7] and later confirmed in Ref. [20] that spatial structure of the saddle-point solution for λ1 mimics the envelope of anomalously localized states described by the saddle point of the σ-model. Therefore the optimal configuration of (d/dr) ln λ1 (r) corresponds to the non-oscillating part of y (r) in the direct optimal fluctuation method. This correspondence allows us to see directly what effect an artificial short-distance cutoff at l∗ would have in our approach. From Eq. (49) we find l∗2 ∼ l|λ|. Introducing such a cut-off at l∗ into Eq. (50) we obtain |λ| ∼ l∗ ln θ and therefore l∗ is estimated as l∗ ∼ l ln θ. If then the integral determining the saddle-point action A0 is also cut off at l∗ , we obtain A0 ∼ l 2

Z

L

l∗

v 2 r2 dr ∼ l

|λ|2 ∼ l2 ln3 θ, l∗

or (pF l) ln3 θ in conventional units, leading to the incorrect answer that was already quoted in the Introduction. It is thus evident that distances shorter than l∗ give an important contribution to the saddle-point action, and this contribution cannot be accounted for by the diffusive nonlinear σ-model. It should be emphasized that excluding the short-distance contribution in this way leads to a significantly smaller estimate for the probability of observing a given (large) value of t. This can be understood by using the following argument. The wave function amplitude grows substantially between r ∼ l∗ and r ∼ 1. Neglecting this growth leads to a need for a faster increase in t between L and l∗ which can only be achieved by means of a “costly” boost in the amplitude of the Bragg mirror. As a result, non-optimal configurations of the random potential are selected. Turning now to the calculation performed with the help of ballistic nonlinear σ-model for the three-dimensional case in Ref. [9] we notice that the distance r∗ which separates the “reaction” and “run-out” zones is of the same order as l∗ . The contribution to the saddle point action from the “run-out” (diffusive) zone leads to the π 2 already quoted √ (pF l) value for C3 , while the “re9 3 action” zone produces a contribution that has one less power of the large logarithm ln θ and is thus neglected. In contrast, in the calculation presented in Section III C of this paper, the contribution from distances of the order of r∗ ∼ ln θ ≪ l∗ dominates the saddle point action and leads to a larger estimate (Eq. (5a)) for the probability of observing anomalously high values of t in the 3dimensional case than the one obtained in Ref. [9]. Thus, the ballistic generalization of the nonlinear σ-model is also not capable of detecting the existence of the scale r∗ at which the quasiclassical approximation breaks down. Moreover, since the calculation within the framework of 13

choice of the cut-off scale made in Ref. [7]. The region of the stretched exponential decay of the wave function envelope described by Eq. (54a) is missed in the σ-model calculation entirely. As in the two-dimensional case, the extra powers of r in the denominators of Eqs. (54) are a consequence of the artificial rotational symmetry. It was conjectured in Ref. [7] that these high-amplitude states have a complicated “snake-like” structure at short distances. The conjecture was based on the fact that the method employed in Ref. [7] was applicable to amplitudes as high as

the ballistic σ-model does not involve any ultraviolet divergencies that would necessitate a short-distance cut-off as in the case of the diffusive σ-model, it appears plausible that the ballistic variant of the model is equivalent to introducing an ultraviolet regularization into the theory. D. Prelocalized wave functions.

To complete the comparison of the results obtained by the direct optimal fluctuation method and those found using the nonlinear σ-models we now turn to the issue of the shape of the envelope of anomalously localized states that are responsible for the large-t tails of the distribution functions P (t). In Ref. [7] it was found that in two dimensions such states are characterized by amplitudes decaying outwards according to a power law: ln t/ ln(L/l) l . |ψ (r)| ∼ r 2

t< ∼

2

ψ 2 (0) r0 ln t/ ln rL0 . r r

(64)

V rather than a na¨ıve t < ∼ ld expected from a cut-off at l. We have found no evidence for such behavior here. Since the solutions of the saddle point equations were assumed to be rotationally invariant from the outset, this cannot be regarded as a conclusive evidence against such a scenario. However, potentially more important is the fact that we did not encounter any limitations on the possible values of t analogous to Eq. (64). Analysis of the fluctuations around the saddle point defined by Eqs. (15) as well as an investigation of the possibility for non-rotationallyinvariant solutions of these equations is needed to settle the issue conclusively.

(61)

This result has to be compared with Eq. (39) which can be rewritten as |ψ (r)| ∼

V pd−1 F , l

(62)

Apart from the difference in the cut-off scale (l vs. r0 ) which was discussed at length above, Eq. (62) contains an extra 1/r factor in the denominator. It is simply a consequence of the idealized model adopted here with its circularly symmetric boundary conditions. In a more realistic model the optimal wave function would become a superposition of different angular momentum eigenstates. √ The 1/ r behavior of the circularly symmetric component would be cancelled in such a superposition. Taking into account fluctuations around the saddle point would √ 2 also have the effect of suppressing the (1/ r) factor in the spatial dependence of the averaged envelope of the optimal solution. Note also that Eq. (61) does indeed describe the averaged envelope of anomalously localized states. In three-dimensional samples the states with anomalously high local amplitudes were found in Ref. [7] to have the following envelope: l 2 |ψ (r)| ∼ exp −A 1 − , (63) r

E. Are the σ-model results universal?

Finally, we would like to make a few remarks concerning the issue of universality. An important consequence of the dominant role played in three dimensions by large local fluctuations of the potential that are responsible for the formation of local resonances is the non-universal character of the distribution functions derived in Section III C. Indeed, such large-amplitude configurations of the potential can only be optimal if they are not too “expensive” – i.e. if their action A is not too high – compared with the low-amplitude “global” fluctuations of the potential (Bragg mirrors). Generalizing the distribution of the random potentials to Z 1 Wn [U ] ∝ exp − U 2n (r) dd r , (65) 2σ where σ is the dispersion and n ≥ 1, we find through a model calculation similar to the one performed in Section II that forming states with large amplitudes at some r by means of local resonances always leads to ln P ∝ − lnd t. For n < d/2 the “cost” of a corresponding Bragg mirror can be estimated from an appropriate generalization of Eqs. (16), and it also leads to a lnd t dependence, as was demonstrated by a detailed study of Gaussian (n = 1) distribution in three dimensions. Thus local large-amplitude fluctuations of the potential are always important for small enough (< d/2) values of n. In the

where A is a constant which in the leading logarithmic approximation is equal to ln t. Comparing this to Eq. (54) we see, in accordance with the discussion above, that the σ-model result gives the correct functional form eConst./r only for the region of large r which corresponds to the region of space where the optimal potential forms a Bragg mirror. The estimate l ln t for the constant in the exponential, however, contains an extra factor of l compared to Eq. (54b) which is a consequence of the 14

extend our method to the quasi-1D and quasi-2D geometries. To conclude, the main results of this work can be summarized as follows. We have demonstrated that the optimal fluctuation method is a useful tool for the investigation of the statistical properties of anomalous electronic eigenstates in disordered two- and three-dimensional conductors. It was shown that in three dimensions this method is preferable to the nonlinear σ-model because the latter does not include effects associated with local resonances that can be formed by the random potential. In the two-dimensional case the results obtained by the optimal fluctuation method essentially coincide with the σ-model results for the unitary ensemble of random potentials which we interprete as an indication that the saddle point of the reduced nonlinear σ-model found in Refs. [4,6–9] corresponds to the same saddle point of the full theory as the one that describes the optimal fluctuation of the potential.

marginal case n = d/2 the weight of the Bragg mirror is given by lnd t d−1 , L ln r0

ln P ∝ −

(66)

making it a much more probable way of achieving a large local wave function intensity t. When n > d/2 the probability of the corresponding Bragg mirror has a faster than lnd t dependence on ln t, but it is compensated by a power of the system size in the denominator: ln P ∝ −

ln2n t . L2n−d

(67)

This situation is realized in the Gaussian case when d = 1, and it was already discussed in detail. Therefore the seemingly universal character of the log-normal behavior of the tails of the distribution functions is tied to the assumption – which is a starting point in the derivation of the σ-model – that the random potential has a Gaussian distribution.

ACKNOWLEDGEMENTS

One of the authors (I. S.) would like to thank N. Wingreen, H. Li, D. Fisher, B. Halperin and A. Andreev for useful discussions. I. S. also acknowledges financial support from NSF grant DMR 9106237.

V. CONCLUSIONS AND OPEN QUESTIONS.

The most important problem arising from the results of the present study is the need to explain the fact that the asymptotics of the distribution functions derived using the direct optimal fluctuation method do not exhibit any dependence on weak magnetic fields. This feature of our answers must be contrasted with all the previous calculations performed in the framework of nonlinear σmodels in which it is quite obvious that when Cooper modes acquire a mass as a consequence of broken invariance with respect to time reversal, the number of independent components of the Q-matrices changes and that has a profound impact on the results. There are two alternative resolutions of this contradiction. One possibility is that our assumption that rotationally invariant solutions of the saddle-point equations dominate the saddle-point action is not valid in the absence of magnetic field. A more detailed study of the properties of Eqs. (15) will be the subject of a future publication. The second possible explanation is that there is some internal inconsistency in the σ-model which causes it in certain circumstances to overcount the number of degrees of freedom. In our view, this possibility must be taken seriously. Whether or not the fluctuations around the saddle point can change the leading order terms in ln t is also one of the questions that were left beyond the scope of this work. A peculiar disagreement in the crossover scale from ln2 t to ln t asymptotics of ln P between the purely one-dimensional and quasi-one-dimensional cases which was noted in Section IV C 1 makes it rather desirable to

APPENDIX:

Introducing an auxiliary field χ as in the main text, we can represent Nψ as Z χ i R drχ(r) pˆ 2 +U(r)−E ψ(r) 2m −1 Nψ = DψD . (A1) e 2π Performing a π/4 rotation in the (χ, ψ) space √ √ χ = (ψ1 + ψ2 ) / 2, ψ = (ψ1 − ψ2 ) / 2,

(A2)

and introducing infinitesimal convergence factors we obtain Z ψ2 ψ1 D √ Nψ−1 = D √ 2π 2π R (+) 1 ˆ (−) ψ2 } ˆ O O ψ −ψ i dr{ψ 1 2 1 × e2 , (A3)

ˆ (±) = pˆ 2 +U (r)−E±i0. Thus, up to a constant, where O 2m Nψ is equal to a symmetrized spectral determinant: r n o n o ˆ + i0 det E − H ˆ − i0 Nψ ∝ det E − H ˆ

= eℜTr ln(E−H+i0) .

(A4)

The variation of the logarithm in the exponent with respect to U (r) is equal to the real part of the Green’s function at coinciding arguments ℜG (r, r; U ). This quantity probes the whole band and is not sensitive to small 15

changes in U . It therefore can be treated as a constant (which we denote as G), so that Nψ becomes R Nψ ∝ eG drU(r) . (A5) R The background potential dr U (r) can be absorbed into a redefinition of energies, which justifies the approximation made in deriving Eqs. (15). It should also be noted that in the sigma model formalism the real part of the Green’s function G (r, r) is effectively set to zero under the assumption of an infinite symmetric band. Therefore any corrections to the distribution function that may arise due to Nψ 6= 1 are at any rate beyond the scope of the σ-model. [16] [17]

[1] Mesoscopic Phenomena in Solids, B. L. Altshuler, P. A. Lee and R. A. Webb, eds., North Holland 1991. [2] B. L. Altshuler, V. E. Kravtsov and I. V. Lerner, in [1], and references therein. [3] F. Wegner, Z. Phys. B 35, 207 (1979). [4] B. A. Muzykantskii and D. E. Khmelnitskii, Phys. Rev. B 51, 5480 (1995). [5] K. B. Efetov, Adv. Phys. 32, 53 (1983). [6] V. I. Falko and K. B. Efetov, Europhys. Lett. 32, 627 (1995). [7] V. I. Falko and K. B. Efetov, Phys. Rev. B52,17413 (1995). [8] A. D. Mirlin, Phys. Rev. B53, 1186 (1996). [9] B. A. Muzykantskii and D. E. Khmelnitskii, unpublished (preprint cond-mat/9601045). [10] B. A. Muzykantskii and D. E. Khmelnitskii, JETP Letters 62, 76 (1995), [Pis’ma Zh. Eksp. Teor. Fiz. vol. 62, p. 68 (1995)]. [11] B. Altshuler and V. Prigodin, JETP Letters 47, 43 (1988), [Pis’ma Zh. Eksp. Teor. Fiz. vol. 47, p.36 (1988)]. [12] B. I. Halperin and Melvin Lax, Phys. Rev. 148, 722 (1966). [13] J. Zittartz and J. S. Langer, Phys. Rev. 148, 741 (1966). [14] Such a formulation of the problem implies perfect rotational symmetry and sidesteps the issue of the role of irregular boundary conditions. However, as will become evident from the answers, the asymptotic form of the three-dimensional distribution is not sensitive to the boundary conditions at all. The two-dimensional case is trickier and may require a separate study of the stability of the solution with respect to the change in boundary conditions. Nevertheless, since the σ-model calculations were performed for similar geometries, the model adopted here is completely adequate for the purposes of comparing different approaches. [15] In the formalism of Falko and Efetov [6] the extra deltafunction in the definition of P was essential in enabling the mapping of the problem onto the nonlinear σ-model. However, the saddle-point solution of the reduced variant of the model which was derived in Ref. [6] to describe de-

[18] [19] [20] [21]

viations from the universal Porter-Thomas statistics [16] was shown to involve only the variables parametrizing the bosonic sector of the theory. This can be viewed as an indication that effects associated with the statistics of energy levels are not important in the approximation corresponding to the saddle-point treatment of the nonlinear σ-model. In the context of the direct optimal fluctuation method this observation receives quite a natural explanation. A given value of energy E can be fine-tuned to become an eigenvalue by a shift in the optimal configuration of U (r) of the order of the mean level spacing 1 ∆= . Such a shift can result in corrections to ln P νd V that are at most ∼ EF τ – much less then the leading contribution when ln t is large. Haake, Quantum signatures of Chaos, Springer-Verlag, Berlin, New York 1991. In several cases the probabilities obtained by our method are smaller than the correponding σ-model results. The relation between the two sets of answers and the likely reasons for the differences between them are discussed in Sections III and IV. E. Kamke, Differential Equations, Chelsea Publ., NY 1971. W. H. Press et.al, Numerical Recipes in C, Cambridge University Press, Cambridge 1988. A. D. Mirlin, unpublished (preprint cond-mat/9512095). I. M. Lifshitz, Soviet Phys. – JETP 17, 1159 (1963).

FIG. 1. The profile of an optimal configuration of the potential in three dimensions as a function of the radial coordinate.

16

0.4 0.2

10

0 -0.2

(ii)

5

-0.4

0

5

10

15

20

0 (iii) (i)

-5

-10

0

5

10

Fig. 1

15

20