# Filling dependence of the Mott transition in the degenerate Hubbard model

## I Introduction

The Hubbard Hamiltonian is a simple model for studying strongly interacting systems. In particular it is used to investigate the Mott-Hubbard metal-insulator transition in half-filled systems. It is clear that for strong correlations such a system should be insulating, since in the atomic limit the states with exactly one electron per lattice site are energetically favored, while all other states are separated from those by a Coulomb gap. For a generalized Hubbard model with degenerate orbitals the same argument implies that for strong correlations not only the half-filled, but all integer filled systems will become Mott-Hubbard insulators. It is then natural to ask how the location of the transition depends on the filling.

As an example we consider a Hamiltonian describing the alkali doped Fullerides. It comprises the three-fold degenerate $`t_{1u}`$ orbital and the Coulomb interaction $`U`$ between the electrons on the same molecule. Using this Hamiltonian, we have recently shown that, although $`U`$ is substantially larger than the band width $`W`$, K<sub>3</sub>C<sub>60</sub> is not a Mott insulator but a (strongly correlated) metal. Prompted by the synthesis of an isostructural family of doped Fullerenes A<sub>n</sub>C<sub>60</sub> with different fillings $`n`$, we now address the question of the Mott transition in integer doped Fullerides. For these systems we have the interesting situation that for fillings $`n=`$ 1, 2, 3, 4, and 5 calculations in the local density approximation predict them all to be metallic, while in Hartree-Fock they all are insulators. Performing quantum Monte Carlo calculations for the degenerate Hubbard model at different fillings and for values of $`U`$ typical for the Fullerides, we find that all the systems are close to a Mott transition, with the critical correlation strength $`U_c`$ at which the transition takes place strongly depending on the filling $`n`$. More generally, our results show how, for an otherwise identical Hamiltonian, the location of the Mott transition $`U_c`$ depends on the filling. $`U_c`$ is largest at half-filling and decreases for fillings smaller or larger than half. We contrast these findings with the results from Hartree-Fock calculations which predict a much too small $`U_c`$ and show almost no doping dependence. We give an interpretation of the results of the quantum Monte Carlo calculations extending the hopping argument introduced in Ref.  to arbitrary integer fillings. Despite the crudeness of the argument it explains the doping dependence found in quantum Monte Carlo. We therefore believe that our simple hopping argument captures the basic physics of the doping dependence of the Mott transition in degenerate systems.

In section II we introduce the model Hamiltonian for doped Fullerenes with a three-fold degenerate $`t_{1u}`$ band. We discuss the fixed-node approximation used in the diffusion Monte Carlo calculations, present the results of our quantum Monte Carlo calculations, and contrast them to the result of Hartree-Fock calculations. Section III gives an interpretation of the results of our calculations in terms of intuitive hopping arguments. We introduce the many-body enhancement of the hopping matrix elements, which explains how orbital degeneracy $`N`$ helps to increase the critical $`U`$ at which the Mott transition takes place and we analyze how frustration leads to an asymmetry of the critical $`U`$ for fillings $`n`$ and $`2Nn`$. A summary in Sec. IV closes the presentation.

## II Model calculations

### A Model Hamiltonian

Solid C<sub>60</sub> is characterized by a very weak inter-molecular interaction. Therefore the molecular levels merely broaden into narrow, well separated bands. The conduction band originates from the lowest unoccupied molecular orbital, the 3-fold degenerate $`t_{1u}`$ orbital. To get a realistic, yet simple description of the electrons in the $`t_{1u}`$ band, we use a Hubbard-like model that describes the interplay between the hopping of the electrons and their mutual Coulomb repulsion:

$`H`$ $`=`$ $`{\displaystyle \underset{ij}{}}{\displaystyle \underset{mm^{}\sigma }{}}t_{im,jm^{}}c_{im\sigma }^{}c_{jm^{}\sigma }^{}`$ (2)   
$`+U{\displaystyle \underset{i}{}}{\displaystyle \underset{(m\sigma )<(m^{}\sigma ^{})}{}}n_{im\sigma }n_{im^{}\sigma ^{}}.`$   
The sum $`ij`$ is over nearest-neighbor sites of an fcc lattice. The hopping matrix elements $`t_{im,jm^{}}`$ between orbital $`m`$ on molecule $`i`$ and orbital $`m^{}`$ on molecule $`j`$ are obtained from a tight-binding parameterization. The molecules are orientationally disordered, and the hopping integrals are chosen such that this orientational disorder is included. The band width for the infinite system is $`W=0.63eV`$. The on-site Coulomb interaction is $`U1.2eV`$. The model neglects multiplet effects, but we remark that these tend to be counteracted by the Jahn-Teller effect, which is also not included in the model.

