# Mesoscopic fluctuations of the ground state spin of a small metal particle

According to Hund’s rule, electrons in a partially filled shell in an atom form a many-body ground state with maximum possible spin. The maximum spin is preferred because it allows a maximally antisymmetric coordinate wavefunction in order to minimize the electrostatic repulsion between the electrons. In recent experiments, Hund’s rule was also observed in a cylindrically-shaped semiconductor quantum dot, or “artificial atom”. The close similarity with real atoms is due to the degeneracy of single-particle levels, caused by the the high degree of symmetry of the device.

In generic ultrasmall systems such as small metal grains, semiconductor quantum dots, or carbon nanotubes there is no systematic degeneracy due to a spherical (or cylindrical) symmetric potential. However, even in the absence of degeneracies, a nonzero value of the ground state spin may occur, as long as the gain in electrostatic energy is larger than the loss in kinetic energy when an antisymmetric coordinate ground state wavefunction is formed. Such a ground state is most likely to be detected in ultrasmall metal and semiconductor devices, since in those systems, unlike in macroscopic samples, the spacing between single-particle energy levels and the typical interaction energies can be larger than the temperature. In fact, the possibility of such a “weakly ferromagnetic” ground state has been suggested as an explanation for some recent experiments, that could not be explained by simple noninteracting models. In addition, a nonzero ground state spin from numerical simulations of a few particles in a chaotic dot, and a theory of spin polarization in larger dots were already mentioned in the literature. The stability of the zero spin ground state in a quantum dot was analyzed for weak interactions in Ref. .

In this paper, we consider small metal grains in the mesoscopic regime, in which fluctuations of wavefunctions and energy levels, caused by, e. g., disorder or an irregular shape, control the behavior of kinetic and interaction energies at the vicinity of the Fermi energy. As a result, the ground state spin becomes subject to sample-to-sample fluctuations. Then, the relevant quantity to consider is the statistical distribution of the ground state spin for an ensemble of small metal grains or chaotic quantum dots, rather than the spin of a specific sample.

Our starting point is a simple toy model that captures the essential mechanisms for mesoscopic fluctuations of the ground state spin. In second-quantized form, our model Hamiltonian $``$ reads

$``$ $`=`$ $`{\displaystyle \underset{n,m,\sigma }{}}c_{n,\sigma }^{}_0(n,m)c_{m,\sigma }^{}`$ (2)   
$`+uM{\displaystyle \underset{n}{}}c_{n,}^{}c_{n,}^{}c_{n,}^{}c_{n,}^{},`$   
where $`c_{n,\sigma }^{}`$ ($`c_{n,\sigma }^{}`$) is the creation (annihilation) operator for an electron with spin $`\sigma `$ at site $`n`$. The indices $`m,n`$ are summed over $`M`$ sites. The Hamiltonian $`_0`$ contains the kinetic energy and the impurity potential. We describe the electron-electron interaction by an on-site (Hubbard) interaction, $`uM`$. While the long-range Coulomb interaction can be trivially included via a charging energy, the model (2) does not include the Coulomb interaction at intermediate distances, which leads to a Gaussian level spacing distribution at the Fermi energy. In this work, we report a calculation of the ground state spin of the Hamiltonian (2) using a restricted version of the Hartree-Fock (HF) approximation with a random-matrix assumption for the eigenvalues and eigenvectors of the self-consistent HF Hamiltonian.

We first present our main result. It consists of an equation that relates the candidate ground state energies $`E_\mathrm{G}(s)`$ for different values of the total spin $`s`$ in terms of eigenvalues $`\epsilon _\mu ^0`$ of a hermitian random matrix with level spacing $`\mathrm{\Delta }`$, the interaction parameter $`\lambda =u/\mathrm{\Delta }`$, and a (nonuniversal) numerical constant $`c`$ that describes the density response to a local perturbation of the impurity potential in $`_0`$,

