text:
# Interactions in Chaotic Nanoparticles: Fluctuations in Coulomb Blockade Peak Spacings

## I Introduction

The ground state properties of electrons in nanoparticles or quantum dots—finite systems of fermions interacting via Coulomb forces—remain incompletely understood. One valuable way to probe these systems is via electrical transport (for reviews see Refs. ). The dominant effect is the suppression of the conductance $`G`$ because of the charging energy needed to add an electron to the nanoparticle and so allow current to flow. This Coulomb blockade (CB) is an essentially classical effect, and the charging energy is simply $`E_C=e^2/2C`$ where $`C`$ is the capacitance of the nanoparticle. The blockade can be lifted by adjusting the potential on the particle—with a gate voltage $`V_g`$, for instance—so that the energy for $`N`$ electrons is the same as that for $`N+1`$. The conductance through the nanoparticle has peaks at these degeneracy points; in fact, $`G(V_g)`$ is periodic with a spacing of $`e^2/C`$ between the peaks .

At low temperature the electron remains coherent in the dot and so quantum interference effects must be included in the Coulomb blockade . Two main effects occur: (1) the peak heights fluctuate from peak to peak because of varying coupling between the wavefunction in the nanoparticle and the leads, and (2) the spacings between the peaks fluctuate because of interference contributions to the ground state energy. Peak heights have been studied first theoretically  and then experimentally ; agreement is good once effects of temperature  and classical dynamics  are taken into account.

In this paper we focus on quantum effects in the peak spacings. These are directly related to the interference contributions to the ground state energy. One obvious contribution is the single-particle energy levels caused by the confinement in the quantum dot. But there may also be contributions coming from the interactions in the dot. We will show that for typical parameters characterizing experimental quantum dots or nanoparticles the contributions coming from the interactions must be included.

The simplest treatment of Coulomb blockade peak spacing fluctuations results from writing the ground state energy as the sum of the classical electrostatic energy and the energies of the occupied single particle states . Such an expression for the ground state energy of a many Fermion system has a long and venerable history in atomic and nuclear physics . Using $``$ for the energy in the presence of the gate voltage while reserving $`E`$ for $`V_g=0`$, we have for $`N`$ particles

$`_N(V_g)`$ $`=`$ $`E_N(C_g/C)eNV_g`$ (1)
$`E_N`$ $`=`$ $`(e^2/2C)N^2+{\displaystyle \underset{\mathrm{occupied}i\sigma }{}}ϵ_{i\sigma }`$ (2)
where $`C_g`$ is the capacitance of the gate to the dot and $`\{ϵ_{i\sigma }\}`$ are the single-particle energy levels. This treatment, known as the constant-interaction model (CI), is still an important point of reference in the field. The gate voltage $`V_g^{}`$ at which a peak in the conductance occurs is given by the condition $`_N(V_g^{})=_{N+1}(V_g^{})`$. The spacing between two peaks is, then, proportional to the (discrete) inverse compressibility $`\mathrm{\Lambda }`$

$`(V_g^{})_{NN+1}(V_g^{})_{N1N}`$ (4)
$`E_{N+1}+E_{N1}2E_N\mathrm{\Lambda }_N.`$
In the case of the CI model, one gets the simple prediction

$`\mathrm{\Lambda }_N=e^2/C`$ for odd $`N`$ (5)
$`\mathrm{\Lambda }_N=e^2/C+(ϵ_{N/2+1}ϵ_{N/2})`$ $`\text{for even }N.`$ (6)
Note the drastic odd/even effect in the quantum correction to the classical result $`e^2/C`$ because of the spin of the electron.

A number of experiments have studied the distribution of CB peak spacings in semiconductor quantum dots , a particularly flexible and convenient type of nanoparticle. The overall result has been striking disagreement with the CI model. The main discrepancy is that there is no sign of the strong odd/even effect predicted by Eq. (5), though there is evidence of a weak odd/even effect . In addition, in some of the experiments the widths of the observed distributions are significantly larger than what is predicted using the CI model and assuming the single-particle energy levels $`ϵ_i`$ are distributed according to the appropriate Gaussian random matrix ensemble. These discrepancies indicate unambiguously that interaction effects play a significant role in the fluctuating part of the Coulomb Blockade peak position.

Theoretically accounting for this disagreement is difficult because of two intrinsic characteristics of the problem: (1) the fairly large number of electrons contained in the dot (typically $`N100`$-$`1000`$), and (2) the long range character of the Coulomb interaction which requires that screening has to be correctly taken into account. The first of these points makes it delicate to extend to realistic dots the results obtained by exact multiparticle diagonalization , which can only treat problems with a few electrons ($`\stackrel{<}{}15`$). To reach the experimental number of particles in a numerical approach, one has to rely on some self-consistent technique. One approach is Hartree-Fock which has been recently used in quantum dots . Then, however, the second point comes into play. Indeed, screening is a higher-order interaction effect, and it is well known for the Coulomb interaction that the exchange terms can be significantly over evaluated in self-consistent Hartree-Fock (for the bulk see Ref.  but this applies also for restricted geometries). This makes it rather difficult to interpret numerical self-consistent Hartree-Fock calculations using the bare Coulomb interaction as input. Another possible self-consistent approach is, of course, density functional theory in the local density or local spin density approximation . Here the screening is properly handled but approximations are made in the exchange term. Since, as we explain in detail in this paper, the subtle interplay between exchange, residual direct interactions, and single-particle level fluctuations is critical in describing quantum dot physics, results from this method are also difficult to generalize.

An alternative theoretical approach is to formally carry out an RPA screening calculation and then use a short-range screened interaction in the subsequent evaluation of the magnitude of various effects . The nature of the screened interaction used is, then, important. Since experimentally the system size is much larger than the screening length, the most natural choice for the interaction is the bulk RPA screened interaction. For momentum smaller than $`2k_F`$ this is correctly approximated by the Thomas-Fermi expression. Throughout this paper we consider two dimensional systems, in which case

$`V_{\mathrm{TF}}(𝐫)`$ $`=`$ $`{\displaystyle \frac{d𝐪}{(2\pi )^2}\widehat{V}_{\mathrm{TF}}(𝐪)\mathrm{exp}[i𝐪𝐫]},`$ (7)
$`\widehat{V}_{\mathrm{TF}}(𝐪)`$ $`=`$ $`{\displaystyle \frac{\nu ^1}{1+r_s^1(|𝐪|/k_F)/\sqrt{2}}}`$ (8)
where $`\nu `$ is the mean density of states (including spin degeneracy), and $`r_s=r_0/a_0`$ with $`\pi r_0^2`$ the average area per electron and $`a_0`$ the Bohr radius in the semiconductor. It should be born in mind that although $`\widehat{V}_{\mathrm{TF}}`$ is proportional to $`r_s`$ for small $`r_s`$, it becomes largely independent of $`r_s`$ when $`r_s\stackrel{>}{}1`$, which is the case experimentally. In fact, within any RPA-like approach the interaction cannot become stronger than

$$V_{\text{zero-range}}(𝐫)=\nu ^1\delta (𝐫)$$
(9)
which we refer to as the zero-range limit of the potential. Once screening is included, the problem is that the interaction is small, in a sense that we shall make more explicit below, even in the strongest zero-range limit. Thus the self-consistent approach should give roughly the same answer as the CI model with first-order perturbation in the screened interaction added. This is exactly the picture one has of weakly interacting Landau quasi-particles, with the obvious drawback that it is not straightforward to extract large contributions using this approach (see e.g. Refs. ).

In order to obtain a large interaction effect, one possibility is to give up the Fermi-liquid-like approach and to assume, for instance, that for some reason screening is not as efficient in a nanoparticle as in the bulk, and so conclude that large zero-range interactions, not amenable to a perturbative treatment, should be used . Before taking this drastic step, it seems reasonable, however, to see how far one can go within the more standard Fermi liquid description. Indeed, even in this framework, we know that the CI model is inadequate.

An accurate treatment of the spacing or inverse compressibility certainly requires that the average effect of the residual interactions be added to the CI model. By “residual interactions” we mean interactions beyond the simple classical effect taken into account in the constant interaction model and, in particular, interactions between the quasiparticles after the bare electrons are dressed through screening. While they make a smaller contribution to the total energy than those considered in Eq. (1), residual interactions make a contribution of order the mean level separation to the second difference $`\mathrm{\Lambda }`$, comparable to the single-particle contribution in Eq. (5) . In the limit of a large nanoparticle—one whose typical dimension $`L`$ is many times the electron wavelength, $`k_FL1`$—only the average effect of residual interactions needs to be added to the CI model . We shall below refer to this approximation as the constant exchange model. The intuitive argument for this is that the effect of residual interactions involves integrating wavefunctions over the entire volume of the nanoparticle and so self-averages: the particular characteristics of any given wavefunction are simply integrated out. Even if one neglects the fluctuations of the residual interaction terms, because of fluctuations of the single-particle levels, the “constant exchange” contribution can modify the total spin of the nanoparticle by favoring non-trivial occupation of the single-particle orbitals and, as a consequence, will significantly modify the fluctuations of the peak spacings . The ground state spin of a large nanoparticle in this regime has been examined in detail , but curiously the corresponding distribution of Coulomb blockade peak spacings has not appeared explicitly (see however Ref. ). Therefore, before turning to effects caused by fluctuations in the interactions—our main interest in this paper—we present below the distribution in the constant exchange limit ($`k_FL1`$) to use as a reference point.

In the experiments, however, $`k_FL`$ is not very large because of the need to have the mean level separation larger than the temperature. This suggests that the variations in interactions caused by properties of individual wavefunctions are important. The main goal of this paper is, therefore, to study the effect of fluctuations in the electron-electron interactions on the distribution of CB peak spacings.

We limit our study here to the case of chaotic dynamics within the dot, for which a random matrix description of the non-interacting limit can be used . Much as in the case of diffusive transport, completely chaotic dynamics introduces an important additional energy scale (and only one): the Thouless energy, $`E_{\mathrm{Th}}=\mathrm{}/\tau `$, is the inverse of the transit time across the dot. It is generally larger than the mean level separation $`\mathrm{\Delta }`$ but smaller than the Fermi energy; we will assume $`\mathrm{\Delta }E_{\mathrm{Th}}E_F`$. On energy scales less than $`E_{\mathrm{Th}}`$, a universal random matrix description holds while for larger energies (shorter times) the individual dynamics of the system comes into play. We will see that taking account of this additional energy scale is critical in describing CB peak spacings in nanoparticles.

We focus on two effects in particular: (1) the fluctuation in the residual Hartree and Fock contributions to the ground state energy and (2) the change in the single-particle energies because of changes in the mean field potential as electrons are added to the dot , an effect we refer to as “scrambling”. For two dimensional quantum dots, the changes in energy are evaluated to leading order in the smallness of the dot, the constant exchange model being used as the starting and reference point. We find that both effects contribute terms which are of order $`\mathrm{\Delta }/\sqrt{k_FL}`$, up to logarithmic corrections. The resulting distribution of CB peak spacings looks more like the experiment—it has a more Gaussian and symmetric form, for instance—but still deviates from the measurements in significant ways.

The paper is organized as follows. Our starting point is a semiclassical expansion of the ground state energy described in Section II—an expansion in $`1/k_FL`$. The important large dot limit is presented in Section III. The next two sections discuss the two main issues—fluctuations of the residual interaction contributions (Section IV) and scrambling (Section V). Section VI takes a first step beyond the Gaussian model introduced in Section IV. Finally, our conclusions, including a discussion of the experimental results, appear in Section VII.

## II Approach: Semiclassical Corrections