We will investigate the above Hamiltonian for different integer fillings $`n`$ of the $`t_{1u}`$ band. The corresponding Hamiltonians describe a hypothetical family of doped Fullerides A<sub>n</sub>C<sub>60</sub> with space group Fm$`\overline{3}`$m, i.e. an fcc lattice with orientationally disordered C<sub>60</sub> molecules. In the calculations we use the on-site Coulomb interaction $`U`$ as a parameter to drive the system across the Mott transition.

### B Quantum Monte Carlo method

As the criterion for determining the metal-insulator transition we use the opening of the gap

$$E_g=E(N+1)2E(N)+E(N1),$$
(3)   
where $`E(N)`$ denotes the total energy of a cluster of $`N_{mol}`$ molecules with $`N`$ electrons in the $`t_{1u}`$ band. Since we are interested in integer filled systems, $`N=nN_{mol}`$, $`n`$ an integer. For calculating the energy gap (3) we then have to determine ground-state energies for the Hamiltonian (2). This is done using quantum Monte Carlo. Starting from a trial function $`|\mathrm{\Psi }_T`$ we calculate

$$|\mathrm{\Psi }^{(n)}=[1\tau (Hw)]^n|\mathrm{\Psi }_T,$$
(4)   
where $`w`$ is an estimate of the ground-state energy. The $`|\mathrm{\Psi }^{(n)}`$ are guaranteed to converge to the ground state $`|\mathrm{\Psi }_0`$ of $`H`$, if $`\tau `$ is sufficiently small and $`|\mathrm{\Psi }_T`$ is not orthogonal to $`|\mathrm{\Psi }_0`$. Since we are dealing with Fermions, the Monte Carlo realization of the projection (4) suffers from the sign-problem. To avoid the exponential decay of the signal-to-noise ratio we use the fixed-node approximation. For lattice models this involves defining an effective Hamiltonian $`H_{\mathrm{eff}}`$ by deleting from $`H`$ all nondiagonal terms that would introduce a sign-flip. Thus, by construction, $`H_{\mathrm{eff}}`$ is free of the sign-problem. To ensure that the ground-state energy of $`H_{\mathrm{eff}}`$ is an upper bound of the ground state of the original Hamiltonian $`H`$, for each deleted hopping term, an on-site energy is added in the diagonal of $`H_{\mathrm{eff}}`$. Since $`|\mathrm{\Psi }_T`$ is used for importance sampling, $`H_{\mathrm{eff}}`$ depends on the trial function. Thus, in a fixed-node diffusion Monte Carlo calculation for a lattice Hamiltonian, we choose a trial function and construct the corresponding effective Hamiltonian, for which the ground-state energy $`E_{\mathrm{FNDMC}}`$ can then be determined by diffusion Monte Carlo without a sign-problem.

For the trial function we make the Gutzwiller Ansatz

$$|\mathrm{\Psi }(U_0,g)=g^D|\mathrm{\Phi }(U_0),$$
(5)   
where the Gutzwiller factor reflects the Coulomb term $`UD=Un_{im\sigma }n_{im^{}\sigma ^{}}`$ in the Hamiltonian (2). $`|\mathrm{\Phi }(U_0)`$ is a Slater determinant that is constructed by solving the Hamiltonian in the Hartree-Fock approximation, replacing $`U`$ by a variational parameter $`U_0`$. Details on the character of such trial functions and the optimization of Gutzwiller parameters can be found in Ref. .

To check the accuracy of the fixed-node approximation, we have determined the exact ground-state energies for a (small) cluster of four C<sub>60</sub> molecules using the Lanczos method. For systems with different on-site Coulomb interaction (Table I) and varying number of electrons (Table II), we consistently find that the results of fixed-node diffusion Monte Carlo are only a few $`meV`$ above the exact energies.

