Interaction matrix element fluctuations in ballistic quantum dots: random wave model

Interaction matrix element fluctuations in ballistic quantum dots: random wave model


We study matrix element fluctuations of the two-body screened Coulomb interaction and of the one-body surface charge potential in ballistic quantum dots. For chaotic dots, we use a normalized random wave model to obtain analytic expansions for matrix element variances and covariances in the limit of large (where is the Fermi wave number and the linear size of the dot). These leading-order analytical results are compared with exact numerical results. Both two-body and one-body matrix elements are shown to follow strongly non-Gaussian distributions, despite the Gaussian random nature of the single-electron wave functions.

73.23.Hk, 05.45.Mt, 73.63.Kv, 73.23.-b

I Introduction

There has been much interest in the properties of quantum dots whose single-particle dynamics are chaotic (1). The generic fluctuation properties of the single-particle spectrum and wave functions in such dots are usually described by random matrix theory (RMT) (2). In open dots that are strongly coupled to leads, electrons can often be treated as non-interacting quasi-particles, and the mesoscopic fluctuations of the conductance have been explained using RMT.

However, in almost-isolated dots, electron-electron interactions are important and must be taken into account. The simplest model of such dots is the constant interaction (CI) model, in which the interaction is taken to be the classical charging energy. Charging energy leads to Coulomb blockade peaks in the conductance versus gate voltage. Each peak occurs as the gate voltage is tuned to compensate for the Coulomb repulsion and an additional electron tunnels into the dot. For a fixed number of electrons, the CI model is essentially a single-particle model, and RMT can be used to derive the statistical properties of the conductance peak heights (3). While the CI plus RMT model has explained (at least qualitatively) (3); (4); (5) several observed features of the peak height fluctuations (6); (7); (8), there have been significant discrepancies with experimental data, in particular regarding the peak spacing statistics (9); (10); (11); (12). Such discrepancies indicate the importance of interactions beyond charging energy.

A more systematic way of treating electron-electron interactions in chaotic ballistic dots is to expand the interaction in a small parameter, the inverse of the Thouless conductance , where is the Fermi wave number and is the linear size of the dot. The Thouless conductance measures the number of single-particle levels within an energy window determined by the time it takes the electron to cross the dot, and increases as the square root of the number of electrons. It can be shown that, in the limit of large Thouless conductance, only a few interaction terms survive, constituting the interacting part of the universal Hamiltonian (13); (14). These universal interaction terms include, in addition to charging energy, a constant exchange interaction. The inclusion of an exchange interaction has explained the statistics of peak heights at low and moderate temperatures as well as the suppression of the peak spacing fluctuations (15); (16). However, at low temperatures, the peak spacing distribution remains bimodal even when the exchange interaction is included, while none of the experimental distributions are bimodal (9); (11); (10); (12).

For finite Thouless conductance, residual interactions beyond the universal Hamiltonian must be taken into account. The randomness of the single-particle wave functions induces randomness in the two-body screened Coulomb interaction matrix elements (17). The possible induced two-body ensembles have been classified according to their underlying space-time symmetries and features of the two-body interaction (18). In a Hartree-Fock-Koopmans (19) approach (assuming the Hartree-Fock single-particle wave functions do not change as electrons are added to the dot), the peak spacing can be expressed directly in terms of certain diagonal interaction matrix elements (20). Sufficiently large fluctuations of these interaction matrix elements can explain the absence of bimodality in the peak spacing distribution (20); (21). The variance of these fluctuations is determined by the spatial correlations of the single-particle wave functions. In a diffusive dot these correlations have been derived to leading order in , and the variance of the matrix elements of the screened Coulomb interaction was shown to behave as  (17); (22), where is the mean level spacing of the single-particle spectrum. However, dots studied in the experiments are usually ballistic. Wave function correlations in ballistic quantum dots are much less understood. Berry’s conjecture (23) regarding the Gaussian nature of wave function fluctuations in chaotic systems provides the leading order behavior of the correlations at short distances, but finite-size contributions at distances that are comparable to the size of the dot can have important effects on the matrix element fluctuations.

