# The Addition Spectrum and Koopmans’ Theorem for Disordered Quantum Dots

## I Introduction

We consider the response of the ground state of spinless fermions to the addition or removal of a particle. To this end, we apply the self-consistent Hartree-Fock (SCHF) approximation: a non-perturbative effective single particle theory.

Koopmans’ theorem  states that the single particle SCHF energy levels describe the affinity and ionisation energy spectra for the unoccupied and occupied states respectively. The approximation involved is that all the other particles do not react to this process. This approximation is generally considered to be good when the single particle wavefunctions are extended: corrections to each wavefunction due the rearrangement of the system following the addition or removal of a particle are expected to be $`𝒪(1/N)`$, where $`N`$ is the number of particles in the system . Moreover, these corrections to the wavefunctions are expected to disappear in the limit of vanishing disorder (having e.g. periodic boundary conditions), where the single particle wavefunctions are free waves, as well as in the limit of vanishing interaction. It is then loosely assumed that in the thermodynamic limit, Koopmans’ theorem becomes exact even for disordered systems, and should be sufficiently accurate for mesoscopic samples. It is evident however, that if the physical quantities at hand require an energy resolution of order the mean single particle level spacing, $`\mathrm{\Delta }`$, the validity of Koopmans’ theorem should be reconsidered.

Our analysis of Koopmans’ theorem for a quantum dot is closely related to the addition spectrum of the latter: the spectrum of energy differences between states with total particle number different by unity. We consider only the energy differences between ground states, which is experimentally accessible through resonant tunnelling , and capacitance  measurements at low bias and temperature. We stress that while our analysis here pertains to some aspects of the experiments, a direct comparison is not feasible: firstly we consider spinless electrons, and secondly we consider disordered systems in the diffusive regime; in Ref. the mean free path is of order the sample size, whereas in Ref. the mean free path is much larger than the system size, and ergodicity is ensured by a chaotic boundary shape.

The position of the observed resonant tunnelling (RT) conductance peaks can be related to the ground state energy difference $`\mu _N=E_G(N,V_G)E_G(N1,V_G)`$ where $`E_G(N,V_G)`$ is the ground state energy of the dot with $`N`$ electrons, at gate voltage $`V_G`$. The spacing between consecutive peaks is thus related to

$$\mathrm{\Delta }_2(N)=E_G(N+1,V_G)E_G(N,V_G)E_G(N,V_G^{})+E_G(N1,V_G^{}).$$
(1)   
Within the constant interaction (CI) model the ground state energy is simply the sum of filled single particle energies $`e(n)`$ plus $`N(N1)V_0/2`$, where $`V_0`$ is the constant interaction. Taking $`V_G=V_G^{}`$, the peak spacing trivially reduces to

$$\mathrm{\Delta }_2(N)=e(N+1)e(N)+V_0$$
(2)   
and so, in the diffusive regime, displays shifted Wigner-Dyson (WD) statistics  up to corrections in one over the dimensionless conductance, $`g`$ : $`P(s)=\pi s/2\mathrm{exp}(\pi s^2/4)`$ for zero magnetic field, the case that we consider here; $`s=(\mathrm{\Delta }_2V_0)/\mathrm{\Delta }`$.

Recent experiments on quantum dots  have shown that whilst the mean peak spacings are well described by the CI model, the fluctuations are not described by Wigner-Dyson statistics. It is found that the distribution of $`\mathrm{\Delta }_2`$ is roughly Gaussian , with broader non-Gaussian tails seen in Ref.. In Refs. the variance of the fluctuations was found to be considerably larger than that given by the WD distribution. Further experimental observations, including correlations of peak heights  and the sensitivity to an Aharanov-Bohm flux  are not consistent with random matrix results, and suggest a breakdown of the naive single particle picture. To investigate this, one needs information on the ground state wavefunction.

Blanter et al  have evaluated the fluctuations within a Hartree-Fock framework, neglecting effects due to the change in gate voltage $`V_G`$ to $`V_G^{}`$. To this end they applied the Random Phase Approximation (RPA) to generate the screened interaction in the confined geometry, and assumed that all HF level spacings, except for the Coulomb gap, are described by WD statistics. Implicitly assuming Koopmans’ theorem  to be valid, and using wavefunction statistics established for non-interacting electrons in a random potential, they calculated the fluctuations of $`\mathrm{\Delta }_2`$ beyond the CI model. These additional fluctuations were found to be parametrically small (in $`1/g`$) and proportional to $`\mathrm{\Delta }`$. Hence, the total fluctuations in $`\mathrm{\Delta }_2`$ were found to be proportional to $`\mathrm{\Delta }`$. The analysis of Ref. is consistent in the limits $`g1`$, $`r_s1`$. The parameter $`r_s`$ characterises the relative importance of interactions in the electronic system, and is defined as the mean electron separation in units of the effective Bohr radius.