### C Quantum Monte Carlo results

Since the quantum Monte Carlo calculations are for finite clusters of $`N_{mol}`$ molecules, we have to extrapolate the calculated energy gaps to infinite system size. An obvious finite-size effect is the fact that the one-particle spectrum is discrete, hence there can be a gap, even for $`U=0`$. Furthermore, in evaluating (3), we add and subtract one electron to a finite system. Even if we distribute the extra charge uniformly over all molecules, there will be an electrostatic contribution of $`U/N_{mol}`$ to the gap. We therefore introduce

$$E_G=E_gE_g(U=0)\frac{U}{N_{mol}}.$$
(6)   
These corrections are expected to improve the finite-size extrapolation. In practice they turn out to be quite small. For a cluster of 32 C<sub>60</sub> molecules, e.g., $`E_g(U=0)`$ is typically already less than $`10meV`$. In the thermodynamic limit both correction terms vanish, as they should.

The results of the quantum Monte Carlo calculations are shown in Fig. 1. Plotting the finite-size corrected gap $`E_G`$ for different values of the Coulomb interaction $`U`$ versus the inverse system size $`1/N_{mol}`$, we read off where the gap starts to open. For the system with one electron per molecule the gap opens around $`U_c0.75\mathrm{}1.00eV`$. At filling 2 the transition takes place later, at $`U_c1.25\mathrm{}1.50eV`$. For both, filling 3 and 4 we find the largest critical $`U`$: $`U_c1.50\mathrm{}1.75eV`$. For the system with 5 electrons per molecule the gap opens around $`U_c1.00\mathrm{}1.25eV`$. The results are summarized in Fig. 2. Thus we find that for an otherwise identical Hamiltonian the critical $`U`$ for the Mott transition depends strongly on the filling. $`U_c`$ is largest at half-filling and decreases away from half-filling. The decrease in $`U_c`$ is, however, not symmetric around half-filling. It is more pronounced for fillings $`<3`$ than for fillings $`>3`$.

We note that the opening of the gap is accompanied by a change in the character of that trial function which yields the lowest energy in the fixed-node approximation. For small $`U`$, where the system is still in the metallic regime, paramagnetic trial functions with small $`U_0`$ (see the discussion after eqn. (5)) are best. When the gap starts to open, trial functions with larger $`U_0`$, which have antiferromagnetic character, give lower energies. The corresponding Slater determinants $`|\mathrm{\Phi }(U_0)`$ describe a Mott insulator in Hartree-Fock approximation.

### D Hartree-Fock calculations

It is instructive to compare the results of the quantum Monte Carlo calculations with the predictions of Hartree-Fock theory. Figure 3 shows the the gap $`E_g`$ calculated for the Hamiltonian (2) within the Hartree-Fock approximation for the different integer fillings. Compared with quantum Monte Carlo, the gap opens much too early, around $`U0.4eV`$ ($`U/W0.65`$). Furthermore, there is only a very weak doping dependence: $`U_c`$ somewhat increases with the filling — in qualitative disagreement with the quantum Monte Carlo results. This failure is a direct consequence of the mean-field approximation. In Hartree-Fock the only way to avoid multiple occupancies of the molecules, in order to reduce the Coulomb repulsion, is to renormalize the on-site energy for the orbitals, thereby localizing the electrons in certain orbitals. For the Hamiltonian (2) this on-site energy is, apart from a trivial offset, given by $`\epsilon _{im\sigma }=U_{m^{}\sigma ^{}}n_{im^{}\sigma ^{}}n_{im\sigma }`$. Lowering the Coulomb energy in this way will, however, increase the kinetic energy. For small changes in the on-site energies this increase will scale like the inverse of the density of states at the Fermi level. This suggests that the critical $`U`$ should be the larger, the smaller the density of states at the Fermi level. Inspecting the density of states $`N(\epsilon )`$ for the non-interacting Hamiltonian (see e.g. Fig. 3 of Ref. ), we find that this is indeed the case: $`N(\epsilon )`$ slightly decreases with filling, explaining the corresponding increase in $`U_c`$. Hence the weak, but qualitatively wrong, doping dependence in Hartree-Fock can be understood as an effect of the small variation in the density of states of the non-interacting system.