An additional contribution to the peak spacing fluctuations originates in surface charge effects (17). In a finite-size system, screening leads to the accumulation of charge on the surface of the dot. The confining one-body potential is then modified upon the addition of an electron to the dot. In a diffusive dot, the variance of one-body matrix elements behaves as .

Another interesting phenomenon, in which interaction matrix element fluctuations play an important role is spectral scrambling as electrons are added to a chaotic or diffusive dot (24); (25). Both the two-body interaction and one-body surface charge effects are responsible for scrambling.

Here we investigate fluctuations of the two-body interaction matrix elements and of the surface charge potential matrix elements in ballistic dots. We use a normalized version of Berry’s random wave model to derive analytically leading-order contributions for an arbitrary dot geometry, and compare with exact numerical simulations.

The outline of this paper is as follows. In Sec. II we review the random wave model for chaotic billiards. The spatial correlator of wave function intensity obtained from this model is geometry-independent but not consistent with the normalization requirement of the wave functions (26). We compute normalization corrections to the correlator in Sec. III. The variances of diagonal, double-diagonal, and off-diagonal two-body interaction matrix elements are calculated in Sec. IV and are all found to be strongly affected by the normalization correction. As a result, ratios of these variances are shown to remain far from their asymptotic values for the range that is relevant to experiments. In Sec. IV.3 we study the covariance of interaction matrix elements, relevant for understanding spectral scrambling when several electrons are added to the dot (24), and supplement the full random wave analysis with a schematic random matrix model. One-body matrix element variances are treated in Sec. V. In all of these studies we compare numerical results of the random wave model with leading-order analytical estimates. In Sec. VI we study the interaction matrix element distributions and show that, in the experimentally accessible range of , they deviate significantly from the Gaussian limit implied by the central limit theorem. Finally, in Sec. VII we address additional implications of the present work, including the necessity to supplement the normalized random wave model with dynamical effects for quantitative comparison with experiment.

Ii Random Wave Model

In a chaotic system without symmetries, a typical classical trajectory will uniformly explore an entire energy hypersurface in phase space. A well-established and extensively tested conjecture by Berry (23) holds that in the quantization of such a system, a typical wave function will spread uniformly over an energy hypersurface, up to the inevitable Gaussian random amplitude fluctuations that arise whenever a random vector is expanded in a generic basis. In a semiclassical picture, the Gaussian fluctuations are the result of exponentially many wave fronts visiting any given region of phase space with quasi-random phases.

For a two-dimensional billiard system, the random wave model implies that a typical chaotic wave function may be written locally as a random superposition of plane waves at fixed energy :


where is distributed as a -correlated Gaussian random variable: and . The normalization is fixed by , where is the area of the billiard. Equivalently, the wave function may be expanded in circular waves with good angular momentum :




The discrete variables () are taken to be uncorrelated Gaussian random variables with


The normalization integral of a wave function (2) is given by


where are the basis state overlaps


Defining to be a column vector with components , the normalization integral of the wave function is simply the norm of with the matrix playing the role of a metric


The random coefficients may be chosen to either satisfy or not satisfy the real wave function condition , corresponding to the absence or presence of an external magnetic field; these two situations are conventionally denoted by , respectively.

Using the addition theorem for Bessel functions we have


It follows from the completeness property (8) of the Bessel basis that


We note that formally the model requires an infinite set of basis states; in practice the effective number of basis states that have appreciable magnitude inside the area scales as , where is a typical linear size of the billiard and is the ballistic Thouless conductance. A more accurate estimate for can be obtained by considering a disk of radius . In a disk, the wave functions are orthogonal because of rotational symmetry and . The effective dimension is then obtained as the “participation number” of the exact :


The approximate completeness of a basis of plane waves at a fixed energy was confirmed in studies of a billiard system (27).

Iii Intensity Correlator

The random wave model of Eqs. (2) and (4) together with Eq. (8) leads to the amplitude correlator


and intensity correlator


Eq. (12) is obtained from Eq. (11) by contracting with and with and noting that there are two equivalent ways of performing the contraction when (i.e., for ).