$`E_\mathrm{G}(s)E_\mathrm{G}(s_0)`$ $`=`$ $`{\displaystyle \underset{\mu =1}{\overset{s}{}}}(\epsilon _{N+\mu +2s_0}^0\epsilon _{N+1\mu }^0)`$ (4)   
$`\lambda \mathrm{\Delta }\left[s^2s_0^2+{\displaystyle \frac{2(ss_0)}{\beta (1\lambda ^2c^2)}}\right].`$   
The total number of electrons is $`2(N+s_0)`$, $`s_0`$ being $`0`$ or $`1/2`$. The spin of the true ground state is found by minimizing Eq. (4) with respect to $`s`$. The parameter $`\beta =1`$ ($`2`$) if time-reversal symmetry is present (absent). The effect of spin orbit coupling and Zeeman splitting is not included here. (The case $`\beta =2`$ is only relevant for semiconductor quantum dots in a weak magnetic field, that affects orbital motion, but causes no Zeeman splitting. It is not relevant for small metal grains, as laboratory magnetic fields do not affect orbital motion in this case.) Equation (4) reflects the competition between kinetic energy (first term on the r.h.s.), which favors small $`s`$ and the on-site interaction (second term), which favors finite $`s`$. The interaction term, in turn, consists of two parts: A term quadratic in $`s`$, which describes the exchange interaction, and a term linear in $`s`$, which describes the additional “dressed” Coulomb repulsion of two particles with the same spatial wavefunction. For large $`s`$, the contribution from the kinetic term is approximately $`s^2\mathrm{\Delta }`$, so that for $`u\mathrm{\Delta }`$ a finite fraction of the total number of spins will align, rather than a small number of spins as in the case $`u<\mathrm{\Delta }`$. The instability at $`u=\mathrm{\Delta }`$ is known as the Stoner instability. In Ref. , a result similar to Eq. (4) was obtained for $`s=1`$, but with a different and fluctuating interaction term. The difference is due to the absence of a self-consistent approximation scheme in Ref. .

As a consequence of the additional dressed Coulomb repulsion of particles with the same wavefunction, we find that already for interaction strengths considerably below the Stoner instability $`u=\mathrm{\Delta }`$, there is an appreciable probability of nonzero ground-state spin. This is illustrated in Fig. 1, where the distribution of the ground state spin at three different values of the interaction parameter $`\lambda `$ is shown: Already at the quite modest interaction strength $`u0.4\mathrm{\Delta }`$ a ground state spin $`s=1`$ is more likely than $`s=0`$.

The effect of a weak magnetic field is twofold: First, it changes the statistics of the $`\epsilon _\mu ^0`$, and, second, it suppresses the interference in the “Cooper channel”, leading to a factor of two reduction of the interaction term linear in $`s`$ \[last term in Eq. (4)\]. Both effects favor lower $`s`$ than without a magnetic field. However, even in the absence of a magnetic field we expect that, similar to 3D metals, inclusion of the electron–electron interaction beyond the HF approximation will also lead to a suppression of the interference in the Cooper channel (logarithmically in the system size $`M`$), and hence to a prefactor in that term that is smaller than 2.

Let us now turn to the details of our calculation. To find the ground state of the Hamiltonian (2) we use a simplified version of the HF approximation: We assume that the ground state has the form of a Slater determinant of single-particle wave functions $`\psi _{\mu ,}`$ and $`\psi _{\mu ,}`$ of particles which have either spin up or spin down. In this case the self-consistent HF equations read

$`_{\mathrm{HF},\sigma }\psi _{\mu ,\sigma }`$ $`=`$ $`\epsilon _{\mu ,\sigma }\psi _{\mu ,\sigma };`$ (5)   
$`_{\mathrm{HF},\sigma }(n,m)`$ $`=`$ $`_0(n,m)+uM\rho _\sigma (n)\delta (n,m),`$ (6)   
$`\rho _\sigma (n)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}f_{\mu ,\sigma }|\psi _\mu (n)|^2.`$ (7)   
The occupation number $`f_{\mu ,\sigma }`$ is $`1`$ ($`0`$) if the level $`\mu ,\sigma `$ is occupied (unoccupied) and $`\delta (n,m)`$ is the Kronecker delta function. The ground state energy $`E_\mathrm{G}`$ is given by

$$E_\mathrm{G}=\underset{\mu ,\sigma }{}f_{\mu ,\sigma }\epsilon _{\mu ,\sigma }uM\underset{n}{}\rho _{}(n)\rho _{}(n).$$
(8)   
Our strategy is as follows: We start from a reference state with zero spin, in which $`N`$ particles of each spin are placed in the same levels $`\epsilon _{\mu ,}=\epsilon _{\mu ,}`$ and with the same wavefunctions $`\psi _{\mu ,}=\psi _{\mu ,}`$. We assume, that for this symmetric case the eigenvectors and the eigenvalues of $`_{\mathrm{HF}}`$ are distributed like those of a random matrix, except that the energy levels below $`E_F`$ are shifted upwards, by a small constant amount, relative to the levels above $`E_F`$, see Eq. (22) below. (If we would have included long-range Coulomb interactions via a charging energy, the shift would have been much larger and in the opposite direction. Omission of the charging energy has no consequence in our case, as we compare ground states with the same number of particles.) The assumption that the single-particle eigenvalues and wavefunctions in a self-consistent potential for a mesoscopic system below and above $`E_F`$ obey random matrix statistics, even though they may be quite different from their counterparts in the noninteracting system, was checked numerically for short range interaction models somewhat similar to ours. The energy shifts in our case result from the spin-degeneracy, which was not present in these calculations.