## III Interpretation

### A Hopping enhancement

To find a simple interpretation for the doping dependence of the Mott transition we consider the limit of large Coulomb interaction $`U`$. In that limit the Coulomb energy dominates and we can estimate the energies entering the gap equation (3) by considering electron configurations in real space. According to the Hamiltonian (2) the contribution to the Coulomb energy from a molecule that is occupied by $`m`$ electrons is $`Um(m1)/2`$. Thus the energy of a system with filling $`n`$ is minimized for configurations with exactly $`n`$ electrons per molecule. The hopping of an electron to a neighboring molecule would cost the Coulomb energy $`U`$ and is therefore strongly suppressed in the large-$`U`$ limit. The energy for a cluster of $`N_{mol}`$ molecules with $`N=nN_{mol}`$ electrons (filling $`n`$) is then given by

$$E(N)=\frac{n(n1)}{2}N_{mol}U+𝒪(t^2/U),$$
(7)   
where $`t`$ is a typical hopping matrix element. Adding an extra electron increases the Coulomb energy by $`nU`$, removing an electron reduces it by $`(n1)U`$. But there will also be a kinetic contribution to the energy $`E(N\pm 1)`$, since the extra charge can hop without any additional cost in Coulomb energy. To estimate the kinetic energy we calculate the matrix element for the hopping of the extra charge to a neighboring molecule. This matrix element will of course depend on the arrangement of the other $`N`$ electrons. It is well known that for the non-degenerate Hubbard model a ferromagnetic arrangement of the spins is energetically favored (Nagaoka’s theorem), allowing the extra charge to hop without disturbing the background spins. For a degenerate Hubbard model, however, the hopping matrix element is larger e.g. for an antiferromagnetic arrangement of the background spins. This is illustrated in Fig. 4 a) for an extra electron in a system with filling 2. Now, instead of only the extra electron, any one out of the three equivalent electrons can hop to the neighboring molecule. Denoting the state with the extra electron on molecule $`i`$ by $`|i`$, we find that the second moment of the Hamiltonian $`i|H^2|i`$ is given by the number of hopping channels $`k`$ (in the present case $`k=3`$) times the number of (equivalent) nearest neighbors $`Z`$ times the single-electron hopping matrix element $`t`$ squared. Thus by inserting the identity in the form $`_j|jj|`$, where $`|j`$ denotes the state where any one of the electrons has hopped form molecule $`i`$ to the neighboring molecule $`j`$, we find

$$i|H|j=\sqrt{k}t,$$
(8)   
i.e. the hopping matrix element is enhanced by a factor of $`\sqrt{k}`$ over the one-particle hopping matrix element $`t`$. In a similar way we find for the system with an extra hole (Fig. 4 b) a hopping enhancement of $`\sqrt{k}`$ with $`k=2`$. The hopping enhancements for other fillings are listed in table III, where $`k_{}`$ denotes the enhancement for a system with an extra hole, and $`k_+`$ is for a system with an extra electron.

For a single electron the kinetic energy is of the order of $`W/2`$, where $`W`$ is the one-electron band width. The enhancement factor $`\sqrt{k}`$ in the many-body case then suggests that the kinetic energy for the extra charge is correspondingly enhanced, implying

$`E(N+1)`$ $``$ $`E(N)+nU\sqrt{k_+}W/2`$   
$`E(N1)`$ $``$ $`E(N)(n1)U\sqrt{k_{}}W/2.`$   
Combining these results we find