Similar correlators can be derived starting from the plane-wave expansion (1). The intensity correlator (12) is valid to leading order in , but becomes problematic when applied to to all , in the finite area . Indeed wave function normalization requires the correlator to vanish on average (22),


and similarly . However, the random wave intensity correlator in (12) is non-negative everywhere and does not satisfy the condition (13). The reason for this failure is that in the random wave model normalization is satisfied only on average, i.e., , but not for each individual wave function.

This deficiency can be corrected by introducing the normalized random wave model, in which each “random” wave function (2) is normalized in area , i.e.,


with . The normalized random wave model is easy to implement numerically by normalizing each random wave. Analytically, the intensity correlator of this model can be written as


This is similar to the situation in random matrix theory (RMT), where the naive guess for ( and are two components of an RMT eigenvector of length with normalization ) must be replaced by the exact expression to obtain correct normalization for finite . In our case, however, the normalized intensity correlator (15) depends on both positions and as well as the system geometry.

We define at a spatial point an excess wave-function intensity by and similarly an excess normalized wave-function intensity . We then have


Eq. (III) is exact but unwieldy. In the semiclassical limit of large , a given superposition of random waves will be almost normalized, i.e., . Thus, to leading order in , we may expand Eq. (III) to obtain


where we have omitted all terms involving three-point and higher-order correlations of the excess intensity . Eq. (III) can be simplified to obtain




The leading-order normalized correlator was derived in Ref. (26) by adding a weak smooth disorder and using the non-linear supersymmetric sigma model. More recently, the same leading correction has been obtained (28) starting from a density matrix and the principle of maximum entropy (29). Here we we have shown that this form can be obtained quite generally by applying a perturbative approach to the unnormalized correlator , assuming only that correlations are weak in the limit of large system size .

We note, however, that is merely a leading order approximation in to the true normalized random wave correlator , as it involves only terms depending on the unnormalized two-point correlator . The complete expression (III) for the normalized two-point correlator involves all unnormalized -point correlators . In particular, the unnormalized three-point correlator gives rise to the correction in Eq. (18). In principle, such higher-order corrections to Eq. (18) may be computed systematically, but in practice their effect on matrix element variances is small, as we will confirm below. The -point correlators for will be important when we discuss the deviation of the matrix element distribution from a Gaussian (see Sec. VI). We also note that our definition of and the approach of Refs. (28); (29) may produce different corrections to at higher orders in a expansion. Although the problem of constructing the “best” random wave ensemble that exactly satisfies normalization constraints is an interesting one in its own right, the effect of these higher-order ambiguities on matrix element statistics is negligible compared to dynamical contributions in real chaotic systems.

In contrast to the unnormalized correlator , the normalized correlator satisfies


implying the normalization of each individual wave function . To show this, we define . We note that by construction, and calculate the variance in wave function normalization:


Eq. (20) implies that the variance of vanishes, and since the average wave function is normalized to unity, every wave function in the ensemble must be normalized to unity.

A more schematic random matrix approximation is obtained if the exact metric is assumed to be diagonal and the diagonal components are replaced by


The distribution (22) satisfies and and therefore reproduces correctly the first two moments (9) and (10) of the exact distribution. In the approximation (22), the normalization of the wave function reads [see Eq. (7)]. The normalized -component vector can then be viewed as an eigenvector of a Gaussian random matrix of the corresponding symmetry class .

In a disk of radius , asymptotic expressions may be used to verify that the exact self-overlaps of the Bessel functions approach a constant for and fall off exponentially as for . In Fig. 1, we show the exact random-wave metric for a disk with and , as well as the schematic random matrix approximation (22).

Figure 1: The exact Bessel function self-overlaps defined by Eq. (6) are computed for a disk of radius , and scaled by the effective Hilbert space dimension . The solid curve indicates , while the dashed curve corresponds to . The dotted rectangle represents the schematic random matrix approximation, where all are taken to be constant for and otherwise. Quantities plotted in this and all subsequent figures are dimensionless.