Starting from this reference state, we build other states by the subsequent addition and removal of electrons. We first discuss the addition of a single up spin in the $`(N+1)`$st level. The first question that needs to be answered is how this addition affects the self-consistent density $`\rho _\sigma (n)`$. The density change $`\delta \rho _{}(n)`$ consists of a direct and an induced contribution, while $`\delta \rho _{}(n)`$ has an induced density shift only,

$`\delta \rho _{}`$ $`=`$ $`\delta \rho _{,\mathrm{dir}}+\delta \rho _{,\mathrm{ind}},\delta \rho _{,\mathrm{dir}}(n)=|\psi _{N+1,}(n)|^2,`$ (9)   
$`\delta \rho _{}`$ $`=`$ $`\delta \rho _{,\mathrm{ind}}.`$ (10)   
Since the density shifts change the HF Hamiltonians by an amount $`\delta _{\mathrm{HF},\sigma }(n,m)=uM\delta \rho _\sigma (n)\delta (n,m)`$, we obtain the following self-consistency equations for $`\delta \rho _\sigma `$,

$`\delta \rho _{\sigma ,\mathrm{ind}}(n)`$ $`=`$ $`2uM\text{Re}{\displaystyle \underset{\mu ,\nu ,m}{}}f_{\mu ,\sigma }(1f_{\nu ,\sigma })\delta \rho _\sigma (m)`$ (12)   
$`\times {\displaystyle \frac{\psi _{\mu ,\sigma }^{}(n)\psi _{\mu ,\sigma }^{}(m)\psi _{\nu ,\sigma }^{}(m)\psi _{\nu ,\sigma }^{}(n)}{\epsilon _{\mu ,\sigma }\epsilon _{\nu ,\sigma }}}.`$   
Both $`\delta \rho _{}`$ and $`\delta \rho _{}`$ are of order $`1/M`$. In Eq. (12) we have computed the induced density change to first order in $`\delta _{\mathrm{HF}}`$. Higher order terms do not contribute to $`\delta \rho _{\mathrm{ind},\sigma }`$ to order $`1/M`$ and are neglected. To evaluate Eq. (12), we first sum the r.h.s. over the space index $`m`$ and then over the energy levels $`\epsilon _\mu `$ and $`\epsilon _\nu `$. Because the eigenfunction has a random sign, a single term in the latter summation is of order $`\lambda M^{3/2}\mathrm{\Delta }/(\epsilon _\mu \epsilon _\nu )`$, which is not relevant in the limit $`M1`$, even if $`\epsilon _\mu `$ and $`\epsilon _\nu `$ are both close to the Fermi level $`E_F`$. For the summation over all levels we may perform an average over the wave functions (since the denominator is a slowly varying function of $`\mu `$ and $`\nu `$ away from the Fermi level). This average is done using that for general $`\mu \nu `$ and in the limit $`M1`$ one has,

$$\psi _{\mu ,\sigma }^{}(n)\psi _{\mu ,\sigma }^{}(m)\psi _{\nu ,\sigma }^{}(m)\psi _{\nu ,\sigma }^{}(n)=\frac{\delta (m,n)}{M^2}\frac{1}{M^3}.$$

Putting everything together, we find the following solution of the self-consistency equations (10) and (12),

$`\delta \rho _{}(n)`$ $`=`$ $`{\displaystyle \frac{1}{1(c\lambda )^2}}\left(|\psi _{N+1,}(n)|^2{\displaystyle \frac{1}{M}}\right)+{\displaystyle \frac{1}{M}},`$ (13)   
$`\delta \rho _{}(n)`$ $`=`$ $`{\displaystyle \frac{c\lambda }{1(c\lambda )^2}}\left(|\psi _{N+1,}(n)|^2{\displaystyle \frac{1}{M}}\right).`$ (14)   
where $`c`$ is a numerical constant of order unity defined by