In the Coulomb blockade through quantum dots we are faced with a classical theory which works remarkably well to which we want to add quantum effects to leading order. The small parameter is the standard semiclassical one, $`\mathrm{}`$. In our context, this is equivalent to an expansion in $`1/k_FL`$ where $`k_F`$ is the Fermi wavevector and $`L`$ is the typical size of the nanoparticle. Another useful way to view the corrections is in terms of the dimensionless conductance, $`g`$, given as the ratio of the transit rate $`\mathrm{}v_F/L`$ to the mean level separation $`\mathrm{\Delta }`$; in this case the expansion is in terms of $`1/g`$, as in the diffusive mesoscopic regime . For 2-dimensional systems, such as the ones we consider here, these two parameters are proportional ($`g=k_FL/2\pi `$).

We proceed by using the method of V. M. Strutinsky  in which the dependence of many-body ground state quantities on particle number can be decomposed into an average and a fluctuating part. While the average part varies smoothly with particle number, the fluctuating part reflects the shell structure of the system. The smooth part is the bulk energy per unit volume integrated over the finite-size system, and the oscillating contributions come from quantum interference effects explicitly caused by the confinement. By supposing that the smooth part is known while the unknown oscillatory contribution is a correction, Strutinsky introduced a physically motivated systematic approach to obtain the shell corrections . This shell correction method is essentially a semiclassical approximation. It rests on the fact that the number of particles in the system considered is large, rather than on the interaction between the particles being weak. (One must, of course, work in a regime where the smooth starting point is basically valid.) Since the quantum dots in which we are interested contain on the order of 100 electrons, they are a perfect place to apply the Strutinsky method. We use the formulation recently developed specifically for quantum dots in Ref. .

Density functional theory guarantees that the ground state energy of a nanoparticle can be written as a functional of its density, $`E_{\mathrm{DFT}}=_{\mathrm{DFT}}[n_{\mathrm{DFT}}]`$. Neglecting quantum interference effects in this functional corresponds to a generalized Thomas-Fermi approximation (generalized because bulk local exchange and correlation can be included) which can also be written as a functional of a (smooth) density, $`E_{\mathrm{GTF}}=_{\mathrm{GTF}}[n_{\mathrm{GTF}}]`$. The quantum interference part of the energy, $`E_{\mathrm{DFT}}E_{\mathrm{GTF}}`$, can be obtained approximately by expanding $`n_{\mathrm{DFT}}`$ about $`n_{\mathrm{GTF}}`$ and solving the resulting density functional equations order by order in the oscillating part of the density $`n_{\mathrm{DFT}}n_{\mathrm{GTF}}`$. To carry out this procedure an explicit form for the exchange-correlation functional is required. Based on the local density approximation, the result derived in Ref.  is

$$E_{\mathrm{DFT}}E_{\mathrm{GTF}}+\mathrm{\Delta }E^{(1)}+\mathrm{\Delta }E^{(2)}$$
(10)
where the first and second order correction terms are

$`\mathrm{\Delta }E^{(1)}`$ $`=`$ $`_{1\mathrm{p}}^{\mathrm{osc}}\left[V_{\mathrm{eff}}[n_{\mathrm{GTF}}]\right]{\displaystyle \underset{i=1}{\overset{N}{}}}\delta ϵ_i`$ (11)
$`\mathrm{\Delta }E^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\stackrel{~}{n}^{\mathrm{osc}}(𝐫)V_{\mathrm{sc}}(𝐫,𝐫^{})\stackrel{~}{n}^{\mathrm{osc}}(𝐫^{})}.`$ (12)
The first-order correction is the oscillatory part of the single-particle energy, $`_{1\mathrm{p}}^{\mathrm{osc}}`$, calculated in the smooth Thomas-Fermi potential, $`V_{\mathrm{eff}}[n_{\mathrm{GTF}}]`$; that is, it is the sum of the deviations $`\delta ϵ_i`$ of the single-particle levels $`ϵ_i`$ from their mean values. In the second-order correction, $`\stackrel{~}{n}^{\mathrm{osc}}`$ is the deviation of the quantum mechanical density calculated in the Thomas-Fermi potential from the smooth starting point, $`\stackrel{~}{n}^{\mathrm{osc}}n\left[V_{\mathrm{eff}}[n_{\mathrm{GTF}}]\right](𝐫)n_{\mathrm{GTF}}`$, and $`V_{\mathrm{sc}}(𝐫,𝐫^{})`$ is the screened interaction in the nanoparticle. To the order that we are working, using $`\stackrel{~}{n}^{\mathrm{osc}}`$ here is equivalent to using $`n_{\mathrm{DFT}}n_{\mathrm{GTF}}`$ .

In this approach, the ground state energy is, then, the sum of a classical contribution—the generalized-Thomas-Fermi result $`E_{\mathrm{GTF}}`$—and two quantum contributions—$`\mathrm{\Delta }E^{(1)}`$ and $`\mathrm{\Delta }E^{(2)}`$. With the inclusion of only the first-order contribution, Eq. (10) is simply the CI model for the ground state energy. Because we have the corrections to this model, we can see exactly when it is applicable. The second-order correction, Eq. (12), has a natural interpretation: the ripples in the density caused by interference (the Friedel oscillations in this setting) interact with each other via an interaction screened by the smooth Thomas-Fermi fluid already present. In fact, the bulk Thomas-Fermi interaction Eq. (8) can often be used for $`V_{\mathrm{sc}}`$. Taking the $`r_s\mathrm{}`$ limit formally yields the zero-range interaction limit (9). This is equivalent to using a Hubbard Hamiltonian on a square lattice with hopping term $`t`$ and interaction $`U`$ with $`U/t=\pi `$.

The result Eq. (10) is, however, inadequate for our purposes here because of the local density approximation used. In particular, the spin degree of freedom is handled poorly in this approximation: while Eq. (12) gives the “direct” interaction between the density ripples, an analogous “exchange” term is missed. A better treatment of the interactions in density functional theory could presumably yield this additional term, by including off-diagonal terms in the density matrix, for instance. In fact, an RPA treatment of the interaction in the diffusive case yields an energy with exactly this structure . Motivated by these physical considerations, we simply add the exchange term by hand to the expression for the energy $`E_{\mathrm{DFT}}`$.

For the purposes of deriving CB peak spacing distributions, we thus consider 3 contributions to the ground state energy: (1) the Thomas-Fermi energy, (2) the deviation of the single-particle energy from its mean, and (3) a Hartree-Fock like treatment of the residual screened interaction. Without the residual interactions, a $`N`$-particle state is, of course, a Slater determinant in the basis of single-particle eigenstates $`\psi _i(𝐫)`$; it is characterized by the occupation numbers $`(n_{0\sigma },n_{1\sigma },\mathrm{})`$ where $`n_{i\sigma }=0`$ or $`1`$, $`\sigma =+1`$ or $`1`$ labels the spin degree of freedom, and $`n_{i\sigma }=N`$. In this occupation number representation, the expression for the energy is, then,