$$E_gU\frac{\sqrt{k_+}+\sqrt{k_{}}}{2}W,$$
(9)   
i.e. the hopping enhancement leads to a reduction of the gap described by the factor $`c=(\sqrt{k_+}+\sqrt{k_{}})/2`$. This reduction is largest ($`1.73`$) for n=3, and becomes smaller away from half-filling: $`c1.57`$ for $`n=2,\mathrm{\hspace{0.17em}4}`$, and $`c1.21`$ for fillings 1 and 5. Extrapolating (9) to intermediate $`U`$ we find that the gap opens for U larger than $`U_c=cW`$. Therefore the above argument predicts that the critical $`U`$ for the Mott transition depends strongly on the filling, with $`U_c`$ being largest at half-filling and decreasing away from half-filling. This is qualitatively the same behavior as we have found in the Monte Carlo calculations. We note, however, that the argument we have presented is not exact. First, the hopping of an extra charge against an antiferromagnetically ordered background will leave behind a trace of flipped spins. Therefore the analogy with the one-electron case for determining the kinetic energy in the large-$`U`$ limit is not exact. Second, using (9) for determining $`U_c`$ involves extrapolating the results obtained in the limit of large $`U`$ to intermediate values of the Coulomb interaction, where the Mott transition takes place. Finally, considering only one nearest neighbor in the hopping argument (cf. Fig. 4) implicitly assumes that we are dealing with a bipartite lattice, where all nearest neighbors are equivalent.

### B Origin of the asymmetry

To analyze the asymmetry in the gaps around half-filling we use the following exact relation for the kinetic energy in the limit of infinite $`U`$, which follows from an electron-hole transformation

$$T_{\genfrac{}{}{0pt}{}{min}{max}}(nN_{mol}\pm 1)=T_{\genfrac{}{}{0pt}{}{max}{min}}((2Nn)N_{mol}1).$$
(10)   
(Note how this symmetry is reflected in the hopping enhancements shown in Table III.) Since the gap for filling $`n`$ is given by

$$E_g(n)=U+T_{min}(nN_{mol}1)+T_{min}(nN_{mol}+1),$$

the asymmetry $`\mathrm{\Delta }=E_g(n)E_g(2Nn)`$ in the gaps can be written entirely in terms of energies for systems with an extra electron:

$$\mathrm{\Delta }=\begin{array}{ccc}T_{max}((2Nn)N_{mol}+1)& +& T_{min}(nN_{mol}+1)\\ T_{min}((2Nn)N_{mol}+1)& +& T_{max}(nN_{mol}+1).\end{array}$$

For a bipartite system the spectrum for a given filling will be symmetric, in particular $`T_{min}+T_{max}=0`$, and thus there will be no asymmetry in the gaps: $`\mathrm{\Delta }=0`$. Frustration breaks this symmetry. To study the effect of frustration we perform a Lanczos calculation in the large-$`U`$ limit, starting from a configuration $`|v_0`$ of the type shown in Fig. 4. The leading effect of frustration is given by the third moment, which already enters after the first Lanczos step. Diagonalizing the Lanczos matrix and expressing everything in terms of the moments of the Hamiltonian, the extreme eigenvalues are given by

$$\epsilon _{\genfrac{}{}{0pt}{}{max}{min}}=\frac{A_3\pm \sqrt{4A_2^3+A_3^2}}{2A_2},$$
(11)   
where $`A_k=v_0|H^k|v_0`$ denotes the $`k^{th}`$ moment of $`H`$, and $`A_1=v_0|H|v_0=0`$ for a state like in Fig. 4. From this expression it is clear that the “band width” $`\epsilon _{max}\epsilon _{min}`$ is essentially given by the second moment, and that an enhancement of $`A_2`$ by a factor of $`k`$ leads to an increase in the band width by a factor of $`\sqrt{k}`$, as already described above. The main effect of the third moment (i.e. of frustration) is to shift the extremal eigenvalues, where the shift is determined by the third moment.

To get a contribution to the third moment the initial state $`|v_0`$ must be recovered after three hops. This is only possible if the extra charge hops around a triangle, without changing spins along its path. For a state with an extra electron this means that one and the same electron has to perform the triangular hop. Therefore, even in the many-body case, for each such electron we get the same contribution to the third moment as in the single electron case. It therefore makes sense to write the third moment $`A_3(n)`$ for a system with $`nN_{mol}+1`$ electrons in terms of the third moment $`A_3^s`$ of the single electron problem: $`A_3(n)=\kappa _+(n)A_3^s`$, where $`\kappa _+(n)`$ describes the many-body effects, just like we introduced $`k_+(n)`$ to describe the many-body enhancement of the second moment. Using these definitions we find that the size of the asymmetry $`\mathrm{\Delta }`$ in the gaps can be estimated by the doping dependence of the (positive) enhancement factors $`\kappa _+(n)`$ of the third moment, while the overall sign is determined by the single-electron moments:

$$\mathrm{\Delta }\left(\frac{\kappa _+(n)}{k_+(n)}\frac{\kappa _+(2Nn)}{k_+(2Nn)}\right)\frac{A_3^s}{A_2^s}.$$
(12)   
To understand the doping dependence of $`\kappa _+/k_+`$ we proceed in two steps. First we observe that the upper limit for the number of different electrons that can perform a triangular hop is given by the number $`k_+`$ of electrons that can hop to a nearest neighbor. Hence, if frustration is not suppressed, $`\kappa _+/k_+=1`$. For filling $`n=1`$, $`N2`$ this upper limit can always be achieved without compromising large next-neighbor hopping by arranging the electrons in such a way as to avoid each other. This is shown in Fig. 5. For the corresponding filling $`2N1`$ the electrons can no longer be completely separated in that way. Thus the channel for triangular hops will be blocked by the the Pauli principle, reducing $`\kappa _+/k_+`$. In that way for the larger fillings frustration is reduced.

This reduction of frustration can already be seen in the simple model of a triangle with orbital degeneracy $`N=2`$ (cf. Fig. 5). Choosing matrix elements $`t=1`$ for hopping between like orbitals we find for filling $`n=1`$ a strong asymmetry $`T_{min}(3n+1)=2`$ and $`T_{max}(3n+1)=+4`$, while at filling $`n^{}=2Nn=3`$ there is no asymmetry in the extremal eigenvalues: $`T_{\genfrac{}{}{0pt}{}{max}{min}}(3n^{}+1)=\pm 2`$. We note that flipping one spin in the configuration for filling $`n=3`$ would allow for a triangular hop. In a Lanczos calculation this spin polarized configuration gives, however, only extremal eigenvalues $`T_{min}=2`$ and $`T_{max}=+1`$. The states described here for a triangle can be easily adapted to the situation in an fcc lattice, where the third moment involves hopping to the nearest neighbor sites, which form connected triangles.

From the non-interacting density of states for our model of the doped Fullerenes (cf. e.g. Fig. 3 of Ref. ) we see that both $`\epsilon _{min}`$ and $`\epsilon _{max}`$ are shifted upwards, compared to the center of the band, hence, looking at (11) we find that for a single electron the third moment is positive: $`A_3^s>0`$. Together with the reduction of the frustration for larger filling, we therefore expect from (12) that for the alkali doped Fullerenes $`E_g(n)>E_g(2Nn)`$; i.e. $`U_c(n)<U_c(2Nn)`$, as is found in the Monte Carlo calculations.

## IV Summary

Using quantum Monte Carlo, we have analyzed a model of alkali-doped Fullerenes and found that the Mott transition strongly depends on the (integer) filling $`n`$. $`U_c`$ is largest for $`n=3`$ and decreases away from half-filling. This result is qualitatively different from both, the results of density functional calculations in the local density approximation, and the results of Hartree-Fock calculations. The doping dependence of the Mott transition can be understood in terms of a simple hopping argument. The key observation is that, due to the orbital degeneracy, there are more hopping channels in the many-body than in the single-body case, thus leading to the degeneracy enhancement $`\sqrt{k}`$ discussed above. In addition, due to frustration, the gaps are not symmetric around half-filling.

The Gutzwiller approximation for a paramagnetic state also predicts a degeneracy enhancement. For a half-filled system, the predicted enhancement is, however, linear in the degeneracy $`(N+1)`$ instead of $`\sqrt{N}`$ as suggested by the hopping argument of Sec. III and as also found in infinite dimensions. The results of the Gutzwiller approximation are reproduced by a slave-boson calculation in the saddle-point approximation. In dynamical mean-field theory a degeneracy enhancement and a reduction of $`U_c`$ away from half-filling, similar to our result, is found.

## Acknowledgments

This work has been supported by the Alexander-von-Humboldt-Stiftung under the Feodor-Lynen-Program and the Max-Planck-Forschungspreis, and by the Department of Energy, grant DEFG 02-96ER45439.