$$c=\underset{M\mathrm{}}{lim}\frac{2\mathrm{\Delta }}{M}_{\mathrm{}}^{E_F}𝑑\epsilon _1_{E_F}^{\mathrm{}}𝑑\epsilon _2\frac{\rho (\epsilon _1)\rho (\epsilon _2)}{\epsilon _2\epsilon _1}.$$
(15)   
In this equation, $`\rho (\epsilon )`$ is the mean density of HF energy levels. \[The mean level spacing $`\mathrm{\Delta }`$ is taken at the Fermi energy, $`\mathrm{\Delta }=1/\rho (\epsilon _F)`$.\] The constant $`c`$ gives the linear density response $`\delta \rho (n)`$ to a shift of the impurity potential $`_0(n,n)`$ at that same site, $`\delta \rho (n)=c\delta _0(n,n)/M\mathrm{\Delta }`$, which can be verified using first order perturbation theory with respect to $`\delta _0`$. Notice that $`c`$ is not a universal constant, but depends on an integration of the density of states over the entire bandwidth. For example, for the Wigner semicircular density of states we find $`c=4/3`$ if the Fermi energy $`E_F`$ is at the band center and $`c0`$ if $`E_F`$ is at a band edge. Equation (14) expresses that the interaction enhances the fluctuations of the spin density: if $`|\psi _{N+1,}(n)|^2`$ is larger than average, the on-site repulsion reduces $`\rho _{}(n)`$, which in turn causes an increase of $`\rho _{}(n)`$, and so on. At the same time, the interaction reduces fluctuations of the charge density $`\rho _{}+\rho _{}`$.

At $`\lambda c=1`$, which may occur before the Stoner instability $`\lambda =1`$ if $`c>1`$, the density changes diverge. Although this instability signals a breakdown of our approach, it is not clear whether it will also lead to a true macroscopic ground state spin. Below, we restrict our discussion to the case $`\lambda c<1`$.

Next we address the HF energy levels $`\epsilon _{\mu ,\sigma }`$ and find

$`\delta \epsilon _{\mu ,}`$ $`=`$ $`\lambda \mathrm{\Delta }{\displaystyle \frac{2}{\beta }}{\displaystyle \frac{c\lambda }{1(c\lambda )^2}}\delta _{\mu ,N+1},`$ (16)   
$`\delta \epsilon _{\mu ,}`$ $`=`$ $`\lambda \mathrm{\Delta }\left(1+{\displaystyle \frac{2}{\beta }}{\displaystyle \frac{1}{1(c\lambda )^2}}\delta _{\mu ,N+1}\right).`$ (17)   
The shift of $`\epsilon _{N+1,\sigma }`$ is extra large, since for that level the interaction effects are enhanced by the spatial fluctuations of the wavefunction. Equation (17) is the result of first order perturbation theory in $`\delta \rho _\sigma `$; second order perturbation theory $`\delta \epsilon _{\mu ,\sigma }`$ gives a correction of order $`\lambda ^2\mathrm{\Delta }\mathrm{ln}M/M`$, which we may neglect in the limit $`M1`$. In the same way, one finds that the changes in each individual wave functions is not significant for $`M1`$.

Finally, we consider the change in the ground state energy $`E_\mathrm{G}`$. Because the summation over $`\mu `$ in Eq. (8) extends over $`𝒪(M)`$ levels, it is important to follow the shifts in the HF levels to second order perturbation theory, although this level of accuracy was not needed for the shift of each level individually, cf. Eq. (17). Putting everything together, we find that

$$\delta E_\mathrm{G}=\epsilon _{N+1,}c\lambda ^2\mathrm{\Delta }[\beta (1c^2\lambda ^2)]^1.$$
(18)   
(No terms proportional to $`\mathrm{log}M`$ appear here since they cancel in the summation over the energy levels.) With the help of Eq. (17), we can rewrite Eq. (18) as $`\delta E_\mathrm{G}=\epsilon _{N+1,}+\frac{1}{2}\delta \epsilon _{N+1,}`$, which is the average of the energy for the newly occupied level $`\epsilon _{N+1,}`$ before and after its occupation. This may be interpreted as a simple extension of Koopmans’ theorem  to the present case, where the modification in each one-electron wavefunction is small (of relative order $`M^{1/2}`$), but the resulting contribution to $`\delta E_\mathrm{G}`$ cannot be neglected to the order we are interested in. In the usual form of Koopmans’ theorem, where one ignores any change in the one-particle wavefunctions, the HF energies of the lowest unoccupied state before addition of the electron, and of the highest occupied state after addition are identical. (The usual Koopmans’ theorem is correct for an infinite system, in general, or for a translationally invariant finite system, as the one-electron states are trivially plane waves in that case.) The simple extension of Koopmans’ theorem also works, in our model, for the addition of several electrons.

We have repeated these calculations for the addition of two electrons with opposite spin in the $`(N+1)`$th level,