Exact numerical calculations on small disordered dots  did produce large Gaussian distributed fluctuations at experimental densities. It was claimed that for strong enough interaction $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2`$ is $`\mathrm{𝑢𝑛𝑖𝑣𝑒𝑟𝑠𝑎𝑙}`$, independent of the interaction strength and disorder, where $`\delta \mathrm{\Delta }_2`$ denotes the typical (RMS) size of the fluctuations, and the angle brackets denote disorder averaging. This $`\mathrm{𝑢𝑛𝑖𝑣𝑒𝑟𝑠𝑎𝑙}`$ constant was found to be approximately $`.10.17`$, to be compared with the WD result for the CI model: $`.52\mathrm{\Delta }/(\mathrm{\Delta }+V_0)`$. We note that the typical experimental value for the charging energy $`V_0`$, is much larger than $`\mathrm{\Delta }`$. The scaling with $`\mathrm{\Delta }_2`$ suggested by this analysis  is in stark contrast to the scaling with $`\mathrm{\Delta }`$ obtained in Ref..

Stopa has considered ballistic chaotic billiards numerically, using local density functional theory . In this case, it was claimed that the fluctuations arise due to strongly scarred wavefunctions in the self consistent potential. As a result of these scars, an asymmetric distribution of $`\mathrm{\Delta }_2(N)`$ was found, including strong correlations over $`N`$. It was then further noted that what is actually measured (i.e. the change in the gate voltage between resonant tunnelling peaks) is not simply related to $`\mathrm{\Delta }_2(N)`$ when the dependence of $`\mathrm{\Delta }_2(N)`$ on the gate voltage is strong. It was then claimed that a self-consistent calculation of $`\mathrm{\Delta }V_G`$ with $`\mathrm{\Delta }_2(N)`$ retrieves a symmetric distribution of $`\mathrm{\Delta }_2(N)`$, and reduces peak to peak correlations.

A further suggestion that the coupling to the gate is important in understanding the fluctuations has been made with reference to the CI model, with WD statistics for the single particle levels . The authors claim that the required distribution of $`\mathrm{\Delta }_2`$ can be generated, except for the non-Gaussian tails, through the de-correlation of neighbouring levels under a parametric change in the Hamiltonian (mediated by $`V_G`$). However, the degree of de-correlation induced by $`\mathrm{\Delta }V_G`$ is left as a fitting parameter.

In this paper we present numerical calculations within the SCHF approximation, considering larger samples than is feasible by exact diagonalisation . This approximation has been seen to be quite good for the calculation of persistent currents in similar systems . We show that fluctuations large compared to the single particle level spacing can arise without recourse to varying the sample shape, size or gate to dot coupling, supposing these to be additional effects. We further demonstrate that approximating the addition spectrum spacings by applying Koopmans’ theorem can lead to large errors in the calculation of the spacing statistics.

We consider separately both a long range (Coulomb) bare interaction and a short range (nearest neighbour) bare interaction. In section II we introduce our model in detail; in section III we present a short discussion of the implications of Koopmans’ theorem; in section IV we present and discuss our numerical results, which are then summarised in the final section.

## II The Model

We address the following tight binding Hamiltonian for spinless fermions

$$H=\underset{i}{}w_ic_i^+c_it\underset{i,\eta }{}c_{i+\eta }^+c_i+\frac{U_0}{2}\underset{ij}{}M_{ijij}c_i^+c_j^+c_jc_i$$
(3)   
where $`i`$ is the site index, $`\eta `$ describes the set of nearest neighbours, $`w_i`$ is the random on site energy in the range $`[W/2,W/2]`$, and $`t`$ the hopping matrix element, henceforth taken as unity. We study separately, both a Coulomb interaction potential,

$$M_{ijij}=1/|𝐫_i𝐫_j|$$
(4)   
and a short range potential plus a constant term $`M_c`$ (see below),

$$M_{ijij}=(\delta _{i,i+\eta }+M_c).$$
(5)   
We consider a 2D system with periodic boundary conditions, and choose to define

$$|𝐫_i𝐫_j|^2[L_x^2\mathrm{sin}^2(\pi n_x/L_x)+L_y^2\mathrm{sin}^2(\pi n_y/L_y)]/\pi ^2,$$
(6)   
where $`(n_x,n_y)𝐫_i𝐫_j`$.

All distances are measured in units of the lattice spacing $`a`$; the physical parameters are therefore $`U_0=e^2/a`$, $`t=\mathrm{}^2/2ma^2`$. The standard definition for $`r_s`$ is given, for low filling, by $`r_s=U_0/(t\sqrt{4\pi \nu })`$, where $`\nu =N/A`$ is the filling factor on the tight binding lattice with $`A`$ sites. The dimensionless conductance $`g`$ can be approximated, again for low filling, using the Born approximation. We find $`g=96\pi \nu (t/W)^2`$, which is valid for $`A,Ng1`$. Here $`\nu 1/4`$ throughout, so that $`r_s.56U_0/t`$, and $`g75(t/W)^2`$.

Having identified the parameters of our model with the standard ones employed in the theory of a continuous electron gas, we note that in the limit of small $`r_s`$ and $`1/g`$ the leading order term for the typical interaction dependent fluctuations predicted by Blanter et al  is $`U_0\mathrm{\Delta }/t\sqrt{g}`$. For the torus geometry considered here, this contribution, being a surface term, vanishes identically. Their prediction then reduces to typical fluctuations in addition to those of the CI model to be of order $`U_0\mathrm{\Delta }/tg`$.

The torus geometry has the advantage over geometries with hard walls whereby in the former, the compensating background charge provides a trivial shift in all the site energies, and can be removed. In a bounded dot, with an overall charge, the excess charge may build up near the boundary, depending on the position of nearby metallic plates and gates. These effects are geometry specific . Upon adding an electron, the average charge configuration may change considerably (the configuration is strongly geometry dependent). As the gate voltage is varied to allow for the next electron addition, the background potential could have changed causing further charge rearrangement. Whilst it is of great interest to analyse this issue (which may play an important role in the peak spacing fluctuations as well as undermining the naive single particle picture by further reducing the accuracy of Koopmans’ theorem ), we concentrate here on effects due entirely to the intrinsic rearrangement of the dot. From this point of view our analysis may be taken as an attempt to establish an upper bound criterion for the breakdown of Koopmans’ theorem. In reality it may break down earlier due to other non-universal factors. During the completion of this work very recent experimental evidence for significant rearrangement has been produced . It is argued that rearrangements due to adding an electron are far greater than rearrangements due merely to a change in shape.

When considering the short range interaction, the mean charging energy $`V_0`$ in Eq.(2) must be put in by hand through $`M_c`$ of equation (5). The way in which this is done depends on the physical situation being modelled, and is highly geometry dependent (vis-à-vis the gates). We stress that the value of $`M_c`$ does not affect the physical results. We choose to insert

$$M_c=V_0/U_04/A,V_0=\underset{𝐫,𝐫^{}}{}^{}\frac{U_0}{|𝐫𝐫^{}|}.$$
(7)   
This value for $`M_c`$ is defined such that if the charge is uniformly spread over the dot, the average charging energies in the Coulomb and nearest-neighbour cases roughly coincide. This choice has been made for simplicity, but corresponds to the premise that the interactions of the $`N`$-electron gas with the positive background and with itself is the same for both models considered. Exchange contributions, which tend to reduce the total charging energy, are included insofar as to cancel both the on-site contributions to the energy, and the unphysical self-interaction of electrons, but are otherwise neglected . The energy associated with charging the system uniformly is $`U_0N(N1)/(2A^2)_{𝐫,𝐫^{}}^{}|𝐫𝐫^{}|^1`$ in the Coulomb case, and $`U_0N(N1)/(2A^2)(4A+M_cA^2)`$ in the nearest-neighbour case. Equation (7) follows from equating these energies. This estimate can be systematically improved if the above premise is taken as the definition, not only by correctly accounting for the exchange contributions, but also by considering single particle wavefunction statistics in the diffusive regime. In this case wavefunction correlation functions such as $`|\psi _i(𝐫)|^2|\psi _j(𝐫^{})|^2`$ are required in order the calculate the average electrostatic energy, where here and after $`\mathrm{}`$ denotes averaging over the disorder ensemble.

In Ref. it was assumed that the (RPA) screening can be taken into account before constructing the Slater determinant ground state, and therefore their result corresponds to a short-ranged effective interaction. It is not clear that this remains a consistent procedure when calculating the ground state energy self-consistently. The reason for the inconsistency is that many of the diagrams generated by the SCHF approximation are already included in the RPA calculation of the screening, resulting in double counting. On the other hand, if the screening is generated externally (e.g. by close metallic gates), then it is consistent to insert a short ranged bare interaction, and this is the point of view taken here.

In some sense, the Coulomb interaction results can be considered as the opposite limit of the nearest-neighbour interaction, and is of interest in this context. However, it is more difficult to physically motivate the use of a Coulombic bare interaction unless one is considering very low electron densities. Screening is indeed weak in a $`2d`$ electron gas in a vacuum, even at high density, but the SCHF procedure cannot correctly generate screening by itself; whilst it can screen the Hartree contributions (as discussed above), it does not screen the exchange (Fock) term. However, we have verified that for the range of parameters considered here, fluctuations in the Hartree energy are larger than the typical fluctuations of the exchange energy. This suggests that the error made in not screening the exchange term correctly is not overly important.

## III Implications of Koopmans’ Theorem

Let us now consider the form of $`\mathrm{\Delta }_2`$, and approximations to it given by applying Koopmans’ theorem. We denote the diagonal matrix elements of the one-body operators in (3) by $`T_i^N`$, and the antisymmetrised Hartree-Fock interaction by $`V_{ij}^N`$, where hereafter the subscripts denote single particle states, and the superscript $`N`$ denotes the number of particles present and identifies the self-consistent basis of single particle wavefunctions being employed, $`\psi _i^N`$. For the torus geometry, where the gate voltage and background potential represent a trivial shift that can be omitted, the SCHF ground state energy is given by

$$E_G(N)=\underset{j}{\overset{N}{}}ϵ_j^N\frac{1}{2}\underset{ij}{\overset{N}{}}V_{ij}^N=\underset{j}{\overset{N}{}}T_j^N+\frac{1}{2}\underset{ij}{\overset{N}{}}V_{ij}^N,$$
(8)   
where $`ϵ_l^m`$ is the $`l`$th SCHF single particle energy for a system of $`m`$ particles in the ground state:

$$ϵ_l^m=T_l^m+\underset{j}{\overset{m}{}}V_{lj}^m.$$
(9)   
Using (8), we find (c.f. Eq.(1))

$$\mathrm{\Delta }_2(N)=T_{N+1}^{N+1}T_N^{N1}+\underset{j}{\overset{N}{}}\left(T_j^{N+1}2T_j^N+T_j^{N1}\right)$$

$$+\underset{j}{\overset{N}{}}\left(V_{N+1j}^{N+1}V_{Nj}^{N1}\right)+\frac{1}{2}\underset{ij}{\overset{N}{}}\left(V_{ij}^{N+1}2V_{ij}^N+V_{ij}^{N1}\right).$$
(10)   
Applying Koopmans’ approximation corresponds to dropping the superscripts and employing an appropriate fixed basis. The theorem implies that the effective single particle states do not depend on the occupation of these states. In particular, Koopmans’ theorem yields $`ϵ_{N+1}^N`$ for the minimum energy required to add a particle to a system of $`N`$ particles, and $`ϵ_N^N`$ for the maximum energy gained by removing a particle from the same system; in both cases the final state is a ground state. Clearly $`ϵ_l^m`$ as well as the ground state energy depend on $`m`$, even in Koopmans’ approximation, through the number of terms in the sum in Eqs.(9) and (8) respectively. It is then easy to see that Koopmans’ approximation yields

$$\mathrm{\Delta }_2^{k_1}(N)=ϵ_{N+1}^Nϵ_N^N.$$
(11)   
We also consider two other approximations to $`\mathrm{\Delta }_2`$ that involve calculating two self-consistent bases rather than just one:

$$\mathrm{\Delta }_2^{k_2}(N)=ϵ_{N+1}^Nϵ_N^{N1}$$
(12)   
$$\mathrm{\Delta }_2^{k_3}(N)=ϵ_{N+1}^{N+1}ϵ_N^N.$$
(13)   
All three estimates (11-13) coincide with $`\mathrm{\Delta }_2(N)`$ of Eq.(10) if Koopmans’ theorem holds. To connect with the notation of Ref., and to demonstrate the difference between the above three approximations and the fully self-consistent result, we provide a schematic diagram of the SCHF spectra in Fig. 1.

Since the self-consistent basis of $`N`$ particles provides the lowest energy for $`N`$ occupied levels, and similarly for $`N1`$ particles, the following relations are clear:

$$\underset{i}{\overset{N1}{}}T_i^N+\frac{1}{2}\underset{ij}{\overset{N1}{}}V_{ij}^N\underset{i}{\overset{N1}{}}T_i^{N1}+\frac{1}{2}\underset{ij}{\overset{N1}{}}V_{ij}^{N1}$$

$$\underset{i}{\overset{N}{}}T_i^N+\frac{1}{2}\underset{ij}{\overset{N}{}}V_{ij}^N\underset{i}{\overset{N}{}}T_i^{N1}+\frac{1}{2}\underset{ij}{\overset{N}{}}V_{ij}^{N1}$$
(14)   
Combining these equations, we find that $`\mathrm{\Delta }ϵ(N)\mathrm{\Delta }_2^{k_1}(N)\mathrm{\Delta }_2^{k_2}(N)0`$, or equivalently

$$\mathrm{\Delta }ϵ(N)ϵ_N^{N1}ϵ_N^N0.$$
(15)   
The equalities in (14),(15) only hold when no modification of the effective single particle wavefunctions occurs following the addition of an electron. In a disordered dot, in which there are no spatial symmetries, such a modification will always take place, and so $`\mathrm{\Delta }ϵ`$ can be considered strictly positive.

The difference $`\mathrm{\Delta }ϵ`$ provides a measure of the effectiveness of Koopmans’ theorem. To demonstrate this we present, in Fig.2, a schematic diagram of the surface of expectation values of the many-body Hamiltonian in the space of Slater determinants of $`N1`$, $`N`$, and $`N+1`$ particles. The SCHF ground states correspond to minima in these surfaces. From the diagram, it is clear that the energies $`ϵ_{N+1}^N`$ and $`ϵ_N^{N1}`$ are upper bounds to the respective addition energies, and the energies $`ϵ_{N+1}^{N+1}`$ and $`ϵ_N^N`$ are lower bounds to the addition energies. The approximation $`\mathrm{\Delta }_2^{k_1}`$, is therefore obtained by subtracting a lower bound ($`ϵ_N^N`$) from an upper bound ($`ϵ_{N+1}^N`$). As a result, the average value contains the average difference between the two bounds in addition to the correct mean $`\mathrm{\Delta }_2`$. It is generally assumed that the difference between the two bounds vanishes in the thermodynamic limit, and therefore so does $`\mathrm{\Delta }ϵ`$. We shall see that our results do not show any indication that this is the case. On the other hand, $`\mathrm{\Delta }_2^{k_2}`$ corresponds to the difference of two upper bounds to the two relevant addition energies. Regardless of the quality of the upper bound, so long as it is not strongly dependent on the number of particles present, both the particle number and disorder averaged results are good. The third approximation to $`\mathrm{\Delta }_2^{k_3}`$ corresponds to the difference of two lower bounds, and like $`\mathrm{\Delta }_2^{k_2}`$ is good in the mean. It is for this reason that we introduce these alternative approximations. It is easy to see that $`\mathrm{\Delta }_2^{k_1}(N)\mathrm{\Delta }_2^{k_3}(N)=\mathrm{\Delta }ϵ(N+1)`$ and therefore provides no further information. On the other hand, the fluctuations of $`\mathrm{\Delta }_2^{k_3}`$ can be different from those of $`\mathrm{\Delta }_2^{k_2}`$, and so are investigated separately. We note that in a clean system at $`r_s`$ below the Wigner crystal transition , the minima would align in Fig. 2, reflecting the validity of Koopmans’ theorem in that limit.

Let us briefly discuss the non-self-consistent single particle picture, for which the Koopmans’ approximations (11-13) and Eq.(10) all coincide:

$$\mathrm{\Delta }_2(N)=T_{N+1}T_N+\underset{j}{\overset{N1}{}}\left(V_{N+1j}V_{Nj}\right)+V_{N+1N}.$$
(16)   
Here, the term non-self-consistent approximation refers to a scheme where a set of effective single-particle states is given (e.g. by solving the N-electron SCHF problem), and utilised for any number of particles present in the system. The nearest neighbour spacings between levels that are both occupied or unoccupied has a similar form:

$$ϵ_{m+1}^Nϵ_m^N=T_{m+1}T_m+\underset{j}{\overset{N}{}}\left(V_{m+1j}V_{mj}\right),$$
(17)   
the major difference between (16) and (17) is the additional unbalanced matrix element $`V_{N+1N}`$ appearing in (16). Let us also suppose that in this simple single particle scheme the electrons interact with a short-ranged effective interaction. Blanter et al  introduce the hypothesis that the (normalised) spacings (17) and $`\mathrm{\Delta }_2V_{N+1N}`$ obey WD statistics up to corrections in $`1/g`$. Further assuming that the wavefunction correlations are still close to those of non-interacting particles leads, for the short-ranged effective interaction, to the result $`\mathrm{Var}(V_{ij})(U_0\mathrm{\Delta }/tg)^2`$ , so that the interaction dependent contribution to $`\delta \mathrm{\Delta }_2`$ scales like $`U_0\mathrm{\Delta }/tg`$ . This analysis is valid in the regime $`r_s1`$ and $`g1`$, implying that Koopmans’ theorem is a good approximation in that regime.

## IV Results and Discussion

In this section we present and discuss the results of the numerical simulations for both the nearest-neighbour and the Coulomb bare potentials. To make each subsection self-contained there is some repetition.

### A Short Range Interactions

We consider first the case of a nearest neighbour bare interaction potential as defined in Eq.(5). We begin by plotting the distributions of both the level spacings (17) and the gap (16) of the SCHF spectrum for finite $`r_s`$, $`g`$. In this case (16) and (17) are calculated in the self-consistent basis of $`N`$ particles. In Fig. 3 it is seen that the normalised level spacings between occupied states show an increasing deviation from WD to Poisson statistics as $`U_0`$ is increased. This is also true for the unoccupied states, but to a much greater extent. The difference between occupied and unoccupied states in the SCHF approximation will be discussed in greater detail later. We interpret the tendency towards Poisson statistics as a signature of the incipient localisation of the effective one particle states.

The normalised gap ($`\mathrm{\Delta }_2^{k_1}`$) distribution tends towards a more symmetric distribution that is approximately Gaussian as $`U_0`$ is increased.

We have also investigated the gap ($`\mathrm{\Delta }_2`$) distribution obtained within the fully self-consistent scheme, which we show in Fig. 4. We find that as $`U_0`$ is increased, the distribution evolves from a WD form to a more symmetric distribution similar to a Gaussian.

We shall concentrate first on the mean values of these distributions. A typical dependence of $`\mathrm{\Delta }_2`$ on the interaction parameter $`U_0`$ is plotted in figure 5. Whilst $`\mathrm{\Delta }_2^{k_2}`$ and $`\mathrm{\Delta }_2^{k_3}`$ provide a good approximation, we see a strong deviation of $`\mathrm{\Delta }_2^{k_1}`$. For all the system sizes considered, this effect occurs at $`r_s𝒪(1)`$. Results for the CI model, evaluated as described above, are plotted for comparison. Deviations of order $`𝒪(U_0/A)`$ from the CI model appear above $`U_02`$ ($`r_s1`$).

Elsewhere , we show that the ground state develops large density modulations as $`U_0`$ is increased beyond $`r_s𝒪(1).`$ These ground state charge density modulations (CDMs) explain the deviations of $`\mathrm{\Delta }_2`$ from the CI model prediction: they reduce the average addition energy by up to $`4U_0/At`$ when $`\nu <1/2`$ for a commensurate lattice. When $`\nu >1/2`$ the mean charging energy can be correspondingly increased.

We are also now in a position to understand why the level spacing statistics between unoccupied states show an increased tendency towards a Poisson distribution: the density modulations that appear in the ground state alter the potential felt by the unoccupied states. These modulations are not spatially ordered, as is demonstrated in Ref.. Hence, as $`U_0`$ is increased, the unoccupied states see an effective potential with increasingly strong modulations and tend to localise, whence the tendency towards a Poissonian distribution.

Let us consider the error in $`\mathrm{\Delta }_2^{k_1}`$ in more detail. As can be seen in Fig. 6, $`\mathrm{\Delta }ϵ/\mathrm{\Delta }`$ increases with the system size for lattices up to about 7\*8; for larger systems ($`A50`$) it seems that the error becomes proportional to $`\mathrm{\Delta }`$, and is of order $`\mathrm{\Delta }`$ when $`r_s`$ is of order unity.

We find that the nature of the disorder dependence of $`\mathrm{\Delta }ϵ`$ depends on the interaction strength, as seen in Fig. 7. The change in dependence occurs at interactions strengths corresponding to $`r_s1`$ for all the sample sizes considered. One might be surprised that deviations from Koopmans’ theorem do not smoothly decrease as $`W0`$ since, in the limit of vanishing disorder, Koopmans’ theorem becomes exact on the torus due to the restoration of translational symmetry. However, when $`W0`$, the spectrum develops many near degeneracies such that the effective perturbation due to $`U_0`$ is magnified as $`W0`$. In the limit $`r_s0`$, $`g1`$, the typical size of the matrix elements $`V_{ijkl}`$ which drive the rearrangement scale like $`\delta V_{ijkl}r_s\mathrm{\Delta }/g`$ , thus one expects that in this regime $`\mathrm{\Delta }ϵ`$ should increase with disorder. We find only a weak increase with disorder for $`r_s1`$.

In figure 8 we plot the interaction dependence of $`\mathrm{\Delta }ϵ/\mathrm{\Delta }`$. We find that at small $`U_0`$, $`\mathrm{\Delta }ϵ/\mathrm{\Delta }(U_0/t)^2`$, with deviations for larger $`U_0`$. In fact, this quadratic behaviour can be understood using second order perturbation theory. To see this we refer back to the schematic diagram of Fig. 2. A shift  occurs in the ground state configuration when a particle is added, which is represented by a misalignment of the minima. This shift is, to leading order, linear in $`U_0`$. Since the SCHF ground state energy is a minimum in the expectation value of the Hamiltonian, the difference of ground state energies for the two configurations will be quadratic in this shift. Furthermore, the local curvature tensor is independent of $`U_0`$ when the interaction matrix elements are small compared to the mean level spacing. Thus, $`\mathrm{\Delta }ϵ`$ scales like $`U_0^2`$ in the perturbative regime. The indication is that second order perturbation theory is qualitatively good even for $`r_s1`$. We note that since both the shift and the local curvature tensor depend on the disorder, there is no such simple $`W`$ dependence.

To summarise the results for the mean Coulomb gap, we find that Koopmans’ approximation (11) makes an error in the mean charging energy which for small $`r_s`$, $`A50`$, and fixed disorder scales like $`\mathrm{\Delta }ϵr_s^2`$. There is also evidence that for the larger sizes ($`A50`$), far beyond that accessible by exact diagonalisation, that $`\mathrm{\Delta }ϵr_s^2\mathrm{\Delta }`$. The latter dependence is consistent with the expectation that for sufficiently small $`r_s`$, $`1/g`$, perturbation theory is valid when the effective interaction is short-ranged. We find that $`\mathrm{\Delta }ϵ𝒪(\mathrm{\Delta })`$ when $`r_s𝒪(1)`$. To understand this result, we return to (16), (17), and fix the basis to be the self-consistent one for $`N`$ particles, so that (16) now describes $`\mathrm{\Delta }_2^{k1}`$. Since we have verified that the level spacings (17) show nearest neighbour separation statistics which are close to WD for all $`mN`$, with an approximately constant density of states, we are led to conclude that $`\mathrm{\Delta }ϵ`$ arises due to the fundamental difference between occupied and unoccupied levels in the SCHF approximation. In short, whilst $`V_{m+1j}^NV_{mj}^N`$ for $`mN`$ vanishes as expected, $`V_{N+1j}^NV_{Nj}^N`$ does not. Indeed $`\mathrm{\Delta }ϵ_j^{N1}V_{N+1j}^NV_{Nj}^N`$.

Turning now to the fluctuations in $`\mathrm{\Delta }_2`$, we plot an example result in Fig. 9. The most striking behaviour is the asymptotic saturation of the fluctuations (verified but not shown for even stronger interactions). As with the deviations of $`\mathrm{\Delta }_2`$ from the CI model, this occurs over the same range of interactions for all the system sizes considered, and is associated with the appearance of CDMs. Over the range of interaction strengths shown the fluctuations have not completely saturated, but the ground state density modulations are already present , such that, at low filling, the short-ranged contribution to the interaction energy is reduced. In the limit of strong interactions the charge segregates at a kinetic energy cost of order $`𝒪(t)`$, and $`U_0`$ plays no further role in the ground state energy fluctuations. The fluctuations therefore become sub-linear in $`U_0`$, and eventually saturate to an interaction independent value. Moreover, the results for the fluctuations become strongly geometry and filling factor dependent . We note that the observed saturation is in fact an artifact of the sharp cut-off in the interaction range: with a longer-ranged interaction, charge segregation cannot eliminate contributions due to the interaction (although it may significantly reduce them), and the fluctuations would no longer be bounded simply by kinetic energy considerations.

In Fig. 9 it can also be seen that for strong interactions, the fluctuations are overestimated by $`\mathrm{\Delta }_2^{k_1}`$ and $`\mathrm{\Delta }_2^{k_2}`$, and underestimated by $`\mathrm{\Delta }_2^{k_3}`$. This can be understood within the picture given above of charge density modulations. In this case, an occupied state that is removed non-self-consistently will yield less energy than can be gained when the system is allowed to reorganise, but the typical size of this error saturates to an interaction independent value for the same reason that the SCHF fluctuations do. If on the other hand an unoccupied state is occupied non-self-consistently, it is not possible to avoid contributions from the short-ranged part of the potential (we do not consider strongly Anderson localised states at very low filling), and the typical error increases indefinitely with the interaction strength. As a result $`\mathrm{\Delta }_2^{k_3}`$ underestimates the fluctuations by an amount that saturates to an interaction independent value, whereas fluctuations in the charging energy predicted by $`\mathrm{\Delta }_2^{k_1}`$ and $`\mathrm{\Delta }_2^{k_2}`$ grow with $`U_0`$ indefinitely; the errors made in employing the latter approximations diverge with the interaction strength.

Concentrating now on the fully self consistent results, the fluctuations in the charging energy are plotted against the sample size in figure 10. As expected, for very weak interactions, the typical fluctuations vanish like $`1/A`$, being dominated by kinetic energy fluctuations. For stronger interactions, this dependence no longer holds: for $`r_s1`$ our results are in broad agreement with Ref., but do not agree with their suggestion that the typical fluctuations remain proportional to $`\mathrm{\Delta }`$ for $`r_s>𝒪(1)`$. This appears to conflict with a simple single-parameter scaling argument . The appearance of fluctuations that do not scale with $`\mathrm{\Delta }`$ coincides with the appearance of density modulations. We stress that in this model there can be no physical connection between the amplitude of the constant interaction and the amplitude of the fluctuations.

For strong interactions the dominant disorder dependence appears to develop only for $`W4`$, where it is consistent with the emergence of a linear dependence to be expected from spatial rearrangements in the disorder potential. An example is plotted in Fig. 11. We reiterate that we also find strong geometry and filling factor dependences. It is extremely difficult to extract disorder scalings in such small systems because $`W/t`$ is required to be fairly large to generate diffusive motion, which in turn stretches the spectrum in the tails. This can be seen at weak interaction, where one would have hoped to see a disorder independent plateau in $`\mathrm{\Delta }`$ (i.e $`\delta \mathrm{\Delta }_2`$ at $`U_0=0`$).

To summarise then, we find $`\delta \mathrm{\Delta }_20.52\mathrm{\Delta }+ar_s\mathrm{\Delta }+𝒪(r_s^2)`$, where $`a`$ is an undetermined constant or function of disorder strength. We note that the disorder scaling is not clear because of the residual dependence of $`\mathrm{\Delta }`$ on $`W`$.

### B Long Range Interactions

We consider here the results for the Coulombic bare potential.

We first study the distributions of both the level spacings (17) and the gap (16) of the SCHF spectrum at finite $`r_s`$, $`g`$. In this case (16) and (17) are calculated in the self-consistent basis of $`N`$ particles. In Fig. 12 it is seen that the normalised level spacings between occupied states obey statistics very close to WD for all interaction strengths considered. Between occupied states (Fig. 12a) there is a mild deviation towards Poisson statistics for the strongest interaction strengths, indicative of a weak tendency towards localisation. Between unoccupied states (Fig. 12b) the distribution is even closer to WD. On the other hand, the normalised gap distribution clearly tends towards a more symmetric distribution that is approximately Gaussian.

We have also investigated the gap ($`\mathrm{\Delta }_2`$) distribution obtained within the fully self-consistent scheme, which we show in Fig. 13. Again we find that as $`U_0`$ is increased, the distribution evolves from a WD form to a symmetric distribution similar to a Gaussian.

We shall concentrate first on the mean values obtained from these distributions, and will come to the variance later in the section. Figure 14 shows a comparison of $`\mathrm{\Delta }_2`$, with the various approximations to it, plotted against $`U_0`$. Whilst $`\mathrm{\Delta }_2^{k_2}`$ and $`\mathrm{\Delta }_2^{k_3}`$ provide a good approximation to $`\mathrm{\Delta }_2`$, we see a clear deviation of $`\mathrm{\Delta }_2^{k_1}`$. Results for the CI model, evaluated as described above, are plotted for comparison. That the CI model is good in the mean indicates that the single particle wavefunctions remain roughly uniformly distributed over the dot for all $`r_s`$ considered.

In figure 15 we plot $`\mathrm{\Delta }ϵ/\mathrm{\Delta }_2`$ against the sample area $`A`$, for an intermediate disorder strength ($`W=4`$). Since we know from Fig.14 that $`\mathrm{\Delta }_2^{k_2}\mathrm{\Delta }_2`$, then $`\mathrm{\Delta }ϵ`$ is very close to $`\mathrm{\Delta }_2^{k_1}\mathrm{\Delta }_2`$, the total error made by applying Koopmans’ theorem. For $`A50`$, we find that for a fixed interaction strength $`\mathrm{\Delta }ϵL`$ ($`\mathrm{\Delta }ϵ/\mathrm{\Delta }_2A`$). For larger samples with $`U_02`$ we find a weakening in the dependence, but see no indication that it will vanish relative to $`\mathrm{\Delta }`$. The result that the deviations from Koopmans’ approximation increase with system size (when compared to $`\mathrm{\Delta }`$), showing no sign of saturation, is admittedly strange, and may be an artifact of the specific model considered here. However, the result that Koopmans’ approximation appears to fail even as the system size tends towards the thermodynamic limit, is in line with our findings for the short-ranged case .

It is interesting to see how the error depends on disorder. In figure 16 we plot $`\mathrm{\Delta }ϵ/t`$ for a range of disorder strengths: the disorder dependence as a function of interaction is weak, but not simple. Similarly to the short-ranged case, the deviations from Koopmans’ approximation do not decrease for small disorder. This occurs for the same reasons as for the short-ranged case. Here too the typical size of the matrix elements $`V_{ijkl}`$ which drive the rearrangement scale inversely with $`g`$. One thus expects that in this regime $`\mathrm{\Delta }ϵ`$ should increase with disorder, and this is indeed seen in the figure. For $`U_04t`$, $`\mathrm{\Delta }ϵ`$ decreases with disorder at sufficiently large $`W`$, with evidence of a turning point ($`d\mathrm{\Delta }ϵ/dW=0`$) at $`U_0=4t`$, $`W4t`$.

In figure 17 we plot the interaction dependence of $`\mathrm{\Delta }ϵ/\mathrm{\Delta }`$. We find that at small $`U_0`$, $`\mathrm{\Delta }ϵ/\mathrm{\Delta }(U_0/t)^2`$, with deviations for larger $`U_0`$. This quadratic behaviour has the same origin as that of the short-ranged case: the indication is that second order perturbation theory is qualitatively good even for $`r_s1`$.

To summarise the results for the mean Coulomb gap, we find that Koopmans’ approximation (11) makes an error in the mean charging energy which for small $`r_s`$ and $`L`$, and fixed disorder scales like $`\mathrm{\Delta }ϵr_s^2L`$. There is also evidence that for the larger sizes ($`A50`$), far beyond that accessible by exact diagonalisation, that the size dependence vanishes: $`\mathrm{\Delta }ϵr_s^2`$. In contrast to the naive expectation however, we find no sign of this error vanishing relative to $`\mathrm{\Delta }`$ in the thermodynamic limit. This is due to the fundamental difference between occupied and unoccupied SCHF levels already discussed in the short-ranged case.

We now consider the fluctuations in $`\mathrm{\Delta }_2`$. As an example of the interaction dependence of these fluctuations in the various approximation schemes, we plot the results for a fixed size in figure 18. It is seen that applying Koopmans’ theorem in the forms (11-13) results in considerably smaller fluctuations than the fully self-consistent calculation. To quantify this error, we plot $`\delta \mathrm{\Delta }_2^{k_2}/\delta \mathrm{\Delta }_2`$ in Fig.19, which shows that the relative error initially increases with interaction strength, but shows signs of saturating. The value of the saturation appears to increase towards unity as the system size is increased.

We now concentrate on the fluctuations of the fully self-consistent peak spacing $`\mathrm{\Delta }_2`$. For comparison with Ref. it is useful to plot $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2`$ against the interaction strength for a range of sample sizes. This is done in figure 20. In the inset we plot $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }`$ which shows that the peak spacing fluctuations are not proportional to $`\mathrm{\Delta }`$ for $`r_sL𝒪(1)`$ . From Fig. 20 it can be seen that the curves $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2`$ do not saturate to a constant as suggested in Ref., although to see this clearly one has to consider larger sample sizes than are accessible by exact calculations. The curve $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2`$ appears to take on the approximate form of a constant term plus a linear term for $`r_sL𝒪(1)`$. The constant contribution identified by Sivan et el. , is here, contrary to their claim, non-universal (i.e. it is disorder dependent).

In figure 21 we plot $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2`$ against disorder for the $`10\times 11`$ lattice with a range of interaction strengths. At $`U_0=0`$ it is seen that for $`W6`$ the systems obeys WD statistics quite well. For the sample size considered we find that in the regime $`0.5U_06.0`$ ($`0.25r_s3.0`$) $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2W`$, and at stronger interactions this dependence weakens. The intermediate dependence, $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2W`$, is consistent with the dependence $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_21/\sqrt{g}`$ recently observed independently by Bonci and Berkovits  for the Buminovich Stadium billiard. Analysis of Fig. 20 leads to the conclusion that the quadratic contribution (in $`U_0`$) to $`\delta \mathrm{\Delta }_2`$ is independent of disorder, which is consistent with Fig. 21.

To identify the system size scaling of the various contributions we plot $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }`$ against $`A`$ in figure 22. For $`U_0=0`$ the system obeys WD statistics and $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }`$ is independent of size. The regime over which the fluctuations are approximately proportional to the mean charging energy (the constant contribution to $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }_2`$ alluded to above, which corresponds to a $`\sqrt{A}`$ dependence in the figure), depends on the system size. As could be seen in Fig. 20 another term begins to dominate the the fluctuations at larger $`U_0`$ which is quadratic in $`U_0`$, this term increases more rapidly with the system size, and so dominates at lower $`U_0`$ in larger systems. Over the range of sizes considered, this term appears to scale like $`L^2`$. The cross-over in dominance therefore occurs at $`U_01/L`$ for fixed disorder strength; clearly the quadratic term will dominate in large samples. The increase in $`\delta \mathrm{\Delta }_2/\mathrm{\Delta }`$ with system size appears to be an artifact of using the unscreened Coulomb interaction in the Hamiltonian.