In the following sections we use the random wave model to estimate the variances of various matrix elements of the screened Coulomb interaction. Normalized random waves can be generated numerically to calculate the spatial wave function correlator and the corresponding variances, which is equivalent to using the exact normalized intensity correlator . Analytical expressions can be obtained as follows. (i) The asymptotic behavior of the unnormalized random wave correlator in Eq. (12) will be sufficient to produce analytically the leading geometry-independent behavior of two-body interaction matrix element variances. (ii) A direct integration of the normalized correlator in Eq. (19) produces geometry-dependent terms for the above matrix elements. This correlator can also be used to calculate to leading order in the size of one-body surface charge matrix element fluctuations.

Iv Two-body interaction matrix elements

Here we model the screened two-body Coulomb interaction in 2D quantum dots as a contact interaction , where the single-particle mean level spacing serves to set the energy scale (30); (31).

iv.1 Fluctuation of diagonal matrix elements

We first discuss the variance of a diagonal two-body matrix element . In the contact interaction model


To leading order in , the dominant contribution to the variance arises from correlations between the intensities of a single wave function at different points (24):


Since this leading-order effect comes from correlations within a single wave function, the only restriction on the energy difference between the two random wave functions is , or equivalently , allowing the variance to be a function of a single parameter . The subleading terms omitted in Eq. (24) involve the two-eigenstate intensity correlator


The correlator (25) produces the leading effect for the covariance and will be discussed in detail in Section IV.3.

The leading-order contribution to Eq. (24) is obtained by using the unnormalized correlator instead of . Changing integration variables in Eq. (24) to and , inserting the unnormalized correlator (12), which depends only on , and substituting the asymptotic form , we obtain


where in the last line we have used .

In Eq. (26) we have not properly included the short distance contribution from and the shape-dependent long-distance contribution from . Both of these contributions to Eq. (24) scale as ; in the first case because defines an volume in -space, and in the second case because for . To obtain the correct result at this subleading order we need to use the normalized correlator of Eq. (19). We then obtain


The leading term, discussed in Ref. (31), depends only on the symmetry class, while the shape-dependent coefficient can be easily evaluated by numerical integration of Eq. (28).

Figure 2: The variance of the two-body matrix element versus in the real () random wave model: (a) the solid curve is the result of exact numerical simulations; (b) the long-dashed line is the result of integrating the square of the normalized correlator , as in Eq. (28); (c) the short-dashed line is the result of using the unnormalized correlator . In the inset, the solid line is the difference between the full result (a) and the approximation (b); the dotted line indicates that omitted terms scale as . Here and in all following figures, the level spacing , which sets the overall energy scale, is set to unity.

In Fig. 2 we show the results for a disk geometry in the range , corresponding roughly to the parameter range relevant for experiments ( electrons in the dot). The numerical simulation of (solid line) is compared with the leading result of Eq. (28) (long-dashed line). For the disk, we find . The difference between the full numerical calculation and the integral of the normalized correlator is plotted in the inset, and is seen to be in good agreement with the scaling of Eq. (28). If the unnormalized correlator (12) is used in (28) we find (for the disk geometry) . The corresponding variance is shown as a short-dashed line in Fig. 2.

Figure 3: The change in the variance (upper two curves) and the variance of the one-body matrix element (lower two curves, see Sec. V) for normalized real random waves are plotted as functions of the aspect ratio of the elliptical geometry. Here is kept fixed at (solid curves) and (dashed curves).

The dependence of the fluctuations on the dot’s geometry is demonstrated in the upper two curves Fig. 3. Here we plot the result of Eq. (28) for elliptical shapes as a function of the aspect ratio between the major and minor axes, while keeping the area and the wave vector fixed. We find that as the aspect ratio changes by a (physically unrealistic) factor of , the shape-dependent parameter only changes from to , resulting in a change of only in the variance for the experimental range of . Similar results are found if elliptical shapes are replaced by rectangles or other geometries. Thus, for all practical purposes, shape effects on the variance can be ignored (at least within the normalized random wave model).