$`\delta \rho _\sigma (n)`$ $`=`$ $`{\displaystyle \frac{1}{M}}+{\displaystyle \frac{1}{1+c\lambda }}\left(|\psi _{N+1,}(n)|^2{\displaystyle \frac{1}{M}}\right),`$ (19)   
$`\delta \epsilon _{\mu ,\sigma }`$ $`=`$ $`\lambda \mathrm{\Delta }\left(1+{\displaystyle \frac{2}{\beta }}{\displaystyle \frac{1}{1+c\lambda }}\delta _{\mu ,N+1}\right),`$ (20)   
$`\delta E_\mathrm{G}`$ $`=`$ $`2\epsilon _{N+1}+{\displaystyle \frac{2}{\beta }}\lambda \mathrm{\Delta }{\displaystyle \frac{1}{1+c\lambda }}.`$ (21)   
As in the case of the addition of a single particle, the individual wavefunctions do not change to order $`M^{1/2}`$. Equation (20) allows us to find the statistics of the HF energy levels $`\epsilon _{\mu ,\sigma }`$ in our reference system with $`N`$ electrons of each spin: The only distribution that is consistent both with the assumption that the $`\epsilon _{\mu ,\sigma }`$ obey random-matrix statistics away from the Fermi level and with the shifts of Eq. (20) is one where the $`\epsilon _{\mu ,\sigma }`$ have the form

$$\epsilon _{\mu ,\sigma }=\epsilon _\mu ^0+\frac{2}{\beta }f_\mu \frac{\lambda \mathrm{\Delta }}{1+\lambda c},$$
(22)   
where the $`\epsilon _\mu ^0`$ have random matrix statistics and $`f_\mu =1`$ ($`0`$) if the level $`\mu `$ is (un)occupied. In other words, the distribution of the $`\epsilon _\mu `$ is the same as that of the eigenvalues of a random matrix, with all occupied levels shifted upwards by an amount $`(2/\beta )\lambda \mathrm{\Delta }/(1+c\lambda )`$.

With the knowledge we have gained above, there is little work left for the calculation of our main result, Eq. (4). Some remarks about the validity of this result are appropriate. First, to make a connection between our random matrix toy model and a laboratory made quantum dot we must identify $`M=(L/\lambda _F)^3`$ as we expect the length scale for wavefunction correlations and the range of the screened Coulomb interaction to be of order of the Fermi wavelength $`\lambda _F`$. Second, while our solution is complete within the HF approximation, one must bear in mind that this approximation scheme does not include correlation effects, such as the Cooper channel renormalization, as we discussed below Eq. (4). Those correlation effects are not expected to affect our result to first order in $`\lambda `$, but it can not be excluded that they are important in the higher order terms in Eq. (4), which involve the factor $`c\lambda `$.

We close this paper with a discussion of the physical consequences of a ground state with spin $`s>1/2`$ and of the experimental situations in which it can be observed.

First, the temperature $`T`$ needs to be smaller than the separation of the ground states for different spin $`s`$. This separation, which is typically smaller than the single particle level spacing $`\mathrm{\Delta }`$, is a fluctuating quantity. Very small values are possible, because, unlike in noninteracting random systems, there is no level repulsion if states of different spin are involved.

For sufficiently low $`T`$, the magnetization of the grain is proportional to the spin $`s`$ of the ground state. However, $`s`$ will also affect other properties which are more easily accessible in an experiment, like current–voltage characteristics: A nonzero ground state spin can serve to explain the absence of an even-odd structure in the addition spectra of Coulomb-blockaded quantum dots, or the presence of kinks in the parametric dependence of Coulomb blockade peak positions, as was noted in Ref. . Spin is also relevant for conductance measurements at a finite bias voltage, which allow for a “spectroscopy” of the quantum dot or metal grain. In the presence of a magnetic field, the ground state is split by the Zeeman energy, and the differential conductance will show two peaks, whose relative intensity differs by a factor $`2s_{N+1}+1`$ or $`2s_N+1`$, depending on whether the tunneling onto or from the grain is the faster process. Even without an external magnetic field the ground state may be split, e.g., by spin-orbit coupling or magnetic impurities, and thus give rise to a multiplet of peaks in the differential conductance. The peak separation within a multiplet is controlled by the strength of the splitting mechanism and may be much smaller than $`\mathrm{\Delta }`$.

We thank D. Davidović and M. Tinkham for stimulating discussions and for sharing their experimental results with us. When this work was nearly completed, we learned of related work by H. U. Baranger, D. Ullmo, and L. I. Glazman (Ref. ), in which some similar results were obtained. We thank H. U. Baranger for discussions on these points. This work was supported in part by the NSF through the Harvard MRSEC (grant DMR 98-09363), and by grants DMR 94-16910, DMR 96-30064, DMR 97-14725.