$`_{\mathrm{tot}}(\{n_{i\sigma }\})=E_{\mathrm{GTF}}+{\displaystyle \underset{i\sigma }{}}n_{i\sigma }\delta ϵ_i`$ (14)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i\sigma ,j\sigma ^{}}{}}n_{i\sigma }M_{ij}n_{j\sigma ^{}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j;\sigma }{}}n_{i\sigma }N_{ij}n_{j\sigma }`$
where

$`M_{ij}`$ $``$ $`{\displaystyle 𝑑𝐫𝑑𝐫^{}\left|\psi _i(𝐫)\right|^2V_{\mathrm{TF}}(𝐫𝐫^{})\left|\psi _j(𝐫^{})\right|^2}`$ (15)
$`N_{ij}`$ $``$ $`{\displaystyle 𝑑𝐫𝑑𝐫^{}\psi _i(𝐫)\psi _j^{}(𝐫)V_{\mathrm{TF}}(𝐫𝐫^{})\psi _j(𝐫^{})\psi _i^{}(𝐫^{})}`$ (16)
are the direct and exchange contributions, respectively. This is the starting point for our study of CB peak spacings.

In order to derive the statistical properties of the CB peak spacings from Eq. (14), the statistical properties of the single-particle eigenvalues and eigenfunctions must be known. For this purpose we assume that the single-particle classical dynamics in the nanoparticle is completely chaotic. In this case it is well-known that the single-particle quantum mechanics is accurately described by random matrix theory (RMT) on an energy scale smaller than the Thouless energy $`E_{\mathrm{Th}}`$ (which, again, is the inverse time of flight across the system) . We consider only the case of no symmetry here, so that the energy levels $`\{ϵ_i\}`$ obey Gaussian Unitary Ensemble (GUE) statistics. Because of the spatial integrals in the expressions for $`M_{ij}`$ and $`N_{ij}`$, the correlations of the wavefunctions are also needed. For the classical Gaussian ensembles, there is no spatial correlation for a given wavefunction, except for that arising from the normalization. This is a direct consequence of the fact that no basis should play a particular role for these matrix ensembles. Here, however, because of the particular role played by the kinetic energy, the plane wave basis is special. More precisely, we should implement that each eigenfunction is localized in this basis on the energy scale $`E_{\mathrm{Th}}`$. We shall therefore build our random matrix ensemble with the following requirements (a practical implementation is given in section VI):

* On energy scales smaller than the Thouless energy $`E_{\mathrm{Th}}`$, one should recover standard GUE properties.
* On scales larger than $`E_{\mathrm{Th}}`$, an eigenstate of energy $`ϵ`$ should appear localized on the energy surface $`|𝐤|=\sqrt{2mϵ}/\mathrm{}`$. More precisely, we require that in a plane wave basis an eigenstate has an envelope of width

$$\delta k=1/L$$
(17)
where $`L`$ is the typical size of the system.
* Finally, we have to implement that the fluctuations of energy levels, and also importantly of local density, saturate at the Thouless energy. That is, if one considers a quantity $`A`$ expressed as a sum over the energy levels in an energy window $`\delta E`$ about the Fermi energy, as the window widens the fluctuations of $`A`$ increase until $`\delta EE_{\mathrm{Th}}`$ beyond which point the fluctuations of $`A`$ do not further increase.

Except for the fact that the Fermi surface is given an explicit width $`E_{\mathrm{Th}}`$ due to the finite size $`L`$ of the system, items i) and ii) are very similar to Berry’s modeling of chaotic wave functions  which assumes that a random superposition of plane waves describes the correlation on scales smaller than the size of the system $`L`$. Item iii) is motivated by the known saturation of, for instance, fluctuations in the density of states of chaotic systems at an energy scale of order $`\mathrm{}`$ divided by the period of the shortest periodic orbit . Note that these requirements will yield a random matrix approach noticeably different from the ones discussed in Refs. .

In this model, the wave function statistics that we find are as follows: for a two dimensional dot of area $`A`$ (the case that we consider throughout this paper)

$$A\psi _i(𝐫)\psi _j^{}(𝐫^{})=\delta _{ij}J_0(k_F|𝐫𝐫^{}|)$$
(18)
for $`|𝐫𝐫^{}|L`$. (In a maximum entropy approach semiclassically restricted by the classical dynamics, this corresponds to keeping only the direct trajectory between $`𝐫`$ and $`𝐫^{}`$; the approach can be generalized to keep more classical paths .) In addition, it can be shown that within our random matrix modeling, one has, up to $`1/(k_FL)`$ corrections

$`\psi _i(𝐫_1)\psi _i^{}(𝐫_2)\psi _i(𝐫_3)\psi _i^{}(𝐫_4)=`$ (21)
$`\psi _i(𝐫_1)\psi _i^{}(𝐫_2)\psi _i(𝐫_3)\psi _i^{}(𝐫_4)`$
$`+\psi _i(𝐫_1)\psi _i^{}(𝐫_4)\psi _i(𝐫_3)\psi _i^{}(𝐫_2).`$
Alternatively, these statistics can be obtained directly by assuming that the distribution of wavefunctions is Gaussian , though in this case special care must be taken with regard to the normalization constraint.

There are alternate routes to our very natural starting point Eq. (14). One that has been developed recently is the “universal Hamiltonian” approach  which uses RPA to treat the interactions, RMT properties of the single-particle wavefunctions, and the small parameter $`1/g`$ to arrive at an effective Hamiltonian. Treating this effective Hamiltonian in the Hartree-Fock approximation leads to the same expression for the ground state energy that we give above.

We end this section by giving explicit expressions for the inverse compressibility, which is proportional to the peak spacing \[Eq. (4)\], in two simple cases in which the standard up-down filling of the states is assumed. First consider the case of $`N`$ even with $`n=N/2`$. For the standard filling, it is a singlet state and the two neighboring states for $`N1`$ and $`N+1`$ electrons are doublets. The peak spacing is then

$`\mathrm{\Delta }^2E_{N\mathrm{even}}=\mathrm{\Delta }^2E_{\mathrm{GTF}}+\delta ϵ_{n+1}\delta ϵ_n`$ (23)
$`+{\displaystyle \underset{i=1}{\overset{n}{}}}\left[2\left(M_{n+1,i}M_{n,i}\right)\left(N_{n+1,i}N_{n,i}\right)\right].`$
On the other hand if $`N`$ is odd, still assuming the simplest up-down filling, the peak spacing is

$$\mathrm{\Delta }^2E_{N+1\mathrm{odd}}=\mathrm{\Delta }^2E_{\mathrm{GTF}}+M_{n+1,n+1}.$$
(24)
The second difference of $`E_{\mathrm{GTF}}`$ is almost equal to the charging energy $`e^2/C`$ but there is a small correction from the residual interactions; this is discussed in detail below.

## III Large Dot Limit

The main simplifying feature in the large dot limit, $`k_FL1`$, is that the fluctuations of the interactions can be neglected, much as for diffusive nanoparticles  in the large $`g`$ limit. The distribution of the spacings is determined, then, by the fluctuations of the single-particle levels and their interplay with the mean residual interactions. Of course, as has been studied previously , the ground state of the nanoparticle does not necessarily follow the simple filling used in Eqs. (23)-(24), and this effect must be included in finding the distribution.

Let us consider the mean and variance of the $`M_{ij}`$ and $`N_{ij}`$ for levels near the Fermi level. We have, for instance, using Eq. (21)

$`M_{ij}={\displaystyle 𝑑𝐫𝑑𝐫^{}\left|\psi _i(𝐫)\right|^2V_{\mathrm{TF}}(𝐫𝐫^{})\left|\psi _j(𝐫^{})\right|^2}`$ (26)
$`+\delta _{ij}{\displaystyle 𝑑𝐫𝑑𝐫^{}\left|\psi _i(𝐫)\psi _i^{}(𝐫^{})\right|^2V_{\mathrm{TF}}(𝐫𝐫^{})}`$
We will make the simplifying assumption that the nanoparticle is a billiard so that the average density is constant: $`\left|\psi (𝐫)\right|^2=1/A`$ where $`A`$ is the area of the particle and the correlation function is given by Eq. (18). In this case the mean values are

$`M_{ij}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}+\delta _{ij}{\displaystyle \frac{\widehat{V}_{\mathrm{TF}}_{\mathrm{fs}}}{A}}`$ (27)
$`N_{ij}`$ $`=`$ $`{\displaystyle \frac{\widehat{V}_{\mathrm{TF}}_{\mathrm{fs}}}{A}}`$ (28)
where

$$\widehat{V}_{\mathrm{TF}}_{\mathrm{fs}}=_0^{2\pi }\frac{d\theta }{2\pi }\widehat{V}_{\mathrm{TF}}\left(k_F\sqrt{2(1+\mathrm{cos}\theta )}\right)$$
(29)
is the average on the Fermi surface of $`\widehat{V}_{\mathrm{TF}}(𝐤𝐤^{})`$. In the zero-range interaction limit, these expressions simplify to

$$M_{ij}=(1+\delta _{ij})\mathrm{\Delta }/2,N_{ij}=\mathrm{\Delta }/2.$$
(30)
The magnitude of the residual interaction effects is thus of order the mean level spacing $`\mathrm{\Delta }`$. (Note that the $`\mathrm{\Delta }`$ that we use here and throughout this paper is the spacing of the orbital levels alone and so does not take into account the spin degeneracy factor.) Since energies of this order are critical in determining the spacing of CB peaks—note, for instance, the case of standard up-down filling Eq. (23)—the residual interaction terms must be included in any theory. In particular, the reason for the failure of the simple constant interaction model, Eqs. (1) and (5), is now clear: it does not consistently keep all terms of order $`\mathrm{\Delta }`$.

The variance or covariance of any of the $`M_{ij}`$ and $`N_{ij}`$ will involve the correlation function of a wavefunction. For example, consider the variance of $`M_{ij}`$:

$`\mathrm{var}(M_{ij})`$ $`=`$ $`{\displaystyle }d𝐫_1d𝐫_2d𝐫_3d𝐫_4[|\psi _i(𝐫_1)|^2A^1]V_{\mathrm{TF}}(𝐫_1𝐫_2)[|\psi _j(𝐫_2)|^2A^1]`$ (32)
$`\times [|\psi _i(𝐫_3)|^2A^1]V_{\mathrm{TF}}(𝐫_3𝐫_4)[|\psi _j(𝐫_4)|^2A^1]`$
$`=`$ $`{\displaystyle 𝑑𝐫_1𝑑𝐫_2𝑑𝐫_3𝑑𝐫_4V_{\mathrm{TF}}(𝐫_1𝐫_2)V_{\mathrm{TF}}(𝐫_3𝐫_4)\left|\psi _i(𝐫_1)\psi _i^{}(𝐫_3)\right|^2\left|\psi _j(𝐫_2)\psi _j^{}(𝐫_4)\right|^2}.`$ (33)
The first equality is due to the normalization of the wavefunction, and the second derives from Eq. (21). Because of the decay of the autocorrelation of $`\psi `$, Eq. (18), the spatial integrals in (32) are restricted to a region where $`𝐫_1`$ is close to $`𝐫_3`$ and likewise $`𝐫_2`$ to $`𝐫_4`$. These restrictions produce factors of the small parameter $`1/k_FL`$ compared to the mean values. We will evaluate these factors below. However, for the large dot limit here, we see that to zeroth order in $`1/k_FL`$ all the fluctuations in the residual interactions may be neglected , as in the diffusive case . Qualitatively, this is natural because of the averaging implicit in the integral in $`M_{ij}`$ and $`N_{ij}`$: when the size of the system is much larger than the scale of oscillation of $`\psi (𝐫)`$, particular features of $`\psi _i`$ and $`\psi _j`$ become less important compared to the mean behavior.

For the CB peak spacing distribution, we thus consider the expression for the ground state energy (14) with constant $`M_{ij}`$ and $`N_{ij}`$. The problem is still complicated because of the interplay of the fluctuations in the single particle levels with these constants, which may, for instance, lead to $`S=1`$ or higher ground state spins . We thus evaluate the distribution of Coulomb blockade peak spacings numerically: Hamiltonians are taken at random in the GUE ensemble, and the lowest energy states for $`N1`$, $`N`$, and $`N+1`$ particles are obtained by determining in each case the occupation numbers $`\{n_{i\sigma }\}`$ which minimize the expression Eq. (14).

The resulting distribution of CB peak spacings is shown in Fig. 1. We show it for two different strengths of interactions, the case of $`r_s=1`$ and the zero-range interaction limit. The first thing to note is that for both of these values the distribution differs substantially from that in the CI model, Eq. (5), with much less odd/even alternation present. Still, there is a $`\delta `$-function in the distribution, as is clear in the cumulative distribution functions in panel (c), because the constant residual interaction term gives a rigid shift to the spacing for odd $`N`$ and standard filling, Eq. (24). The origin of spacing in these plots is $`\mathrm{\Delta }^2E_{\mathrm{GTF}}+M_{ij}`$. This is the mean spacing for adding a particle into the next available state neglecting all spin/filling effects; it thus corresponds to the classical charging energy $`e^2/C`$.

The fluctuations of the peak spacings are smaller in this constant exchange model than in the CI model. This is natural since the even and odd parts of the distribution are brought closer together. The values are rms$`(\mathrm{\Delta }^2E_N)/\mathrm{\Delta }=0.58`$ in the CI model and $`0.28`$ ($`0.24`$) in the present large dot case for $`r_s=1`$ (zero-range interaction).

## IV Fluctuations of Residual Interactions

In order to evaluate the statistical effect of fluctuating residual interactions, we need to find the distribution of the $`M_{ij}`$ and $`N_{ij}`$. We will show that in the limit of large but not infinite $`k_FL`$, the $`M`$’s and $`N`$’s are Gaussian distributed and uncorrelated.

The first step is to evaluate the variance of the $`M`$’s and $`N`$’s. Using the random matrix ensemble introduced in Section II, we find to a good approximation,

$`\mathrm{var}(M_{ij})`$ $``$ $`{\displaystyle \frac{8}{\pi A^2}}{\displaystyle _{\pi /L}^{2k\pi /L}}{\displaystyle \frac{dq}{q}}{\displaystyle \frac{V_{\mathrm{TF}}(q)^2}{(2kL)^2(qL)^2}}{\displaystyle \frac{1}{\pi A^2}}{\displaystyle \frac{\mathrm{ln}(kL)}{(kL)^2}}[2\widehat{V}_{\mathrm{TF}}^2(0)+\widehat{V}_{\mathrm{TF}}^2(2k)]`$ (34)
$`\mathrm{var}(N_{ij})`$ $``$ $`{\displaystyle \frac{2}{\pi A^2}}{\displaystyle _{\pi /L}^{2k\pi /L}}{\displaystyle \frac{dq}{q}}{\displaystyle \frac{[\widehat{V}_{\mathrm{TF}}(q)+\widehat{V}_{\mathrm{TF}}(2k)]^2}{(2kL)^2(qL)^2}}`$ (35)
$`M_{ij}N_{ij}M_{ij}N_{ij}`$ $``$ $`{\displaystyle \frac{4}{\pi A^2}}{\displaystyle _{\pi /L}^{2k\pi /L}}{\displaystyle \frac{dq}{q}}{\displaystyle \frac{\widehat{V}_{\mathrm{TF}}(q)[\widehat{V}_{\mathrm{TF}}(q)+\widehat{V}_{\mathrm{TF}}(2k)]}{(2kL)^2(qL)^2}};`$ (36)
details for the case of $`\mathrm{var}(M_{ij})`$ are given in Appendix A. The diagonal part of the direct residual interaction has an extra contribution because of the additional way in which the wavefunctions can be paired:

$`\mathrm{var}(M_{ii})2\mathrm{v}\mathrm{a}\mathrm{r}(M_{ij}^2)+{\displaystyle \frac{8}{\pi A^2}}{\displaystyle _\pi ^{2k\pi /L}}{\displaystyle \frac{dq}{q}}{\displaystyle \frac{\widehat{V}_{\mathrm{TF}}(q)[\widehat{V}_{\mathrm{TF}}(0)+\widehat{V}_{\mathrm{TF}}(\sqrt{(2k)^2q^2})]}{(2kL)^2(qL)^2}}.`$ (37)
In the zero-range interaction limit, the expressions for the variance of the $`M`$’s and $`N`$’s simplify considerably, as for the mean values given in Eq. (30). In this case, we find

$$\mathrm{var}(M_{ij})=\mathrm{var}(N_{ij})=\frac{3\mathrm{\Delta }^2}{4\pi }\frac{\mathrm{ln}(kL)}{(kL)^2}(1+3\delta _{ij}).$$
(38)
Note that the decay of the wavefunction correlations appearing in the variance produces a factor of $`1/k_FL`$ in the rms compared to the mean. The $`\mathrm{ln}(kL)`$ factor is special to two dimensions; it comes from the $`1/kL`$ decay of the wavefunction correlator in this case. Thus in the large dot limit, $`k_FL1`$, the distributions of the $`M`$’s and $`N`$’s are narrowly peaked about their mean values. This is the justification for ignoring these fluctuations to a first approximation in the last section.

Higher moments of the distribution involve coupling between more pairs of wavefunctions. The non-Gaussian part of the distribution is described by the cumulants, in which all of the wavefunctions are coupled in ways not present in the lower moments. In this case, the spatial integrals will have more restrictions coming from the wavefunction correlation (18); these restrictions will in turn produce additional factors of the small parameter $`1/k_FL`$. Thus in the large dot limit all higher moments can be neglected. In the case of third and fourth moments this can be easily verified by explicit calculation. Likewise it is straightforward to check that cross correlations among the $`M`$’s and $`N`$’s do not need to be considered, except for the obvious correlation between $`M_{ij}`$ and $`N_{ij}`$ which are equal in the zero-range interaction limit. We shall return to this point in Section VI. In this Section, however, we shall either take $`M_{ij}=N_{ij}`$ in the zero-range interaction limit, or neglect this correlation for $`r_s=1`$, and therefore take the $`M_{ij}`$ and $`N_{ij}`$ to be uncorrelated random Gaussian variables with mean and variance given above.

We now use the statistics of the $`M_{ij}`$, $`N_{ij}`$, and single particle levels (assumed to follow the random matrix GUE ensemble as in Section III) to deduce the statistics of the CB peak spacings. It is clear from the expressions for the spacing in the simple up-down filling case, Eq. (23), that sums such as $`_jM_{ij}`$ will be important. Because of the saturation of the fluctuations, item iii) in the random matrix model of Section II, the variance of such a sum involves only levels within $`E_{\mathrm{Th}}`$ of each other. In calculating $`\mathrm{var}\left(_jM_{ij}\right)`$, then, the sum should be taken over of order $`k_FL`$ independent $`M`$’s:

$$\mathrm{var}\left(\underset{j}{}M_{ij}\right)\frac{\mathrm{ln}(k_FL)}{k_FL}\mathrm{\Delta }^2.$$
(39)
The variance of other sums over the residual interaction terms will have the same dependence on $`k_FL`$. Up to a $`\sqrt{\mathrm{ln}}`$ factor, we thus expect the contribution of the fluctuations of the residual interactions to the CB peak spacings to be of order $`\mathrm{\Delta }/\sqrt{k_FL}`$ in rms.

In order to evaluate the magnitude more accurately, we turn again to numerical evaluation. We consider a collection of $`E_{\mathrm{Th}}/\mathrm{\Delta }`$ levels with GUE energies and Gaussian residual interaction terms. We fill the levels with $`N1`$, $`N`$, and $`N+1`$ particles successively, and find the occupation numbers $`\{n_{i\sigma }\}`$ minimizing the energy Eq. (14) in each case. Fig. 2 shows the resulting distribution of CB peak spacings. The fluctuations clearly act to smooth out all the sharp features of the $`k_FL\mathrm{}`$ case: the $`\delta `$-function becomes a finite peak, and the discontinuity at the origin is removed. The difference between the $`N`$ even and odd cases is quite clear—even for the relatively small $`N=80`$ (for which $`1/\sqrt{k_FL}0.21`$). The fluctuations of the residual interaction terms lead, of course, to an increase in the fluctuations of the spacing; now rms$`(\mathrm{\Delta }^2E_N)/\mathrm{\Delta }=0.31`$ for $`r_s=1`$ and $`N=80`$-$`200`$ compared with $`0.28`$ in the absence of such fluctuations. While the inclusion of the fluctuations has clearly made the distributions more like the experiments, substantial differences remain.

## V Scrambling: Added Electron Changes Confinement

As electrons are added to the nanoparticle, the mean field potential which confines the electrons changes. This causes a change in the single particle levels, whose energies appear in the expression for the ground state energy Eq. (14), which then in turn change the CB peak spacing. Because we consider a chaotic quantum dot, the change in the energy levels will be unpredictable; hence, we call this effect “scrambling”. Note that it is an intrinsic effect connected to the charge of the electron not the geometry or environment of the dot. This effect has been studied previously for diffusive quantum dots in Ref. ; here we give a derivation directly in the context of ballistic chaotic dots.

We know from the Strutinsky approach (Section II) that the change in the single particle levels should be found from the change in the smooth Thomas-Fermi potential, Eq. (11). So we first evaluate this change in potential, and then use perturbation theory to find the change in the energy levels. Perturbation theory is justified since the resulting shift is smaller than the mean level separation $`\mathrm{\Delta }`$ by a factor depending on $`k_FL`$. As a check, one can evaluate the magnitude of the fluctuation in the event of a complete scrambling of all the single particle levels. The properties of sums of random matrix eigenvalues in Ref.  imply a very large fluctuation of the peak spacings—of the order $`2`$-$`4\mathrm{\Delta }`$ which is larger than observed experimentally.

### A Evaluation of $`\delta V_{\mathrm{eff}}`$

The generalized Thomas-Fermi problem is specified in terms of the density functional

$$_{\mathrm{GTF}}[n]=𝒯_{\mathrm{TF}}[n]+_{\mathrm{ext}}[n]+_{\mathrm{coul}}[n]+_{\mathrm{xc}}[n]$$
(40)
comprising four contributions: the kinetic energy

$$𝒯_{\mathrm{TF}}[n]=𝑑𝐫\left[^{n(𝐫)}𝑑n^{}ϵ(n^{})\right]$$
(41)
where $`ϵ(n^{})`$ is the maximum kinetic energy of non-interacting particles with density $`n^{}`$, the confinement potential defining the geometry of the dot

$$_{\mathrm{ext}}[n]=𝑑𝐫V_{\mathrm{ext}}(𝐫)n(𝐫),$$
(42)
the Coulomb energy

$$_{\mathrm{coul}}[n]=\frac{e^2}{2}𝑑𝐫𝑑𝐫^{}\frac{n(𝐫)n(𝐫^{})}{|𝐫𝐫^{}|},$$
(43)
and the exchange-correlation potential $`_{\mathrm{xc}}`$. The Thomas-Fermi density is found by solving

$$\frac{\delta _{\mathrm{GTF}}}{\delta n}[n_{\mathrm{GTF}}]=\mu _{TF},$$
(44)
where $`n(𝐫)`$ is the electronic density and $`\mu _{TF}`$ is found such that $`n(𝐫)𝑑𝐫=N`$. The effective potential is, then, defined by

$$V_{\mathrm{eff}}(𝐫)\frac{\delta (_{\mathrm{ext}}+_{\mathrm{coul}}+_{\mathrm{xc}})}{\delta n}[n_{\mathrm{GTF}}].$$
(45)
Let $`\delta n_{\mathrm{GTF}}=n_{\mathrm{GTF}}(N+1)n_{\mathrm{GTF}}(N)`$ be the change in the Thomas-Fermi density when one electron is added to the system. Bearing in mind that Eq. (44) is a classical-like equation for which $`\delta n_{\mathrm{GTF}}`$ is small, one can write

$`{\displaystyle 𝑑𝐫^{}\frac{\delta ^2𝒯_{\mathrm{TF}}}{\delta n^2}[n_{\mathrm{GTF}}](𝐫,𝐫^{})\delta n_{\mathrm{GTF}}(𝐫^{})}`$ (47)
$`+{\displaystyle 𝑑𝐫^{}\frac{\delta V_{\mathrm{eff}}}{\delta n}(𝐫,𝐫^{})\delta n_{\mathrm{GTF}}(𝐫^{})}=\delta \mu _{TF}`$
with $`\delta \mu _{TF}`$ such that $`\delta n_{\mathrm{GTF}}(𝐫)=1`$. In two dimensions, the Coulomb energy of the added charge will dominate any variations in the kinetic energy. One can therefore write the variation of the density as $`\delta n_{\mathrm{GTF}}=\delta n_{\mathrm{GTF}}^0+\delta n_{\mathrm{GTF}}^1+\mathrm{}`$ with $`\delta n_{\mathrm{GTF}}^0`$ the solution of the electrostatic problem

$`\delta \mu ^0`$ $`=`$ $`{\displaystyle 𝑑𝐫^{}\frac{\delta ^2_{\mathrm{coul}}}{\delta n^2}(𝐫,𝐫^{})\delta n_{\mathrm{GTF}}^0(𝐫^{})}`$ (48)
$`=`$ $`e^2{\displaystyle 𝑑𝐫^{}\frac{\delta n_{\mathrm{GTF}}^0(𝐫^{})}{|𝐫𝐫^{}|}}.`$ (49)
At this level of approximation, the variation of the potential inside the dot is just the constant $`\delta \mu ^0`$, which will just shift the one particle energies $`ϵ_i`$ but not give rise to any fluctuation. One should keep in mind though that Eq. (44) applies only at places where $`n_{\mathrm{GTF}}(𝐫)0`$. Therefore the boundaries of the effective potential can be affected, in a way that depends, of course, largely on the external confining potential, and in particular on the second derivative of this latter at the the boundary. If the curvature is weak, fairly large displacements of the boundary can occur, which could play a significant role.

If, on the other hand, one assumes a billiard-like confining potential, the only source of modification of the effective potential comes from $`\delta n_{\mathrm{GTF}}^1`$, which is obtained through the equation

$$𝑑𝐫^{}\frac{\delta ^2𝒯_{\mathrm{TF}}}{\delta n^2}(𝐫,𝐫^{})\delta n_{\mathrm{GTF}}^0(𝐫^{})+𝑑𝐫^{}\frac{\delta V_{\mathrm{eff}}}{\delta n}(𝐫,𝐫^{})\delta n_{\mathrm{GTF}}^1=\delta \mu ^1.$$
(50)
Noting that $`(\delta V_{\mathrm{eff}}/\delta n)\delta n_{\mathrm{GTF}}^1`$ is, up to an inessential additive constant, the variation $`\delta V_{\mathrm{eff}}`$ of the effective potential we are interested in, and that

$$\frac{\delta ^2𝒯_{\mathrm{TF}}}{\delta n^2}(𝐫,𝐫^{})=\nu ^1\delta _d(𝐫𝐫^{}),$$
(51)
we get

$$\delta V_{\mathrm{eff}}(𝐫)=\mathrm{\Delta }\left(\delta \stackrel{~}{n}(𝐫)\overline{\delta \stackrel{~}{n}}\right)$$
(52)
where $`\delta \stackrel{~}{n}A\delta n_{\mathrm{GTF}}^0(𝐫)`$ is a smooth function of order $`1`$ ($`A`$ is the area of the dot). Note that we choose $`\delta V_{\mathrm{eff}}`$ with mean value zero. This ensures that there is no change in the mean energy levels; any such change in the mean should be incorporated in the charging energy.

### B Consequences for the peak spacing distribution

The variation of the effective potential induces a change of the one particle energies. We let $`\epsilon _N`$ denote the single particle energy of the $`N^{\text{th}}`$ electron, whatever orbital it is in (for $`N`$ even and simple up/down filling, for instance, $`\epsilon _N=ϵ_{N/2}`$). In perturbation theory, $`\delta \epsilon _j=\delta \epsilon _j^{(1)}+\delta \epsilon _j^{(2)}+\mathrm{}`$ with

$`\delta \epsilon _j^{(1)}`$ $`=`$ $`\psi _j|\delta V_{\mathrm{eff}}|\psi _j`$ (53)
$`\delta \epsilon _j^{(2)}`$ $`=`$ $`{\displaystyle \underset{kj}{}}{\displaystyle \frac{|\psi _j|\delta V_{\mathrm{eff}}|\psi _k||^2}{\epsilon _k\epsilon _j}}.`$ (54)
When going back to the inverse compressibility, one should bear in mind, however, that the variation of $`V_{\mathrm{eff}}`$ when going from $`N1`$ to $`N`$ electron is the same as when going from $`N`$ to $`N+1`$. As a consequence, it can be checked easily that the linear variation $`\delta \epsilon _j^{(1)}`$ cancels for all levels except the two top ones. At this level of approximation one therefore gets

$$\mathrm{\Delta }^2E_N=\mathrm{\Delta }^2E_N^{\mathrm{CI}}+\mathrm{\Delta }^2E_N^{\mathrm{scramb}}$$
(55)
where $`\mathrm{\Delta }^2E_N^{\mathrm{CI}}`$ is the constant interaction model result given in Eq. (5), and

$$\mathrm{\Delta }^2E_N^{\mathrm{scramb}}=\delta \epsilon _{N+1}+\delta \epsilon _N+\underset{j=1}{\overset{N1}{}}\delta ^2\epsilon _j$$
(56)
is the correction due to scrambling. The levels used for $`\mathrm{\Delta }^2E_N^{\mathrm{CI}}`$ are those for the $`N`$-electron $`V_{\mathrm{eff}}`$. Throughout this scrambling argument we assume that the filling of the levels does not change.

#### First order terms

The mean of the first order terms is zero by construction. The variance of $`\delta \epsilon _i^{(1)}`$ is

$$\mathrm{var}\left(\delta \epsilon _i^{(1)}\right)=\mathrm{\Delta }^2𝑑𝐫𝑑𝐫^{}|\psi _i(𝐫)\psi _i^{}(𝐫^{})|^2\delta \stackrel{~}{n}(𝐫)\delta \stackrel{~}{n}(𝐫^{}).$$
(57)
Using the wavefunction correlations (18) and approximating $`J_0^2(k_F|𝐫𝐫^{}|)1/(\pi k_F|𝐫𝐫^{}|)`$ plus oscillating terms ($`k_F|𝐫𝐫^{}|>1`$), we find

$$\mathrm{var}\left(\delta \epsilon _i\right)\alpha \frac{\mathrm{\Delta }^2}{\pi k_FL}.$$
(58)
The prefactor depends on the geometry of the dot under consideration:

$$\alpha =d^2𝐮d^2𝐮^{}\frac{1}{|𝐮𝐮^{}|}\delta \stackrel{~}{n}(𝐮)\delta \stackrel{~}{n}(𝐮^{})$$
(59)
is a dimensionless quantity of order 1 (note $`𝐮=𝐫/L`$). For a circular dot, for instance, with diffuse boundary scattering to make it chaotic, $`\delta \stackrel{~}{n}=\{[1(r/R)^2]^{1/2}2\}/2`$ where $`R`$ is the radius; in this case, $`\alpha 1.2`$.

Bearing in mind that $`k_FL\sqrt{2\pi N}`$ (taking into account the spin degeneracy), we see that the typical amplitude of the variation of the $`\delta \epsilon _i`$ scales as $`N^{1/4}`$ times the mean level spacing. This may, in fact, be numerically significant for experimental dots—for a circular dot with 100 electrons, for instance, the magnitude would be $`0.12\mathrm{\Delta }`$ .

#### Second order terms

The variance of the sum $`_1^{N1}\delta ^2\epsilon _i`$ can be evaluated along the same lines. Since it turns out to be parametrically smaller, we will be less careful with the constant.

First we note that the fluctuation of an individual matrix element is

$$\mathrm{var}\left(i|\delta V_{\mathrm{eff}}|j|\right)=\frac{\mathrm{\Delta }^2}{A^2}d𝐫d𝐫^{}J_0(k_i\delta 𝐫)J_0(k_j\delta 𝐫)\delta \stackrel{~}{n}(𝐫)\delta \stackrel{~}{n}(𝐫^{})$$
(60)
where $`\delta 𝐫=𝐫𝐫^{}`$. Following the same reasoning as above, this implies that

$$\mathrm{var}\left(i|\delta V_{\mathrm{eff}}|j|\right)\{\begin{array}{ccc}0& \mathrm{if}& \delta kL>1\\ \alpha \mathrm{\Delta }^2/\pi k_FL& \mathrm{if}& \delta kL<1\end{array}$$
(61)
where $`\delta k=k_ik_j`$.

This result is now inserted in the expression (54) for the energy shift $`\delta \epsilon _j^{(2)}`$ and then for the inverse compressibility Eq. (56). Assuming the matrix elements $`i|\delta V_{\mathrm{eff}}|j|`$ are Gaussian distributed independent variables, we find

$`\sigma _2^2`$ $``$ $`\mathrm{var}\left({\displaystyle \underset{i=1}{\overset{N1}{}}}\delta ^2\epsilon _i\right)`$ (62)
$``$ $`{\displaystyle \underset{i=NN_c}{\overset{N1}{}}}{\displaystyle \underset{j=N}{\overset{N+N_c}{}}}{\displaystyle \frac{1}{(\overline{\epsilon }_i\overline{\epsilon }_j)^2}}\mathrm{var}(|\psi _i|\delta V_{\mathrm{eff}}|\psi _j||^2).`$ (63)
Because of the saturation of the fluctuations of the levels, point (iii) in the random matrix criteria of Section II, the sum over levels is cut off at $`N_c=k_FL`$. The final result is

$$\sigma _2^2\mathrm{ln}k_FL\left(\frac{\mathrm{\Delta }}{k_FL}\right)^2$$
(64)
which is parametrically smaller than the first order result (58). We therefore see that the the effect of scrambling is dominated by the variation of the top levels, and that this adds contributions to the peak spacing fluctuations which while parametrically small on the scale of $`\mathrm{\Delta }`$ could be numerically important for dots which are not too large.

### C Numerics

It is a simple matter to incorporate the scrambling of the top two levels into our numerical investigation of the ground state energy specified by Eq. (14). In random matrix theory, the change in an energy level upon varying a parameter is uncorrelated with its value , so we simply choose random $`\delta \epsilon _N`$ and $`\delta \epsilon _{N+1}`$ from a Gaussian distribution with variance Eq. (58). The ground state configuration is then found as before.

Fig. 3 shows the resulting probability density of CB peak spacings. In the $`r_s=1`$ case there is still a substantial difference between the distributions of even and odd spacings: for $`N`$ even, the large spacing tail coming from the level-spacing distribution (Wigner surmise) survives even for less than 100 electrons. When the interactions are stronger, as in the zero-range interaction case shown, the even and odd distributions are more similar, but there is still a discernible difference. The magnitude of the fluctuations increases further: we now find rms$`(\mathrm{\Delta }^2E_N)/\mathrm{\Delta }=0.33`$-$`0.32`$ for $`r_s=1`$ and $`N=80`$-$`200`$.

## VI Beyond the Gaussian model

Our treatment of the fluctuations of the residual interactions in the previous sections relies on a certain number of assumptions. We have for instance neglected correlations among the $`M_{ij}`$ and $`N_{i^{}j^{}}`$ (except for the obvious $`M_{ij}=N_{ij}`$ for a zero-range interaction) when we know, for instance, that the correlation coefficient of $`M_{ij}`$ with $`N_{ij}`$ is of order 1 even in the infinite size limit. Similarly we have assumed the absence of correlation between these residual interaction terms and the one particle energies $`ϵ_i`$. Furthermore, we have used expressions such as Eqs. (34)-(37) which are derived in the large $`k_FL`$ limit. To check that these assumptions do not drastically modify the resulting distribution of CB peak spacings, we shall in this section implement the random matrix ensemble described in Section II \[paragraph containing Eq. (17)\] by directly performing numerical Monte Carlo calculations.

To implement the constraints of our random matrix model, we start by considering the system to be a square of size $`L`$, for which we can use a plane wave basis

$$|n_x,n_y=L^2\mathrm{exp}[i2\pi n_xx/L]\mathrm{exp}[i2\pi n_yy/L]$$
(65)
ordered by increasing momentum. In this basis, we write the one particle Hamiltonian as $`H=D+V`$ where $`D`$ is diagonal with $`D_{ii}=i`$. $`V`$ is taken to be a banded random Hamiltonian with independent complex matrix elements which are Gaussian distributed with variance $`v_{ij}^2=a(i)exp(b(i)|ij|)`$. For $`a(i)`$ in the correct range  such an ensemble yields (1) eigenvalues with a mean spacing close to one, (2) local fluctuations which are typical of the GUE, and (3) eigenfunctions which are localized in energy. The localization in energy can be characterized by the inverse participation ratio in the plane wave basis: $`I_j=_{n_x,n_y}|n_x,n_y|ϵ_j|^4`$. Our requirement Eq. (17) implies the constraint

$$I_j^1=(3/4)(k_jL).$$
(66)
In practice, we use a simple algorithm to tune the value of $`a(i)`$ and $`b(i)`$ in such a way that this constraint holds. Diagonalizing the random matrix specified leads to eigenvalues and eigenvectors which are then used to evaluate ground state energies for different number of particles via Eq. (14), and hence arrive at the CB peak spacings in this model.

In this approach, deviations from the wavefunction correlations (18) required by normalization are naturally present. Furthermore, correlation between eigenvalues and eigenfunctions is present as well. Of course, in a real physical system, these effects would be automatically included because of the scattering at the boundary of the system—the transit time $`\tau `$, for instance, gives both the scale over which fluctuations saturate ($`\mathrm{}/\tau `$) and the length scale at which the Bessel function correlations become suspect ($`v_F\tau `$). We believe that our random matrix model includes the main features of these effects in a generic way.

Fig. 4 shows numerical results for this model. We consider only the smaller value of $`N`$, $`N=80`$, because the effects are expected to be only significant for smaller quantum dots. On the other hand, we present results both with and without the scrambling effect of Section V, and for different interaction strengths—$`r_s=1`$ corresponds to a typical value for a GaAs quantum dot while the zero-range interaction case corresponds to the maximum effect possible within a perturbative approach to the interactions. When these results are compared to those of Figs. 2 and 3, we see that there is a small but noticeable effect. This is highlighted in the quantile-quantile plots in the insets to Fig. 4 where an ordered set of RMT spacings is plotted against a similar set for the Gaussian model. The present RMT yields a more pronounced large spacing tail than the Gaussian model.

## VII Conclusions

In this paper we have studied the Coulomb blockade peak spacings for large but not infinite quantum dots in the case that the interactions can be treated perturbatively and the single-particle properties described by random matrix theory. Motivated by the Strutinsky approach to finite Fermi systems , we use an expression for the ground state energy \[Eq. (14)\] which includes the charging energy, the single particle quantization effects, and the residual (screened) interactions treated in first order (Hartree-Fock).

We presented results in four cases. First, in the limit $`k_FL\mathrm{}`$, fluctuations of the residual interactions are negligible. The CB peak spacing distribution has a discontinuity at $`e^2/C`$ and a $`\delta `$-function because of the odd $`N`$ case (Fig. 1). The distribution is very different from that predicted by the simple CI model which neglects some terms of order $`\mathrm{\Delta }`$ while keeping others.

Second, fluctuations in the residual interactions contribute a term typically of order $`\mathrm{\Delta }/\sqrt{k_FL}`$ to the CB peak spacings. The sharp features of the infinite dot limit are rounded out (Fig. 2).

Third, the addition of an electron in the CB process changes the confinement potential in the quantum dot and so changes the single particle properties. This scrambling effect also contributes a term typically of order $`\mathrm{\Delta }/\sqrt{k_FL}`$ to the CB peak spacings, and the distribution is rounded further. However, the difference between $`N`$ even and $`N`$ odd remains substantial at large spacings where a tail is present for $`N`$ even (Fig. 3) which is a remnant of the Wigner surmise distribution for level spacings.

Finally, we took a first step towards investigating corrections to the Gaussian uncorrelated model for the fluctuations of the residual interactions. We found that the large spacing part of the distribution is somewhat underestimated in the Gaussian model.

Before ending this paper, we compare our theoretical distributions to those obtained experimentally . The observed distribution has little structure, being Gaussian near its peak but with longer tails for both larger and smaller spacings. The width of the distribution is contentious: early work claimed that the width scaled with the charging energy  but this has not been verified by more recent experiments . In the more recent experiments, the width is clearly related to the mean level separation $`\mathrm{\Delta }`$ and in the range $`0.5`$-$`1.5\mathrm{\Delta }`$, somewhat larger than given by, e.g., the CI model.

The theoretical results presented here certainly look more like the experiments than the benchmark CI model: the distribution is more rounded and so less structured, and the odd/even alternation is much less pronounced. We suspect that further development of this approach will produce qualitative agreement with the experiments.

In the end, however, the theory presented here does not provide an adequate quantitative description of the experiments. To start with, the magnitude of the fluctuations is actually less than in the CI model (0.58 versus 0.32 for $`N=200`$) and so differs from experiment. With regard to the odd/even effect, the situation is somewhat less clear, and so more interesting. Looking at Fig. 3, for instance, one observes that a substantial odd/even effect is still present for $`r_s=1`$, while it has been largely washed out in the zero-range limit of the screened interaction. It should be stressed that this drastic modification of the distributions can be traced back mainly to the relatively modest change in the mean values $`M_{ii}\mathrm{\Delta }/2`$ and $`N_{ij}`$, Eqs. (27)-(30), which change from about $`0.3\mathrm{\Delta }`$ for $`r_s=1`$ to $`0.5\mathrm{\Delta }`$ in the zero-range limit. The sensitivity to the numerical value of these quantities should not come as a surprise since we know that the Stoner instability, for which a spontaneous polarization of the spin of many electrons takes place, occurs for $`N_{ij}=\mathrm{\Delta }`$. These observations suggest an experiment: placing a good conductor close to the 2DEG in such a way that it can screen the interaction between electrons within the quantum dot should significantly affect the peak spacing distribution. If the auxiliary screening is strong enough, a strong odd/even effect should be visible. A double quantum well structure seems like a possible setting for such an experiment.

We can divide the potential explanations for the quantitative discrepancy between theory and experiment into basically three categories. The first one, which we cannot absolutely exclude, is that noise in the experiments is corrupting the results . The second is that effects not included here, but interpretable within the kind of Fermi liquid theory on which we have relied throughout this paper, dominate. Some candidates are (a) the effect of the changing gate potential on the single particle levels —this would, however, somewhat contradict the remarkable stability observed in the wavefunction —(b) boundary effects in the scrambling which need to be considered using a realistic description of the dot, and (c) consequences of non-chaotic dynamics within the dots . Finally, a third completely different kind of interpretation is that the perturbative approach to interactions taken here is not accurate even though $`r_s1`$-$`2`$ . Further experimental and theoretical work is needed to discern between these options.

## Acknowledgements

We thank G. Usaj for many valuable conversations, and appreciate helpful discussions with I. L. Aleiner, L. I. Glazman, and C. M. Marcus. The LPTMS is a “Unité de recherche de l’Université Paris 11 associée au C.N.R.S.”

## A Variance of the $`M_{ij}`$

In this appendix we sketch the argument which leads to Eqs. (34)-(37). As in Section VI, we assume that the configuration space is a square of size $`L`$, for which the plane wave basis $`|𝐤=|n_x,n_y`$ introduced in Eq. (65) can be used.

Expanding the eigenfunctions in this basis,

$$|\phi _i=\underset{𝐤}{}v_{i𝐤}|𝐤,$$
(A1)
we can write the $`M_{ij}`$’s as

$`M_{ij}`$ $`=`$ $`{\displaystyle \underset{𝐤_1,𝐤_2,𝐤_3,𝐤_4}{}}v_{i𝐤_1}v_{i𝐤_2}^{}v_{j𝐤_3}v_{j𝐤_4}^{}\delta _{𝐤_1𝐤_2,𝐤_3+𝐤_4}\widehat{V}_{\mathrm{TF}}(𝐤_1𝐤_2)`$ (A2)
$`=`$ $`{\displaystyle \underset{𝐪}{}}\widehat{V}_{\mathrm{TF}}(𝐪)W_{i𝐪}W_{j𝐪}^{}`$ (A3)
where we define

$$W_{i𝐪}\underset{𝐤_1𝐤_2=𝐪}{}v_{i𝐤_1}v_{i𝐤_2}^{}.$$
(A4)
Our random matrix model (Section II) implies

$`v_{i𝐤_1}v_{j𝐤_2}^{}`$ $`=`$ $`(k_FL)^1\delta _{ij}\delta _{𝐤_1𝐤_2}\text{ if }\delta k<\pi /L`$ (A5)
$`=`$ $`0\text{ if }\delta k>\pi /L`$ (A6)
with $`\delta k=\left||𝐤_1|k_i\right|`$ and $`k_i=\sqrt{2mE_i}/\mathrm{}`$. From this, we deduce

$$\mathrm{var}(M_{ij})\frac{1}{A^2}\underset{𝐪0}{}\widehat{V}_{\mathrm{TF}}^2(𝐪)|W_{i𝐪}|^2|W_{j𝐪}|^2$$
(A7)
$`|W_{j𝐪}|^2`$ can be interpreted as $`(2\pi k_i)^2`$ times the area of the intersection of two rings of diameter $`k_i`$ and width $`2\pi /L`$, centered at a distance $`|𝐪|=q`$. Simple geometry therefore gives, for $`2\pi /L|𝐪|2k_i2\pi /L`$

$$|W_{i𝐪}|^2\frac{4}{(qL)\sqrt{(2k_iL)^2(qL)^2}}.$$
(A8)
We obtain for $`ij`$

$$\mathrm{var}(M_{ij})\frac{8}{\pi A^2}_{\pi /L}^{2k\pi /L}\frac{dq}{q}\frac{V_{\mathrm{TF}}(q)^2}{(2kL)^2(qL)^2}.$$
(A9)
The variance of $`M_{ii}`$ and $`N_{ij}`$ and the covariance between $`M_{ij}`$ and $`N_{ij}`$ can be computed along the same lines.


format: text

text:
# Phase Transitions in Load Transfer Models of Fracture.

## 1 Introduction

Fracture phenomena have attracted a lot of interest in the last several years, both experimentally and theoretically. In the lab, a disordered material subjected to an increasing external load can be experimentally studied by measuring the acoustic emissions before the global rupture. It has been shown  that this intense precursory activity in the form of bursts of different microscopic sizes follows a well-defined power law. From the theory side, the understanding of fracture in heterogeneous materials has progressed due to the use of lattice models and large scale numerical simulations . Recently, the introduction of models of material failure has lead to the evidence that rupture can be viewed as a kind of critical phenomenon . Nevertheless, the question of whether rupture exhibits the properties of a first-order or a second-order phase transition remains under discussion as well as what is the order parameter that determines the type of transition.

In this field, it is important to use models able to describe the complexity of the rupture process, although they should be simple enough to permit analytical insights. To this class of models belong the well-known Fiber Bundle Models (FBM) widely used since their introduction more than forty years ago . In static FBM, a set of fibers (elements) is located on a supporting lattice and one assigns to its elements a random strength threshold sampled from a probability distribution. The lattice is loaded and fibers break if their loads exceed their threshold values. Now, one can assume different load transfer rules to mimic the range of interactions among the fibers in the set. The Global Load Sharing rule is the simplest theoretical approach one can adopt to make the problem analytically tractable, which implies that the load carried by failed elements is equally distributed among the surviving elements of the system, representing in this way a long-range interaction among the constituent parts of the system. This is a kind of mean field approach to the more extreme Local Load Sharing scheme, where the interaction among the elements is short-range and the load supported by the failing elements is transferred only to a neighborhood.

Very recently , we have developed a probabilistic approach suited to smooth the fluctuations around the point of final collapse. In this paper, we extend the method to analyze the local load sharing scheme and contrast the two limiting cases of range interactions: the long range and the short range. While for the case of global load sharing, we obtain different scaling relations that point out that the fracture of a fiber set with long range interactions may be a continuous transition, for the case of short range interactions, the failure exhibits no significant precursory activity and the system breaks abruptly showing, in this case, the properties of a first order transition where no scaling laws can be identified. The rest of the paper is organized as follow. In Section 2, we present the standard approach to the FBM. In Section 3, we introduce the new probabilistic method and the results derived from it and from Monte Carlo simulations. Finally, Section 4 is devoted to discussion and conclusions.

## 2 Static Fiber Bundle Models

Let us first recall the basic ingredients of the static FBM and how one proceeds in numerical simulations. The term static means that time plays no role in the model . The system under consideration is a set of $`N_0`$ elements (fibers) located in a supporting lattice each one having at the initial state a zero load and a fixed strength threshold value (quenched disorder) sampled randomly from a probability distribution $`P(\sigma )`$. The system is then subjected to an external force $`F`$ which is shared democratically among the fibers, that is, the distribution of the total force produces that each element increases its load, $`\sigma `$, in the same amount. This individual stress acts as the control parameter. The loading process is done quasistatically; i.e., the external force is increased at a sufficiently slow rate as to produce a single breaking event when the stress in the weakest element equals its threshold value. Then, the increase of external force stops and the load of the broken fiber is transferred according to the transfer scheme assumed. In the global load sharing FBM, this implies that the load on fiber $`i`$ is given by $`\sigma _i=F/n_s(F)`$, where $`n_s`$ is the number of surviving elements for a given external load $`F`$. In the local load sharing FBM, the above relation does not hold any more because there the redistribution of load is performed among the nearest neighbors of the failed fiber giving rise to the appearance of regions with different stress concentrations.

In both schemes, the rupture of a fiber may induce secondary failures which in turn may trigger more failures. This process of induced failure at constant external load, termed avalanche, stops when all surviving elements carry a load lower than their thresholds. The system is then loaded again and the process is repeated until the final catastrophic avalanche provokes the total rupture of the material, which occurs at a critical load $`\sigma _c`$ that depends on the probability distribution from where the individual strengths were drawn as well as on the system size. At this point, it is worth recalling that the exact value of $`\sigma _c`$ can be analytically obtained in the thermodynamic limit for the global load sharing case while for local load sharing schemes there is no theoretical approach leading to $`\sigma _c`$. It is also known that the critical load $`\sigma _c`$ for the global load sharing FMB is non-zero and independent of $`N_0`$ for $`N_0\mathrm{}`$, whereas $`\sigma _c`$ tends to zero as $`N_0`$ goes to infinity in the case of local load sharing models . Fiber bundle models have also been recently used in self-organized criticality (SOC), a theoretical framework widely used for the study of avalanche phenomena in disordered systems. It has been shown using these models that systems with plastic behavior can reach a SOC state just before the global rupture . A second case of self-organization with power law distributions in several quantities corresponds to the situation in which the fracture process coexists with a healing process .

In numerical simulations, the cycle of complete breakdown of the material is performed many times in order to average out the effect of fluctuations. Nevertheless, the stress history of a particular element is made of steps, so that the fluctuations around the critical load can never be completely avoided. Let fiber $`k`$ be supporting a stress $`\sigma _k^0`$ at a given step of the failure process. It will continue to support $`\sigma _k^0`$ until a global driving or a stress transfer from one (or more) of the failed elements occurs. At this moment, element $`k`$ instantaneously changes its load to $`\sigma _k^1=\sigma _k^0+\sigma _{dist}`$, where $`\sigma _{dist}`$ is the load coming from the failed fibers or from the global driving. So, in a later step, element $`k`$ could receive again load, suffering a second step-like increase in stress. This step-like stress history continues until fiber $`k`$ fails. As we are interested in studying the behavior of the system as the critical point is approached, it is of utmost importance to find a simple method able to capture the evolution of the system near the critical point avoiding as much as possible the fluctuations.

## 3 Probabilistic Approach to Static Fiber Bundle Models

Let there be a set of $`N_0`$ elements located on a supporting lattice. Suppose that each element carries a given load $`\sigma `$, which is set to zero at the initial state. Fibers break depending on their strengths which are distributed according to a probability distribution $`P(\sigma )`$. Different probability distributions can be considered. In materials science the Weibull distribution is widely used,

$$P(\sigma )=1e^{(\frac{\sigma }{\sigma _0})^\rho },$$
(1)
where $`\rho `$ is the so-called Weibull index, which controls the degree of disorder in the system (the bigger the Weibull index, the smaller the disorder), and $`\sigma _0`$ is a load of reference. In the following, we will assume $`\sigma _0=1`$, and therefore the loads are dimensionless. Although the results we show hereafter are for the Weibull distribution to gain in definiteness, they have been also obtained for a wide class of distributions. Besides, for continuous distributions decaying fast enough, the approach to the critical state does not depend on the detailed form of the disorder .

Equation (1) represents the probability that an element fails under the individual load $`\sigma `$. Now, consider the case in which an element drawn from Eq. (1) supports a load $`\sigma _1`$, but breaks under a new load $`\sigma _2`$. The probability that this happens is given by

$$p(\sigma _1,\sigma _2)=\frac{P(\sigma _2)P(\sigma _1)}{1P(\sigma _1)},$$
(2)
which is equal to

$$p(\sigma _1,\sigma _2)=1e^{(\sigma _2^\rho \sigma _1^\rho )}.$$
(3)
So, the probability $`q(\sigma _1,\sigma _2)`$ that an element that has survived to the load $`\sigma _1`$ also survives to the load $`\sigma _2`$ will be given by

$$q(\sigma _1,\sigma _2)=1p(\sigma _1,\sigma _2)=e^{(\sigma _2^\rho \sigma _1^\rho )}.$$
(4)
All these probabilities depend on the state of stress of each element, which is a complex function of both the control parameter and the stress redistributions due to fiber failures. To mimic the quasistatic increase in load on the system as is applied in MC simulations, we impose the condition that under an external force $`F`$, the next breaking event consists of one single failure. Now, we will explore how the driving process controls the approach to global failure in the global load sharing case and in a particular version of the local load sharing schemes.

### 3.1 Global Load Sharing case

Let suppose that after the latest avalanche, there are $`N_k`$ surviving elements each one bearing a load $`\sigma _k`$. The assumption of quasistatic increase in load on the system implies that the load $`\sigma _l`$ needed to provoke the failure of only one element is ,

$$\sigma _l=\left[\sigma _k^\rho \mathrm{ln}\left(1\frac{1}{N_k}\right)\right]^{\frac{1}{\rho }},$$
(5)
where in Eq. (5) $`N_k=N_0`$ and $`\sigma _k=0`$ at the initial state. Elevating the external force up to the $`N_k\sigma _l`$ level a first element breaks. As we are dealing with a global load sharing set, the choice of the broken element is irrelevant because all of them are equivalent. Once the first element fails, the redistribution of its stress takes place which may induce secondary failures and so on, until the end of the avalanche.

Now, the number of surviving elements in the set after an avalanche has come to an end can be calculated by the recursive relations ,

$$N_{j+1}=N_jq(\sigma _{j1},\sigma _j),$$
(6)
$$N_j\sigma _j=N_{j1}\sigma _{j1}$$
(7)
$$N_j=N_{j+1}$$
(8)
Equation (6) relates the number of surviving elements between two successive avalanche steps, whereas Eq. (7) allows the computation of the new load acting on the $`N_j`$ surviving elements. This is possible because we are under a global load sharing scheme and thus all surviving elements support the same load. Finally, Eq. (8) determines the end of the ongoing avalanche.

The dynamics of the system is completely determined by Eq. (5), (6), (7). In this way, the size of an avalanche is given by the number of elements that break between two successive steps of external loading. The total stress accumulated in the system can be calculated multiplying the number of intact fibers before an avalanche starts by the load given by Eq. (5). The critical load, defined as the load needed to provoke the total collapse of the system, is equal to the load on the intact fibers just before the final catastrophic avalanche. Note that in this probabilistic approach, in contrast to Monte Carlo simulations, we need to store only the information concerning the loads of the intact elements, that is, the threshold dynamics is omitted with the subsequent advantages of saving computer resources and the possibility of exploring systems of larger size.

For Eq. (8) to be satisfied, we can proceed in two different ways in order to determine when an avalanche ends, to which we will refer as the continuous and the discrete cases. For the continuous case, the number $`N_{j+1}`$ of surviving elements is considered a real number. This means that condition (8) is never fulfilled before the final avalanche. So, condition (8) is replaced in numerical simulations by a factor $`\nu 1`$ that determines the end of an avalanche, i.e., if $`N_jN_{j+1}\nu `$ the avalanche stops; otherwise it continues. In the discrete case, $`N_{j+1}`$ is considered to be an integer number, so that after each iteration of Eq. (6), $`N_{j+1}`$ has to be rounded up. This is done comparing the remainder of $`N_{j+1}`$, $`\lambda `$, with a random number $`\alpha `$ uniformly distributed in the interval $`[0,1[`$. Thus, if $`\alpha \lambda `$, $`N_{j+1}`$ is equal to its whole part, and if not, $`N_{j+1}`$ is equal to its whole part plus one . Next, we check whether condition (8) is satisfied for the rounded value of $`N_{j+1}`$ or if a new iteration of Eq. (6) has to be performed. The continuous approach has the great advantage that the fluctuations are ruled out, whereas for the discrete case the results are similar to those obtained by Monte Carlo simulations where it is necessary to average over many realizations in order to get accurate mean values. This is because in the global fiber bundle model the central limit theorem applies .

Now, we proceed to explore the behavior of the system near the critical point. In particular, we are interested in inspecting the evolution of some quantities as the critical point is reached. In order to avoid unnecessary fluctuations around the critical point, the results shown below have been obtained for the continuous case of the probabilistic approach ($`\rho =2`$). It turns out  that the avalanche size near to the critical point diverges with an exponent $`\gamma =\frac{1}{2}`$ as $`s(\sigma _c\sigma )^\gamma `$. A similar behavior, through a mapping of a fuse network model to the global fiber bundle model used here, has been recently reported . This mapping between fiber bundle models and fuse networks with strong disorder was first noted  a few years ago; nevertheless, it has also been pointed out that 3D fuse networks apparently do not follow the FBM picture . Note, additionally, that for the probabilistic continuous version, this relation is obtained for a single realization avoiding, in this way, the large number of iterations performed in MC simulations. In Fig. 1 we have depicted the derivative of the number of broken fibers, ($`dN/d\sigma `$), as a function of the distance to the critical point $`\sigma _c\sigma `$ in a log-log plot. This rate $`dN/d\sigma `$ diverges as $`(\sigma _c\sigma )^\gamma `$ also with $`\gamma =\frac{1}{2}`$, thus qualifying a critical mean field behavior as was already shown  by means of analytical analysis of fiber bundle models. In Ref. , a similar scaling behavior is addressed for the derivative of the strain carried by the fibers with respect to the driving force.

Another way to shed light on the critical behavior of this type of system is to define a branching ratio $`\zeta `$ for each avalanche. This magnitude represents the probability to trigger future breaking events given an initial individual failure , and is related to the number of broken fibers by

$$\zeta =\frac{<z>1}{<z>}.$$
(9)
The above relation can be obtained by thinking of the evolution of fracture as a kind of branching process . In this process, each node gives rise to a number $`n`$ of new branches in the next time step. The average number $`<n>`$ of new branches is called the branching ratio. Let us denote by $`n_t`$ the number of branches at a given step $`t`$ of the branching process, and by $`t_{max}`$ the total number of time steps before it dies. Then,

$$\zeta =<n>=\frac{\underset{t=0}{\overset{t_{max}1}{}}n_{t+1}}{\underset{t=0}{\overset{t_{max}}{}}n_t},$$

and

$$\zeta =1\frac{n_0}{\underset{t=0}{\overset{t_{max}}{}}n_t}.$$

As $`n_0=1`$, $`\zeta =1\frac{1}{n_{tot}}`$ where $`n_{tot}`$ is the total number of nodes developed in the branching process. For a fracture process, $`n_{tot}`$ is equal to the average number of failure events. So, Eq. (9) defines the branching ratio. We represent by $`<z>`$ the average number of elements that fail in one avalanche, which is a function of the control parameter $`\sigma `$ and coincides with $`s`$. This analogy between fracture and branching processes has been previously used to study the criticality in the process of fragmentation of Hg drops . The branching ratio will then act as the order parameter. It takes the value 1 when the system is critical thereby representing a measure of the distance of the system from the critical state . We would like to remark here two characteristics of the branching ratio defined as above. First, it coincides with its general definition for a branching process. Here, the values of the branching ratio are always less than or equal to one because we are analyzing an irreversible breaking process that cannot continue forever, so that a value of the branching ratio greater than one would imply a physically unreachable situation. Secondly, we have introduced a definition in terms of the avalanche size and not in terms of the cluster size. The avalanche size in fiber bundle models is a measure of causally connected broken sites while the cluster size is a measure of spatially connected broken sites. It is clear that in the global fiber bundle model, the spatial correlations are ruled out and then the distribution of avalanche sizes does not coincide with the distribution of cluster sizes.

We have numerically computed $`\zeta `$ for the continuous version of the probabilistic model. The results obtained for a system of $`N_0=50000`$ elements and $`\rho =2`$ have been plotted in Fig. 2. It can be seen in this figure that the branching ratio approaches the unity as the critical load is reached. In the figure, the values of $`\zeta `$ are collected for all the avalanches except for that which provokes the catastrophic event leading to the collapse of the system. So, this dependency of $`\zeta `$ with $`\sigma `$ means that very close to the critical point, the initial failure provokes the breaking of another fiber which in turn induces tertiary ruptures and so on. This is the result that should be expected from the divergence of the avalanche size at the critical point. Figure 3 shows the behavior of the order parameter near the critical point. As can be seen, near the critical point, the relation $`1\zeta (\sigma _c\sigma )^\beta `$, where $`\beta =\frac{1}{2}`$ applies. Note the similarity of the behavior with those obtained for the magnetization in known magnetic systems with second-order phase transitions. It is of additional interest to note that the branching ratio also captures the feature that the precursory activity is significant only for strong disorder (for the Weibull distribution, small values of $`\rho `$); whereas for small disorder the system behaves more similar to the breakdown of homogeneous materials. On the other hand, the branching ratio does not depend on the size of the system for large systems, in contrast to previous results in other fracturing systems . In Fig. 4 we illustrate this behavior by plotting the value of the branching ratio for the last non-catastrophic avalanche, $`\zeta ^{}`$, as a function of the system size. In all cases, the numerical simulations were performed for $`\rho =2`$. It is clear from the figure that $`\zeta ^{}`$ reaches the unity as $`N_0`$ goes to infinity.

### 3.2 Local Load Sharing Case

It is well-known that the local load sharing FBM is much more complicated than the global load sharing case. The complexity of the fracture problem increases because the load borne by failed elements is transferred to nearest neighbors and then there appear regions of stress concentration throughout the system. The distribution of load is now not homogeneous, and we have to carry on several lists to record the individual load of the fibers in the system. However, with the probabilistic approach, we are able to study a few things about the behavior of the system, which are sufficient to remark the great differences between the long range and the short range interactions schemes. One important difference between the two schemes is that contrary to the global load sharing case, in the local scheme there is no significant precursory activity signaling the approach to the final collapse and the system undergoes an abrupt catastrophic avalanche. For maximum simplicity, we consider next a one-dimensional periodic local load system where the stress transfer is done by adding the load of the failed element always to the element on its right. Despite of its simplicity, it has been shown  that the general properties of this particular model are identical to those of more complex local schemes.

The probability that the system fails in only one avalanche is given by:

$$p(\sigma _1,\sigma _2)p(\sigma _1,\sigma _3)p(\sigma _1,\sigma _4)\mathrm{}p(\sigma _1,\sigma _{N_0}),$$
(10)
where $`\sigma _i=i\sigma _1`$ and $`\sigma _1`$ is

$$\sigma _1=\left(0\mathrm{ln}\left(1\frac{1}{N_0}\right)\right)^{\frac{1}{\rho }}.$$

That is, the probability of having a one-step failure is given by the probability that an element fails under a load $`\sigma _2`$ having survived to the load $`\sigma _1`$, multiplied by the probability that the next element in the lattice also fails under the new load $`\sigma _3`$ having survived to the load $`\sigma _1`$ and so on.

On the other hand, Eq. (4) takes the form

$$q(i,j)=e^{\sigma _1^\rho (j^\rho i^\rho )},$$
(11)
which gives the probability that an element that has survived to the load $`i\sigma _1`$ also survives to a new load $`j\sigma _1`$, being $`\sigma _1`$ the load that produces the first breaking event in the system at the initial state. The probability that the first avalanche consists of, for example in one failure, is given by $`q(1,2)`$. The first avalanche will have size two with probability $`p(1,2)q(1,3)`$. In general, the probability that the first avalanche does not provoke the rupture of the whole system, i.e., be finite, is given by

$$P_1=\underset{l=1}{\overset{\mathrm{}}{}}\mathrm{\Pi }_lq(1,l+1),$$
(12)
where

$$\mathrm{\Pi }_l=p(1,2)p(1,3)\mathrm{}p(1,l)=\underset{i=2}{\overset{l}{}}p(1,i),l<N_0,$$

with $`\mathrm{\Pi }_1=1`$. Thus, the average size of the first avalanche, derived from Eq. (12) is

$$S_1=\underset{l=1}{\overset{\mathrm{}}{}}l\mathrm{\Pi }_lq(1,l+1).$$
(13)
Numerical simulations of the continuous case, the discrete approximation (both defined as above for the global load sharing scheme) and Monte Carlo method confirm the validity of Eq. (13). For example, for $`\rho =2`$ and $`N_0=1000`$ we obtain $`S_1=1.0030`$, $`S_1=1.0022`$ and $`S_1=1.0039`$ respectively. This result is valid only for the first avalanche in the system.

We have also verified by means of Monte Carlo simulations that for this local transfer rule, there is no power law in the distribution of avalanche sizes. Besides, no scaling relations appear and the system is more sensitive to parameters such as the degree of disorder ($`\rho `$) and the size of the system $`N_0`$. Catastrophic avalanches arise very often in the first stages of the rupture process.

Finally, we illustrate how one can describe the precursory activity when it is limited to a few steps before the global rupture. Let us solve the particular situation in which the first avalanche is finite and the second provokes the final breakdown of the system. The load needed to break the first element at the initial state is

$$\sigma _1=\left(\mathrm{ln}\left(1\frac{1}{N_0}\right)\right)^{\frac{1}{\rho }}\left(\frac{1}{N_0}\right)^{\frac{1}{\rho }},N_01.$$

So $`\sigma _1^\rho \frac{1}{N_0}`$, and

$$q(1,2)=e^{(2^\rho 1)\frac{1}{N_0}}.$$
(14)
If the value of $`q(1,2)`$ in Eq. (14) is bigger than $`0.5`$, then the fracture will likely consist of one single event of rupture and the avalanche will not progress beyond this first failure. As shown above, for $`\rho =2`$ a system size of $`N_0=1000`$ elements is enough to produce first avalanches of size $`1`$ on average. Next, the driving force is increased and the new load acting on the surviving $`N_01`$ fibers is

$$\sigma _2=\left(\sigma _1^\rho \mathrm{ln}\left(1\frac{1}{N_01}\right)\right)^{\frac{1}{\rho }}\left(\frac{2}{N_0}\right)^{\frac{1}{\rho }},N_011$$

that is, $`\sigma _2^\rho =\frac{2}{N_0}`$. In the derivation of $`\sigma _2`$ it has been assumed that the second avalanche starts with high probability in a fiber different of that which is located just in the crack tip of the first avalanche. This approximation is well-justified by numerical simulations of this local load model for which it is known that the first breaking events are randomly distributed in the sample, leading to the appearance of several small cracks which much later coalesce and grow provoking the final breakdown. Now, there are $`(N_02)`$ elements bearing a load $`\sigma _2`$ and one element supporting a load $`\sigma _2+\sigma _1`$. The system will experience on average the final breakdown in this second breaking event if the probability

$$q^{}(1,2)=e^{(2^\rho 1)\frac{2}{N_0}}$$
(15)
is lower than $`0.5`$. Nevertheless, $`q(1,2)`$ and $`q^{}(1,2)`$ are not independent since the relation

$$q^{}(1,2)=\left(q(1,2)\right)^2$$
(16)
has to be verified. This dependence, together with the restrictions $`q(1,2)>0.5`$ and $`q^{}(1,2)<0.5`$ leads to $`0.5<q(1,2)<0.7`$ and $`0.25<q^{}(1,2)<0.5`$ as acceptable values for $`q(1,2)`$ and $`q^{}(1,2)`$ in order to have a two-step system collapse. This approximation seems reasonable taking into account the way in which numerical simulations of the continuous and the discrete versions of the probabilistic approach proceed. So, the condition for the first avalanche to be finite and that the second breaks the system is reduced to

$$N_02^{\rho +1}2$$
(17)
for $`q(1,2)=0.6`$, a value in the middle of its range. This relationship between $`N_0`$ and $`\rho `$ is confirmed in numerical simulations. The situation considered above corresponds to a very brittle rupture since the ratio of the size of the typical local damage that induces the system failure to the system size is very close to one $`(N_01)/N_0`$. Expression (17), although very simple, gives a rough estimate of the qualitative relation between the size of the system and the amount of disorder when the degree of brittleness of the fracture is high. The situation analyzed becomes important and more accurate as we go to larger system sizes because for the local load transfer rule and a larger number $`N`$ of fibers, the probability to find a weak region somewhere in the system is high and due to the local redistribution rule the crack will propagate fast enough to break the whole system in a few steps (strictly speaking, in the thermodynamic limit the critical load is zero).

## 4 Discussion and Conclusions

We have proposed a probabilistic approach to fiber bundle models of fracture. The cases of long range interactions among the fibers of the bundle and a local load sharing transfer scheme (short range interactions) were considered in order to investigate the behavior of the system near the global breakdown of the material. The results obtained for the local scheme indicate that the system undergoes a kind of first-order phase transition in agreement with previous reports . The system fails with no significant precursors announcing the incipient rupture, and no scaling relations can be found even for large system sizes and strong disorder. We illustrated this behavior by considering the situation of a two-step global failure. In this particular case, several quantities have a discrete jump at the critical point, like the branching ratio $`\zeta `$ which goes from zero to $`\frac{N_02}{N_01}`$, i.e., from zero to one in the limit of $`N_0\mathrm{}`$.

The scenario is quite different for the global load sharing FMB. In this case, the type of phase transition is not very clear. There are basically two possible ways in which such a transition may occur. In a first-order phase transition, we find discontinuous behavior in various quantities as we pass through the critical point. Contrary to this case, in a second-order (or continuous) phase transition, the fluctuations are correlated over all length scales, therefore, not only does the correlation length diverges continuously as the critical point is approached, but also other quantities show scaling . A particular case is the first-order phase transition close to a spinodal-like instability, for which scaling relations can be obtained despite that some macroscopic quantities, like the elastic modulus, have a discrete, finite jump at the critical point . The question then is whether the global load sharing FBM is a case of a second-order transition or a first-order spinodal transition.

The results obtained with the probabilistic approach seem not to support the claim of Ref. that the failure of a fiber bundle model under a global load sharing transfer scheme behaves as a first-order phase transition close to a spinodal point. There, by simulating models of electric breakdown and fracture, the authors presented ample numerical and theoretical evidence of several scaling relations and of the discrete jump of the macroscopic properties. We have obtained the same scaling relation for the rate of fiber failures as the critical point is reached, as well as for the avalanche sizes, which also diverge at the transition. Note that the scaling exponent derived from numerical simulations of the continuous version of the probabilistic approach fits very well the mean field result $`\gamma =\frac{1}{2}`$. Besides, the fraction of unbroken fibers just before the global rupture has a discontinuity, which is not size dependent for big systems as can be observed in Fig. 5. It has also been shown by analytic means that the distribution of avalanche sizes in the case of global load sharing transfer rule follows a universal power law with an exponent $`5/2`$ . However, we should notice that in driven disordered systems, the concepts related to spinodal nucleation are not sufficiently well established. In particular, it seems that the finite size of the system makes it possible that the last catastrophic avalanche provokes the elastic modulus to drop to zero with a finite, and very likely, N-dependent jump. Thus, from our point of view, that is not enough to set the conclusion that fracture can be described as a first-order phase transition.

Our alternative point of view is to consider the branching ratio defined previously as a measure of the distance of the system from the critical point. According to the results obtained, the branching ratio goes continuously from zero to one. The branching ratio measures the probability that a fiber that has failed triggers none, one or more breaking events. Therefore, a branching ratio equal to one at the critical point implies that, on average, one failure will induce another failure so that the system reaches a state where any perturbation propagates across the entire system, which is, in essence, a critical phase. This invokes a continuous phase transition as claimed in other analysis of fracture models . Note, additionally, that what changes discontinuously at $`\sigma _c`$ is the rate of change of $`\zeta `$ rather than $`\zeta `$ itself.

In summary, we have introduced a probabilistic approach to fiber bundle models which allows to avoid the fluctuations near the critical point. The different rupture behaviors of the system for the global load sharing scheme and the local one were shown. For the short-range interactions case, the rupture is no-doubt of the first-order phase transition type. For the global transfer rule, several scaling relations were obtained with mean field critical exponents. The branching ratio was defined as an appropriate order parameter. According to the results obtained, the branching ratio goes continuously from zero to one. This suggests that fracture in heterogeneous systems with long range interactions can be described as a phase transition of the second-order type, at least within the FBM picture.

## 5 ACKNOWLEDGMENTS

Y.M thanks A. Vespignani and H. J. Jensen for very useful and stimulating discussions. Y. M also thanks The Abdus Salam International Centre for Theoretical Physics for hospitality and financial support. This work was supported in part by the Spanish DGICYT under Project PB98-1594.