At this subleading order, we must in principle also take into account the energy difference , if this energy difference is larger than the ballistic Thouless energy (but still small compared with the average ). Then , and the factor in Eq. (26) must be replaced by , which has the average value only for and is reduced to an average of for . After performing the integration over in Eq. (26), we find that in the final result must be replaced with . Assuming remains large but fixed while , this corresponds merely to a modification of the geometry-dependent coefficient in Eq. (28): for . In practice, for reasonable energy windows, e.g., , the consequent reduction in the variance is at most and may be safely ignored compared to the much larger dynamical effects present in real chaotic systems.

iv.2 Fluctuation of and

The preceding section studied fluctuations in the intensity overlap between two random wave functions and . Similar techniques may be applied to the self-overlap of a single wave function , also known as the inverse participation ratio in coordinate space, and to the “off-diagonal” four-wave-function overlap .

To leading order, , all the results arise from integrating the unnormalized correlator (12) and differ only by combinatorial factors. For example, the factor in Eq. (26) may be understood from the fact that there are four distinct ways to contract pairs of same-wave-function amplitudes between and if the wave functions are real (), but only one way to perform this contraction if the wave functions are complex (). An analogous counting argument for the variance of leads to


with and . We note that such combinatorial factors can also be derived from the invariance of the second moments of the matrix elements under a change of the single-particle basis (18).

Finally, if , , , and are all different, there is only one way to perform the contraction in either the real or complex case, leading to


While the leading behavior in Eqs. (29) nd (30) is identical to that for up to geometry-independent combinatorial factors, the coefficients of the subleading geometry-dependent terms are not so simply related, as indicated by in the above expressions. That is because the normalization-related subtraction, which enters at , works differently for the three matrix elements. For example, the variance of involves the product of two normalized intensity correlators . On the other hand, the variance of involves the product of four amplitude correlators (i.e., ), which do not require subtraction. As indicated earlier, the absence of subtraction in the integral results in for a disk geometry, in contrast with . For fluctuations, we find for a disk (32). Thus, although and are related to by geometry-independent universal factors in the limit, the convergence to this universal limit is logarithmically slow. Specifically, in the presence of time-reversal symmetry (),


to leading order in .

Fig. 4 shows these ratios versus in the physical range of interest. Note in particular the very slow convergence of Eq. (31) to the asymptotic value of because of the large difference between the coefficients of the subleading terms .

Figure 4: Ratios of matrix element variances are plotted for real random waves inside a disk as functions of : (thick solid line) and (thick dashed line). For comparison, the leading-order analytic results of Eqs. (31) and (32) are shown as thin solid and dashed lines, respectively.

iv.3 Matrix element covariance

Another quantity of interest is the covariance , where , , and are three distinct wave functions of a single dynamical system. Such matrix element covariances are important in quantitative estimates of scrambling of the Hartree-Fock single-particle spectrum as electrons are added to the dot (24). As in Eq. (24), the leading contribution for large can be written as a double integral of a product of two intensity correlators:


where is the intensity correlator for wave function and is the intensity correlator (25) between two distinct wave functions and . In a diffusive dot, expression (33) leads to a covariance , where is the diffusive Thouless conductance (24).

We proceed to evaluate in the normalized random wave approximation. Clearly, if , are chosen to be independent random wave functions, the correlator vanishes by construction. Therefore, we need to impose the orthogonality condition


which may be done through the Gram-Schmidt procedure. We begin by generating two independent random vectors and at the same energy : and , project onto the subspace orthogonal to , and normalize the result to unit length to find :


Using the basis states in Eq. (3), we find


Eq. (36) can be evaluated numerically for a given geometry in order to obtain the covariance of diagonal matrix elements via Eq. (33). However, we obtain insight into the qualitative behavior of the correlator (36) by using the schematic random matrix model (22) in a disk geometry. We find


where is the effective dimension defined by (10). Substituting the random matrix model relation (37) in (33) and using Eq. (24), we have


This simple relation between the variance and covariance of the interaction matrix elements could also be obtained directly from an effective completeness condition (valid for any )


where . Eq. (39) follows from the completeness of the eigenstate basis