Summarising the results presented in figures 20 to 22, and the above discussion, we find for $`r_sL𝒪(1)`$ an approximate form: $`\delta \mathrm{\Delta }_2.52\mathrm{\Delta }+a\mathrm{\Delta }_2/\sqrt{g}+br_s^2`$ where $`a,b`$ are constants. One would normally expect the fluctuations to be linear in the interaction strength, (i.e. $`b=0`$). A possible source for such a quadratic interaction dependence in the typical fluctuations is the development of correlations that grow like $`r_s^2`$ in products of eight wavefunctions. Elsewhere  we present evidence for increased fluctuations in the ground state density in this regime, as compared to a non-interacting system. It is not yet clear whether this result is an artifact of the SCHF approximation (which has also very recently been observed in $`1d`$ systems  using a similar approximation scheme), or a genuine physical effect.

Finally, it is worth noting that since the exchange interaction is not correctly screened, that errors in the SCHF scheme might be expected to diverge with respect to $`\mathrm{\Delta }`$ as the system size is increased. We can neither confirm nor counter this argument, but have verified that for $`r_s5`$ the fluctuations in the exchange contribution are smaller than those of the direct contribution.

## V Summary

We have investigated the addition spectra of disordered quantum dots employing an effective single particle approximation, both using a fully self consistent analysis, and by invoking Koopmans’ theorem. We were able to consider system sizes with up to 144 sites, and 37 particles, compared to the latest exact calculations on samples with 24 sites and 6 particles . The larger sample size also allows us to consider smaller values of $`g`$ than exact calculations whilst retaining an ergodic non-interacting limit, and therefore approaches the experimental parameters more closely. The inclusion of spin in a consistent manner is left for a future project.

Our SCHF results for the typical fluctuations of the peak spacings for particles possessing short-ranged bare interactions are entirely different from the results for long-ranged bare interactions. In the short-ranged case we find the same scaling as Ref. for very weak interactions ($`r_s1`$), but for $`r_s1`$ deviations from this behaviour become significant and coincide with the appearance of interaction induced density modulations. We find no size dependence in the onset of these effects. We find that strong filling factor and geometry dependences arise due to these density fluctuations, and therefore do not expect that the disorder ensemble statistics can be mapped to statistics over the ensemble of filling factors: ergodicity is lost. We suggest that employing a short ranged bare interaction in a self-consistent scheme is not an appropriate model for the quantum dots of Refs. for which $`r_s>1`$, but may be a useful model for dot geometries sandwiched between very close metallic gates which provide a good external source of screening . In this respect we identify some experiments on the addition spectrum which possess a metallic source (heavily doped $`n^+`$ GaAs) and drain (Cr/Au), at separations of the order of the average inter-particle spacing in the dot .

In the Coulomb case the SCHF approximation to $`\mathrm{\Delta }_2`$ yields typical fluctuations that do not scale with $`\mathrm{\Delta }`$ for $`r_s1/L`$. In Ref. it is claimed that $`\delta \mathrm{\Delta }_2`$ is universally proportional to $`\mathrm{\Delta }_2`$ for strong interactions (but still far from the accepted Wigner Crystal transition point). In contrast, we find, in addition to the small interaction independent contribution, a contribution to $`\delta \mathrm{\Delta }_2`$ that is proportional to $`\mathrm{\Delta }_2/\sqrt{g}`$ (i.e. non-universal), and a further contribution that scales like $`r_s^2`$, which is independent of disorder, and appears to be due to the development of charge density modulations . The latter is not detectable in the small systems examined numerically in , and so our results are not numerically inconsistent with exact calculations . Whilst we do not include spin, the observed decrease in the fluctuations with $`g`$ is consistent with the experimental indications  that in cleaner samples the fluctuations are smaller.

We show that a direct application of Koopmans’ theorem overestimates $`\mathrm{\Delta }_2`$. This overestimate, a manifestation of the breakdown of Koopmans’ approximation, does not vanish on the scale of $`\mathrm{\Delta }`$ in the thermodynamic limit. The error seems to scale differently with sample size for sample areas above or below $`A50`$. In the nearest neighbour case, with $`A50`$, this error scales with $`\mathrm{\Delta }`$, but in smaller systems, accessible by exact methods, it is independent of system size. In the Coulomb case the error grows with the system size as $`L`$ for $`A50`$. For larger sizes the error appears to tend towards a $`1/L`$ scaling, i.e. in proportion with the charging energy, and therefore diverges with respect to the mean effective single-particle level spacing. This result for the Coulomb interaction case appears to be non-physical, and may be an artifact of the model considered. However, the result that Koopmans’ theorem is not recovered in the thermodynamic limit also occurs in the short-ranged interaction case. In both cases we find that initially this error grows in proportion to $`U_0^2`$, to be expected since the lowest order contribution is second order, but for strong interactions it grows more slowly in $`U_0`$, and that the disorder dependence of this error is weak and non-monotonic. We identify the source of the error $`\mathrm{\Delta }ϵ`$ to be the fundamental difference between occupied and unoccupied states that is inherent in the SCHF approximation. We introduce two improved applications of Koopmans’ theorem, $`\mathrm{\Delta }_2^{k_2}`$, $`\mathrm{\Delta }_2^{k_3}`$, which provide a good approximation to $`\mathrm{\Delta }_2`$, but not to $`\delta \mathrm{\Delta }_2`$.

Whilst preparing the manuscript, two related works appeared that confirm some of the points discussed above .

In both cases, fluctuations in the ground state density develop with $`r_s`$ , and have significant effects of the addition spectrum statistics. It remains to be seen whether these density modulations are an artifact of the SCHF approximation (i.e. due the neglect of dynamical correlations), or in fact interesting results on the continuous transition to a Wigner-type solid in disordered samples with short- and long-ranged bare interactions.

###### Acknowledgements.

We acknowledge discussions with H. Orland and F. von Oppen in the early stages of this project, as well as with Ya. Blanter, S. Levit, A. Mirlin and D. Orgad. We acknowledge support from the EU TMR fellowship ERBFMICT961202, the German-Israeli Foundation, the U.S.-Israel Binational-Science Foundation and the Minerva Foundation. One of us (YG) would also like to acknowledge support from an EPSRC senior professorial fellowship, grant number GR/L67103. Much of the numerical work was performed using IDRIS facilities.